{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SA13 Selection Functions\n",
    "## Depth maps and selection functions\n",
    "\n",
    "The simplest selection function available is the field MOC which specifies the area for which there is Herschel data. Each pristine catalogue also has a MOC defining the area for which that data is available.\n",
    "\n",
    "The next stage is to provide mean flux standard deviations which act as a proxy for the catalogue's 5$\\sigma$ depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This notebook was run with herschelhelp_internal version: \n",
      "708e28f (Tue May 8 18:05:21 2018 +0100)\n",
      "This notebook was executed on: \n",
      "2018-06-07 13:18:53.426428\n"
     ]
    }
   ],
   "source": [
    "from herschelhelp_internal import git_version\n",
    "print(\"This notebook was run with herschelhelp_internal version: \\n{}\".format(git_version()))\n",
    "import datetime\n",
    "print(\"This notebook was executed on: \\n{}\".format(datetime.datetime.now()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#%config InlineBackend.figure_format = 'svg'\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rc('figure', figsize=(10, 6))\n",
    "\n",
    "import os\n",
    "import time\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.table import Column, Table, join\n",
    "import numpy as np\n",
    "from pymoc import MOC\n",
    "import healpy as hp\n",
    "#import pandas as pd #Astropy has group_by function so apandas isn't required.\n",
    "import seaborn as sns\n",
    "\n",
    "import warnings\n",
    "#We ignore warnings - this is a little dangerous but a huge number of warnings are generated by empty cells later\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "from herschelhelp_internal.utils import inMoc, coords_to_hpidx, flux_to_mag\n",
    "from herschelhelp_internal.masterlist import find_last_ml_suffix, nb_ccplots\n",
    "\n",
    "from astropy.io.votable import parse_single_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "FIELD = 'SA13'\n",
    "#FILTERS_DIR = \"/Users/rs548/GitHub/herschelhelp_python/database_builder/filters/\"\n",
    "FILTERS_DIR = \"/opt/herschelhelp_python/database_builder/filters/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Depth maps produced using: master_catalogue_sa13_20180501.fits\n"
     ]
    }
   ],
   "source": [
    "TMP_DIR = os.environ.get('TMP_DIR', \"./data_tmp\")\n",
    "OUT_DIR = os.environ.get('OUT_DIR', \"./data\")\n",
    "SUFFIX = find_last_ml_suffix()\n",
    "#SUFFIX = \"20171016\"\n",
    "\n",
    "master_catalogue_filename = \"master_catalogue_{}_{}.fits\".format(FIELD.lower(), SUFFIX)\n",
    "master_catalogue = Table.read(\"{}/{}\".format(OUT_DIR, master_catalogue_filename))\n",
    "\n",
    "print(\"Depth maps produced using: {}\".format(master_catalogue_filename))\n",
    "\n",
    "ORDER = 10\n",
    "#TODO write code to decide on appropriate order\n",
    "\n",
    "field_moc = MOC(filename=\"../../dmu2/dmu2_field_coverages/{}_MOC.fits\".format(FIELD))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I - Group masterlist objects by healpix cell and calculate depths\n",
    "We add a column to the masterlist catalogue for the target order healpix cell <i>per object</i>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Add a column to the catalogue with the order=ORDER hp_idx\n",
    "master_catalogue.add_column(Column(data=coords_to_hpidx(master_catalogue['ra'],\n",
    "                                                       master_catalogue['dec'],\n",
    "                                                       ORDER), \n",
    "                                   name=\"hp_idx_O_{}\".format(str(ORDER))\n",
    "                                  )\n",
    "                           )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Convert catalogue to pandas and group by the order=ORDER pixel\n",
    "\n",
    "group = master_catalogue.group_by([\"hp_idx_O_{}\".format(str(ORDER))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Downgrade the groups from order=ORDER to order=13 and then fill out the appropriate cells\n",
    "#hp.pixelfunc.ud_grade([2599293, 2599294], nside_out=hp.order2nside(13))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II Create a table of all Order=13 healpix cells in the field and populate it\n",
    "We create a table with every order=13 healpix cell in the field MOC. We then calculate the healpix cell at lower order that the order=13 cell is in. We then fill in the depth at every order=13 cell as calculated for the lower order cell that that the order=13 cell is inside."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "depths = Table()\n",
    "depths['hp_idx_O_13'] = list(field_moc.flattened(13))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140074613482888-892934\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>idx</th><th>hp_idx_O_13</th></tr></thead>\n",
       "<tr><td>0</td><td>177665210</td></tr>\n",
       "<tr><td>1</td><td>177656680</td></tr>\n",
       "<tr><td>2</td><td>177665538</td></tr>\n",
       "<tr><td>3</td><td>177665204</td></tr>\n",
       "<tr><td>4</td><td>177656718</td></tr>\n",
       "<tr><td>5</td><td>177665539</td></tr>\n",
       "<tr><td>6</td><td>177656666</td></tr>\n",
       "<tr><td>7</td><td>177665205</td></tr>\n",
       "<tr><td>8</td><td>177665206</td></tr>\n",
       "<tr><td>9</td><td>177665207</td></tr>\n",
       "</table><style>table.dataTable {clear: both; width: auto !important; margin: 0 !important;}\n",
       ".dataTables_info, .dataTables_length, .dataTables_filter, .dataTables_paginate{\n",
       "display: inline-block; margin-right: 1em; }\n",
       ".paginate_button { margin-right: 5px; }\n",
       "</style>\n",
       "<script>\n",
       "\n",
       "var astropy_sort_num = function(a, b) {\n",
       "    var a_num = parseFloat(a);\n",
       "    var b_num = parseFloat(b);\n",
       "\n",
       "    if (isNaN(a_num) && isNaN(b_num))\n",
       "        return ((a < b) ? -1 : ((a > b) ? 1 : 0));\n",
       "    else if (!isNaN(a_num) && !isNaN(b_num))\n",
       "        return ((a_num < b_num) ? -1 : ((a_num > b_num) ? 1 : 0));\n",
       "    else\n",
       "        return isNaN(a_num) ? -1 : 1;\n",
       "}\n",
       "\n",
       "require.config({paths: {\n",
       "    datatables: 'https://cdn.datatables.net/1.10.12/js/jquery.dataTables.min'\n",
       "}});\n",
       "require([\"datatables\"], function(){\n",
       "    console.log(\"$('#table140074613482888-892934').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",
       "    $('#table140074613482888-892934').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "depths.add_column(Column(data=hp.pixelfunc.ang2pix(2**ORDER,\n",
    "                     hp.pixelfunc.pix2ang(2**13, depths['hp_idx_O_13'], nest=True)[0],\n",
    "                     hp.pixelfunc.pix2ang(2**13, depths['hp_idx_O_13'], nest=True)[1],\n",
    "                     nest = True),\n",
    "                     name=\"hp_idx_O_{}\".format(str(ORDER))\n",
    "                        )\n",
    "                 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140074613484064-312113\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>idx</th><th>hp_idx_O_13</th><th>hp_idx_O_10</th></tr></thead>\n",
       "<tr><td>0</td><td>177665210</td><td>2776018</td></tr>\n",
       "<tr><td>1</td><td>177656680</td><td>2775885</td></tr>\n",
       "<tr><td>2</td><td>177665538</td><td>2776024</td></tr>\n",
       "<tr><td>3</td><td>177665204</td><td>2776018</td></tr>\n",
       "<tr><td>4</td><td>177656718</td><td>2775886</td></tr>\n",
       "<tr><td>5</td><td>177665539</td><td>2776024</td></tr>\n",
       "<tr><td>6</td><td>177656666</td><td>2775885</td></tr>\n",
       "<tr><td>7</td><td>177665205</td><td>2776018</td></tr>\n",
       "<tr><td>8</td><td>177665206</td><td>2776018</td></tr>\n",
       "<tr><td>9</td><td>177665207</td><td>2776018</td></tr>\n",
       "</table><style>table.dataTable {clear: both; width: auto !important; margin: 0 !important;}\n",
       ".dataTables_info, .dataTables_length, .dataTables_filter, .dataTables_paginate{\n",
       "display: inline-block; margin-right: 1em; }\n",
       ".paginate_button { margin-right: 5px; }\n",
       "</style>\n",
       "<script>\n",
       "\n",
       "var astropy_sort_num = function(a, b) {\n",
       "    var a_num = parseFloat(a);\n",
       "    var b_num = parseFloat(b);\n",
       "\n",
       "    if (isNaN(a_num) && isNaN(b_num))\n",
       "        return ((a < b) ? -1 : ((a > b) ? 1 : 0));\n",
       "    else if (!isNaN(a_num) && !isNaN(b_num))\n",
       "        return ((a_num < b_num) ? -1 : ((a_num > b_num) ? 1 : 0));\n",
       "    else\n",
       "        return isNaN(a_num) ? -1 : 1;\n",
       "}\n",
       "\n",
       "require.config({paths: {\n",
       "    datatables: 'https://cdn.datatables.net/1.10.12/js/jquery.dataTables.min'\n",
       "}});\n",
       "require([\"datatables\"], function(){\n",
       "    console.log(\"$('#table140074613484064-312113').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",
       "    $('#table140074613484064-312113').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1, 2], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140074613766912-127783\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>idx</th><th>hp_idx_O_13</th><th>hp_idx_O_10</th><th>ferr_ap_ukidss_j_mean</th><th>f_ap_ukidss_j_p90</th><th>ferr_ukidss_j_mean</th><th>f_ukidss_j_p90</th><th>ferr_ap_90prime_g_mean</th><th>f_ap_90prime_g_p90</th><th>ferr_90prime_g_mean</th><th>f_90prime_g_p90</th><th>ferr_ap_90prime_r_mean</th><th>f_ap_90prime_r_p90</th><th>ferr_90prime_r_mean</th><th>f_90prime_r_p90</th><th>ferr_ap_mosaic_z_mean</th><th>f_ap_mosaic_z_p90</th><th>ferr_mosaic_z_mean</th><th>f_mosaic_z_p90</th></tr></thead>\n",
       "<thead><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th>uJy</th><th>uJy</th><th></th><th></th><th>uJy</th><th>uJy</th><th></th><th></th><th>uJy</th><th>uJy</th></tr></thead>\n",
       "<tr><td>0</td><td>177655159</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>12.153136</td><td>597.1377803802492</td><td>0.21581419</td><td>2.9271871089935306</td><td>0.2552266</td><td>2.4469470024108886</td><td>0.3040327</td><td>7.038074016571047</td><td>0.3783155</td><td>4.802640533447265</td><td>0.8010596</td><td>21.58711814880371</td><td>0.47714946</td><td>21.053738117218</td></tr>\n",
       "<tr><td>1</td><td>177655134</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>12.153136</td><td>597.1377803802492</td><td>0.21581419</td><td>2.9271871089935306</td><td>0.2552266</td><td>2.4469470024108886</td><td>0.3040327</td><td>7.038074016571047</td><td>0.3783155</td><td>4.802640533447265</td><td>0.8010596</td><td>21.58711814880371</td><td>0.47714946</td><td>21.053738117218</td></tr>\n",
       "<tr><td>2</td><td>177655135</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>12.153136</td><td>597.1377803802492</td><td>0.21581419</td><td>2.9271871089935306</td><td>0.2552266</td><td>2.4469470024108886</td><td>0.3040327</td><td>7.038074016571047</td><td>0.3783155</td><td>4.802640533447265</td><td>0.8010596</td><td>21.58711814880371</td><td>0.47714946</td><td>21.053738117218</td></tr>\n",
       "<tr><td>3</td><td>177655165</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>12.153136</td><td>597.1377803802492</td><td>0.21581419</td><td>2.9271871089935306</td><td>0.2552266</td><td>2.4469470024108886</td><td>0.3040327</td><td>7.038074016571047</td><td>0.3783155</td><td>4.802640533447265</td><td>0.8010596</td><td>21.58711814880371</td><td>0.47714946</td><td>21.053738117218</td></tr>\n",
       "<tr><td>4</td><td>177655164</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>12.153136</td><td>597.1377803802492</td><td>0.21581419</td><td>2.9271871089935306</td><td>0.2552266</td><td>2.4469470024108886</td><td>0.3040327</td><td>7.038074016571047</td><td>0.3783155</td><td>4.802640533447265</td><td>0.8010596</td><td>21.58711814880371</td><td>0.47714946</td><td>21.053738117218</td></tr>\n",
       "<tr><td>5</td><td>177655158</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>12.153136</td><td>597.1377803802492</td><td>0.21581419</td><td>2.9271871089935306</td><td>0.2552266</td><td>2.4469470024108886</td><td>0.3040327</td><td>7.038074016571047</td><td>0.3783155</td><td>4.802640533447265</td><td>0.8010596</td><td>21.58711814880371</td><td>0.47714946</td><td>21.053738117218</td></tr>\n",
       "<tr><td>6</td><td>177655133</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>12.153136</td><td>597.1377803802492</td><td>0.21581419</td><td>2.9271871089935306</td><td>0.2552266</td><td>2.4469470024108886</td><td>0.3040327</td><td>7.038074016571047</td><td>0.3783155</td><td>4.802640533447265</td><td>0.8010596</td><td>21.58711814880371</td><td>0.47714946</td><td>21.053738117218</td></tr>\n",
       "<tr><td>7</td><td>177655166</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>12.153136</td><td>597.1377803802492</td><td>0.21581419</td><td>2.9271871089935306</td><td>0.2552266</td><td>2.4469470024108886</td><td>0.3040327</td><td>7.038074016571047</td><td>0.3783155</td><td>4.802640533447265</td><td>0.8010596</td><td>21.58711814880371</td><td>0.47714946</td><td>21.053738117218</td></tr>\n",
       "<tr><td>8</td><td>177655132</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>12.153136</td><td>597.1377803802492</td><td>0.21581419</td><td>2.9271871089935306</td><td>0.2552266</td><td>2.4469470024108886</td><td>0.3040327</td><td>7.038074016571047</td><td>0.3783155</td><td>4.802640533447265</td><td>0.8010596</td><td>21.58711814880371</td><td>0.47714946</td><td>21.053738117218</td></tr>\n",
       "<tr><td>9</td><td>177655167</td><td>2775861</td><td>6.322884</td><td>607.8445449829104</td><td>12.153136</td><td>597.1377803802492</td><td>0.21581419</td><td>2.9271871089935306</td><td>0.2552266</td><td>2.4469470024108886</td><td>0.3040327</td><td>7.038074016571047</td><td>0.3783155</td><td>4.802640533447265</td><td>0.8010596</td><td>21.58711814880371</td><td>0.47714946</td><td>21.053738117218</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(\"$('#table140074613766912-127783').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",
       "    $('#table140074613766912-127783').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for col in master_catalogue.colnames:\n",
    "    if col.startswith(\"f_\"):\n",
    "        errcol = \"ferr{}\".format(col[1:])\n",
    "        depths = join(depths, \n",
    "                      group[\"hp_idx_O_{}\".format(str(ORDER)), errcol].groups.aggregate(np.nanmean),\n",
    "                     join_type='left')\n",
    "        depths[errcol].name = errcol + \"_mean\"\n",
    "        depths = join(depths, \n",
    "                      group[\"hp_idx_O_{}\".format(str(ORDER)), col].groups.aggregate(lambda x: np.nanpercentile(x, 90.)),\n",
    "                     join_type='left')\n",
    "        depths[col].name = col + \"_p90\"\n",
    "\n",
    "depths[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III - Save the depth map table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "depths.write(\"{}/depths_{}_{}.fits\".format(OUT_DIR, FIELD.lower(), SUFFIX))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IV - Overview plots\n",
    "\n",
    "### IV.a - Filters\n",
    "First we simply plot all the filters available on this field to give an overview of coverage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'90prime_g', '90prime_r', 'mosaic_z', 'ukidss_j'}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tot_bands = [column[2:] for column in master_catalogue.colnames \n",
    "             if (column.startswith('f_') & ~column.startswith('f_ap_'))]\n",
    "ap_bands = [column[5:] for column in master_catalogue.colnames \n",
    "            if column.startswith('f_ap_') ]\n",
    "bands = set(tot_bands) | set(ap_bands)\n",
    "bands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Passbands on SA13')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAEgCAYAAACU3FvWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl8VNX5+PHPM1t2tpCwE3ZlX0RQpJWiRbAWl+Jeq1YrbaVSQa3falFB/Vat2PqtWqhbtVq3Vn8UFJdaEUVklTUgiwRICIQA2ZdZzu+PmYkhJJMJZDJ3Ms/79eI1y71z7zMJk2eec849R4wxKKWUUsq6bNEOQCmllFKhabJWSimlLE6TtVJKKWVxmqyVUkopi9NkrZRSSlmcJmullFLK4jRZq1ZHRO4Xkb+3tnMppeKXJmvVIkRkj4hUiEipiBwUkRdEJDXacbUGInKTiGwTkZLAz3aJiKTV2ed+ETEiMqbO811EZJGI5AW296qz/VER2ScixSKSIyL3RP4dKaXq0mStWtIPjTGpwCjgTODeKMcT80TkXOBh4GpjTBowEHijzj4CXAccAa6vcwgfsBT4UQOneA443RjTBhgHXCMilzXfO1BKhUOTtWpxxphc4D1gCICI3Cgi2YHKcLeITA/uKyIdRWSxiBwTkSMislxEbIFtvxGR3MDrtovIebVOkygirwe2rROR4bWOebeI7Aps2yoil9badoOIfCYifxCRoyLyjYhMqbW9t4gsC7z2Q6BjrW2JIvJ3ESkMxLtaRDrV9zMQkYEi8klgvy0iMrXWthdF5KlAhVwiIl+KSN8GfpxnAl8YY9YHfrZHjDF/M8aU1NrnO0BXYCZwlYi4av0uDhpjngZW13dwY8x2Y0xZrad8QL8GYlFKRYgma9XiRKQHcCGwPvDUIeAioA1wI/CEiIwKbJsN7AcygE7AbwEjIqcBM4AzAxXlBcCeWqe5GHgT6AC8CrwjIs7Atl34E1hb4AHg7yLSpdZrxwLb8SfiR4HnAtUpgWOtDWybx/GV6vWBY/YA0oGfAxX1vH8n8G/gAyAT+BXwSuA9BV0diK09sBN4qO5xAr4ELhCRB0TkHBFJqGef6wPnez3w+KIGjlWvwJebUvy/hxT8PwOlVAvSZK1a0jsicgz4DFiGv/kWY8wSY8wu47cMfxL7TuA1bqALkGWMcRtjlhv/hPZeIAEYJCJOY8weY8yuWudaa4x5yxjjBuYDicBZgfO9aYzJM8b4jDGvAzuA2n25OcaYvxpjvMDfAufvJCI98VeyvzPGVBljPsWfBIPc+JN0P2OM1xiz1hhTXM/P4SwgFfi9MabaGPMxsBh/gg76lzFmlTHGA7wCjKjvB2qMWQ5chr9rYQlQKCLzRcQOICLJwOXAq4GfxVuc2BQekjHm90Ba4BwvA0VNeb1S6tRpslYt6RJjTDtjTJYx5pfGmAoAEZkiIisDzdzH8Ffdweblx/BXlh8EmsjvBjDG7AR+DdwPHBKR10Ska61z7QveMcb48FeFXQPn+4mIfBVogj6Gvzm+Y63X5td6bXngbmrg9UfrNAvn1Lr/MvA+8FpgwNajtar52roC+wJx1T5Ot/piAMoD56+XMeY9Y8wP8bciXAzcANwc2Hwp4AHeDTx+BZgiIhkNHa+Bc5hAU3sF/opfKdWCNFmrqAo02/4T+APQyRjTDn9iEQBjTIkxZrYxpg/wQ2BWsG/aGPOqMWY8kAUY4JFah+5R6xw2oDuQJyJZwF/xN6GnB863OXi+RhwA2otISq3negbvBCr/B4wxg/APxroI+Ek9x8kDegT73msdJzeMGBoUaCn4D/AxgfEA+KvoVGCviOTj7xpwcnwV3xQOoKH+c6VUhGiyVtHmwt+cXQB4AoO5JgU3ishFItIv0GdcjL/52ysip4nIxECyr8Rf8XlrHfcMEblMRBz4K/AqYCX+PlcTOB8iciPfJraQjDE5wBrgARFxich4/F8ggrF+T0SGBpqgi/E3i3vrOdSXQBlwl4g4RWRC4DivhRNHbSJysYhcJSLtxW8McC6wUkS6Aefh/9IwIvBvOP4vNdfXOkYi/t8BQELgMSJiE5HpdY59K/CfpsaplDo1mqxVVAVGLd+G/3Kjo8A1wKJau/QHPgJKgS+Ap40xn+BPLr8HDuNvMs7EP/gs6P8BVwaOeR1wWaDy3Qo8HjjWQWAo8HkTQr4G/wC0I8B9wEu1tnXG3ydcDGTj75c/YcIUY0w1MBWYEoj/aeAnxphtTYgj6CjwM/z97sWB8z1mjHkF//v+yhjzgTEmP/gPeBIYJiLBLykV+H++ANs4flDcpfgH5JUEjv1/gX9KqRYk/rE6SimllLIqrayVUkopi9NkrZRSSlmcJmullFLK4jRZK6WUUhanyVoppZSyOEe0A2iqjh07ml69ekU7DKWUiilr1649bIxp0sx19Rwj0+FwPIt/bgIt9pqXD9js8XhuPuOMMw7V3RhzybpXr16sWbMm2mEopVRMEZGcxvcKzeFwPNu5c+eBGRkZR202m17324x8Pp8UFBQMys/Pfxb/PAzH0W9GSimlwjUkIyOjWBN187PZbCYjI6OIBmZU1GStlFIqXDZN1JET+NnWm5cjlqxF5HkROSQimxvYLiLypIjsFJGNtdYvVkoppVQtkaysXwQmh9g+Bf+8z/2BW4BnIhiLUkqpVmDevHmZ/fv3H9yvX7/Bc+fOzQQ4ePCgfdy4cf2zsrKGjBs3rn9BQYG9qccdOXLk6c0fbfOJWLI2xnyKf7GDhlwMvBRYJ3cl0E5EukQqHqWUUrFt9erViS+99FLGunXrsrOzs7csXbq03aZNmxLuu+++LhMmTCjJycnZPGHChJI5c+Z0DveYHo8HgPXr15/MQjotJpqjwbsB+2o93h947kB0wlFKKRWuO9/a0OPr/JLk5jzmgM5p5Y9NG76voe2bNm1KGjVqVGlaWpoP4Jxzzil5/fXX2y1durTdsmXLtgNMnz698Nxzzz0NyJ01a1bX3bt3J+Tn5zsPHDjguu222/Jnz559ePHixWnz5s3rkpmZ6d66dWvyrl27tiQnJ48sLy9fv3jx4rQHHniga0ZGhnvr1q3JF1544dGhQ4dWPP30052qqqrk7bff3jV48OCqvLw8x4033piVm5vrApg/f/7eSZMmldUXd15enmPatGm9jx075hgxYkT5J5980mbt2rXZXbp08YT7s4nmADOp57l6By6IyC0iskZE1hQUFEQ4LKXinM8LRblwZDcc3ArrX4GS/GhHpRQjRoyo+PLLL9Py8/PtJSUltg8//LDtvn37XIWFhY6srCw3QFZWlvvIkSM1hWh2dnbSRx99tGPlypXbHnvssa579uxxAmzcuDHlsccey921a9eWuufZtm1b0jPPPLMvOzt7y1tvvZX+9ddfJ27atCn7uuuuO/z4449nAkyfPr3HrFmzDm7evDn77bff3vXzn/+8V0Nx33333V3PPffckq1bt2ZfdtllRw8cOOBq6nuPZmW9H+hR63F3IK++HY0xC4GFAKNHj9aRiEo1B58XPFXgChRH1WXw4X3w1avgrqdAGHMLnH8/uFLqP1bRPkhsC0ntIxm1sohQFXCkjBo1qnLmzJn5EydOHJCcnOwbNGhQucMROo1NmTLlWGpqqklNTfWcffbZxcuXL09p3769d9iwYWWnn356dX2vGTp0aFkw+ffs2bNqypQpRQDDhw+vWLZsWRrA559/3mbHjh1JwdeUlpbajx49amvfvr2v7vFWrVqV+s477+wEmDZtWnGbNm28TX3v0UzWi4AZIvIaMBYoMsZoE7hSLeHYPlh4rj9Bj50O7XvBp49D8X5I6wqdh4CnEnp/F/auhP2rYdVC/787d0FKR/9xctfC4tshfzOYwN+fLsPhB09A9zP8jz1VcHQPiA069o/Gu1WtyO2333749ttvPwwwY8aMbt27d69OT0/35OTkOLOystw5OTnODh061DQvixzfiBt8nJycfEJSDUpISKgpCm02G4mJiSZ43+v1CoAxhjVr1mSnpqY2WkAac+o1ZsSStYj8A5gAdBSR/cB9gBPAGPMX4F3gQmAnUA7cGKlYlFJ1rHsJyguh0xD4/E/+59r3ghvehV7n1P+aLxfAe3fBY31hyDTY/QmUH/ZvO+0HkJIOlcXwzTJ4diJ0HeWvwnPXgrvcv99374KJ90T63alWLDc319GtWzfPjh07XEuWLGm3atWqbd98803CggUL0h9++OH8BQsWpE+ePPlYcP/33nuv3UMPPXSguLjYtnLlyrQnnngid/PmzYmnGsf48eOLH3nkkcx58+YdBFixYkXSuHHjKurbd8yYMaUvv/xyh4ceeij/X//6V5vi4uImj1aPWLI2xlzdyHYD3Bqp8yulQvh6KfQcBz99z18VVx6DHmeBPcSfhLHTIb0vvHI5bH7L/9wZN8D4WdA+69v9qkpg+eOw7V0ozoUhl0Gnof5zfvoo9BgL/c+P6NtTrdfUqVP7Hjt2zOFwOMwf//jHvRkZGd4HHnjgwKWXXto3KyurY9euXavfeeedXcH9R44cWXbeeef1z8vLc91xxx0HevXq5W6OZL1w4cJ9N998c88BAwYM8nq9Mnbs2JJx48btrW/f3//+93nTpk3rM2jQoPZnn312aUZGhrtdu3ZNagqX5ijPW9Lo0aONzg2u1CnwuuGhzjDuNjj/vpM/BhI6udflqYanx4IzBX6+HKS+MaYqUkRkrTFm9KkcY8OGDXuGDx9+uLliirRZs2Z1TU1N9c6dO/dgNOOoqKgQh8NhnE4nH330UcqMGTOytm3btrW+fTds2NBx+PDhveo+H3MLeSilTtGR3eDzQMZpJ38Mu7Ppr3G44Kxfwrt3wMHN0HnoyZ9fqRiyc+dO1xVXXNHX5/PhdDrNggUL9jT1GJqslYo3Rfv9t+2yQu8XCYMu8fd7r38Fpvy+5c+v4sr8+fPrvcIoUv70pz+lP/PMM51qP3fmmWeWvvzyy3uzs7PrraTDpclaqXhTGmgRTOsUer9ISM2APhNg2xKY/L/aFK5alZkzZxbOnDmzMBLH1lW3lIo3wWSdkhmd8w+ZBkV7Ye8X4K6MTgxKxRhN1krFm5KD4EqFhNTonH/QxZDQBl6YAvMHQsWxxl+jVJzTZK1UvCk7BKlRqqrB/yXhsoXQoQ9UHIGcFdGLRakYoclaqXhTWQSJ7aIbw2lT4KaP/PeP1XtpqlKqFk3WSsWbqhJISIt2FJDcAWxOKNFZhpX1WG19a03WSsUbqyRrEf8c4+UxM8eGiiNWW99aL91SKt5UlfgHeFmBK9W/mIiKPe/c2oNDW5t1PWsyB5VzyVMNrua1fft21+TJk/uPGTOmdN26dakDBw4s/+lPf3p47ty53QoLCx0vvvji7kGDBlVde+21vfbu3ZuQlJTkW7hwYc7YsWMrlixZkjp79uye4F/MY8WKFdtsNhuTJ0/uV1RUZPd4PDJnzpy8H//4x8cAgutbA9x7772d3njjjXQR4bzzzit6+umnc+vGtmfPHufkyZNrVqrZsWNHUnZ29qYBAwbUu7JXU2myVireVBVDokWSdUIqVJVGOwoVQ/bt25f4+uuv7z7jjDNyhg0bNvCVV15JX7NmzbZXX3213UMPPdSlW7du1cOHDy//6KOPdi1atCjt+uuv771t27atjz/+eOcnn3wyZ9KkSWVFRUW24KpbS5Ys2dmhQwffgQMHHGPHjj39mmuuOWazfdvo/MYbb7RZsmRJ+7Vr125LS0vzHTx4sN5FOHr16uUOTiH6v//7vxnLly9Pa65EDZqslYovxlinGRy0so5lISrgSOrWrVvVmDFjKgAGDBhQMXHixGKbzcaoUaPKH3zwwa65ubkJ//znP3cCTJ06teSWW25xFBYW2s8666zSO+64o8cVV1xx5Oqrrz7at29fX1VVlfz617/uvnLlylSbzcahQ4dc+/fvd/Ts2bNmic0PP/ywzY9//OPDaWlpPoBOnTqFXIDjgw8+SHnppZcyVq5c2azN6NpnrVQ8cZeD8VkoWadAdUm0o1AxxOVy1bvWtN1ux+v11rs4lYiYhx9+OP/ZZ5/NqaiosI0bN27g+vXrExcsWNChsLDQsWnTpuxt27ZtTU9Pd1dUVByXF40xJ6yJ3ZCcnBzn9OnTe73++uu72rZt2+B62SdDk7VS8SR4mZRlkrVW1qp5nXXWWSUvvPBCOsDixYvT2rdv7+nQoYNvy5YtCWPGjKl46KGH8ocOHVq2efPmxKKiInvHjh3dCQkJ5t///ndaXl6eq+7xJk+eXPzyyy93LCkpsQE01AxeVVUll112WZ958+blDhs2rKq535c2gysVT56f7L91pkQ3jiBXMlSXRzsK1Yo88sgjeddcc02vAQMGDEpKSvK9+OKL3wA8+uijmStWrGhjs9nMgAEDKqZNm1Z07Ngx+5QpU/oNGTJk4ODBg8t79+59wvy306ZNK163bl3yiBEjBjqdTnP++ecX/fnPfz5hgNlHH32Usnnz5pQHH3yw64MPPtgVYOnSpTt69erlbo73petZKxVP7m/rv/3RczB0WnRjAVhyB2x+C36zJ9qRtHrxuJ51LGpoPWttBlcqHrmiNC94XY4E8DTbgFmlWi1tBlcqHjmToh2Bn90F3mbv3lMqoq677rqeq1evPu4b7y9+8YuDkVoeEzRZKxWfrJKsHQng84DPBzZt6FOx4eWXX27xCe3106FUPErvF+0I/OyBwbdaXSsVkiZrpeLJ0Mv9S1Mmd4h2JH6ORP+tR5O1UqFoslYqnnir/StdWYUjWFnrIDOlQtFkrVQ88bq/bXq2AnuC/1Yra6VC0mStVDzxusFupco6kKy1slZhmjdvXmb//v0H9+vXb/DcuXMzwT+r2Lhx4/pnZWUNGTduXP+CgoJ6ZxkLxWrrV9elyVqpeOKttlhlHYhFK2sVhtWrVye+9NJLGevWrcvOzs7esnTp0nabNm1KuO+++7pMmDChJCcnZ/OECRNK5syZ0zncY3o8/jU7Wmr9arf75CY000u3lIonVqusg8na1ywzMqoW9LvPf9dj59Gdzbqedb/2/crnnTOvwdW8Nm3alDRq1KjS4ApY55xzTsnrr7/ebunSpe2WLVu2HWD69OmF55577mlA7qxZs7ru3r07IT8/33ngwAHXbbfdlj979uzDixcvTps3b16XzMxM99atW5N37dq1Jbh+9eLFi9MeeOCBrhkZGe6tW7cmX3jhhUeHDh1a8fTTT3eqqqqSt99+e9fgwYOr8vLyHDfeeGNWbm6uC2D+/Pl7J02aVO9E97Nmzep64MAB5969e10dOnTw/Pvf//6mqT8bTdZKxRNvNTgtspY1gC3wJ8jrCb2fUsCIESMq5s6d2y0/P9+ekpJiPvzww7bDhw8vKywsdGRlZbkBsrKy3EeOHKnJbdnZ2Ulr167NLikpsY8cOXLQj370oyKAjRs3pqxfv37L6aeffkIfzLZt25Leeuut3ZmZmZ6srKyhCQkJhzdt2pQ9b968zMcffzzz+eef3zd9+vQes2bNOnjBBReU7tixw3XBBRf0371795aGYt+4cWPyl19+uS01NfWk5vjWZK1UPPFZbICZLdC16NNkHWtCVcCRMmrUqMqZM2fmT5w4cUBycrJv0KBB5Q5H6DQ2ZcqUY6mpqSY1NdVz9tlnFy9fvjylffv23mHDhpXVl6gBhg4dWhZM/j179qyaMmVKEcDw4cMrli1blgbw+eeft9mxY0fN7EKlpaX2o0eP2tq3b1/v0piTJ08+drKJGjRZKxVfLNcMHohFk7UK0+2333749ttvPwwwY8aMbt27d69OT0/35OTkOLOystw5OTnODh061PyHqrsWdfBxcnJyg+tNJyQk1Ltmts1mw+v1CvjXuV6zZk12uAk4JSXllNa31gFmSsUTq11nHWwG1z5rFabc3FwHwI4dO1xLlixpd9NNNx254IILji1YsCAdYMGCBemTJ08+Ftz/vffea1deXi75+fn2lStXpo0fP75ZFlAfP3588SOPPJIZfLxixYqIzuEb0cpaRCYDfwLswLPGmN/X2d4T+BvQLrDP3caYdyMZk1JxzWqVdU2y1spahWfq1Kl9jx075nA4HOaPf/zj3oyMDO8DDzxw4NJLL+2blZXVsWvXrtXvvPPOruD+I0eOLDvvvPP65+Xlue64444DvXr1cm/evDnxVONYuHDhvptvvrnngAEDBnm9Xhk7dmzJuHHjIjZneMTWsxYRO/A18H1gP7AauNoYs7XWPguB9caYZ0RkEPCuMaZXqOPqetZKnYL5g6HPBLjkqWhH4pe7Dv76PbjmDRhwQbSjadXicT3rWbNmdU1NTfXOnTv3YLRjCVc01rMeA+w0xuw2xlQDrwEX19nHAMGhqW2BvAjGo5Tyub8d1GUFNaPBtRlcqVAi2QzeDag9WnA/MLbOPvcDH4jIr4AU4PwIxqOU8nm0GVzFjfnz57doAfinP/0p/ZlnnulU+7kzzzyztDmW1IxkspZ6nqvb5n418KIx5nERORt4WUSGGGOOGzUnIrcAtwD07NkzIsEqFRe8nm8TpBVoslatyMyZMwtnzpxZGIljR7IZfD/Qo9bj7pzYzH0T8AaAMeYLIBHoWPdAxpiFxpjRxpjRGRkZEQpXqTjgs1iytmuyViockUzWq4H+ItJbRFzAVcCiOvvsBc4DEJGB+JN1QQRjUiq++dzWStZaWSsVlogla2OMB5gBvA9kA28YY7aIyFwRmRrYbTbwMxHZAPwDuMFEani6Ukov3VIqRkX0K3bgmul36zw3p9b9rcA5kYxBKRXg8wHGYpV14IuDjgZXKiSdwUypeBGcJcxSyTo4N7g3unGomDZmzJjTPv300xNWADv33HP7HT58+IRrFWfNmtV1zpw5neo+31R79uxxTp48uc+pHiccFvrUKqUiKtjUbKlkrc3gKnKWLVu2M5LH79Wrl3vp0qW7I3mOIAt9apVSERVsarZkn7U2g8eavN/e06Nqx45mXc86oX//8q4PP9Tgal7bt293XXTRRf137NixBWDOnDmdSktLaypnr9fL5Zdf3qt79+7VTz75ZF63bt2GrlmzJrtLly6e3/zmN51ff/31jl27dq1OT093jxw5shzgwQcfzHzhhRcy7Ha7GTBgQOXixYt3L1myJHX27Nk9wb/wx4oVK7bVt5pW3XgiSZO1UvEi2NRspcpaV91SzcTtdssll1zSe9CgQRWPPPJIfu1ty5cvT3777bc7bNq0aavb7WbEiBGDgsn6ySef7JyTk7MpKSnJBJvMH3/88c5PPvlkzqRJk8qKiopsoVboaikW+tQqpSLKin3Won3WsSpUBRwNv/zlL7MuueSSI3UTNcB///vf1AsvvPBYWlqaD2DSpEk1q3KddtppFZdeemnvqVOnHrv22muPAZx11lmld9xxR48rrrjiyNVXX320b9++UU/WOsBMqXgRrF4t1QxuA7HpaHAVFofDYXy+b/NmZWVlTQ4bPXp06fLly9uUl5fXN3vmCetaB/33v//dceuttxasXbs2Zfjw4YPcbjcPP/xw/rPPPptTUVFhGzdu3MD169ef8ipdp0qTtVLxwmvByhr8l29pM7gKQ/fu3T1Hjhxx5Ofn2ysqKuT9999vG9w2ffr0w5MmTSq66KKL+rrdx3/5mzhxYumSJUvalZaWytGjR20ffvhhO/D3ce/atcv1wx/+sOTpp5/eX1JSYi8qKrJv2bIlYcyYMRUPPfRQ/tChQ8uaY0nNU2WxT61SKmJq+qwtVFmD/8uDJmsVhoSEBDN79uwDY8aMGdi9e/eqfv36Vdbefv/99x+8/fbb7Zdddlnvd95555vg8+PHjy+/9NJLjwwZMmRwt27dqsaMGVMK4PF45JprruldUlJiN8bI9OnTD3bs2NE7e/bsritWrGhjs9nMgAEDKqZNm1bUUEwi0iITeWmyVipe1PRZW2iJTAgka+2zVuG59957D917772HGtr+xBNP1KxBkZubuyl4/5FHHsmvrz977dq12+s+97e//S2s/vhDhw452rZt2yL/ebUZXKl4YcU+a/Av5qGXbqkY8+mnnyZfd911fWbMmHGwJc6nlbVS8cKyfdbaDK6sbdWqVUk/+clPetd+zuVy+fbs2bO5pWKw2KdWKRUx2metTp3P5/OJzWaLqwWXxowZU7Ft27atkT6Pz+cToN7LxLQZXKl4Ydk+azt4NVnHiM0FBQVtA0lFNSOfzycFBQVtgXqrda2slYoXVu2ztjnA6ACzWODxeG7Oz89/Nj8/fwha7DU3H7DZ4/HcXN9GTdZKxQvts1an6IwzzjgETI12HPFIvxkpFS+s2mctdk3WSjVCk7VS8cKyfdZ6nbVSjdFkrVS8sGyftV2TtVKN0GStVLzQPmulYpYma6XihRXXswZN1kqFQZO1UvHCiutZgzaDKxUGTdZKxYtgM7gV+6z1OmulQtJkrVS8CDY1W66y1mZwpRqjyVqpeKHJWqmYFdanVkQ6Az1r72+MWRGpoJRSEWDVZK2ToijVqEY/tSLyMPBjYBsQ7FgywIURjEsp1dys3Gftq3ehIaVUQDhfsX8EDDDGVEY6GKVUBOmlW0rFrHD6rL8Jcz+llJVZ9tItTdZKNSacT20JsF5EPgKqgk8aY2ZFLCqlVPPzefyJUSy2FLFN+6yVakw4yXpp4J9SKpZ53darqkEX8lAqDI1+co0xz4mIA+gXeGqnMUa/BisVa3xe6y2PCTopilJhaLQvWkS+A+wEngOeB74WkXPCObiITBaR7SKyU0TubmCfK0Rkq4hsEZFXmxK8UqoJfG7rLY8J2metVBjCaRN7ArjQGLMVQEQGAi8Do0O9SETswFPA94H9wGoRWRQ8TmCf/sD/AOcYY46KSObJvQ2lVKN8HutdtgV6nbVSYQhnlLerdoI1xmQDrjBeNwZ/k/luY0w18BpwcZ19fgY8ZYw5Gjj2ofDCVko1mfZZKxWzwknW60RkgYiMD/x7Blgfxuu6AftqPd4feK62AcAAEflcRFaKyOT6DiQit4jIGhFZU1BQEMaplVInsHKftSZrpUIKJ1n/HNgF3AX8BtgNTA/jdfVdH2LqPHYA/YEJwNVg/xv0AAAgAElEQVTAsyLS7oQXGbPQGDPaGDM6IyMjjFMrpU6gfdZKxaxwRoNXAo8G/jXFfqBHrcfdgbx69llpjHED34jIdvzJe3UTz6WUaoxV+6w1WSvVqAYraxH5R+B2vYisq/svjGOvBvqLSG8RcQFXAYvq7PMO8L3AeTribxbffTJvRCnVCMv2WQcu3TJ1G96UUkGhPrl3Bm6nncyBjTEeEZkBvA/YgeeNMVtEZC6wxhizKLBtkohsxb9IyJ3GmMKTOZ9SqhE+r0WTdSAm4/OPDFdKnaDBT64xZn/gbh5QaYwxItIXOA34IJyDG2PeBd6t89ycWvcNMCvwTykVST4LV9YQmA5Vk7VS9QlngNlyIElEugDLgF/gnxxFKRVLrNxnDdpvrVQI4SRrmzGmHP9SmX82xvwQGBbZsJRSzc7rsWZlLbUqa6VUvcJK1iJyJnANsDjwnLZVKRVrfBZN1jWVtV5rrVRDwknWs4AHgCXGmM0i0gd/07hSKpZYvs9ak7VSDQnnOuuPgY8BRESAg8aYX0Y6MKVUM/O6tc9aqRgVzqpbL4lIGxFJBrbgn7xER28rFWsse+mW9lkr1ZhwmsGHGmOKgUvwX7LVHbghkkEppSLAss3gWlkr1ZiwVt0SEQf+FbPeCayg5YtsWEqpZmf1AWZG/6wo1ZBwPrnPAnuBzcAyEekJlEY0KqWiYF/JPj7P/ZwEewJTek8h0ZF43PZVB1bhsrsYkTkiShGeIq9Vr7PWZnClGhPOALMngCeCj0VkHzAxkkEp1RI+3vsxy/Yv42jlUdw+N+sOrqPcUw7A7qLdzB49G4D3vnmPbqnduOmDmwD4/OrPaeNqE7W4T5pVZwjTZnClGtVgshaRq40x/xCR2xrY5ckIxaRUxKzIXUGHpA6UVJfw6//+GqfNSWZyJinOFCb2nMh1g65j5n9nklfqXyCuuLqYuz6967hjPL/peX59xq+jEf6p8bmtuZ61ToqiVKNCVdbtA7e6gLSKecYYXtr6En9Y84ea57LaZPHKha/QNqHtcfu2dbWl2lcNwPt73q95/vye5yMivPH1G/xixC9IsCe0TPDNxfLTjep11ko1JNRCHk8Hbn/XcuEo1by2FG7h6yNfsyJvBUv3LGVk5kiGZwzHLnZuHHLjCYkawGV34fa5eX/P+8z9Yi4DOwzknrPuYWjHoazMW8mHOR/yRd4XTOgxoeXf0Kmw6nSjmqyValSjn9zAgLIZQK/a+xtjLotcWEqdvOzCbErdpQxOH8xVi68CwCY2ZoyYwc+G/QybhL4Iwmlz4va6eeqrpwC47+z7GNxxMAB92/UFoKCiIILvIEIsOxpcm8GVakw4n9xFwEvAh+glW8rCytxlPLH2CV7f/voJ2+4/+34u7X9pWMdx2p3kleaxr2Qfd46+syZRA6Q4U/znqi5rnqBbkmWvs9ZkrVRjwvnkVhtj5kc8EqVOQZW3iukfTmdDwQauGHAFPdv05G9b/kZWmyxmjJzBqMxRYR/LaXNyoPQAAF1Tux63LdmZDECZJ8aStTHW77M22gyuVEPCSdb/JyL3Au8DVcEnjTEbIxaVUk1Q7i7nrk/vYkPBBh777mNM7j0ZgOsGXYcxBnsTL1dy2Vx4jL/KS09KP26bTWwkO5Ipc8dYsg72B1uystZLt5RqTDif3AHAzcAUvm0GN8B3IxWUUo0pd5dzuOIw3VK7cdend7E8dzn3jr23JlGDP7EiTT+2s1b1mZ6YfsL2VGcqa/LXsP3Idk7rcNpJxd/igonQkslaV91SqjHhfHKvAHoZY6oa3VOpCNpXso83tr/BhoINbCzYiLdWs+lvx/6WK0+/slnO47K5au53TOp4wvZkZzLZR7KZ9u9pbLp+U7OcM+J8bv+tFZO1XmetVKPCmRt8I5AW6UCUaojX5+UPq//AJe9cwivZr1DtrebK067EEUg8P+jzA6467apmO9/ekr0194N91LX5V4qNMcFEaOU+a03WSjUonK/Z6cA2EfmS4/us9dItFXH5Zfl8/63vAzCl9xRmnzGbTimdAFiRt4I9xXuY2mdqsybQxvqjS6pLmu1cLcZr5WZwvc5aqcaE88l9KOJRKFXHqgOreH/P+3yZ/yUA1w68lrvH3H3cPsMzhrOneA/92vdr1nMHk/E/fvCPerffMfoO7l5+Nx0SOzTreSOqps/aynODa7JWqiHhJOsVQKUxxohIX+A0/OtaK9XsjDH8ZeNfePqrp0l1pjKg/QBuHXErU3pPOWHfe866hytPu5LM5MxmjaHU7V9UrlNyp3q3/6DPD1h/aD0f7Imhj0FNn7UVm8EDvXHaDK5Ug8JJ1suB74pIW2AZsB64CvhJJANT8edA6QEeX/s47+95n6l9pzLn7Dkh599OciQxNGNos8cxJH0IX+Z/SbvEdiHPXeGpaPZzR4w3kKy1z1qpmBROsrYZY8pF5KfAn40xvxeRryIdmGq9ytxlvLn9Tdw+Nxf3u5hfffwrthZurdk+64xZ3DD4hqgN5Hrie0+wr2QfzhBVaJIjiUpvJV6ft8nXcUeFpS/d0klRlGpMWMlaRM4ErgFuCTwXA3+dlBWVu8u5Zsk17C7aDcCT679daTXJkcT8CfMZ3218tMIDIM2VxqD0QSH3SXIkAVDprSTFltISYZ2aWEjWWlkr1aBwPrmzgAeAJcaYzSLSB3/TuFJN9uymZ9ldtJunz3uarYVb+fNXfwZg0SWL6N22d5SjC18wWVd4KmrmC7e0mGgG18paqYY0mqyNMR8DH9d6vBv4ZSSDUq3T4t2LeXbTs0ztO5XvdP8OZ3U5iypvFeO6joupRA21krW7ApKiHEw4aiprCyZr0QFmSjUmnCUy++Gvrntx/BKZkyIXlmptPt3/Kfd+di9ndj6T353lXyLdaXdy26jbohzZyQlOllLuKY9yJGGqqayt3AyulbVSDQnnk/sW8Bzwd0A/TarJFu9ezJzP5zCg/QCenPgkiY7EaId0ymo3g8cE7bNWKqaF88n1GWP+L+KRqFZp/aH13PPZPZzR6Qz+cO4fYqN/Nwyxl6ytfJ21VtZKNSacucH/n4jcIiIZItIm+C+cg4vIZBHZLiI7ReTuEPtNExEjIqPDjlzFhGc3PUt6Yjp/nvjn2JrxqxExl6y9Vp4bXBfyUKox4VTWNwduf1frOQP0DPUiEbEDTwHfB/YDq0VkkTFma5390oDbgC/DDVrFhgpPBSvzVnLFaVfUuyBGLAsm65jps7ZyM7iIf5CZXmetVIPCGQ3e4ySPPQbYGRg9joi8BlwMbK2z3zzgUeCOkzyPsqiDZQep9lU3es1yLIq5ytrKS2SCPy6trJVqUFifXBE5HRgE1IwMMsa82sjLugH7aj3eD4ytc9yRQA9jzGIRaTBZi8gtBCZk6dkzZEGvLCS4elWqMzXKkTS/RLv/o1DtrY5yJGGy8nXWoMlaqUaEc+nWvcAk4HTgfeAC4DOgsWRd31yRptZxbcATwA2NxWCMWQgsBBg9erRpZHdlESVu/+pVqa7Wl6yDa2l7YiXBWLkZHEDsOsBMqRDCGWB2JfA94IAx5jpgOOFV5PuB2k3o3YG8Wo/TgCHAJyKyBzgLWKSDzFqPsurWW1kHk7U72LxsdT4LDzAD/yCzWPnio1QUhJOsK4wxXsATGAyWD/QJ43Wrgf4i0ltEXPhX6loU3GiMKTLGdDTG9DLG9AJWAlONMWua/C6UJQWXmtTK2gK8sdBnrZW1Ug0JJ1mvF5F2wPPAGmAVsK6xFxljPMAM/E3n2cAbxpgtIjJXRKaeQswqinYVlPLG6n2N70itZN0KK2u7+C83iplkbeXpRkH7rJVqRMiv2eJfo/B+Y8wx4CkReR9oY4xpNFkDGGPeBd6t89ycBvadEFbEqsUVllYB0CHFxV1vbWT93qNcPrp7o0tYlla33mQtIjhtzthpBrfydKMQaAbXylqphoT85BpjjIgsBs4IPN7ZIlEpS3B7fSz8dDfzP/waAb4/qBNrc44CUOn2keQKvVJqmbuMBHsCTqv2k54ih80Rg5W1lZN1jPwslYqCcD65q0RkVLjVtIpdPp9h2Y4C9h+t4FhZNX//MoeDxVVMHtyZnQWlvLc5v2bfCre30WRd4i5pNdOL1ie2krWFpxsF/5cInRRFqQY1mKxFxBHodx4P/ExEdgFl+C/JMsaYUS0Uo2oBxZVubnlpDSt3H6l5rm9GCnOuGcyFQztT5fHx4oo9bM8v4e31uZRXe+iQ4gp5zLLqslbZBB7ktDljJ1lbebpR0D5rpRoRqrJeBYwCLmmhWFSUbM4tYvYbG9hVUMpDlw7h7D7pdG2XRKLz28o50Wnn5+f25d8b8nh7fS4V1Y1XQaXu0lY5EjzIYXPETp91TWUdujUkajRZKxVSqGQtAMaYXS0Ui2phb6zZx4uf72HrgWLaJjl58cYxjO/fMeRrkgNN3+XhJmutrK3B6wZ7QrSjaJhOiqJUSKGSdYaIzGpoozFmfgTiURHk9Rn+tmIPz332DXabsPdIOQO7tOH+Hw7i0pHdaZvceBNpkrNpybpH6slOLW99MdVn7XWDPXS3RVTpaHClQgqVrO1AKvVPG6piyMHiSpZsPMA7X+WycX8Ro3q2w2m38aNR3fnl9/ritIdzub1fcFBZpbvxP6xl1WWtuxlcHHhMrCTrauv2V4M2gyvViFDJ+oAxZm6LRaKa3ab9Rbz0xR7e25xPaZWHXunJ/OmqEUwd3rXRa6Qbkuzy/5cJp7KOh9HgMdNn7a22eGWtyVqpUBrts1axpdrjw2B4+YscHl26nWqvj0mDOjF70mmc1jntlI8f7LNem3OUHwzrEnLfCndFq07WMTcpiqUra20GVyqUUMn6vBaLQp2yKo+XP3+8k78s24Xb61+Y7Dv9O/L4FcPJTEts5NXhCzaDP//5N8z5YcPrVPuMD4/x4LJZuJo7RbHVZ231ZnA7eGJkuVGloqDBZG2MOdLQNmUdB4oqKKn08OvXvmLrgWLOOz2TkT3b0aND8ik1dzckvZFrq4OC6zy31tnLIBaTtYW/ONkcYCqiHYVSlmXRuQdVY/YdKedP/9nBW2v3A5CW6ODP14zkwiFdsNki14MhItwwrhdvr88NuV+wedhp1RmzmoHD5qDCEyMJxuexeGWtfdZKhaLJOoYYY/hP9iFW7TnCv9blcqy8mhvP6UWXtol8f1Bnendsmf5hp11we30h9wlW1i4rV3OnSCvrZiQ6N7hSoWiyjgFen6Gows2DS7byr3W5uBw2hnVry7ybxjCwS5sWj8flsFHtCZ2sg5V1a+6zjq0BZhZP1jrATKmQNFlbnM9nuP75VXy28zB2m3Dr9/oy87wBuBzhXxvd3Fx2Ox6fweczDTa5uwNLMmqftUVYfjS4Q5O1UiFosrYgYwzLvi7gL8t2sWl/EWXVXvpnpvLElSMY0q1ttMPD6fAn6Gqvj8QG5pqu9gWawVtxZR1byboaXBa+jE77rJUKSZO1xVS6vdz2j/V8sPUgHVNd/GBYF0b1bM+VZ/Zo9pHdJ8sVmPGs2us7brGP2mpGg7fiAWaxNTe41ZvBNVkrFYomawvZnl/CXf/cyMb9x/jthadz/bheJDist0pSsAneHaLfumY0uJWbXk9RbPVZW70ZXPuslQpFk7UFbNh3jL+vzOHdTQdIcNr501UjmTq8a7TDalDtyrohOhrcYixfWdvBaLJWqiGarKPI7fVx36ItvPrlXgAuGNyJOT8cTLd2SVGOLLTgwh9uj2lwn3i5zjp2krUbrPy70GZwpULSZB0lhaVVzH5zA59sL+CK0d356fjenN655S/DOhnBZvBqb8OVUDxcuqWrbjUjTdZKhaTJuoUZY/jnulweXbqNY+Vu7v3BQG7+Tp9oh9Ukwcq6OkRlHS/N4LHTZ23xZnDRPmulQtFk3YIOlVTy239t4qPsQwzr3pb5V4xgfP+O0Q6ryRIcjfdZx0MzuNPuHw1ujLHMSP0Ged3WTtY6wEypkDRZtwCfz/DKlzn84YOvqXR7+d1Fg7hxXK+IzuEdSTV91mEMMGvNo8Ed4v/4eIwHp1j8fVp+NLg2gysViibrCPP6DNc99yUrdhUyrm86cy8eQr/M1GiHdUpq+qxDXLoVL5OiAHh8Hmu3IBhj/WZwTdZKhaTJOoLKqz388P8+Y1dBGb8+vz+3Tewfs9V0bU77tzOYNSRephsFrD8i3OcFjMWTdeDSLWPA6l0KSkWBJusIeuq/O9lVUMY9Fw7kpvG9W0WihvAq67gYDR4ryTrQJWH5ZnDwf7Gw658lperST0UEVLq9/PbtTfxrXS4XDevCz74bW6O9G1MzKUqoZvA4GA0ebPq2/IjwmmRt4d9FcI5540X/LCl1Iv1UNLOCkir+518b+Sj7EDeM68Wvz+8f7ZCaXc10o2GMBg9Wn61RMFlbv7IOfJmIicraAyRENRSlrKj1/iWNgs93HubWV9fF7PXT4QprgJm3GofNgU2it5RnpMVeM7iFK2sJVNZW/1kqFSWarJvJ5twirn9+Fd3aJ/Hc9aM5I6tDtEOKmLAu3fJVt+r+atBk3axq91nHAeP14i0uxtG+fbRDUTEiomWPiEwWke0islNE7q5n+ywR2SoiG0XkPyKSFcl4ImXnoVJufHE1HVJcvPPLc1p1ooZvK+uqUAPMvO5WPRIcYqnPOhaaweOnsj7y0stsGzyEHWeP4+vx38FXXR3tkFQMiFiyFhE78BQwBRgEXC0ig+rsth4YbYwZBrwFPBqpeCJl56ESrlq4EmPg1Z+NpX2KhauXZuKqqaxDL+ShlbVFxMJocEegn9pTFd04mlnpsmXs+O65VGzYAID74EEOPvxwzXbv4cN888OpVGzZQvbpA9kxcSLG1/CXYBW/IllZjwF2GmN2G2OqgdeAi2vvYIz5rzGmPPBwJdA9gvE0uzV7jnDpUysQgdduOYt+mWnRDqlFOMMYDe72ua09UUgzCCZr61fWMdAM7kj033pbT5V59M032Tf953gOHWLPlVfhzstj57kTarafnr2VTr+7l+qcHPb8aBoAnrwDbBs0OEoRKyuLZLLuBuyr9Xh/4LmG3AS8V98GEblFRNaIyJqCgoJmDPHkHSmr5jf/3EibJCf/+sW4mJ+VrCnsNsFuk0anG23Nl21BLFXWMdAMHvy/4qmMbhwnKf/Bh9h7yy1UbNwIQNmKFeTf/8Bx++yceF7N/cw770RE6HDttfUez52fH7lgVUyK5ACz+mYAqbfdVER+DIwGzq1vuzFmIbAQYPTo0Q23vbaAY+XVvPxFDm+t28+BokpeuOFMenRIjmZIUeGy20LOYFbtrW71fda15wa3tGC1auWWjmBlHYPJ+utzxuMtLASg7NPltLlwCiX//YSEPn3IevUV9lx+BdV79hz3mvSbflpzf+C2bA4++hhHX3kFU+XvBqjasQNn584t9h6U9UWyst4P9Kj1uDuQV3cnETkfuAeYaoyxdIfVil2HOX/+p8z/6GvaJDp55eaxnNMv9lbNag5Ou8R9M3jwy0hwalXLCvYDO5OiG0coNX3WsdUMfuyf/6pJ1Ok/u5nUc8+lbMUXJI85kx7PPYs9LY2+S49vMOz7/tITjtPprjs5fcNXNY+9R49GNnAVcyJZWa8G+otIbyAXuAq4pvYOIjISWABMNsYcimAsp2zfkXJufGE1PTsk87efnsngrm2jHVJUuRz20JW1XrplHe7AsJCYSNbWrqxNdTW+igrsbdviKSzkwD33ANBv2TKcnTLDOoYtreGxLV0fe5S8O+/Cc+RIs8SrWo+IVdbGGA8wA3gfyAbeMMZsEZG5IjI1sNtjQCrwpoh8JSKLIhXPqThYXMmMf6xHBF6+aWzcJ2oAV2OVtdfd+vusY6UZ3F3hv3VauLvG4qPBg6O1tw0bztdjz/KP3D5nfM32xhJ1tyfm19wXR8M1UpsLLwSgcOFfTzFi1dpEdFIUY8y7wLt1nptT6/75kTx/c1j1zRF++co6yqs9/PHKEXRumxjtkCzB5bA1Ot1oksPClVwziJnpRmOhsrYHkrXXesn6yMt/5+BDDzW4/fStWxo9Rup53w4uE2fD3UNi919v7tXKWtXReueCPEXGGF74/Buu+etK0hIdvHPrOUwe0iXaYVmG025rdLrRVj/ALGaawYOVtYWTdc0AM2sla2PMCYm6w/U/qbnf49lnEVvjf0Ztrm9bmUJV1ko1RP/X1MPrM9z9z428uXY/5w/sxPwrh9MmsXUnnqZqrLKOhz7rmJnBrKay1mbwpvIcOHDc44HbsgHIuO02KjZsIGXcuKYfNMxkbXy+sL4IqPig/xPq8PoMs9/4ijfX7udXE/ux8LozNFHXw2m3hZxutNJTSaKjdXcZxFRlLTaLT4pizQFmhxcsrLmfeeedNfdtKSknl6gBkdDr2mfe/RsAfMXFJ3V81TppZV3Lmj1HeGTpNlbvOcodkwYwY2LrW96yuTRWWVd4Klp9n3XMzGDmqfRX1Y0kiaiyaGVtb9cOgP6ff4YjPb1Fzhlc3MN77FjN+ZXSyjrg420HuXzBF+QUlvP7y4Zqom5EgiN0n3Wlp7LVJ+uYGmBm5f5qqDXdqLWSdeGCBQDYOzTD4jyBwWON7hZI0O46TfAqvmllHfDXT7+ha9skPrj9u6Qk6I+lMc4QM5j5jI9Kb+tvBq+ZFMXqlbW7Aqz+xalmulFrJeugxpquw9F36XsnzGRWn2CyLn73PVLOPvuUz6taB81KwMb9x/hidyF3XnCaJuowuew23J76Z36tClRHrb2yDl5nbf1kHQOVtYj/8i2L9VnbUlNJnfi9ZjmWq0cPXD16NLpf4tChADgyMprlvKp10GZw4JlPdtE+2clPzo7J5bSjwulouLKu8PgvFUq0t+7K2m6zYxe79acbdVdYP1mDv9/aQtON+ioq8JWWktCnb4ueNzgCvHz9uhY9r7K2uE/WlW4vy74u4AfDupCmo77D5gpxnXVloDpq7ZU1gMvuioHKusLal20FOaxVWQdXvnJ2ic6CGuVfrorKeZU1xX2b7+ur91Fe7eX7g3SFm6ZwOaTByjqekrXD5rB+sq4uheSWGcl8ShyJlkrWnoMHAXB0bvnJkJJGjQIT1QUGlcXEfWX91+W7GdO7A9/tH5+rZ50sl73hS7dqmsFb+QAz8I8It3wzeGURJMbAJUCuVP8XC4sIJutwF+hoTo5OmXiPHWvx8yrriutkbYzhUEkVI3q0a5bRnvEk1HSjcZesrV5ZVxZBYgwsPpOQClUl0Y6ihvugfyFAR2YUknX79jo/uDpOXCfrsmov1R4fHVIsPLOTRYWaFKXSGz/N4C67i2qfdQZFncCYGErWaVBlrcralpaGLbnl+/vt7drjLS7GeCx+Db9qMXGdrI+W+f/IarJuOqfdhttr8PlO7FeLl9HgEAPN4NVl4PPERrK2WjN4QUHULp+yd+gAxuDVKUdVQFwn68JAsk7XZN1kLof/v059g8ziaYCZ0+ak2luNzzQ8m1tUVRb5b2MhWSe0sVQzuOdIYYtNMVqXvb1/jIE2hauguE7WR8r8k3e012TdZC67/79OfU3hwco6XpL1J/s/4axXz4p2KPWLqWSdaqlmcO/hQuwdo5Osa+YHP3o0KudX1hPXybqwVCvrk1VTWdczyCyeBpi5AtNkBt+z5QSTdVKsjAYvscwlS57CQhzp0blKJDgXuUeTtQqI62R9RPusT5qzprI+8Q9rsBk8HpJ1cDEPy4qpyjoNjO/b9bejyFdVha+kBEd6MyzgcRLswcr6iCZr5Rf3ydplt5Gq84E3WajKutJbicPmsH4iawaWb+qvSdYxUFknpPpvLdAUHuwrtketzzq4TKYma+UX18m6sKyaDikuvcb6JIQaYFbhqSDJbvEk1ky6p3WPdgihlRf6b5PaRzeOcCS08d9WRX8EtOew/+fm6BidZnCby4UtJUX7rFWNuE7WRwLJWjWdy+7/grPj4Imjdys9rX95zKAeaY2vohRVxbn+aTxjIVkHp0QtOxzdOADPQf+84NFc+UoSE6ncsjVq51fWEtfJurCsmvRUTdYnI9FpB+Chd7NP2FbhqbB+83Az6dmmZ7RDCK04F9p08y9BaXWpnfy3ZYeiGwdQvX8/AM7u0Ws58RYWUr5mTdTOr6wlrpP1kbIqraxP0ll9/FVQ9/YnJuUKT0XcVNbdUrtFO4TQivOgTddoRxGe1MC0nqXRT9bu3DxsycnY20W/r99XHv0Bdyr64jpZHy1za7I+SYlOO2N6d6CeCcwoKC8gPTEGVnlqBt1TLd5nXZQLbS0eY1ByOogNSg9GOxLceXk4u3WzxHiWY2+9Fe0QlAXEbbKu8ngprfLoNdanINllp6Lae8LzuaW5dEuzeMXZTJx2JzcPvRmHWPCKAp8XSg7ETmVts0NKhjUq67w8HF1bfmnM2no+/xwAkhAfrVQqtLhN1sEJUTqkJEQ5ktiV7LJTXn38QgNl7jKOVh21fsXZjBw2Bx7jwVhkMo8aRfvBeGOnsgZ/v3VxXrSj8FfWXaP7JcfVuzcA+ffdF9U4lDXEbbIuKPFPNZqRpsn6ZCU5HVS6j790a3+Jf2BOvFTWQE1V7TEWWyEpb53/tsvw6MbRFOn9oHBHVENwHzyIr6iIhN59ohpH8FprpSCOk/WhQLLO1GR90uqrrHNLcwHokWrxS5qakcMWSNY+iyXr/WvAngCdhkY7kvB1HABHc8BdGbUQqr/5BoCE/v2iFgOALfHb5m9TbeFlWFWLiONk7f9jkNlGk/XJ8ifrb/us/7P3P8z870wgBkZJNyPLJuvctdBlGDhiaFxGxgDAwOGvoxZC5ebNwLfN0NGUfsstABS//0GUI1HRFr/JurgKEeiYqsn6ZCW57FR5fHgDQ8Lf2flOzba2CTEwF3UzsWSyri6DvPXQ/cxoR9I03Ub7b/d+EbUQDv3hcQCcnTtHLYagpFEjAci7884oR9p3vAIAAAunSURBVKKiLX6TdUkVHZJdNQtSqKZLdvknRqlwe/EZHxsObeB7Pb7HoksWWeKSl5YSnAPdUsl66yLwVMLpF0U7kqZpnwUd+sD2d6NyevdB/2VjiYMHR+X8daVNmFBz31cZva4BFX0RzVQiMllEtovIThG5u57tCSLyemD7lyLSK5Lx1HagqEIHl52iZJe/oiyucPNZ7mccrTrKeT3Po3fb6DcftiTLVdblR+Cj+6HTEOh5drSjabqhV8DuZf5m/BZ26LE/ANBt/uMtfu6GJJx+OgDbR4yMciQqmiKWrEXEDjwFTAEGAVeLyKA6u90EHDXG9AOeAB6JVDy1HSqp5ItdhYztHZ3l72KdMYbCikI8CZtI6Pw21753Obf+51YykzP5TvfvRDu8Fpdo9w8E+mT/J9ENBKDga3j9x1B+GC55Gmwx2HI0drr/crN/XAOb/wm+ExeLaW6+ykoO//WvFC9ejKtXL1xZWRE/Z7iyXvpbzf3STz+NYiQqmiI5k8MYYKcxZjeAiLwGXAzUnpn+YuD+wP23gD+LiJgIXLC6dPMBXvlyL3sKy9h/tAKAH50RQ9efRtELm19g7cG1FFUVUVRdRH5ZPhUe/8/Q1U44VNSTi/v9knvPvZ5kZ3KUo215E3pMYGyXsTz85cMs37+c/xnzP/Ro04Kj4avLYMlsOJQNB77yL9xxyV9i65Kt2pI7wLVvwt9/BG/9FJwz/O9l2JVwxg0RmefcV1FBwePzAWtV1QD2Nm2wd+yI9/BhytesJfW73412SCoKIpmsuwH7aj3eD4xtaB9jjEdEioB04Lhld0TkFuAWgJ49T27hhLIqL8UVbkb1bM8VZ/TgOwMyGNY9+vP+xoL8snwOlR+iTUIb+iX345yu59A9rTt92/Wle+JA7nl7Oz/q3z8uEzVAsjOZBecv4LnNz/HG9jda/ufgSIJ9q6BtNzj/fhjxY0iN3mpRzSJzINy6yt+cX34YDm6FPZ/5k3UEONq3p9c/3yKhTx9sSdZbhKb/8k/xHjuGQ6+9jlsSqVmXRORy4AJjzM2Bx9cBY4wxv6q1z5bAPvsDj3cF9ils6LijR482a3QlGmVRxpi4GlzXYnw+8Hli6zI0ixGRtcaY0dGOQ52cSHZo7QdqtwV2B+rOI1izj4g4gLbAkQjGpFREaaKOEJtNE7WKa5FM1quB/iLSW0RcwFXAojr7LAKuD9yfBnwcif5qpZRSKpZFrM860Ac9A3gfsAPPG2O2iMhcYI0xZhHwHPCyiOzEX1FfFal4lFJKqVgV0XX9jDHvAu/WeW5OrfuVwOWRjEEppZSKdTF4EaZSSikVXzRZK6WUUhanyVoppZSyOE3WSimllMVFbFKUSBGRAiAnAoduCxRF4LhWEGvvzarxRjuuaJy/pc7ZkTozF6pml2WMifGp7eJXzCXrSBGRhcaYW6IdRyTE2nuzarzRjisa52+pc4rIGp1dS6mGaTP4t/4d7QAiKNbem1XjjXZc0Th/tN+zUgqtrJVSFqCVtVKhaWWtlLKChdEOQCkr08paKaWUsjitrJVSSimL02StlFJKWZwm61MgIgNF5C8i8paI/CLa8cQz/V0opVqzVpOsRcQuIutFZPEpHON5ETkkIpvr2TZZRLaLyE4RuRvAGJNtjPk5cAWgI1kBEWkXSJjbRCRbRM4+yePo7yLOiUiKiKwVkYuiHYtS0dZqkjUwE8iub4OIZIpIWp3n+tWz64vA5HpebweeAqYAg4CrRWRQYNtU4DPgP6cSfCvyJ2CpMeZ0YDh1fif6u4hfDX0Bq+/LV8BvgDdaNkqlrKlVJGsR6Q78AHi2gV3OBf6fiCQG9v8Z8GTdnYwxnwJH6nn9GGCnMWa3MaYaeA24OPCaRcaYccC1p/xGYpyItAG+CzwHYIypNsYcq7Ob/i7i14vU+QLW0JcvETkf2AocbOkglbIiR7QDaCZ/BO4C0urbaIx5U0R6A6+JyJvAT4HvN+H43YB9tR7vB8aKyATgMiABePck4m5t+gAFwAsiMhxYC8w0xpQFd9DfRfwyxnwqIr3qPF3z5QtARIJfvlKBFPwJvEJE3jXG+FowXKUsJeaTdaA/65AxZm3gD3a9jDGPBv4QPAP0NcaUNuU09R/SfAJ80oTjtHb/v707i7VrCsA4/v9UaJWgqjG7iJijNSRaJKgnU2JugpgSEkHEkGhCiCnEFMQDRUiIIU21JELFGEUVpVrDg5CoeaiqmZvPw9q32U6O3nvacvc5vt/T3uustdfe9zx8e629zt1rAnsA59ieK+kW4GLg0nqlfBdR0/bmy/bZAJJOAb5JUMf/XS9Mg+8LHCHpY8qU6EGS7m+tJGl/YFfgUeCyDvtYDGxZ298C+Gylzra3LQYW255b7U+nhPff5LuImrY3X8s37Httr/Si0Yhe0fVhbXuq7S1s9wFTgGdtn1ivI2kCMI0yvXYqMEbSVR10Mw/YXtI2ktaq+nlstVxAD7H9BfCJpB2qosmU547L5buIFrn5ihiCrg/rIVoHONb2h9V02sm0eSe2pAeBV4AdJC2WdDqA7T+Bs4GnKKubH7G96D87++5yDvCApAXAeOCals/zXURdbr4ihiD/Gzwi/hPVDdgBwFjKKu/LbN8t6RDKItERwD22rx6+s4xopoR1REREw/1fpsEjIiK6VsI6IiKi4RLWERERDZewjoiIaLiEdURERMMlrCMiIhouYR0REdFwCeuIiIiGS1hH9AhJt0l6U9Lew30uEbF6JawjeoCk0cA44EzgsGE+nYhYzRLW0XUk3SzpvNr+U5Luqu3fKOn81dxnJ+/cHsrxNpB0Vm2/T9LCIbYdJekFSSMGymz/BGxKeaf3rZLWkvSipK5/Z31EJKyjO70MTAKQtAblxRC71D6fBMwZhvPqxAbAWYPWau80YIbt/oECSRtR3mi2DOi3/TvwDHD8qp5oRAy/hHV0ozlUYU0J6YXAMkkbSlob2AmYL2mmpDckLZJ0xkBjSde1jGovl3SBpBMlvSbpLUl31Eeutbpt61Qj4/ckTav6my1pVPXZpZLel/S0pAclXQhcC2xXHef66vAj2rVv4wRgVkvZJcANwCJg56psZlU3Irpcwjq6ju3PgD8lbUUJ7VeAucBEYC9gQTWyPM32nlXZudXoE+Ah/j7iPA54vSrb1/Z4oJ+WoJO00yB1tgdut70L8D1wtKS9gKOBCcBR1bkAXAx8aHu87Yv+qX3rtVfvfN7W9se1sr7q7/Aw5R3fA7MMC4EsNovoAXmeFd1qYHQ9CbgJ2LzaXkqZJocS0EdW21tSwvBb2/MljZO0GbAxsATYDdgTmCcJYBTwVUufkwep85Htt6rtN4A+yhT9LNu/AEh6fAXX1K59q7GUIK+7CrjCtiUtD2vb/ZJ+l7Se7WUr6DciGi5hHd1q4Ln1bpQR5CfABcAPwD2SDgAOBiba/lnS88DIWvvpwDHAJpSRtoD7bE9dQZ+D1fmttt1PCXN1cE3t2rf6hdp1SBpPGbHvJ+n26rN3avXXBn7t4BwiooEyDR7dag7lJ0rf2e63/R1l0dZEyrT4+sCSKqh3BPZpaf8QMIUS2NMpi7GOkTQOQNIYSVu3tBlKnVYvAYdLGilpXeDQqnwZsF6nF217CeXZ9kBgXwccbrvPdh+wO9XIupr2/9r2H532ExHNkrCObvUOZUr41Zaypba/AZ4E1pS0ALiypR62F1HC8lPbn9t+l7JIa3bV5mnKT6HqbQat08r2POAx4G1gBuXZ+FLb3wJzJC2sLTAbqtmUkfRBwGjbz9T6+xIYLWkMcCDwRIfHjogGku3hPoeIniZpXds/SloHeBE4w/abq3C8CcD5tk8apN4MYKrtD1a2r4hohjyzjvj33SlpZ8rz5PtWJagBqgVyz0kaUf+tdV21anxmgjqiN2RkHRER0XB5Zh0REdFwCeuIiIiGS1hHREQ0XMI6IiKi4RLWERERDZewjoiIaLiEdURERMMlrCMiIhruL3tkLWTKOg0yAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for b in bands:\n",
    "    plt.plot(Table(data = parse_single_table(FILTERS_DIR + b + '.xml').array.data)['Wavelength']\n",
    "            ,Table(data = parse_single_table(FILTERS_DIR + b + '.xml').array.data)['Transmission']\n",
    "            , label=b)\n",
    "plt.xlabel('Wavelength ($\\AA$)')\n",
    "plt.ylabel('Transmission')\n",
    "plt.xscale('log')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.title('Passbands on {}'.format(FIELD))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### IV.a - Depth overview\n",
    "Then we plot the mean depths available across the area a given band is available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ukidss_j: mean flux error: 6.922429084777832, 3sigma in AB mag (Aperture): 20.60655057467303\n",
      "90prime_g: mean flux error: 0.2022828757762909, 3sigma in AB mag (Aperture): 24.44229906520416\n",
      "90prime_r: mean flux error: 0.31267088651657104, 3sigma in AB mag (Aperture): 23.969478250706167\n",
      "mosaic_z: mean flux error: 1.0103346109390259, 3sigma in AB mag (Aperture): 22.69603378612627\n",
      "ukidss_j: mean flux error: 13.648526191711426, 3sigma in AB mag (Total): 19.86948246944096\n",
      "90prime_g: mean flux error: inf, 3sigma in AB mag (Total): -inf\n",
      "90prime_r: mean flux error: inf, 3sigma in AB mag (Total): -inf\n",
      "mosaic_z: mean flux error: inf, 3sigma in AB mag (Total): -inf\n"
     ]
    }
   ],
   "source": [
    "average_depths = []\n",
    "for b in ap_bands:\n",
    "    \n",
    "    mean_err = np.nanmean(depths['ferr_ap_{}_mean'.format(b)])\n",
    "    print(\"{}: mean flux error: {}, 3sigma in AB mag (Aperture): {}\".format(b, mean_err, flux_to_mag(3.0*mean_err*1.e-6)[0]))\n",
    "    average_depths += [('ap_' + b, flux_to_mag(1.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(3.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(5.0*mean_err*1.e-6)[0])]\n",
    "    \n",
    "for b in tot_bands:\n",
    "    \n",
    "    mean_err = np.nanmean(depths['ferr_{}_mean'.format(b)])\n",
    "    print(\"{}: mean flux error: {}, 3sigma in AB mag (Total): {}\".format(b, mean_err, flux_to_mag(3.0*mean_err*1.e-6)[0]))\n",
    "    average_depths += [(b, flux_to_mag(1.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(3.0*mean_err*1.e-6)[0], \n",
    "                        flux_to_mag(5.0*mean_err*1.e-6)[0])]\n",
    "    \n",
    "average_depths = np.array(average_depths,  dtype=[('band', \"<U16\"), ('1s', float), ('3s', float), ('5s', float)])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def FWHM(X,Y):\n",
    "    \n",
    "    half_max = max(Y) / 2.\n",
    "    #find when function crosses line half_max (when sign of diff flips)\n",
    "    #take the 'derivative' of signum(half_max - Y[])\n",
    "    d = half_max - Y\n",
    "    #plot(X,d) #if you are interested\n",
    "    #find the left and right most indexes\n",
    "    low_end = X[np.where(d < 0)[0][0]]\n",
    "    high_end = X[np.where(d < 0)[0][-1]]\n",
    "    return low_end, high_end, (high_end - low_end)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = []\n",
    "\n",
    "for b in ap_bands:\n",
    "    data += [['ap_' + b, Table(data = parse_single_table(FILTERS_DIR + b +'.xml').array.data)]]\n",
    "    \n",
    "for b in tot_bands:\n",
    "    data += [[b, Table(data = parse_single_table(FILTERS_DIR + b +'.xml').array.data)]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ap_ukidss_j (11695.0, 13280.0, 1585.0)\n",
      "ap_90prime_g (4180.0, 5470.0, 1290.0)\n",
      "ap_90prime_r (5680.0, 7150.0, 1470.0)\n",
      "ap_mosaic_z (8552.0, 10018.0, 1466.0)\n",
      "ukidss_j (11695.0, 13280.0, 1585.0)\n",
      "90prime_g (4180.0, 5470.0, 1290.0)\n",
      "90prime_r (5680.0, 7150.0, 1470.0)\n",
      "mosaic_z (8552.0, 10018.0, 1466.0)\n"
     ]
    }
   ],
   "source": [
    "wav_range = []\n",
    "for dat in data:\n",
    "    print(dat[0], FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission'])))\n",
    "    wav_range += [[dat[0], FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission']))]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Depths on SA13')"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAEgCAYAAABcsLEAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XlYlWX6B/DvfdgRXFAQRQQ3VBQVMyytn46VYpm5ZYtaP21xmaZMzSxnmsmyckpLpyxtHR37ZbnlkprOuKVTCa6sooYrKIogCCJw7t8f52BHBAHlwOHl+7kur+uc933e97nPoav7fZ73Oe8tqgoiIiIyJlN1B0BERET2w0RPRERkYEz0REREBsZET0REZGBM9ERERAbGRE9ERGRgTPRUK4lIsojcW91xEBHZGxM9VTtr0s0VkSwRyRCRXSIyTkQq5b9PEflKRN6sjHPZk4i8KiK/iUi2iJwUkaUltPlKRApEpGmx7R1FZKOInBOR6x6OISL/EpEUEbkoIodE5Gl7fhYichxM9OQoHlRVbwBBAN4B8DKAz6s3pKojIk8CGAXgXlX1AtANwL+LtakDYCiATAAjip0iH8C3AJ4qpYu3AQSral0AAwG8KSK3Vd4nICJHxURPDkVVM1V1NYBHADwpIh0BQETcROQ9ETkuImdE5BMR8bDu620dAb9qHdEmi8gI675nYUmKU60j5TU23XURkQMikikiS0XE3XpMIxFZa51dSBeRHaXNLohIDxHZbT3HbhHpYbNvq4i8ISI7rbMVP4pIo1I++u0ANqrqEev3kKqqC4u1GQogA8AMAE8W+94SVfVzALGlfK+xqppX9Nb6r1UpsRCRgTDRk0NS1V8BnARwt3XTLAAhALoAaA0gAMBrNof4A2hk3f4kgIUi0taaLJcA+LuqeqnqgzbHDAcQCaAFgE4A/te6fbK1b18AjQG8CktivIaI+ABYB2AegIYA5gBYJyINbZo9DmA0AD8ArgCmlPKRfwbwhIi8JCLdRMSphDZPAvg/AN8AaCciXUs5V4lEZL6I5ABIAJAC4IeKHE9ENRMTPTmy0wB8REQAPAPgRVVNV9UsAG8BeLRY+7+oap6qboMlAQ8v4/zzVPW0qqYDWAPLRQRgmQZvAiBIVfNVdYeWXBTiAQBJqrpYVQtU9f9gSaK2FxNfquohVc2FZWq9Swnngar+C8CfAPQDsA3AWRGZVrRfRJoD+AOAr1X1DCzT+k+WdK7SqOoEAN6wXDytAJB34yOIyAiY6MmRBQBIh2Vk7Qkg2jqdngFgg3V7kQuqesnm/TEA1yxYK0GqzescAF7W1+8COAzgRxE5aptwi2lq7cfWMWvcZfVxHVVdoqr3AqgPYByAGSLSz7p7FIB4Vd1nfb8EwOMi4lLa+Urpo1BVfwLQDMD4ihxLRDUTEz05JBG5HZaE+ROAcwByAXRQ1frWf/Wsi9aKNLAuVivSHJYZAaCEafcbUdUsVZ2sqi1hGZ1PEpF7Smh6GpbFg7aaAzhVkf5K6D9fVb8DcABAR+vmJwC0FJFUEUmF5TZBIwD9b7IbZ/AePVGtwERPDkVE6orIAFjuQ/9LVQ+qqhnApwDeFxE/a7sAm9FukddFxFVE7gYwAMB31u1nALSsQAwDRKS19ZbBRQCF1n/F/QAgREQeFxFnEXkEQCiAteX/xFf7/F8ReUBEvEXEJCL9AXQA8IuI3AlLUo6AZeq/CywXAF/DOn0vFu6wrAOAiLiLiJv1tZ+IPCoiXiLiZP3eHgPwn4rGSUQ1j3N1B0BktUZECgCYAcTBMmL9xGb/y7AsvvvZunL9FICPAWy07k8FcAGWUXYOgHGqmmDd9zmA76xT/ltVdVAZsbQB8CEstwYuAJivqluLN1LV89aLkrnWWA4DGKCq5yrywa0uwrLo718AnGC5BTBeVX8SkU8AfK+qB20PEJG5AHZYFwXWBfCbze5c6zmCYZnRGA/L92mybp+oqt/fRJxEVMNIyWuMiGoOEekNy+i/WXXHQkTkaDh1T0REZGBM9ERERAbGqXsiIiID44ieiIjIwJjoiYiIDKxG/LyuUaNGGhwcXN1hEBHVKNHR0edU1bfslqUe7+fs7PwZLM9t4MDQMZkBxBQUFDx92223nS2pgd0SvYgEAlgES7ERM4CFqjrXZv8UWB416lvW746Dg4MRFRVlr1CJiAxJRIo/orlCnJ2dP/P392/v6+t7wWQycUGXAzKbzZKWlhaampr6GSwlqK9jzyu0AgCTVbU9gDsA/FFEQoGrFwH3AThux/6JiOjWdPT19b3IJO+4TCaT+vr6ZuL3x2Vf38Zenatqiqrusb7OAhCP34t9vA9gKir4DHIiIqpSJiZ5x2f9G5Waz6vknouIBAMIh+W53QMBnFLV/VXRNxERUW1m90QvIl4AlgOYCMt0/nRYnlle1nHPikiUiESlpaXZOUoiIqpt5s2b1/CJJ55oXnz73//+d98PP/ywYfHtiYmJrm3atOlQGX0/8sgjQdHR0e6Vca6y2HXVvbVW9nIAS1R1hYiEAWgBYL+lMBiaAdgjIhGqalu3G6q6EMBCAOjWrRunjoiIqEpMnTrV7qPLpUuX3tJCyYqw56p7gaVqWLyqzgEAa/UtP5s2yQC63WS1LyIiqiIvLdsfeCg1y7Myzxni753z7rDOJ8pqd++997ZKSUlxzcvLM40bN+7MlClTznl6eoaPGDEibefOnd716tUrXL58+dGmTZsWlHR8RERE2/fee+/E//zP/+SkpKQ4d+vWrf2pU6euqQb5zTff1HvnnXearF+//vC7777r5+XlVThjxowzO3bs8Hz66aeDPTw8zN27d88uah8VFeU+evToFvn5+WI2m7F8+fIjQUFB+QMHDmyZkpLiajabZerUqaefeeaZC2XFVNHvraLsOaLvCWAUgIMiss+67VVV/cGOfV5r/TQg9WDZ7Wor/zCg/zvVHQUR0Q0tWbIkuXHjxoXZ2dkSHh4eOnLkyAu5ubmmrl275nz66acnp0yZ0mTatGlNFy1adFO/5Fq0aFH9uXPnNt60aVOSr69voe2+p556Kvj9998//sADD2SPHTv2aoXMf/zjH74TJkw4M378+PTLly9LQUEBli1bVs/f3z9/69athwHg/PnzTrf2ySuH3RK9qv4EQMpoE2yv/omIqPKUZ+RtL7NmzWq8bt26+gCQmprqEhsb624ymfD000+nA8CYMWPODxkypPXNnHvXrl3e+/fv99yyZcshHx8fs+2+8+fPO2VlZTk98MAD2UX9/Oc//6kHAHfeeeel9957r8nJkyddH3300QthYWF5Xbt2zZ0+fXrg+PHjAx566KHMyMjI7JL6rGo14sl4N42jVSKiGm3t2rXe27Zt846Kikrw9vY2R0REtM3Nzb1uIbl13VeJnJ2dtbDQMlDPycm5pmHz5s3zjh8/7hYTE+NefBpdVUs977hx49LvvvvuSytXrqzXv3//kPnz5ycPHDgwa8+ePXHLly+vN3369IDNmzdffO+991Ju4mNXKj7SkIiIHFZGRoZTvXr1Cr29vc179+51379/fx0AMJvN+PLLLxsAwFdffdUwIiIiq7RzBAYG5v366691AGDJkiUNbPc1a9bsyvLlyw+PHj26RVRU1DWr4Bs1alTo5eVVuHHjRi9rPz5F++Li4lzbt2+f9+c///ls3759M/bt2+eRnJzs4u3tbZ4wYUL6xIkTz+zbt69S1zTcLGOP6ImIqEYbOnRo5sKFC31DQkJCW7Vqdblz586XAMDDw8McGxvr0aFDB39vb+/CFStWHC3tHNOmTTvzyCOPtPzmm28a3n333ReL7+/cuXPeokWLjj7yyCOtVq9efdh23+eff55ctBivT58+V49dvHixz3fffdfQ2dlZfX19899+++3TP/30U51XXnmlmclkgrOzs86fP/+GK+tvNAtRmWpEPfpu3bopn3VPRFQxIhKtqt1u9vj9+/cnd+7c2SF/FeXp6Rmek5Ozt7rjuFkhISGhq1evPtyuXbsrlXG+/fv3N+rcuXNwSfs4dU9ERFSFevTo0aZt27a5lZXky8KpeyIiqnFKGs2PGjWq+e7du71st40fP/7MCy+8cL7qIrvWfffd1+rEiRNutttmzpx5cujQodfdQrAXJnoiIjKExYsXO1xF1E2bNh2p7hg4dU9ERGRgTPREREQGxkRPRERkYEz0REREBsZET0REhvbf//7Xo0uXLu1CQkJC+/Tp0zo9Pf1q7nvllVf8mzdv3jE4OLjj8uXL61b03BMnTmy6atUq78qNuHIx0RMRkaE988wzwTNnzjx56NChuIEDB154/fXX/QEgOjrafcWKFT6JiYmxGzZsODRx4sTmBQUlVrotUUFBAT744IPTgwYNKvXxu46AP68jIqKyrfpjIM7GVe6z2/1CczDoI7vXo09OTnbv379/NgAMGDDgYr9+/ULmzp17etmyZfWHDBmS7uHhoe3atbsSFBSUt3Xr1jqBgYH5kZGRbcLDwy/FxMR4tmzZ8vJ3332X7O3tbQ4ICAh77LHHzm3ZsqXu2LFjz27cuLHegAEDMkePHn0hICAgbPDgwek//fSTd0FBgXzyySfHpk2bFnDs2DG3P/3pT2emTp2aBgB/+ctfGq9cudLnypUr8sADD2S8//77p0v77C+99FKTZcuW+TRp0uRKw4YNC8LDw3NmzJhxpiJfM0f0RETk0JYsWZIcGxsbv2/fvrgFCxY0Tk1NdSqqRx8XFxffs2fPrGnTpjUt7fg2bdrkfv311/UB4F//+pdPamqqKwCcOnXKNTAw8OrT6Zo2bXrlxIkTroDl4mDcuHFphw4divP29ja/++67vkXt3N3dzdHR0YnPPvvsheJ9BQYGXtm3b19C9+7ds8eMGRO8Zs2aI7/88kvCO++80xQAVqxYUffw4cPuBw4ciI+Pj4/bt2+f5/r1672KnwcAtm/f7rlmzZoGBw8ejFu3bt2RAwcO1LmZ748jeiIiKls5Rt72cqv16L/44ovk5557LvDtt99uEhkZmeHi4qKApQxtcSKiAODv73+lb9++lwBg1KhR5+fNm+cH4AwAPPHEE9cl+CLDhw/PAICwsLCcS5cumRo0aGBu0KCB2c3NzXzu3DmnDRs21N2+fXvd0NDQUADIyckxJSQkXJ1xsLV161av/v37Z3h5eSkAve+++zLK+ZVdg4meiIgcVmXUow8PD7+8c+fOJAA4cOCA248//lgfsJSoLRrBA8Dp06ddmzVrll/S+Wzfe3t7m0vry93dXQHAZDLB1dX16pWEyWRCfn6+qComTpyY8tJLL5VZLKiyis5x6p6IiBxWZdSjP3XqlDMAFBYW4q9//WuTp5566iwADB06NGPFihU+ubm5kpCQ4JqcnOzeu3fvSwCQkpLiunnz5joA8PXXX/v06NHjuhH3zejfv//FxYsXN8rMzDQBwG+//eZSFF9xvXv3zt64cWO9nJwcyczMNG3evLn+zfTJET0RETmsyqhH/8UXX/h8/vnnfgBw//33X3j++efPA0C3bt0uDxo0KD0kJKSDk5MT5syZc8zZ2ZIWW7ZsefmLL75oOGHChKAWLVrkTZkyJa0yPs+QIUMuxsbGut9+++3tAMDT09O8ZMmS3wICAq5bSNirV6+cyMjIzNDQ0A4BAQF5nTp1ulSvXr3CivbJevR0Q7N+nYWE9ITqDqPKtfNph5cjXq7uMIhuCevR35zExETXAQMGtElKSoq1x/krIjMz01SvXj1zVlaW6c4772z7ySefHLvrrrtyire7UT16juiJiIgc1MiRI4OSkpI88vLy5NFHHz1fUpIvCxM93RBHtUTkiOxZj75t27ZXqnI0n5qa6tS7d++2xbdv3bo1cc2aNb/d6vmZ6ImIyBAcsR59efj7+xcmJCTE2ev8XHVPRERkYEz0REREBsZET0REZGBM9ERERAbGRE9ERIbGevREREQGZoR69GazGYWFFX4oHgD+vI6IiMrhLzv/Enj4wuFKrUffukHrnDd6vsF69KXUo09MTHTt379/mx49emRFR0d7ff/994dDQkKulNT2RjiiJyIih1Zb69EXxTF69Ojz8fHxcTeT5AGO6ImIqBzKM/K2l9pajx4AmjRpcuWee+65VJ7vqTRM9ERE5LBqcz16wFLdrjztboRT90RE5LBqcz36ymK3k4tIIIBFAPwBmAEsVNW5IvIGgIes284C+F9VLXEhAhER1W61uR59ZbFbPXoRaQKgiaruERFvANEABgE4qaoXrW2eBxCqquNudC7WoyciqjjWo785jlSPvrxuVI/eblP3qpqiqnusr7MAxAMIKEryVnUA2OdKg4iIiKpmMZ6IBAMIB/CL9f1MAE8AyATwh6qIgYiIjKO21KP39/e/uafk2LB7ohcRLwDLAUwsGs2r6nQA00XkFQDPAfhrCcc9C+BZAGjevLm9wyQiohqO9ehLZtdV9yLiAkuSX6KqK0po8jWAoSUdq6oLVbWbqnbz9fUtqQkRERGVwW6JXiw/OvwcQLyqzrHZ3sam2UAACfaKgYiIqLaz59R9TwCjABwUkX3Wba8CeEpE2sLy87pjAG644p6IiIhunt0Svar+BKCkRxX9YK8+iYiI6Fp8Mh4REVElccT69HzWPRERUSX54IMPHO5Jr0z0RERUptOvTg/MS0qq1Hr0bm3a5DR9a6bd69FHRES0DQsLy9m/f79nenq685dffvnbzJkzmyQmJno89NBD6fPmzTsNAH/7298aL1mypBEAjBo1Ku211147e/HiRdPAgQNbpqSkuJrNZpk6derpZ5555sKUKVOabNiwoX5eXp6pW7du2UuWLDlmMpkwdOjQ4KL69Nu2bfOcOHFi85ycHJOrq6tu3749sUGDBtcVqXnkkUeCip7hf+bMGZcxY8acnT17dsqtfbu/49Q9ERE5tFutRw8Arq6u5qioqMTRo0enPfzww60//fTT4wkJCbFLly5tlJqa6rRjxw7Pr7/+umF0dHR8VFRU/KJFi3x37tzpsWLFirr+/v75iYmJcUlJSbFDhgy5CAAvvfTS2ZiYmPikpKTY3Nxc0zfffFPPtr/Lly/LiBEjWn3wwQfHExMT47Zt25bo5eVVYiW6pUuXHktISIhbvXr14fr16xeMHTu2Qg/4KQtH9EREVKbyjLzt5Vbr0QPA4MGDMwCgc+fOua1bt84NCgrKB4DAwMC8o0ePum7dutXr/vvvz6hbt64ZAB544IELW7Zs8R44cGDm9OnTA8ePHx/w0EMPZUZGRmYDwPr1673nzJnjf/nyZVNGRoZzaGhoLixPewUAHDhwwN3Pzy+/V69eOQDg4+Nzw3KzOTk5MnTo0Fbvv//+8ZCQkCs3/WWVgCN6IiJyWLb16BMTE+Pat2+fW9F69MC1deLd3NyuqRNfUFBQaoG3Tp065e3ZsycuLCwsd/r06QFTpkxpkpOTI5MnTw5asWLFkUOHDsWNHDny3OXLl6+JSVUhIuWu5TJq1KigBx988MKgQYNKLbd7s5joiYjIYVVGPfry6NOnT/YPP/xQPysry3Tx4kXTDz/80OAPf/hDVnJysou3t7d5woQJ6RMnTjyzb98+z5ycHBMA+Pv7F2RmZprWrFnToPj5OnfufPnMmTOu27Zt8wSACxcumPLz80vs++233/bNzs52euutt1Jv5TOUhlP3RETksCqjHn153HXXXTmPP/74+a5du7YHLIvxevbsmbt8+fK6r7zySjOTyQRnZ2edP3/+sUaNGhWOGDEiLTQ0tEOzZs2uFMVky93dXZcsWXLk+eefb3758mWTu7u7efv27Yfq1at33RT+hx9+6O/i4qLt2rULBYAxY8akTZ06Ne1WPo8tu9Wjr0ysR09EVHGsR197VEs9eiIiIqp+nLonIqIax5716O1l+fLldadPn97MdltgYGDepk2bjtizXyZ6IiIyBEevRz906NCLQ4cOtVvd+dJw6p6IiMjAmOiJiIgMjImeiIjIwJjoiYiIDIyJnoiIapyIiIi227dvv66aXq9evVqfO3fOqfj2SZMmNX3ttdca32q/ycnJLpGRkS1v9TxViavuiYjIMLZt23bYnucPDg7O37Bhwy09ha+qMdETEVGZ/r0oPjD9VHal1qP3CfDKueeJ9jesipeYmOg6YMCANklJSbEA8NprrzXOzs6+OmIvLCzEww8/HNysWbMr8+bNOx0QEBAWFRUV36RJk4KXX37Zf+nSpY2aNm16pWHDhvnh4eE5APDmm2/6ffnll75OTk4aEhJyee3atUfXrVvnNXny5OaApUDOrl27EkqqHV88npqAiZ6IiGqk/Px8GTRoUIvQ0NDcWbNmXVMQZseOHZ4rV670OXjwYFx+fj66dOkSWpTo582b53/s2LGDHh4eWjTNP3v2bP958+Yd69u376XMzEyTp6fnDcvK1iRM9EREVKayRt7VYcKECUGDBg1KL57kAWDLli1e999/f4a3t7cZAPr27ZtRtK9t27a5gwcPbjFw4MCMESNGZADAHXfckT1lypTA4cOHpz/22GMXWrVqZZhEz8V4RETksJydndVs/j3n2tZ979atW/aOHTvq5uTklFiMvrQa9Vu2bEn64x//mBYdHV2nc+fOofn5+XjrrbdSP/vss2O5ubmmHj16tN+7d697ZX+W6sJET0REDqtZs2YF6enpzqmpqU65ubmycePGekX7xo4de65v376ZAwYMaFW81nufPn2y161bVz87O1suXLhg2rRpU33Ack//yJEjrg8++GDW/PnzT2ZlZTllZmY6xcbGukVEROTOnDkzNSws7FJMTIxhEj2n7omIyGG5ubnp5MmTUyIiIto3a9Ysr3Xr1pdt9//tb3878+KLLzoNGTKkxapVq34r2n7XXXflDB48OL1jx44dAgIC8iIiIrIBoKCgQB5//PEWWVlZTqoqY8eOPdOoUaPCyZMnN921a1ddk8mkISEhucOGDcssLSYRcfz67jZYj56IyKCMXI++uuzYscNz0qRJgbt3706s7lhssR49ERHRLdq+fbvnqFGjWj733HNnqjuWiuDUPRERUTG//vqrxxNPPNHCdpurq6s5OTk5prpiullM9ERERMVERETkJiQkVHnteHvg1D0REZGBMdETEREZGBM9ERGRgfEePZGBpL71FvLiE6o7jDK5tW8H/1dfre4wiGoFjuiJiMihvfHGG35t2rTp0Lp16w4zZszwA4AzZ8449ejRo01QUFDHHj16tElLS7uuBn1ZwsPD21V+tI6HI3oiA+EomYxm9+7d7osWLfLds2dPvLu7u7lXr14hgwcPzvzoo498e/funfXWW28lvfrqq/6vvfaa/8cff3yqPOcsKCiAs7Mz9u7d6/jTX5WAiZ6IiMq08eMPAs+dOFap9egbBQbl9Bs/8YZV8Q4ePOjRtWvX7KIqdD179sxaunRp/Q0bNtTftm1bIgCMHTv2fK9evdoCODVp0qSmR48edUtNTXVJSUlxff7551MnT558bu3atd5vvPFGEz8/v/y4uDjPI0eOxHp6eobn5OTsXbt2rffrr7/e1NfXNz8uLs7z/vvvvxAWFpY7f/78xnl5ebJy5cojHTp0yDt9+rTz6NGjg06dOuUKAHPmzDnet2/fSyXFffr0aedhw4a1yMjIcO7SpUvO1q1b60ZHR8c3adKkoDK/w/Kw29S9iASKyBYRiReRWBF5wbr9XRFJEJEDIrJSROrbKwYiIqrZunTpkvvLL794p6amOmVlZZk2bdpU78SJE67nz593DgoKygeAoKCg/PT09KsD1/j4eI/Nmzcn/fzzzwnvvvtu0+TkZBcAOHDgQJ1333331JEjR2KL95OQkODx8ccfn4iPj49dtmxZw0OHDrkfPHgwftSoUedmz57tBwBjx44NnDRp0pmYmJj4lStXHhk3blxwaXFPmzataa9evbLi4uLihwwZciElJcW10r+ccrLniL4AwGRV3SMi3gCiRWQTgE0AXlHVAhGZBeAVAC/bMQ4iIrpFZY287aVr166XX3jhhdQ+ffqEeHp6mkNDQ3OcnW+cuvr375/h5eWlXl5eBXfeeefFHTt21GnQoEFhp06dLrVr1+5KSceEhYVdKrpwaN68eV7//v0zAaBz586527Zt8waAnTt31k1KSvIoOiY7O9vpwoULpgYNGlxXu/7XX3/1WrVq1WEAGDZs2MW6desW3vSXcIvsluhVNQVAivV1lojEAwhQ1R9tmv0MYJi9YiAioprvxRdfPPfiiy+eA4DnnnsuoFmzZlcaNmxYcOzYMZegoKD8Y8eOufj4+FydEi9eh77ovaen53UJuYibm9vVCm8mkwnu7u5a9LqwsFAAQFURFRUV7+XlVWY1OEcqGFclq+5FJBhAOIBfiu0aA2B9Kcc8KyJRIhKVlpZm3wCJiMhhnTp1yhkAkpKSXNetW1f/qaeeSu/Xr1/GggULGgLAggULGkZGRmYUtV+/fn39nJwcSU1Ndfr555+977rrrhLvo1fUXXfddXHWrFl+Re937drlUVrbiIiI7MWLF/sAwIoVK+pevHixwr8KqCx2X4wnIl4AlgOYqKoXbbZPh2V6f0lJx6nqQgALAUuZWnvHSUREjmngwIGtMjIynJ2dnfWDDz447uvrW/j666+nDB48uFVQUFCjpk2bXlm1atWRovbh4eGX7rnnnjanT592nTJlSkpwcHB+TEyM+63GsXDhwhNPP/1085CQkNDCwkLp3r17Vo8ePY6X1Padd945PWzYsJahoaEN7rzzzmxfX9/8+vXrV8v0vV3r0YuIC4C1ADaq6hyb7U8CGAfgHlXNKes8rEdPRFRxtbEe/aRJk5p6eXkVzpgxo1pLyebm5oqzs7O6uLhg8+bNdZ577rkgexbJuVE9eruN6MVyU+RzAPHFknwkLIvvepUnyRMREdU0hw8fdh0+fHgrs9kMFxcXXbBgQXJ1xWLPqfueAEYBOCgi+6zbXgUwD4AbgE3WBRI/q+o4O8ZBRES1xJw5c05XZX9z585t+PHHHze23Xb77bdnL168+Hh8fLxDlLktV6IXETcAQwEE2x6jqjNKO0ZVfwIgJez6oWIhEhEROaYXXnjh/AsvvHC+uuO4kfKO6L8HkAn/tG2gAAATk0lEQVQgGkCe/cIhIiKiylTeRN9MVSPtGgkRERFVuvL+jn6XiITZNRIiIiKqdDcc0YvIQQBqbTdaRI7CMnUvAFRVO9k/RCIiIrpZZY3oBwB4EEB/AK0B9LW+L9pORERkV6xHf2tuOKJX1WMAICKLVXWU7T4RWQzLz+eIiMjg0pcdCsxPvVSpZWpd/Ovk+AwLuWGxHCPUo8/Pz4eLi0tVdFWi8t6j72D7RkScANxW+eEQERH9zrYevYuLyzX16MeOHXsesNSjX79+fQPA8mS8QYMGtbjjjjtCgoKCOs6ePbsRAKxdu9a7e/fuIQ8++GCLtm3bdgAAT0/P8KJ9t99+e9v777+/ZXBwcMcJEyYEfPzxxz5hYWHtQ0JCQmNjY90AS435fv36terYsWP7jh07tv/xxx/rlBb3pEmTmj722GNBPXv2bDNkyJAW9v6ebqSse/SvwPKQGw8RuYjffxd/Bdbn0BMRkfGVNfK2ly5duuTOmDEjIDU11alOnTq6adOmep07d75UVj366Ojo+KysLKfw8PDQoUOHZgKWevR79+6NLalUbUJCgseyZcuO+vn5FQQFBYW5ubmdO3jwYPwbb7zhN3v2bL8vvvjiRFE9+n79+mUnJSW59uvXr83Ro0evq21f5MCBA56//PJLQnmq3dlTWVP3bwN4W0TeVtVXqigmIiIiADW3Hj0AREZGZlR3kgfK/zv6V0VkCIC7YFmFv0NVV9kvLCIiIouaWI8eAOrUqVNqf1WpvPfoP4Kl2txBADEAxonIR3aLioiIyKom1qN3JOUd0fcC0FGtNW1F5J+wJH0iIiK7qon16B1JuerRi8gKAC/a/NwuCMA7qvqYneMDwHr0REQ3g/Xoa4/KqEffEEC8iPxqfX87gP+KyGoAUNWBtxwlERERVbryJvrX7BoFERFRJXCkevRVGceNlCvRq+o263R9G1XdLCIeAJxVNcu+4RERETmumlCPvlyr7kXkGQDLACywbmoGgD+vIyIicnDl/XndHwH0BHARAFQ1CYDfDY8gIiKialfeRJ+nqlefJiQizrA8OIeIiIgcWHkT/TYRKXrm/X0AvgOwxn5hERERVS2jlq0tb6KfBiANlofkjAXwA4A/2ysoIiKiqlZVZWurWnlX3ZtFZBWAVaqaZueYiIjIwaxatSrw7NmzlVqP3s/PL2fQoEE3rIqXmJjoGhkZ2SYiIiJ7z549Xu3bt88ZM2bMuRkzZgScP3/e+auvvjoaGhqaN2LEiODjx4+7eXh4mBcuXHise/fuuevWrfOaPHlyc8DyvPtdu3YlmEwmREZGts7MzHQqKCiQ11577fTIkSMzAEvZ2pycnL0A8Oc//7nxt99+21BEcM8992TOnz//ulr3ycnJLpGRkW2K3iclJXnEx8cfDAkJKbFwTnUpq0ytAPgrgOdgKVErIlII4B+qOqMK4iMiolruxIkT7kuXLj162223HevUqVP7JUuWNIyKikr4+uuv68+cObNJQEDAlc6dO+ds3rz5yOrVq72ffPLJFgkJCXGzZ8/2nzdv3rG+ffteyszMNBUVtVm3bt1hHx8fc0pKinP37t3bPf744xkm0+8T3N9++23ddevWNYiOjk7w9vY2nzlzxqmkuIKDg/MTEhLiAODtt9/23bFjh7ejJXmg7BH9RFhW29+uqr8BgIi0BPCxiLyoqu/bO0AiIqp+ZY287SkgICAvIiIiFwBCQkJy+/Tpc9FkMqFr1645b775ZtNTp065LV++/DAADBw4MOvZZ591Pn/+vNMdd9yRPWXKlMDhw4enP/bYYxdatWplzsvLk4kTJzb7+eefvUwmE86ePet68uRJ5+bNm1+tfrdp06a6I0eOPOft7W0GgMaNGxfeKL4ff/yxzqJFi3x//vlnh5z6L+se/RMAHitK8gCgqkcBjLTuIyIisitXV9cSS8g6OTmhsLCwxJotIqJvvfVW6meffXYsNzfX1KNHj/Z79+51X7Bggc/58+edDx48GJ+QkBDXsGHD/Nzc3GtyoapeV+q2NMeOHXMZO3Zs8NKlS4/Uq1fPIcrSFldWondR1esKGljv07vYJyQiIqLyu+OOO7K+/PLLhgCwdu1a7wYNGhT4+PiYY2Nj3SIiInJnzpyZGhYWdikmJsY9MzPTqVGjRvlubm66Zs0a79OnT7sWP19kZOTFxYsXN8rKyjIBQGlT93l5eTJkyJCWb7zxxqlOnTrl2fdT3ryypu5vdK/B4e5DEBFR7TNr1qzTjz/+eHBISEioh4eH+auvvvoNAP7+97/77dq1q67JZNKQkJDcYcOGZWZkZDj179+/dceOHdt36NAhp0WLFpeLn2/YsGEX9+zZ49mlS5f2Li4ueu+992Z++OGH1y3G27x5c52YmJg6b775ZtM333yzKQBs2LAhKTg4ON/+n7r8blim1rrw7lJJuwC4q2qVjOpZppaIqOJqY5na2uqmy9SqaonTFURERFQzlLdMLRERUa01atSo5rt37/ay3TZ+/Pgzjl65DmCiJyIiKpMj1ZevqPI+ApeIiGofs9lsLt/vzKjaWP9Gpf60j4meiIhKE5OWllaPyd5xmc1mSUtLqwcgprQ2nLonIqISFRQUPJ2amvpZampqR3Bg6KjMAGIKCgqeLq2B3RK9iAQCWATA3xrIQlWdKyIPA/gbgPYAIlSVv5sjInJAt91221kAA6s7Dro19hzRFwCYrKp7RMQbQLSIbIJlemEIgAV27JuIiIhgx0SvqikAUqyvs0QkHkCAqm4CUO7nCBMREdHNq5J7LiISDCAcwC8VOOZZEYkSkai0tDR7hUZERGRodk/0IuIFYDmAiap6sbzHqepCVe2mqt18fX3tFyAREZGB2TXRi4gLLEl+iaqusGdfREREdD27JXqx3IT/HEC8qs6xVz9ERERUOnuuuu8JYBSAgyKyz7rtVQBuAP4BwBfAOhHZp6r97BgHERFRrWXPVfc/wVLOtiQr7dUvERER/Y5POiIiIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA3O214lFJBDAIgD+AMwAFqrqXBHxAbAUQDCAZADDVfWCveIgIqqo19fEIu70xeoOAwAQ2rQu/vpgh+oOg2owe47oCwBMVtX2AO4A8EcRCQUwDcC/VbUNgH9b3xMREZEd2G1Er6opAFKsr7NEJB5AAICHAPS2NvsngK0AXrZXHEREFcURNBlJldyjF5FgAOEAfgHQ2HoRUHQx4FfKMc+KSJSIRKWlpVVFmERERIZj90QvIl4AlgOYqKrlvumlqgtVtZuqdvP19bVfgERERAZm10QvIi6wJPklqrrCuvmMiDSx7m8C4Kw9YyAiIqrN7JboRUQAfA4gXlXn2OxaDeBJ6+snAXxvrxiIiIhqO7stxgPQE8AoAAdFZJ9126sA3gHwrYg8BeA4gIftGAMREVGtZs9V9z8BkFJ232OvfomIiOh3fDIeERGRgTHRExERGRgTPRERkYEx0RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExERGRgTPRERkYEx0RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExERGRgTPRERkYEx0RMRERkYEz0REZGBMdETEREZGBM9ERGRgTHRExERGRgTPRERkYEx0RMRERkYEz0REZGBOVd3AEREjmbHt4dw7kR2dYcBAGgU6IW7h4dUdxhUg3FET0REZGAc0RMRFcMRNBkJR/REREQGxkRPRERkYEz0REREBsZET0REZGBM9ERERAbGRE9ERGRgTPREREQGxkRPRERkYKKq1R1DmUQkDcAxm031AGRWUzhVqaZ/zpoQvyPFWJ2xVHXfjQCcq8L+aqsgVfWt7iCoetWIRF+ciCxU1WerOw57q+mfsybE70gxVmcsVd23iESpareq6o+oNqupU/drqjuAKlLTP2dNiN+RYqzOWBzpeyCiSlQjR/REVLNxRE9UdWrqiJ6IaraF1R0AUW3BET0REZGBcURPRERkYEz0REREBlZrEr2ItBeRT0RkmYiMr+54qHT8WxERVZ5qSfQi4iQie0Vk7S2c4wsROSsiMSXsixSRRBE5LCLTAEBV41V1HIDhALjatxxEpL412SaISLyI3HmT5+Hfim5IROqISLSIDKjuWIiMprpG9C8AiC9ph4j4iYh3sW2tS2j6FYDIEo53AvARgP4AQgE8JiKh1n0DAfwE4N+3EnwtMhfABlVtB6Aziv3N+Lei0pR2cVfShZ3VywC+rdooiWqHKk/0ItIMwAMAPiulSS8A34uIu7X9MwDmFW+kqtsBpJdwfASAw6p6VFWvAPgGwEPWY1arag8AI275gxiciNQF8D8APgcAVb2iqhnFmvFvRaX5CsUu7kq7sBORewHEAThT1UES1QbO1dDnBwCmAvAuaaeqficiLQB8IyLfARgD4L4KnD8AwAmb9ycBdBeR3gCGAHAD8MNNxF3btASQBuBLEekMIBrAC6p6qagB/1ZUGlXdLiLBxTZfvbADABEpurDzAlAHluSfKyI/qKq5CsMlMrQqTfTW+29nVTXa+j/zEqnq363/E/gYQCtVza5INyWfUrcC2FqB89R2zgC6AviTqv4iInMBTAPwF9tG/FtRBZR4YaeqzwGAiPwvgHNM8kSVq6qn7nsCGCgiybBM0/YRkX8VbyQidwPoCGAlgL9WsI+TAAJt3jcDcPqmoq3dTgI4qaq/WN8vgyXxX4N/K6qAEi/srr5Q/UpVb3qBLhGVrEoTvaq+oqrNVDUYwKMA/qOqI23biEg4gE9hmdIbDcBHRN6sQDe7AbQRkRYi4mrtZ3WlfIBaRFVTAZwQkbbWTffAch/1Kv6tqIJ4YUdUDRzxd/SeAB5W1SPWKbwncW0tegCAiPwfgP8CaCsiJ0XkKQBQ1QIAzwHYCMsq8W9VNbbKojeWPwFYIiIHAHQB8Fax/fxbUUXwwo6oGvBZ90RU6awXd70BNIJlNf1fVfVzEbkflgW5TgC+UNWZ1RclUe3ARE9ERGRgjjh1T0RERJWEiZ6IiMjAmOiJiIgMjImeiIjIwJjoiYiIDIyJnoiIyMCY6ImIiAyMiZ6IiMjAmOiJajAR+YeI7BGR26s7FiJyTEz0RDWUiNQB4AdgLIAB1RwOETkoJnqqEUTkfRGZaPN+o4h8ZvN+tohMquQ+syv5fPVFZILN+2ARiSnnsR4isk1EnIq2qeolAE0AbAUwT0RcRWS7iDhXZtxEVLMx0VNNsQtADwAQERMsxVI62OzvAWBnNcRVEfUBTCizVcnGAFihqoVFG0SkISwVBLMAFKrqFQD/BvDIrQZKRMbBRE81xU5YEz0sCT4GQJaINBARNwDtAewVkVUiEi0isSLybNHBIjKr2Gj6byIyWURGisivIrJPRBbYjpht2pbYxjoijxeRT639/SgiHtZ9fxGRBBHZJCL/JyJTALwDoJX1PO9aT+9U0vElGAHg+2Lb/gzgPQCxAEKt21ZZ2xIRAWCipxpCVU8DKBCR5rAk/P8C+AXAnQC6AThgHdGOUdXbrNuet456AeAbXDvSHQ4gyrqtp6p2AVCIYklSRNqX0aYNgI9UtQOADABDRaQbgKEAwgEMscYCANMAHFHVLqr6UmnHF//s1trtLVU12WZbsPV7WAogHr/PbsQA4MI8IrqK9/KoJika1fcAMAdAgPV1JixT+4AluQ+2vg6EJZGeV9W9IuInIk0B+AK4ACAMwG0AdosIAHgAOFusz3vKaPObqu6zvo4GEAzLbYXvVTUXAERkzQ0+U0nHF9cIlosAW28CmKGqKiJXE72qForIFRHxVtWsG/RLRLUEEz3VJEX36cNgGbmeADAZwEUAX4hIbwD3ArhTVXNEZCsAd5vjlwEYBsAflhG+APinqr5ygz7LapNn87oQlgsBqcBnKun44nJh8zlEpAssMwV3ichH1n0Hbdq7AbhcgRiIyMA4dU81yU5YfkaWrqqFqpoOywK3O2GZyq8H4II1ybcDcEex478B8CgsyX4ZLAvXhomIHwCIiI+IBBU7pjxtivsJwIMi4i4iXgAesG7PAuBd0Q+tqhdguZdflOxnAXhQVYNVNRhAZ1hH9NZbFWmqml/RfojImJjoqSY5CMs09s/FtmWq6jkAGwA4i8gBAG8UawdVjYUl0Z5S1RRVjYNlQduP1mM2wfJzNdtjymxTnKruBrAawH4AK2BZC5CpqucB7BSRGJvFeOX1Iywj+D4A6qjqv236OwOgjoj4APgDgB8qeG4iMjBR1eqOgchwRMRLVbNFxBPAdgDPquqeWzhfOIBJqjqqjHYrALyiqok32xcRGQvv0RPZx0IRCYXl/vk/byXJA4B1MeEWEXGy/S29Levq/FVM8kRkiyN6IiIiA+M9eiIiIgNjoiciIjIwJnoiIiIDY6InIiIyMCZ6IiIiA2OiJyIiMjAmeiIiIgNjoiciIjKw/wfeI5uIhJQ0/wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for dat in data:\n",
    "    wav_deets = FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission']))\n",
    "    depth = average_depths['5s'][average_depths['band'] == dat[0]]\n",
    "    #print(depth)\n",
    "    plt.plot([wav_deets[0],wav_deets[1]], [depth,depth], label=dat[0])\n",
    "    \n",
    "plt.xlabel('Wavelength ($\\AA$)')\n",
    "plt.ylabel('Depth')\n",
    "plt.xscale('log')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.title('Depths on {}'.format(FIELD))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### IV.c - Depth vs coverage comparison\n",
    "\n",
    "How best to do this? Colour/intensity plot over area? Percentage coverage vs mean depth?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Depths (5 $\\\\sigma$) vs coverage on SA13')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAEYCAYAAAC0mTTAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOX1P/DPyUYSMkCAQEgCWSB7IGwGBCqIymKVClitKFpFBSk/BcSt9GutgjtoqdXiQi0UFJWlIAKCZZUChi2ELIQlyJJA2MJAAoTk/P64d3AICQmZTJbh83695sXMXc8zk3DyPM+de0RVQURERK7JrbYDICIiIudhoiciInJhTPREREQujImeiIjIhTHRExERuTAmeiIiIhfGRE9EROTCmOiJiIhcGBM9lUtEskXk9mo61hsiMrY6jlXGsTeLSLwzjk1EVN8x0dcjZuItFBGriJwWkQ0iMkpEquVzrM7EXuq4AQAeBjDdbtlqETkvImfNR6YDp3gXwKuOxknOISK9zJ/VfBE5KSI/ishNpbZZLSKnRKRBGfuPEZFkEbkgIp+XWvdvEckRkTMisltEHndyc4jqHSb6+uduVbUACAXwJoAXAHxWuyFV6PcAvlPVwlLLx6iqn/mIduD4iwDcKiKtHDhGnSYiHrUdQ1WISCMA3wL4G4CmAIIB/AXABbttwgD8CoACGFTGYY4AmARgRhnr3gAQpqqNzH0niUiX6msBUf3HRF9PqWq+qi4CcD+AR0QkAQBEJEhE5olInojsF5GnbfuYPfaXRCTN7D39U0S8zXWzALQBsNjsYT9v7tZRRFLM3thc2/bmPi+IyGFzhCFTRG4rJ9yBANZUta0i4ikik834i0REzccO8704D2ALgH5l7PuiiHxTatlfRWTadbYBItJaROab7+0JEfnAXB5r9khPi8guERlUmXObryv6vF4QkRQA50TEwzzmXjPeNBEZbLd9ZxHZZq772vy8JlXmXGW0tcw22cU1obyfi1KiAEBVv1DVYlUtVNXvVTXFbpuHAWwE8DmAR0ofQFXnq+pCACfKWLdLVW1/NKj5aFteu4huSKrKRz15AMgGcHsZy38G8BSMP9y2AHgZgBeACAD7APS32z8VQGsYvasfAUwq7/jm680Agszt0wGMMtdFAzgIIMh8HQagbTlx5wG4qdSy1eby42Ycfa7R7rdgJILWABoCWAlgPoAIu22mAZhaxr6hAAoANDJfuwPIAdD9OtvgDmAHgPfMGLwB9ALgCWAPgD+a73lfAFbz2OWe23xdmc9ru9luH3PZb83Pww3GH3nnALQy9z8A4BkzpiEALto+34rOVaqt5bapop+LMo7VCEaC/heMP/j8y9hmD4DRALoAKALQspxjTQLweRnLPzTfZwWwFYBfbf+u8sFHXXqwR+8ajsD4D/cmAAGq+qqqXlTVfQA+AfA7u20/UNWDqnoSwGQAD1Rw7GmqesTcfjGAjubyYgANAMSJiKeqZqvq3nKO0QRGorD3AoxkEwzgYxgjCVf1xETEAuBpAMPNuM8BmAegqdk+G6t5niuo6gEY//nfYy7qC6BAVTdeZxuSYCS251T1nKqeV9X1MP5g8APwpvme/xfGUPUDFZwbqNznNc1sd6HZnq/Nz6NEVecCyDJj6w7Aw9y+SFXnw0jGuI5z2ZTbplJxlfVzcQVVPQPjDyI1z5cnIotEpCVgzN/D+IPoK1XdAmAvgGFlHas8qjoagAXG8P982E0LEBGH7l1FMICTMP7DDDKHW0+LyGkYvbKWdtsetHt+AEbyupZcu+cFMBIAVHUPgLEAXgFwTES+FJHyjnUKxn/El6nqJlW1quoFVf0XjF79nWXsewuAfaqaZbfMv1RcMI9/upzzz8EvSWqY+fp629AawAFVvVRqeRCAg6paYrfsAIzPpNxzm67384KIPCwi2+22TwDQ3IzjsKpqOftW5lyVbRNQzs9FWVQ1XVV/r6ohZrxBAN43Vz8C4HtVPW6+noMyhu8rosa0wHoAITBGt4jIxERfz4lx9XIwgPUw/mPfr6pN7B4WVbVPoK3tnreBMRpgY58kKqSqc1TV1iNTGEPsZUmBOVd7rcMBkDKWB8D4QwEAICICYDCMHqa9WBhD62X5GkAfEQkx972cbK+jDQcBtCnjorgjAFrLld98aAPgcEXnRuU+r8ufiYiEwugVjwHQTFWbwJiKERhTAsHm+2Nj/1lX5lyVbVOVqWoGjLn4BBHxAXAfgN4ikisiuQDGAUgUkcQqnsIDnKMnugITfT0lIo1E5C4AXwL4t6ruhDFUe8a8gMtHRNxFJEGu/CrTH0QkRESawujRzbVbdxTGcHplzh8tIn3F+DrUeQCFMIbCy/IdgN52+zYRkf4i4m1eYPYgjJ778jL2TQXQWUQ6monhDRjJb67d8RrAmN9dUdbJVTUPxjUB/4SR7NKr0IbNMJLpmyLS0Iy9J4BNMObJnxfjosE+AO6G8bmUe267Y1b0edlraLY9z4z/URg9ZAD4nxn7GPM9/Q2MIf2qnOuabboeIhIjIs+af+hARFrDGOHYCGNKoxhAHIyh/44w/mBbB+MCPdsxPMyL/dwBuNv93LQQkd+JiJ/Znv7msf97vXESubTavkiAj8o/YFwEVQhjPjofxn/ufwDgbrdNEIAvYAytnoLxH+rtdvu/BCANxjD3vwD42u37GxgX9p0GMAFXX5z3Cow/KgCgA4zkYYUxbfAtzIvayoi7OYBD+OWCsgAAP5n7njZjvOMa7Z4Io5eZA6M32LzU+t8CmF/BezccRpJ8zm5Zpdtgbt8GgO3q7+Mw5qkBIB7Gtwryzfd2cEXnvo7P6/ZS2082Yz0OYKp53sfNdV1hXLx3FsZIwnwA/1eZc5URV7ltutbPRRnHCQbwFYzRgHPmv9NhXKS3DMCUMva5z4zRw+74WurxivlztMb8GToDYCeAJ2r795QPPuraQ1Sva7SW6jERyYaRFFbWwrlfB3BMVd+vcOPrP/YmACNUNbW6j12fme/LP1T1n7UdCxHVnnp5Ew6qf1T1j048djdnHbs+EZHeADJh9PYfhDFisaxWgyKiWsdET+Q6omEMk/vB+JravaqaU7shEVFt49A9ERGRC+NV90RERC6sXgzdN2/eXMPCwmo7DCKiemXLli3HVTXAgf1beHh4fArja5zsGNZNJQBSL1269HiXLl2OlbVBvUj0YWFhSE5Oru0wiIjqFRE54Mj+Hh4enwYGBsYGBASccnNz4zxvHVRSUiJ5eXlxubm5n6Ls6o/O+wtNjEpfq0Qk3ax+9Uyp9RPEqEDW3FkxEBGRQxICAgLOMMnXXW5ubhoQEJCPX26edRVn9ugvAXhWVbeahUm2iMgKVU0z7451B4ybsxARUd3kxiRf95mfUbkdd6f16FU1R1W3ms+tMEpZ2opivAfgeVznvdWJiIjo+tTIxRUiEgagE4BNIjIIRpWt8gqQ2PZ5UkSSRSQ5Ly+vBqIkIiJyPU5P9CLiB6N++FgYw/kTAbxc0X6q+rGqdlXVrgEBVb5olIiIasC7yzNbrkw/ekU56pXpRy3vLs8sqxRynTBt2rRmDz/8cJvSy99+++2ADz74oFnp5ZmZmV6RkZHx1XHu+++/P3TLli3e1XGsijg10YuIJ4wkP1tV58MoHxkOYId53/UQAFtFJLBaT7z+fWD/2iuX7V9rLCciomrXsU2TgvFfbY+wJfuV6Uct47/aHtGxTZOC2o7tej3//PN5Y8aMOeHMc8ydO/dAly5dzjvzHDbOvOpeAHwGIF1VpwKAqu5U1RaqGqaqYTAqmnVW1dxqPXlwZ+Dr3/+S7PevNV4Hd67W0xARkeH22JbWqfd13Df+q+0Rf1m8K2j8V9sjpt7Xcd/tsS2tDh/79tvbxsfHx7Zr1y7+3XffbQ4Avr6+nZ544omQuLi42JtvvjnqyJEj5V5cnpSUFL127VpfAMjJyfEIDg5uX3qbL7/8snHHjh1jcnJyPMaPHx/08ssvtwSAdevW+UZHR8d17NgxZurUqS1s2ycnJ3u3b98+NiYmJi4qKipu586dDc6cOePWp0+fdtHR0XGRkZHxn3zyiX9lYnI2Z/boe8Ioz9lXRLabjzudeL5fhN8C/PZzI7n/d7Lx728/N5YTEZFT3B7b0jq0c0jeP3/MbjW0c0hedSR5AJg9e3b2rl270rdv3542ffr0lrm5ue6FhYVunTt3LkhLS0vv2bOn9cUXXwyq6vFnzpzZ5J133glcsWJFVqtWrS7ZrxsxYkTY1KlTf96+fXuG/fK//e1vAaNHjz6akZGRlpKSkh4eHn5x/vz5jQIDA4syMzPTsrKydg0ZMuRMVWOqTs686n69qoqqdlDVjubju1LbhKnqcacEEH4L0HUEsPZt418meSIip1qZftQyb+uhgEd7huXM23oooPScfVW99dZbLaOjo+O6dOkSm5ub67lr1y5vNzc3PP744ycB4LHHHjuxefNmv6oce8OGDZYpU6YErlixIisgIKDYft2JEyfcrVar+69//euztvPY1t18883npkyZ0mrixImBWVlZXn5+ftq5c+fCdevWNXrqqaeCly1b5tesWbPi0uerDa57S8P9a4Hkz4Bbnjf+LT1nT0RE1cY2Jz/1vo77/nx3/BHbML6jyf7bb7+1rFmzxpKcnJyRmZmZFhsbW1hYWHhV7jJmi8vm4eGhxcVGzi0oKLhiwzZt2lw4d+6ce2pq6lUXxqlquccdNWrUyf/85z97fHx8SgYOHBi1aNEiS4cOHS5s3bo1rX379oUTJ04MnjBhQqvrbK5TuGait83J//ZzoO/EX4bxmeyJiJxi+8+nfe3n5G1z9tt/Pu3QPPTp06fdGzduXGyxWEq2bdvmvWPHjoYAUFJSgn/+85/+APD55583S0pKKneaoHXr1hc2b97cEABmz559xbx5SEjIxXnz5u159NFHw5OTk69I9s2bNy/28/MrXr58uZ95nqa2dWlpaV6xsbEX/vSnPx3r16/f6e3bt/tkZ2d7WiyWktGjR58cO3bs0e3bt9fIHHxF6sW97q/b4a1Xzsnb5uwPb+UQPhGRE0zoH3209LLbY1taHZ2nHzp0aP7HH38cEBUVFde2bdvziYmJ5wDAx8enZNeuXT7x8fGBFouleP78+fvKO8aLL7549P7774/48ssvm/3qV7+6at48MTHxwsyZM/fdf//9bRctWrTHft1nn32W/fjjj4f5+PiU9O3b9/K+s2bNavr111838/Dw0ICAgKI33njjyPr16xu+9NJLIW5ubvDw8NAPP/zwmrUGrjUKUZ3qRT36rl27KovaEBFdHxHZoqpdq7r/jh07shMTE51zHZWDfH19OxUUFGyr7TiqKioqKm7RokV7YmJiLlbH8Xbs2NE8MTExrKx1rjl0T0REVEf16NEjMjo6urC6knxFXHPonoiIXFpZvfnhw4e3+emnn664+v6pp546+swzzzj15jfXcscdd7Q9ePBgA/tlkydPPjR06NAa++odEz0REbmEWbNm1bmKqCtWrNhb2zFw6J6IiMiFMdETERG5MCZ6IiIiF8ZET0RE5MKY6ImIyHE/vNYSmUuvvN1t5lILfnit1uvR/+9///Pp2LFjTFRUVFzfvn3bnTx58nLue+mllwLbtGmTEBYWljBv3rxG13vssWPHBi1cuLBa7unvLEz0RETkuJCuBVgwKuJyss9casGCUREI6Vrr9eifeOKJsMmTJx/avXt32qBBg0795S9/CQSALVu2eM+fP79pZmbmrmXLlu0eO3Zsm0uXLlV0uMsuXbqE999//8g999xTLVX6nIWJnoiIHBc90IrB/9iHBaMisPTFICwYFYHB/9iH6IG1Xo8+Ozvbe+DAgWcB4K677jrz7bff+gPAN99802TIkCEnfXx8NCYm5mJoaOiF1atXN8zMzPQKDw+PHzJkSFhUVFTcgAEDIqxWqxsABAcHt58wYUKrLl26RM+YMcN/6NChYbZ77gcHB7cfM2ZMcMeOHWMSEhJi169f79urV6/I1q1bJ7z99tsBtnj+7//+r2VCQkJsVFRU3Lhx465ZXve5555rFR4eHt+jR4/Iu+++O/zll1++7hESJnoiIqoe0QOtSHwgD5s+aoXEB/KqI8kDjtejj4yMLJwzZ04TAPj3v//dNDc31wsADh8+7NW6devLd6cLCgq6ePDgQS/A+ONg1KhRebt3706zWCwl77zzzuVE7e3tXbJly5bMJ5988lTpc7Vu3fri9u3bM7p163b2scceC1u8ePHeTZs2Zbz55ptBADB//vxGe/bs8U5JSUlPT09P2759u+/SpUvLLLG7du1a38WLF/vv3LkzbcmSJXtTUlIaVuX9Y6InIqLqkbnUgh1fBKDbUznY8UXAVXP2VeRoPfoZM2Zkf/TRRwHx8fGxVqvVzdPTUwGjDG1pIqIAEBgYeLFfv37nAGD48OEnNmzYcPn4Dz/88FUJ3ua+++47DQDt27cv6Ny58zl/f/+SoKCgSw0aNCg5fvy4+7JlyxqtXbu2UVxcXFx8fHzc3r17vTMyMq4qkQsAq1ev9hs4cOBpPz8/9ff3L7njjjtOV+oNK4V3xiMiIsfZ5uRtw/URva3VMXxvX4/eYrGUJCUlRV9vPfpOnTqd//HHH7MAICUlpcH333/fBDBK1Np68ABw5MgRr5CQkKKyjmf/2mKxlJR3Lm9vbwUANzc3eHl5Xf5Lws3NDUVFRaKqGDt2bM5zzz1XYbGg6io655I9+hmpM7A5Z/MVyzbnbMaM1Bm1FBERkYs7lOx7RVK3zdkfSq71evSHDx/2AIDi4mL8+c9/bjVixIhjADB06NDT8+fPb1pYWCgZGRle2dnZ3n369DkHADk5OV4rV65sCABz5sxp2qNHj7OOtMNm4MCBZ2bNmtU8Pz/fDQD279/vaYuvtD59+pxdvnx544KCAsnPz3dbuXJlk6qc0yV79AnNEjBhzQS82/tdJLVKwuaczZdfExGRE9z2f1fVo0f0QKuj8/TVUY9+xowZTT/77LMWAHDnnXeeevrpp08AQNeuXc/fc889J6OiouLd3d0xderUAx4eRlqMiIg4P2PGjGajR48ODQ8PvzBhwoQ8R9phM2TIkDO7du3yvummm2IAwNfXt2T27Nn7g4ODr7rcv3fv3gUDBgzIj4uLiw8ODr7QoUOHc40bNy6+3nO6bD16W3K/L/o+fJX51eWkT0R0o2A9+qrJzMz0uuuuuyKzsrJ2OeP41yM/P9+tcePGJVar1e3mm2+O/sc//nGgV69eV31l8Vr16F2yRw8ASa2ScF/0fZieMh0jO4xkkicionrnoYceCs3KyvK5cOGC/O53vztRVpKviMsm+s05m/FV5lcY2WEkvsr8CkmBSUz2REQuwpn16KOjoy/WZG8+NzfXvU+fPtGll69evTpz8eLF+x09vksmevs5+aRWSUgKTLriNRERuZ66WI++MgIDA4szMjLSnHV8l7zqPvVE6hVJPalVEt7t/S5ST6TWcmREREQ1yyV79I8lPHbVsqRWHLonIqIbj0v26ImIiMjARE9EROTCmOiJiMhh07ZOa7n64Oor7m2/+uBqy7St01iPvpYx0RMRkcM6BHQomLh+YoQt2a8+uNoycf3EiA4BHViPvhqUlJSguPi6b4oHgImeiIiqQZ/WfayTe03eN3H9xIg3N78ZNHH9xIjJvSbv69O6D+vRV7EefWZmpldERET8Qw891MasdOdV3rbXwkRPRETVok/rPta7296dNzt9dqu7296dVx1JHrhx69Hb4nj00UdPpKenp0VFRV0sb7trYaInIqJqsfrgasvivYsDHox9MGfx3sUBpefsq+pGrUcPAK1atbp42223navwTboGl/wePRER1SzbnLxtuL57q+7W6hi+v5Hr0QNGdbvKbHct7NETEZHDUvJSfO2Tum3OPiUvhfXo7VxPPfrq4rSDi0hrADMBBAIoAfCxqv5VRF4D8Btz2TEAv1fVI86Kg4iInO/pzk9fVY++T+s+Vkfn6W/kevTVxWn16EWkFYBWqrpVRCwAtgC4B8AhVT1jbvM0gDhVHXWtY1WlHj0R0Y2O9eirpi7Vo6+sa9Wjd9rQvarmqOpW87kVQDqAYFuSNzUE4Jy/NIiIiKhmLsYTkTAAnQBsMl9PBvAwgHwAt9ZEDERE5DpulHr0gYGBVbtLjh2nJ3oR8QMwD8BYW29eVScCmCgiLwEYA+DPZez3JIAnAaBNmzbODpOIiOo51qMvm1OvuhcRTxhJfraqzi9jkzkAhpa1r6p+rKpdVbVrQEBAWZsQERFRBZyW6MX40uFnANJVdard8ki7zQYByHBWDERERDc6Zw7d9wQwHMBOEdluLvsjgBEiEg3j63UHAFzzinsiIiKqOqclelVdD6CsWxV956xzEhER0ZV4ZzwiInLYsfffb2ldteqKe9tbV62yHHv//VqvR1+T6mJ9eiZ6IiJymE9iYsGRF16MsCV766pVliMvvBjhk5hY6/Xoa1JN1ae/Hkz0RETkMMutt1qD3npz35EXXozIff31oCMvvBgR9Nab+yy33lrr9eiTkpKiR4wY0bpr167RERER8WvWrPHt169f29DQ0ISnn376cnnbV155pWVkZGR8ZGRk/KuvvtoCAM6cOePWp0+fdtHR0XGRkZHxn3zyiT8ATJgwoVVCQkJsZGRk/AMPPBBaUmLUnrGvT79mzRrfTp06xURHR8e1b98+9tSpU2Xm3Pvvvz80JiYmLiYmJs7f3z/x2WefbeXoe2aPiZ6IiKqF5dZbrY3v+U3eqZmzWjW+5zd51ZHkAcfr0QOAl5dXSXJycuajjz6a99vf/rbdJ5988nNGRsauuXPnNs/NzXVft26d75w5c5pt2bIlPTk5OX3mzJkBP/74o8/8+fMbBQYGFmVmZqZlZWXtGjJkyBkAeO65546lpqamZ2Vl7SosLHT78ssvG9uf7/z58/Lggw+2ff/993/OzMxMW7NmTaafn1+Zlejmzp17ICMjI23RokV7mjRpcmnkyJHXdYOfijDRExFRtbCuWmXJX/ifAP+Hh+fkL/xPQOk5+6pytB49AAwePPg0ACQmJha2a9euMDQ0tMjHx0dbt259Yd++fV6rV6/2u/POO083atSopHHjxiW//vWvT61atcrSuXPnwnXr1jV66qmngpctW+bXrFmzYgBYunSppUOHDjFRUVFxGzZssKSmpvrYny8lJcW7RYsWRb179y4AgKZNm5Z4enqWG19BQYEMHTq07XvvvfdzVFTURQffsisw0RMRkcNsc/JBb725L/CPfzxiG8Z3NNnb16PPzMxMi42NLbzeevTAlXXiGzRocEWd+EuXLpVb4K1Dhw4Xtm7dmta+ffvCiRMnBk+YMKFVQUGBPPvss6Hz58/fu3v37rSHHnro+Pnz56+ISVUhIpWu5TJ8+PDQu++++5Qz5veZ6ImIyGGFO3b42s/J2+bsC3fsqPV69JXRt2/fs999910Tq9XqdubMGbfvvvvO/9Zbb7VmZ2d7WiyWktGjR58cO3bs0e3bt/sWFBS4AUBgYOCl/Px8t8WLF/uXPl5iYuL5o0ePeq1Zs8YXAE6dOuVWVFRU5rnfeOONgLNnz7q//vrruY60oTw1UtSGiIhcW4uxY6+qR2+59Varo/P01VGPvjJ69epVMGzYsBOdO3eOBYDhw4fn9ezZs3DevHmNXnrppRA3Nzd4eHjohx9+eKB58+bFDz74YF5cXFx8SEjIRVtM9ry9vXX27Nl7n3766Tbnz5938/b2Llm7du3uxo0bXzVP/8EHHwR6enpqTExMHAA89thjec8//3yeI+2x57R69NWJ9eiJiK4f69HfOGqlHj0RERHVPg7dExFRvePMevTOMm/evEYTJ04MsV/WunXrCytWrNjrzPMy0RMRkUuo6/Xohw4dembo0KFOqztfHg7dExERuTAmeiIiIhfGRE9EROTCmOiJiMhhG/+zt+X+lONX3AVvf8pxy8b/7L2hytTWRUz0RETksJbhjQt++Dwtwpbs96cct/zweVpEy/DGTilTm5SUFL127dqr7rrXu3fvdsePH3cvvXz8+PFBL7/8ssN/dGRnZ3sOGDAgwtHj1CRedU9ERA4L79Dcetvv4/b98HlaRHT3wLzMjbkBt/0+bl94h+Y1Wpt9zZo1e5x5/LCwsKJly5Y5dBe+msYePRERVYvwDs2t0d0D81L+e6hVdPfAvOpI8pmZmV6RkZHxttcvv/xyy/Hjx18uSVtcXIwhQ4aE2erKBwcHt8/JyfEAgBdeeCEwLCwsoUePHlFZWVkNbPtMmjSpRdu2beOjoqLi7rrrrggAWLJkiZ+tJnxsbGxcebXjS8dTH7BHT0RE1WJ/ynFL5sbcgA59Q3IyN+YGhMQ0tTqzR19UVCT33HNPeFxcXOFbb711RUGYdevW+S5YsKDpzp0704qKitCxY8e4Tp06FQDAtGnTAg8cOLDTx8dHbcP8U6ZMCZw2bdqBfv36ncvPz3fz9fUts3Z8fcQePREROcw2J3/b7+P2/eq+qCO2YfzSF+hVp9GjR4eWleQBYNWqVX533nnnaYvFUtK0adOSfv36nbati46OLhw8eHD4hx9+2NTT01MBoHv37mcnTJjQetKkSS2OHz/ufq3a8fUNEz0RETns6P58X/s5educ/dH9+Q6VqfXw8NCSkl861/Z137t27Xp23bp1jQoKCsosRl9ejfpVq1Zl/eEPf8jbsmVLw8TExLiioiK8/vrruZ9++umBwsJCtx49esRu27bN25G46xImeiIiclj337Q9WnqYPrxDc2v337S9qnzt9QgJCbl08uRJj9zcXPfCwkJZvnx5Y9u6kSNHHu/Xr1/+XXfd1bZ0rfe+ffueXbJkSZOzZ8/KqVOn3FasWNEEMOb09+7d63X33XdbP/zww0NWq9U9Pz/ffdeuXQ2SkpIKJ0+enNu+fftzqampLpPoOUdPRER1VoMGDfTZZ5/NSUpKig0JCbnQrl278/brX3nllaPjxo1zHzJkSPjChQv325b36tWrYPDgwScTEhLig4ODLyQlJZ0FgEuXLsmwYcPCrVaru6rKyJEjjzZv3rz42WefDdqwYUMjNzczfAYQAAAdQUlEQVQ3jYqKKrz33nvzy4tJROp+fXc7rEdPROSiXLkefW1Zt26d7/jx41v/9NNPmbUdiz3WoyciInLQ2rVrfYcPHx4xZswYh6YjahqH7omIiErZvHmzz8MPPxxuv8zLy6skOzs7tbZiqiomeiIiolKSkpIKMzIyarx2vDNw6J6IiMiFuWSiP/Hppzi3cdMVy85t3IQTn35aSxERERHVDpdM9N4J7XF43LjLyf7cxk04PG4cvBPa13JkRERENcslE33D7t0Q/N57ODxuHPKmTcPhceMQ/N57aNi9W22HRkTkktZ/ObPl3i2br7jd7d4tmy3rv5zJevS1zCUTPWAke/8HfofjH34E/wd+xyRPRORErSJjCpb+fUqELdnv3bLZsvTvUyJaRcY4XI/+tddeaxEZGRnfrl27+FdffbUFABw9etS9R48ekaGhoQk9evSIzMvLu6oGfUU6deoU42hs9YHLJvpzGzfh1Bdfovnop3Dqiy+vmrMnIqLq07ZLknXgH57dt/TvUyJWff5x0NK/T4kY+Idn97XtkuRQ9bqffvrJe+bMmQFbt25NT09P37Vs2bImO3fubPDnP/+5VZ8+fawHDhxI7dOnj/Xll18OrOwxL126BADYtm1bhiOx1Rcumehtc/LB772HgKefvjyMz2RPROQ8bbskWeNvuS1v69JFreJvuS3P0SQPADt37vTp3LnzWYvFUuLp6YmePXta586d22TZsmVNRo4ceQIARo4ceWLp0qX+ADB+/Pige+65J7x79+5RoaGhCVOmTGkOAN9++62lW7duUXfffXd4dHR0PAD4+vp2sq276aabou+8886IsLCwhNGjRwd/9NFHTdu3bx8bFRUVt2vXrgYAcOTIEY/+/fu3TUhIiE1ISIj9/vvvG5YX95EjRzx69OgRGRcXFzts2LDQoKCg9jk5ObXylXanJXoRaS0iq0QkXUR2icgz5vJ3RCRDRFJEZIGINKnuc59P3XnFnLxtzv586s7qPhUREZn2btls2bX2h4DOAwfl7Fr7Q0DpOfuq6NixY+GmTZssubm57lar1W3FihWNDx486HXixAmP0NDQIgAIDQ0tOnny5OUkmp6e7rNy5cqsjRs3ZrzzzjtB2dnZngCQkpLS8J133jm8d+/eXaXPk5GR4fPRRx8dTE9P3/XNN9802717t/fOnTvThw8ffnzKlCktAGDkyJGtx48ffzQ1NTV9wYIFe0eNGhVWXtwvvvhiUO/eva1paWnpQ4YMOZWTk+Pl6HtRVc786+ISgGdVdauIWABsEZEVAFYAeElVL4nIWwBeAvBCdZ642eOPX7WsYfdunKcnInIS25y8bbi+TfuO1uoYvu/cufP5Z555Jrdv375Rvr6+JXFxcQUeHtdOXQMHDjzt5+enfn5+l26++eYz69ata+jv71/coUOHczExMRfL2qd9+/bnbH84tGnT5sLAgQPzASAxMbFwzZo1FgD48ccfG2VlZfnY9jl79qz7qVOn3Pz9/UtKH2/z5s1+Cxcu3AMA995775lGjRoVV/U9cJTTevSqmqOqW83nVgDpAIJV9XtVvWRuthFAiLNiICKimpGTleFrn9Rtc/Y5WRkO1aMHgHHjxh1PS0tLT05OzmzatGlxZGTk+WbNml06cOCAJwAcOHDAs2nTpra8clUdettrX1/fqxKyTYMGDS5XeHNzc4O3t7fanhcXFwsAqCqSk5PTMzIy0jIyMtKOHTuWUlaSt21bV9TIHL2IhAHoBKD0JPljAJaWs8+TIpIsIsl5eXnODZCIiBzS63cPHy3dc2/bJcna63cPO1wA5vDhwx4AkJWV5bVkyZImI0aMONm/f//T06dPbwYA06dPbzZgwIDTtu2XLl3apKCgQHJzc903btxo6dWr1zlHYwCAXr16nXnrrbda2F5v2LDBp7xtk5KSzs6aNaspAMyfP7/RmTNnrvtbAdXF6RcGiIgfgHkAxqrqGbvlE2EM788uaz9V/RjAx4BRptbZcRIRUd00aNCgtqdPn/bw8PDQ999//+eAgIDiv/zlLzmDBw9uGxoa2jwoKOjiwoUL99q279Sp07nbbrst8siRI14TJkzICQsLK0pNTfV2NI6PP/744OOPP94mKioqrri4WLp162bt0aPHz2Vt++abbx659957I+Li4vxvvvnmswEBAUVNmjSpleF7p9ajFxFPAN8CWK6qU+2WPwJgFIDbVLXC71iyHj0R0fW7EevRjx8/PsjPz6/41VdfrdVSsoWFheLh4aGenp5YuXJlwzFjxoQ6s0jOterRO61HL8akyGcA0ksl+QEwLr7rXZkkT0REVN/s2bPH67777mtbUlICT09PnT59enZtxeLMofueAIYD2Cki281lfwQwDUADACvMCyQ2quooJ8ZBREQ3iKlTpx6pyfP99a9/bfbRRx9dcZvfm2666eysWbN+Tk9PrxNlbiuV6EWkAYChAMLs91HVV8vbR1XXA5AyVn13fSESERHVTc8888yJZ5555kRtx3Etle3R/wdAPoAtAC44LxwiIiKqTpVN9CGqOsCpkRAREVG1q+z36DeICIu5ExER1TPX7NGLyE4Aam73qIjsgzF0LwBUVTs4P0QiIqrr8pdnt/RqYynwiW12+aY5heknLBd/tvo27h9Wq191u9FV1KO/C8DdAAYCaAegn/natpyIiAhebSwFJ7/aHVGYfsICGEn+5Fe7I7zaWFiPvpZdM9Gr6gFVPQBgku25/bKaCZGIiOo6n9hm1qb3Re07+dXuiNOL9wad/Gp3RNP7ovbZ9/CrwhXq0RcVFdXEacpV2Tn6ePsXIuIOoEv1h0NERPWVT2wza8POLfLO/nikVcPOLfIcTfJA/a1HP378+KAHHnggtGfPnpFDhgwJd/R9cERFc/QvwbjJjY+InMEv34u/CPM+9ERERIAxXH9u67EAv55BOee2Hgto0K6J1dFk37Fjx8JXX301ODc3171hw4a6YsWKxomJiecqqke/ZcuWdKvV6t6pU6e4oUOH5gNGPfpt27btKqtUbUZGhs8333yzr0WLFpdCQ0PbN2jQ4PjOnTvTX3vttRZTpkxpMWPGjIO2evT9+/c/m5WV5dW/f//Iffv2XVXb3iYlJcV306ZNGX5+frVar+WaiV5V3wDwhoi8oaov1VBMRERUz9jm5G3D9Q3aNbFWx/B9fa1HDwADBgw4XdtJHqj89+j/KCJDAPSCcRX+OlVd6LywiIioPrn4s9XXPqnb5uwv/mz1dbRXP27cuOPjxo07DgBjxowJDgkJuWirRx8aGlpU0/XoK5u8GzZsWO75alJl5+j/DqPa3E4AqQBGicjfnRYVERHVK437hx0tndB9YptZq+OrdfWxHn1dUtkefW8ACWrWtBWRf8FI+kRERE5VH+vR1yWVqkcvIvMBjDO/VgcRCQXwpqo+4OT4ALAePRFRVbAe/Y2jOurRNwOQLiKbzdc3AfifiCwCAFUd5HCUREREVO0qm+hfdmoURERE1aAu1aOvyTiupVKJXlXXmMP1kaq6UkR8AHioqsM3QyAiIqqv6kM9+kpddS8iTwD4BsB0c1EIAH69joiIqI6r7Nfr/gCgJ4AzAKCqWQBaXHMPIiIiqnWVTfQXVPXy3YRExAPGjXOIiIioDqtsol8jIrZ73t8B4GsAi50XFhER1Sc//PBDy8zMTIv9sszMTMsPP/zQsrx96hpXLVtb2UT/IoA8GDfJGQngOwB/clZQRERUv4SEhBQsWLAgwpbsMzMzLQsWLIgICQlxuB59TampsrU1rVKJXlVLYFx8N1pV71XVT7Qyd9ohIqIbQnR0tHXw4MH7FixYELF06dKgBQsWRAwePHhfdHS0Q9/OyszM9AoPD4+///77QyMjI+MHDRoUvnDhQkvnzp1jQkNDE1atWuV79OhR99tvv71tVFRUXGJiYsymTZt8AGDJkiV+MTExcTExMXGxsbFxp06dcsvPz3e7+eabo+Li4mKjoqLi/v3vfzexnctWthYA/vSnP7WMioqKi46Ojhs9enRwWbFlZ2d72o4fExMT5+7u3mX37t1ejrTXGSoqUysA/gxgDIwStSIixQD+pqqv1kB8RERUT0RHR1sTExPzNm3a1Kpbt245jiZ5m4MHD3rPnTt3X5cuXQ506NAhdvbs2c2Sk5Mz5syZ02Ty5MmtgoODLyYmJhasXLly76JFiyyPPPJIeEZGRtqUKVMCp02bdqBfv37n8vPz3WxFbZYsWbKnadOmJTk5OR7dunWLGTZs2Gk3t1/6vV999VWjJUuW+G/ZsiXDYrGUHD161L2suMLCwooyMjLSAOCNN94IWLdunSUqKqrM6ni1qaIe/VgYV9vfpKrNVLUpgG4AeorIOKdHR0RE9UZmZqZlx44dAd26dcvZsWNHQOk5+6oKDg6+kJSUVOju7o6oqKjCvn37nnFzc0Pnzp0LDh061GDz5s2WESNGnACAQYMGWU+fPu1x4sQJ9+7du5+dMGFC60mTJrU4fvy4u6enJ0pKSmTs2LEhUVFRcbfeemvUsWPHvA4dOnRFp3fFihWNHnrooeMWi6UEAFq2bFl8rfi+//77hjNnzgz44osvsqujvdWtokT/MIAHVHW/bYGq7gPwkLmOiIjo8pz84MGD9w0cOPCIbRi/OpK9l5dXmSVk3d3dUVxcXGbNFhHR119/PffTTz89UFhY6NajR4/Ybdu2eU+fPr3piRMnPHbu3JmekZGR1qxZs6LCwsIrcqGqXlXqtjwHDhzwHDlyZNjcuXP3Nm7cuE6UpS2tokTvqapXFTRQ1TwAns4JiYiI6ptDhw752s/J2+bsDx065Ovsc3fv3t36z3/+sxkAfPvttxZ/f/9LTZs2Ldm1a1eDpKSkwsmTJ+e2b9/+XGpqqnd+fr578+bNixo0aKCLFy+2HDly5Ko59QEDBpyZNWtWc6vV6gYA5Q3dX7hwQYYMGRLx2muvHe7QocMF57ay6iq6Be615hrq3DwEERHVjttuu+2qanHR0dHW6pqnv5a33nrryLBhw8KioqLifHx8Sj7//PP9APD222+32LBhQyM3NzeNiooqvPfee/NPnz7tPnDgwHYJCQmx8fHxBeHh4edLH+/ee+89s3XrVt+OHTvGenp66u23357/wQcfHC693cqVKxumpqY2nDRpUtCkSZOCAGDZsmVZYWFhRc5u8/W4Zpla88K7c2WtAuCtqjXSq2eZWiKi63cjlqm9UVW5TK2qljlcQURERPVDZcvUEhER3bCGDx/e5qeffvKzX/bUU08dreuV6wAmeiIiogrVpfry16uyt8AlIqIbT0lJSUnlvmdGtcb8jMr9ah8TPRERlSc1Ly+vMZN93VVSUiJ5eXmNAaSWtw2H7omIqEyXLl16PDc399Pc3NwEsGNYV5UASL106dLj5W3gtEQvIq0BzAQQaAbysar+VUR+C+AVALEAklSV35sjIqqDunTpcgzAoNqOgxzjzB79JQDPqupWEbEA2CIiK2AMLwwBMN2J5yYiIiI4MdGrag6AHPO5VUTSAQSr6goAlb6PMBEREVVdjcy5iEgYgE4ANl3HPk+KSLKIJOfl5TkrNCIiIpfm9EQvIn4A5gEYq6pnKrufqn6sql1VtWtAQIDzAiQiInJhTk30IuIJI8nPVtX5zjwXERERXc1piV6MSfjPAKSr6lRnnYeIiIjK58yr7nsCGA5gp4hsN5f9EUADAH8DEABgiYhsV9X+ToyDiIjohuXMq+7XwyhnW5YFzjovERER/YJ3OiIiInJhTPREREQujImeiIjIhTHRExERuTAmeiIiIhfGRE9EROTCmOiJiIhcGBM9ERGRC2OiJyIicmFM9ERERC6MiZ6IiMiFMdETERG5MCZ6IiIiF8ZET0RE5MKY6ImIiFwYEz0REZELY6InIiJyYUz0RERELoyJnoiIyIUx0RMREbkwJnoiIiIXxkRPRETkwpjoiYiIXBgTPRERkQtjoiciInJhTPREREQujImeiIjIhTHRExERuTAmeiIiIhfGRE9EROTCmOiJiIhcGBM9ERGRC2OiJyKqY/6xZi827D1+xbINe4/jH2v21lJEVJ8x0RMR1TEdQhpjzJxtl5P9hr3HMWbONnQIaVzLkVF95LRELyKtRWSViKSLyC4RecZc3lREVohIlvmvv7NiICKqj3q0bY4PhnXCmDnbMPX7TIyZsw0fDOuEHm2b13ZoVA85s0d/CcCzqhoLoDuAP4hIHIAXAfygqpEAfjBfExGRnR5tm+Ohbm0w7b978FC3NkzyVGVOS/SqmqOqW83nVgDpAIIB/AbAv8zN/gXgHmfFQERUX23Yexz/3vQznu7bDv/e9PNVc/ZElVUjc/QiEgagE4BNAFqqag5g/DEAoEU5+zwpIskikpyXl1cTYRIR1Qm2OfkPhnXC+H7Rl4fxmeypKpye6EXED8A8AGNV9Uxl91PVj1W1q6p2DQgIcF6ARER1TMqh/Cvm5G1z9imH8ms5MqqPPJx5cBHxhJHkZ6vqfHPxURFppao5ItIKwDFnxkBEVN+M6t32qmU92jbnPD1ViTOvuhcAnwFIV9WpdqsWAXjEfP4IgP84KwYiIqIbnTN79D0BDAewU0S2m8v+COBNAF+JyAgAPwP4rRNjICIiuqE5LdGr6noAUs7q25x1XiIiIvoF74xHRETkwpjoiYiIXBgTPRERkQtjoiciInJhTPREREQujImeiIjIhTHRExERuTAmeiIiIhfGRE9EROTCmOiJiIhcGBM9ERGRC2OiJyIicmFM9ERERC6MiZ6IiMiFMdETERG5MCZ6IiIiF8ZET0RE5MKY6ImIiFwYEz0REZELY6InIiJyYUz0RERELoyJnoiojtm6/AAOZZ66YtmhzFPYuvxALUVE9RkTPRFRHdMirBGWf5J6OdkfyjyF5Z+kokVYo1qOjOojj9oOgIiIrhQS7Y/+TyRg+SepSLglGKlrD6P/EwkIifav7dCoHmKPnoioDgqJ9kfCLcFI/i4bCbcEM8lTlTHRExHVQYcyTyF17WF0vTMMqWsPXzVnT1RZTPRERHWMbU6+/xMJ6DYo4vIwPpM9VQUTPRFRHXMs+8wVc/K2Oftj2WdqOTKqj3gxHhFRHdO5f+hVy0Ki/TlPT1XCHj0REZELY6InIiJyYUz0RERELoyJnoiIyIUx0RMREbkwUdXajqFCIpIHoKrVHJoDOF6N4dQmtqXucZV2AGxLXeVIW0JVNaA6g6H6p14kekeISLKqdq3tOKoD21L3uEo7ALalrnKltlDt4NA9ERGRC2OiJyIicmE3QqL/uLYDqEZsS93jKu0A2Ja6ypXaQrXA5efoiYiIbmQ3Qo+eiIjohsVET0RE5MLqdaIXkQEikikie0TkxTLWtxGRVSKyTURSROROc7mniPxLRHaKSLqIvFTz0V8RZ0XtCBWRH8w2rBaRELt1j4hIlvl4pGYjv1pV2yIiHUXkfyKyy1x3f81Hf1WsVf5czPWNROSwiHxQc1GXzcGfsTYi8r35u5ImImE1GXupOB1px9vmz1e6iEwTEanZ6K+KdYaIHBOR1HLWixnnHrM9ne3W1anfe6rjVLVePgC4A9gLIAKAF4AdAOJKbfMxgKfM53EAss3nwwB8aT73BZANIKwOt+NrAI+Yz/sCmGU+bwpgn/mvv/ncv45/JuW1JQpApPk8CEAOgCb1sS126/8KYA6AD2qrHdXRFgCrAdxhPvcD4Fvf2gGgB4AfzWO4A/gfgD61/LncAqAzgNRy1t8JYCkAAdAdwCZzeZ36veej7j/qc48+CcAeVd2nqhcBfAngN6W2UQCNzOeNARyxW95QRDwA+AC4COCM80MuU2XaEQfgB/P5Krv1/QGsUNWTqnoKwAoAA2og5vJUuS2qultVs8znRwAcA1Cbd/Ry5HOBiHQB0BLA9zUQa0Wq3BYRiQPgoaorAEBVz6pqQc2EfRVHPhMF4A3jD4QGADwBHHV6xNegqmsBnLzGJr8BMFMNGwE0EZFWqHu/91TH1edEHwzgoN3rQ+Yye68AeEhEDgH4DsD/M5d/A+AcjF7jzwDeVdVr/cI5U2XasQPAUPP5YAAWEWlWyX1rkiNtuUxEkmD8h7zXSXFWRpXbIiJuAKYAeM7pUVaOI59LFIDTIjLfnAJ7R0TcnR5x2arcDlX9H4zEn2M+lqtqupPjdVR57a1rv/dUx9XnRF/W/Frp7wo+AOBzVQ2BMQw2y/xPOAlAMYwh4nAAz4pIhDODvYbKtGMCgN4isg1AbwCHAVyq5L41yZG2GAcweiyzADyqqiXOCrQSHGnLaADfqepB1A2OtMUDwK/M9TfBGDb/vdMivbYqt0NE2gGIBRACIyn2FZFbnBlsNSivvXXt957qOI/aDsABhwC0tnsdgl+G5m1GwBzSUtX/iYg3jAIRwwAsU9UiAMdE5EcAXWHMddW0CtthDmUPAQAR8QMwVFXzzZGKPqX2Xe3MYCtQ5baYrxsBWALgT+ZQZW1y5HO5GcCvRGQ0jDltLxE5q6pXXTxWQxz9GdumqvvMdQthzBd/VhOBl+JIO54EsFFVz5rrlsJox9qaCLyKymtvXfu9pzquPvfofwIQKSLhIuIF4HcAFpXa5mcAtwGAiMTCmKPLM5f3Na9qbQjjFz6jxiK/UoXtEJHm5kgEALwEYIb5fDmAfiLiLyL+APqZy2pLldtibr8Axpzk1zUYc3mq3BZVfVBV26hqGIwe5sxaTPKAYz9jPwHwFxHb9RJ9AaTVQMxlcaQdP8Po6XuIiCeM3n5dH7pfBOBh8/+p7gDyVTUHde/3nuq62r4a0JEHjOH43TDmcieay14FMMh8HgfjStsdALYD6Gcu94Nxde4uGP9pPVfH23EvgCxzm08BNLDb9zEAe8zHo/XgMymzLQAeAlBkfk62R8f62JZSx/g9avmq+2r4GbsDQAqAnQA+B+BV39oB40r76TCSexqAqXXgM/kCxvUCRTB66SMAjAIwylwvAP5utnUngK52+9ap33s+6vaDt8AlIiJyYfV56J6IiIgqwERPRETkwpjoiYiIXBgTPRERkQtjoiciInJhTPR0QxCRQBH5UkT2mhXYvhORqNqOi4jI2ZjoyeWZ5UgXAFitqm1VNQ7AH2EUnanuc9XWfeCJiMrERE83glsBFKnqP2wLVHU7gPVmkZZUEdkpIvcDgIjMFZE7bduKyOciMlRE3M3tfzLrg4801/cRkVUiMgfGjU0gIgtFZIsY9c+ftDvWCBHZLUat9E/ErFUvIgEiMs889k8i0rNG3hkicnn1+V73RJWVAGBLGcuHAOgIIBFGDYSfRGQtjPKn9wP4zrzV6m0AnoJx57J8Vb1JRBoA+FFEbGVokwAkqOp+8/VjqnpSRHzM486DUR71/2DUILcC+C+MuzYCRu3691R1vYi0gXFL09jqewuI6EbFRE83sl4AvlDVYgBHRWQNjAptSwFMM5P5AABrVbVQRPoB6CAi95r7NwYQCeAigM12SR4AnhaRwebz1uZ2gQDWqFkSWUS+hlEGFgBuBxBnzDIAABqJiEVVrdXfbCK6kTDR041gF4x7oJdWVrlPqOp5EVkNoD+Mnv0Xdtv/P1W9ooCIiPQBcK7U69sB3KyqBeaxvMs7n8nN3L6w4uYQEVUe5+jpRvBfAA1E5AnbAhG5CcApAPebc+8BAG4BsNnc5EsAj8KoxW5L7MsBPGVWP4OIRJnVD0trDOCUmeRjYFRHhHns3mbVMQ8AQ+32+R7AGLv4OjrUYiIiE3v05PJUVc1h9PdF5EUA5wFkAxgLo5LhDgAK4HlVzTV3+x7ATACLVPWiuexTAGEAtppX8ucBuKeMUy4DMEpEUgBkAthoxnFYRF4HsAlGXfE0APnmPk8D+Lu5jweMOumjquUNIKIbGqvXEdUgEfFT1bNmj34BgBmquqC24yIi18Whe6Ka9YqIbAeQCmA/gIW1HA8RuTj26ImIiFwYe/REREQujImeiIjIhTHRExERuTAmeiIiIhfGRE9EROTC/j/5hEpqdhmGggAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for dat in data:\n",
    "    wav_deets = FWHM(np.array(dat[1]['Wavelength']), np.array(dat[1]['Transmission']))\n",
    "    depth = average_depths['5s'][average_depths['band'] == dat[0]]\n",
    "    #print(depth)\n",
    "    coverage = np.sum(~np.isnan(depths['ferr_{}_mean'.format(dat[0])]))/len(depths)\n",
    "    plt.plot(coverage, depth, 'x', label=dat[0])\n",
    "    \n",
    "plt.xlabel('Coverage')\n",
    "plt.ylabel('Depth')\n",
    "#plt.xscale('log')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.title('Depths (5 $\\sigma$) vs coverage on {}'.format(FIELD))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (herschelhelp_internal)",
   "language": "python",
   "name": "helpint"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}