{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Some simple examples of interacting with the HELP database\n",
    "\n",
    "In this notebook we will demostrate some ways you can use HELP data. The simplest way to use HELP data is probably to download a fits table from HeDaM and use the table with TOPCAT or astropy. In this example we will look at the ELAIS-N1 field. You can download the file (10Gb) here:\n",
    "\n",
    "http://hedam.lam.fr/HELP/dataproducts/dmu32/dmu32_ELAIS-N1/data/ELAIS-N1_20171016.fits\n",
    "\n",
    "That is a large fiel that would take some time to download. Perhaps you only want to look a small region with the field or only want a subset of columns for you analysis. This is where the Virtual Obersvatory (VO) can be very useful. It allows you to query the database for only the specific information you are interested in.\n",
    "\n",
    "In order to do this we will use PyVO...\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pyvo as vo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If that failed you probably haven't installed PyVO. If you install the herschelhelp_python package and run this notebook with that kernel then you will have PyVO installed.\n",
    "\n",
    "https://github.com/H-E-L-P/herschelhelp_python\n",
    "\n",
    "Alternatively you can install it individually using your preferred method for installing Python packages.\n",
    "\n",
    "Now in order to query the VO we need to write a query. As an example lets get all the ELAIS-N1 objects with a redshift greater than 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# SELECT - columns you want to select\n",
    "# FROM -  which table you want to \n",
    "# WHERE - the criteria the search results must meet\n",
    "\n",
    "example_query=\"\"\"\n",
    "SELECT ra, dec, redshift, cigale_dustlumin\n",
    "FROM herschelhelp.main\n",
    "WHERE (\n",
    "herschelhelp.main.field='ELAIS-N1'\n",
    "AND herschelhelp.main.redshift>4\n",
    ")\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You could copy paste the above query into the html query page:\n",
    "\n",
    "https://herschel-vos.phys.sussex.ac.uk/__system__/adql/query/form\n",
    "\n",
    "But here we will execute it using PyVO:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Then we establish the VO connection to our database\n",
    "service = vo.dal.TAPService(\"https://herschel-vos.phys.sussex.ac.uk/__system__/tap/run/tap\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: W27: None:2:1176: W27: COOSYS deprecated in VOTable 1.2 [astropy.io.votable.tree]\n",
      "WARNING: W06: None:2:2476: W06: Invalid UCD 'phys.dust;phys.luminosity': Unknown word 'phys.dust' [astropy.io.votable.tree]\n"
     ]
    }
   ],
   "source": [
    "#Then we execute the query\n",
    "resultset = service.search(example_query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table masked=True length=85572</i>\n",
       "<table id=\"table4924444968\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>ra</th><th>dec</th><th>redshift</th><th>cigale_dustlumin</th></tr></thead>\n",
       "<thead><tr><th>deg</th><th>deg</th><th></th><th>W</th></tr></thead>\n",
       "<thead><tr><th>float64</th><th>float64</th><th>float32</th><th>float64</th></tr></thead>\n",
       "<tr><td>246.555662513057</td><td>55.003551059499301</td><td>4.5766001</td><td>--</td></tr>\n",
       "<tr><td>246.546474785057</td><td>54.997026013499301</td><td>4.8196998</td><td>--</td></tr>\n",
       "<tr><td>246.397152569057</td><td>54.927381462499298</td><td>4.2312999</td><td>--</td></tr>\n",
       "<tr><td>246.38414028705699</td><td>54.952376668499298</td><td>4.8771</td><td>--</td></tr>\n",
       "<tr><td>246.391144082057</td><td>54.9396872434993</td><td>4.0481</td><td>--</td></tr>\n",
       "<tr><td>246.44518433205701</td><td>55.000396006499301</td><td>5.0402002</td><td>--</td></tr>\n",
       "<tr><td>246.41164064105701</td><td>55.0204406554993</td><td>5.0925999</td><td>--</td></tr>\n",
       "<tr><td>246.413862278057</td><td>55.0054131774993</td><td>4.6952</td><td>--</td></tr>\n",
       "<tr><td>246.401368546057</td><td>55.012446681499299</td><td>4.9246001</td><td>--</td></tr>\n",
       "<tr><td>...</td><td>...</td><td>...</td><td>...</td></tr>\n",
       "<tr><td>239.10977223905701</td><td>54.642360868499303</td><td>4.6788998</td><td>--</td></tr>\n",
       "<tr><td>239.100740939057</td><td>54.648426408499297</td><td>4.5992999</td><td>--</td></tr>\n",
       "<tr><td>239.11018292746701</td><td>54.7110898655787</td><td>4.7849998</td><td>--</td></tr>\n",
       "<tr><td>239.09112570705699</td><td>54.715674113499297</td><td>4.5995002</td><td>--</td></tr>\n",
       "<tr><td>239.083845867924</td><td>54.729565023113899</td><td>5.2941999</td><td>--</td></tr>\n",
       "<tr><td>239.11001283105699</td><td>55.001597679499298</td><td>6.9067998</td><td>--</td></tr>\n",
       "<tr><td>239.10892315305699</td><td>54.995338971499301</td><td>4.3011999</td><td>--</td></tr>\n",
       "<tr><td>239.105466580057</td><td>54.998672856499297</td><td>4.1054001</td><td>--</td></tr>\n",
       "<tr><td>239.095654958057</td><td>55.0079141494993</td><td>4.6278</td><td>--</td></tr>\n",
       "<tr><td>239.09667155205699</td><td>54.9927280284993</td><td>4.7912002</td><td>--</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table masked=True length=85572>\n",
       "        ra                dec          redshift cigale_dustlumin\n",
       "       deg                deg                          W        \n",
       "     float64            float64        float32      float64     \n",
       "------------------ ------------------ --------- ----------------\n",
       "  246.555662513057 55.003551059499301 4.5766001               --\n",
       "  246.546474785057 54.997026013499301 4.8196998               --\n",
       "  246.397152569057 54.927381462499298 4.2312999               --\n",
       "246.38414028705699 54.952376668499298    4.8771               --\n",
       "  246.391144082057   54.9396872434993    4.0481               --\n",
       "246.44518433205701 55.000396006499301 5.0402002               --\n",
       "246.41164064105701   55.0204406554993 5.0925999               --\n",
       "  246.413862278057   55.0054131774993    4.6952               --\n",
       "  246.401368546057 55.012446681499299 4.9246001               --\n",
       "               ...                ...       ...              ...\n",
       "239.10977223905701 54.642360868499303 4.6788998               --\n",
       "  239.100740939057 54.648426408499297 4.5992999               --\n",
       "239.11018292746701   54.7110898655787 4.7849998               --\n",
       "239.09112570705699 54.715674113499297 4.5995002               --\n",
       "  239.083845867924 54.729565023113899 5.2941999               --\n",
       "239.11001283105699 55.001597679499298 6.9067998               --\n",
       "239.10892315305699 54.995338971499301 4.3011999               --\n",
       "  239.105466580057 54.998672856499297 4.1054001               --\n",
       "  239.095654958057   55.0079141494993    4.6278               --\n",
       "239.09667155205699   54.9927280284993 4.7912002               --"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "resultset.table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 85572 objects in ELAIS-N1 with a redshift above 4.\n"
     ]
    }
   ],
   "source": [
    "print(\"There are {} objects in ELAIS-N1 with a redshift above 4.\".format(len(resultset.table)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets just look at the galxies with a CIGALE dust luminosity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/rs548/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/astropy/table/column.py:965: RuntimeWarning: invalid value encountered in greater\n",
      "  return getattr(self.data, op)(other)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<i>Table masked=True length=548</i>\n",
       "<table id=\"table4924389864\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>ra</th><th>dec</th><th>redshift</th><th>cigale_dustlumin</th></tr></thead>\n",
       "<thead><tr><th>deg</th><th>deg</th><th></th><th>W</th></tr></thead>\n",
       "<thead><tr><th>float64</th><th>float64</th><th>float32</th><th>float64</th></tr></thead>\n",
       "<tr><td>246.14931803705699</td><td>54.809186737499303</td><td>5.0386</td><td>2.5676551478449e+39</td></tr>\n",
       "<tr><td>246.142722151057</td><td>54.796410900499303</td><td>5.5542998</td><td>4.9495341508333003e+39</td></tr>\n",
       "<tr><td>246.052311434057</td><td>54.833570216499297</td><td>5.7277002</td><td>1.9968836505858899e+39</td></tr>\n",
       "<tr><td>245.662680998057</td><td>54.803000700499297</td><td>4.9025998</td><td>2.9256870365008103e+39</td></tr>\n",
       "<tr><td>245.53703475905701</td><td>54.636527328499298</td><td>4.9355001</td><td>2.7411183850620899e+39</td></tr>\n",
       "<tr><td>245.427574447057</td><td>54.692954164499298</td><td>4.8910999</td><td>5.4584847835451502e+39</td></tr>\n",
       "<tr><td>245.25461235305701</td><td>54.697848129499299</td><td>5.6761999</td><td>6.2213053662153604e+39</td></tr>\n",
       "<tr><td>245.83274087705701</td><td>55.346899091499303</td><td>4.1290002</td><td>2.42122800018857e+39</td></tr>\n",
       "<tr><td>245.37802747605701</td><td>55.282590311499298</td><td>5.5229001</td><td>4.5451596391438499e+39</td></tr>\n",
       "<tr><td>...</td><td>...</td><td>...</td><td>...</td></tr>\n",
       "<tr><td>240.21271985405701</td><td>54.6885793514993</td><td>4.9288998</td><td>7.6360474232459204e+39</td></tr>\n",
       "<tr><td>240.08311886605699</td><td>54.6845476884993</td><td>4.9776001</td><td>3.6463978077662102e+39</td></tr>\n",
       "<tr><td>240.054497359057</td><td>54.710187106499298</td><td>5.1499</td><td>1.5075099182612501e+39</td></tr>\n",
       "<tr><td>240.27762134005701</td><td>55.2502169034993</td><td>4.9631</td><td>5.22705283434084e+39</td></tr>\n",
       "<tr><td>240.22574852705699</td><td>55.268350270499297</td><td>5.0251002</td><td>4.4490955781957303e+39</td></tr>\n",
       "<tr><td>240.15862158205701</td><td>55.100210505499298</td><td>4.5486999</td><td>2.5656889062731399e+39</td></tr>\n",
       "<tr><td>240.173950076057</td><td>54.9782225514993</td><td>4.4177999</td><td>6.8265540614961194e+39</td></tr>\n",
       "<tr><td>239.56528639605699</td><td>54.808217529499302</td><td>5.0035</td><td>5.9250862731313694e+39</td></tr>\n",
       "<tr><td>239.356651982057</td><td>55.013747431499297</td><td>4.6564999</td><td>1.51565546266832e+39</td></tr>\n",
       "<tr><td>239.598031613057</td><td>55.130811222499297</td><td>4.9376001</td><td>6.2590117397044798e+39</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table masked=True length=548>\n",
       "        ra                dec          redshift    cigale_dustlumin   \n",
       "       deg                deg                             W           \n",
       "     float64            float64        float32         float64        \n",
       "------------------ ------------------ --------- ----------------------\n",
       "246.14931803705699 54.809186737499303    5.0386    2.5676551478449e+39\n",
       "  246.142722151057 54.796410900499303 5.5542998 4.9495341508333003e+39\n",
       "  246.052311434057 54.833570216499297 5.7277002 1.9968836505858899e+39\n",
       "  245.662680998057 54.803000700499297 4.9025998 2.9256870365008103e+39\n",
       "245.53703475905701 54.636527328499298 4.9355001 2.7411183850620899e+39\n",
       "  245.427574447057 54.692954164499298 4.8910999 5.4584847835451502e+39\n",
       "245.25461235305701 54.697848129499299 5.6761999 6.2213053662153604e+39\n",
       "245.83274087705701 55.346899091499303 4.1290002   2.42122800018857e+39\n",
       "245.37802747605701 55.282590311499298 5.5229001 4.5451596391438499e+39\n",
       "               ...                ...       ...                    ...\n",
       "240.21271985405701   54.6885793514993 4.9288998 7.6360474232459204e+39\n",
       "240.08311886605699   54.6845476884993 4.9776001 3.6463978077662102e+39\n",
       "  240.054497359057 54.710187106499298    5.1499 1.5075099182612501e+39\n",
       "240.27762134005701   55.2502169034993    4.9631   5.22705283434084e+39\n",
       "240.22574852705699 55.268350270499297 5.0251002 4.4490955781957303e+39\n",
       "240.15862158205701 55.100210505499298 4.5486999 2.5656889062731399e+39\n",
       "  240.173950076057   54.9782225514993 4.4177999 6.8265540614961194e+39\n",
       "239.56528639605699 54.808217529499302    5.0035 5.9250862731313694e+39\n",
       "  239.356651982057 55.013747431499297 4.6564999   1.51565546266832e+39\n",
       "  239.598031613057 55.130811222499297 4.9376001 6.2590117397044798e+39"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "resultset.table[resultset.table['cigale_dustlumin'] > 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then investigate the data using typical Python commands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/rs548/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/matplotlib/__init__.py:1401: UserWarning:  This call to matplotlib.use() has no effect\n",
      "because the backend has already been chosen;\n",
      "matplotlib.use() must be called *before* pylab, matplotlib.pyplot,\n",
      "or matplotlib.backends is imported for the first time.\n",
      "\n",
      "  warnings.warn(_use_error_msg)\n"
     ]
    }
   ],
   "source": [
    "#You must have matplot lib installed to produce these plots\n",
    "%matplotlib inline\n",
    "#%config InlineBackend.figure_format = 'svg'\n",
    "\n",
    "import matplotlib as mpl\n",
    "mpl.use('pdf')\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/rs548/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/astropy/table/column.py:965: RuntimeWarning: invalid value encountered in greater\n",
      "  return getattr(self.data, op)(other)\n"
     ]
    }
   ],
   "source": [
    "example_table = resultset.table[resultset.table['cigale_dustlumin'] > 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1244c2160>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEDCAYAAAAyZm/jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+QHdV1579nnp7EjLAZARMDD4S0KSIFGaQxU4I1VIxc\nayTAwMSQAMHZbEouFSmzVWid2ci7XpCwd62UKusfazusiqKIKxhkDIyFJVvYBQkJRAmjjASIIC8G\nDHp2wtjS8EMa0Mzo7B/v9ain3723b/++/d75VFFo+sfr27dvnz733PODmBmCIAhC59BVdAMEQRCE\nfBHBLwiC0GGI4BcEQegwRPALgiB0GCL4BUEQOgwR/IIgCB2Gs4KfiO4lojeJ6IUI51xPRExEA75t\nf0RE/6/53x9l01pBEITyQK768RPR7wB4F8C3mfnDFsd/AMAOAHMB3MbMI0R0KoARAAMAGMAeABcx\n8+HsWi4IguA2zmr8zPwUgEP+bUT0m0T0IyLaQ0R/R0RLfbu/CODPAbzn27YawI+Z+VBT2P8YwJqs\n2y4IguAyzgp+DVsB/GdmvgjAnwL4FgAQ0UcAnMPMOwLH1wC84fv7YHObIAhCxzKn6AbYQkQnA/go\ngIeIyNs8j4i6APxvAP+poKYJgiCUitIIfjRmJ+PMvMK/kYhOAfBhAH/T/CCcAWA7EV0LoA7gct/h\nZwP4mzwaKwiC4CqlMfUw89sAXiWi3wMAarCcmd9i5tOZeREzLwKwG8C1zDwCYBeAK4hoAREtAHBF\nc5sgCELH4qzgJ6IHAPwDgCVEdJCI1gK4BcBaItoHYD+A60y/wcyH0Fj0fbb5313NbYIgCB2Ls+6c\ngiAIQjY4q/ELgiAI2eDk4u7pp5/OixYtKroZgiAIpWHPnj2/YuY+m2OdFPyLFi3CyMhI0c0QBEEo\nDUT0c9tjxdQjCILQYYjgFwRB6DBE8AuCIHQYoTZ+IroXwCcBvKnKkklEQ2j413u/99sA+pj5EBG9\nBuAdANMApph5IHi+IAiCkC82Gv99MGS0ZOYtzLyimUrh8wD+NhAktaq5X4S+IAiCA4Rq/Mz8FBEt\nsvy9mwE8kKRBgtBuDI/WsWXXAfxifAJn9XZjaPUSDPZLklihOFKz8RNRDxozg4d9mxnAT5r589eF\nnL+OiEaIaGRsbCytZglCoQyP1vH5R55HfXwCDKA+PoHPP/I8hkfrRTdN6GDSXNy9BsDTATPPZU0T\n0JUAPtusqqWEmbcy8wAzD/T1WcUgCILzbNl1ABOT07O2TUxOY8uuAwW1SBDSFfw3IWDmYeZ68/9v\nAngUwMoUrycIzvOL8YlI2wUhD1IR/M2c+B8D8H3ftvnNOrggovlopES2LpwuCO3AWb3dkbYLQh6E\nCn5VemQiupWIbvUd9rsAHmfmI75tHwLw980Uyv8EYAcz/yjNxguC6wytXoLuamXWtu5qBUOrlxTU\nIkGw8+q52eKY+9Bw+/RvewXA8rgNE4R2wPPeEa8ewSWcTNImCO3EYH9NBL3gFJKyQRAEocMQwS8I\ngtBhiKlHEFJEonSFMiCCXxBSwovS9QK2vChdACL8O4gyfPzF1CMIKSFRukJZUnSI4BeElJAoXaEs\nH38R/IKQEhKlK5Tl4y+CXxBSQqJ0hbJ8/EXwC0JKDPbX8OVPXYBabzcIQK23G1/+1AXOLewJ2VGW\nj7949QhCikiUbmdTlhQdIvgFQRACJHHJLMPHXwS/IAiCj06IxxDBLwg5UIagHqGBySWzXZ6ZCH5B\nyJhO0CDbibK4ZCZBvHoEIWPKEtQjNCiLS2YSRPALQsZ0ggbZTpTFJTMJIvgFIWM6QYNsJzohHkNs\n/IKQMUOrl8yy8QPtp0G2G2VwyUyCCH5ByJi8g3rEg0gII1TwE9G9AD4J4E1m/rBi/+UAvg/g1eam\nR5j5rua+NQC+BqAC4B5m3pxSuwWhVOSlQYoHkWCDjY3/PgBrQo75O2Ze0fzPE/oVAN8EcCWA8wHc\nTETnJ2msIAhmxINIsCFU42fmp4hoUYzfXgngZWZ+BQCI6EEA1wF4McZvCYLTuGJeEQ8iwYa0vHo+\nSkTPEdEPiWhZc1sNwBu+Yw42tykhonVENEJEI2NjYyk1SxCyx6WqS+JBJNiQhuD/ZwALmflCAP8H\nwHCcH2Hmrcw8wMwDfX19KTRL6FSGR+u4dPMTWLxhBy7d/ETmAtgl80on+KALyUks+Jn5bWZ+t/nv\nnQCqRHQ6gDqAc3yHnt3cJgiZUYT27ZJ5pRN80IXkJHbnJKIzAPwbMzMRrUTjY/JrAOMAziOixWgI\n/JsA/EHS6wmCiSISbJ3V2426QsgXZV5pdx90ITmhGj8RPQDgHwAsIaKDRLSWiG4lolubh9wA4AUi\n2gfg6wBu4gZTAG4DsAvAvwD4LjPvz+Y2BKFBEdq3mFeEsmHj1XNzyP5vAPiGZt9OADvjNU0QolOE\n9l2WqkuC4CGRu0JbUVR6BDGvCGVCBL/QVoj2nT6uxCgI6SGCX2g7RPtOD0kB0Z5IWmZBELS4FKMg\npIcIfkEQtLgUoyCkh5h6BEHQ4lqMgquUbR1ENH5BELRIjEI4LuVqskU0fkGISdm0vDiIl1Q4RUSL\nJ0UEvyDEoJO8XcRLykwZ10HE1CMIMRBvF8GjjKmwRfALQgzKqOUJ2VDGdRAR/IIQgzJqeUI2lDEV\nttj4BSEGReUEEtykbOsgIvgFIQbi7SKUGRH8ghCTsml5guAhNn5BEIQOQwS/IAhChyGCXxAEocMQ\nG78gCB1LJ6TdUCGCXxCEFjpBIHZS2o0goaYeIrqXiN4kohc0+28houeI6HkieoaIlvv2vdbcvpeI\nRtJsuCAI2VDGbJNx6OS0GzY2/vsArDHsfxXAx5j5AgBfBLA1sH8VM69g5oF4TRQEIU86RSB2ctqN\nUFMPMz9FRIsM+5/x/bkbwNnJmyUIQlG4KhDTNj91cpGZtL161gL4oe9vBvATItpDROtMJxLROiIa\nIaKRsbGxlJslCIItLuYhysL8VMbkammRmuAnolVoCP4/822+jJlXALgSwGeJ6Hd05zPzVmYeYOaB\nvr6+tJolCKVjeLSOSzc/gcUbduDSzU/kblt3USBmYX4qY3K1tEjFq4eILgRwD4ArmfnX3nZmrjf/\n/yYRPQpgJYCn0rimILQjLniauJiHKCvzU6em3Ugs+IloIYBHAPwhM//Ut30+gC5mfqf57ysA3JX0\neoIAtK+7oStl/FwTiJ1sj88CG3fOBwD8A4AlRHSQiNYS0a1EdGvzkDsAnAbgWwG3zQ8B+Hsi2gfg\nnwDsYOYfZXAPQofRzu6Gri6sFo2L5qcyY+PVc3PI/s8A+Ixi+ysAlreeIQjJcEUrzgLRbNW4aH4q\nMxK5K5SOvLTiIsxJUuBFj2vmpzIjgl8oHXloxUUtsopmK+SBCH6hdOShFRdpThLNVsgaEfxC6chD\nK45iTmpXDyMhHVwcHyL4hVKStVZsa05ywe9ecBdXx4cUYhEEBbbug52S0EyIh6vjQzR+QVBga04S\nv3vBhKvjQwS/IGiwMSfpTEKndFdx6eYnnLLrCvnjalyGmHoEIQEqk1C1i3Dk2FRbRhYL0XA14lgE\nvyAkQJXh8eST5mBymmcd54JdV8gfVzOAEjOHH5UzAwMDPDIilRqFcrJ4ww6o3ioC8Ormq/NujtAh\nENEe20qHovELQsq4WMhEEPyI4BfaDilkIghmxKtHaCtcCJiRfDuC64jgF9qKPHLs2ITgS74dwWVE\n8AttRdYBMy7MKAQhKWLjF9qKrBdWXQ3BF4QoiOAX2oqsF1bjzCiKXmwWhCBi6hHaiqwXVqOG4Itp\nSHAREfxC5uSdjzyNhVVdm6MWgWnn+sBCeQk19RDRvUT0JhG9oNlPRPR1InqZiJ4joo/49q0hogPN\nfRvSbLhQDjyNt0x5a0xtjhqC72p2RqGzsdH47wPwDQDf1uy/EsB5zf8uBvCXAC4mogqAbwL4BICD\nAJ4lou3M/GLSRgvlIarG60K1orA2R5lRuJqdUehsQjV+Zn4KwCHDIdcB+DY32A2gl4jOBLASwMvM\n/AozHwPwYPNYoYOIWsLQhdlBmlq6RPEKLpKGV08NwBu+vw82t+m2KyGidUQ0QkQjY2NjKTRLcIEo\n7pWuuEqm6RLqanZGobNxxp2Tmbcy8wAzD/T19RXdHCElomi8rtjDw9o8PFrHik2PY9GGHVi0YQf6\n73rcOCvxFoXP6u3GL8YnsGXXAafXOIT2Jw2vnjqAc3x/n93cVtVsFzqIKO6VednDw9YRTG0eHq1j\n6KF9mDx+IvHy4aOTGPrevlnnBq8nLp2CS6Qh+LcDuI2IHkRjcfctZv4lEY0BOI+IFqMh8G8C8Acp\nXE8oGbaLoVFdJeNgK4R1bd6y68Asoe8xOc3aBWtx6RRcI1TwE9EDAC4HcDoRHQRwJxraPJj5bgA7\nAVwF4GUARwH8cXPfFBHdBmAXgAqAe5l5fwb3ILQJeWS11Anhjdv3W13XZHaKaqoSl06hKEIFPzPf\nHLKfAXxWs28nGh+GjsYFF8WykHVWS52wHZ+YxPjEJACzKUZnjvL26baLS6fgEs4s7rYrrrgoCg1s\nha3Om2ho9RJUu6hle7VCWpOUuHQKriGCP2NccVEUGqiEsA7V7GCwv4Ytv7ccvd3VmW0LeqrYcsNy\n7UxFXDoF15BcPRkj9l23UK0jHD02hcNHJ1uODc4Ogia7jdcusxbeUpglX8S8akYEf8aIfdc9gkI4\n6OkDtJpixCWzPMR9Vp30sRBTT8aIfdd9bEwxYrIrD3GeVaetxYnGnzFSeLschJlixGRXHuI8q06L\ntRDBnwNltO920rTXBp3J7hTfIq/gBnHMq532YRdTj9BCp017bdC5cR45NhXaL1J6MV/imFdNMRjt\n+PxE8AstFG3PjvuiZfmCDvbXcPJJrRPkyWnG5767T3st+YjmTxz3Wd3HYtXSvrZ8fmLqEVooctqb\nxCMja6+bcYXLJwBMM2uv1Wm2Y1eIal7VrcW16/MTjV9oIc189FGJO9vIY5Ziun/dtTrNdlxmBvtr\neHrDx/Hq5qvx9IaPY7C/1rbPTwS/0EKRLqhxX7Q8XtCh1UtQrbTa+U3XKvIjKiQnyfNzeW1ABL/Q\nQpEpBuK+aHkI2MH+GubP1VtHVddadJr6+rrtglvEVYJcX9sRwS8oUU178yDui5bXLOWtCbWd32tD\nkN2vHFYe+/TPDjmpCXY6QS0dQCwlqGgHiTBkcVfIjDixAHED3vIKlNP5iC/oqc5U6PK3YZpbi7Z4\n+DVB/z2UlbLHfugcBL78qQvw9IaPR/ot19cGiA0DsygGBgZ4ZGSk6GYICdDlvyl7VkrTfQFo2WdL\nrbc7snBxiTI+7+CHSpesL86zuXTzE0oFIcvnTER7mHnA5lgx9QiZ4PpUNy6m9Q/VPdviiiYYl7I9\nb5UNXiX0gXjPxvUcXWLqETLB9aluEnQ+4qZ7qxAZzT5l9/Ip2/OO8pGO82xcz9Elgl9IFW/6rBNx\nZRdwJnT2f//0/gvDz+P+3a/P6h+XNMG4uJp+XLfuYPtBSvJsXM7RJYJfSA2VnddPOwi4IH7B0ttT\nRbWLMHn8hFgnNMwIl25+AquW9uHhPfVZQp8AXH+RuwLClqHVS0JrGuSNKZrbVDt5QU8V40cnndPS\n08RK8BPRGgBfA1ABcA8zbw7sHwJwi+83fxtAHzMfIqLXALwDYBrAlO3ig1A+TNPnWhu+REHBorIR\ne0K+Pj7Roul7+598aSzTduaBi6YN07rD0OolWL9tr3Jm2jN3DkbvuKJle1yvJRe9nUIFPxFVAHwT\nwCcAHATwLBFtZ+YXvWOYeQuALc3jrwGwnpkP+X5mFTP/KtWWC86hmz4TUGqPFR1RF3N15q8odnAX\nhYiHa6YN07rDYH8Nt2/ba32eyzmk4mDj1bMSwMvM/AozHwPwIIDrDMffDOCBNBonlItOS0+Q1sKl\nbf+4Hg3qGmHjsRZhvLqcQyoONoK/BuAN398Hm9taIKIeAGsAPOzbzAB+QkR7iGid7iJEtI6IRoho\nZGysHFNfl3NxFIHrLmxp4T33NCJgovSPq0LEVcLGY5Tx6nIOqTikvbh7DYCnA2aey5i5TkS/AeDH\nRPQSMz8VPJGZtwLYCjQCuFJuV+q4OoUrEhftvGkTtoBtQ623O1b/uCpEXCVsPEYZr3G9llz1drIR\n/HUA5/j+Pru5TcVNCJh5mLne/P+bRPQoGqajFsFfNto1T3dSXLPzpk2SIC0gWeSmq0LEZcLGo+14\njeu1FHZeUWs2NqaeZwGcR0SLiWguGsJ9e/AgIjoFwMcAfN+3bT4RfcD7N4ArALyQRsOLwpvm61zB\nRPtqb5I836Rmr04xpblI3Iy1pvOKXLMJ1fiZeYqIbgOwCw13znuZeT8R3drcf3fz0N8F8DgzH/Gd\n/iEAjxKRd63vMPOP0rwBE2l/TW2m+aJ9tTcm/28TFaLEeWs6wZTmMnFns7rzirQaWNn4mXkngJ2B\nbXcH/r4PwH2Bba8AWJ6ohTHJwgYfNs13WfvK4iPYiQJoaPUSrRugDgLwF7+/PJX+aXdTWiehmz3G\nUSyi0rZJ2rLwgDBN8/MsVhKVtKeUnexWONhfw4KeaqRzPE8F8QAT/OisAwRkPj7aVvBn4QGhe1De\ngp2LQh9I/yPY6W6Fd16zDNUufQnGIAt6qh37oRT0rFrap9zOQObvUtsK/iyCicqyuBaML0h7IbrT\n3QoH+2vY8nvL0V0Nf30IjVQOnfyhFNSYUnVk/S61reDPQkgXWYvWFpUZRqebxv0IdkqErilAb7C/\nhn/54pXo7TabfUwBKZ3yoWwHsgjWNNnys36X2jY7Z1YeEC4trqkWWFVmGEZD89SlAo66UOtiJsa0\nsXUOMNXgDaPdPpTtShaOIsOj9ZZ30oOgrt+cJlJ6saToSt2ZvI5UEaNxS+a1u1fPik2PY1wh1CtE\nOM48c88bt+9XHheG62UJhRPozKULeqromTsn1jtgMsF++pKF+NLgBZHbGaX0Yttq/O2OboFVV+lJ\nFzEa15fYpZlP2gyP1rXC3Ovb+vhEZLdOPxHWhgUDeSggOpPc4aOTM6m4o84CTGa+OEI/KiL4S4pu\n4Ewzt2j+JjNMpy/UegyP1rHpsf3auqtpc+TYNIa+tw9A5+Z2Skpe+bJsg/aiBF+ZqrXlQdsu7rY7\nJtfSKAvQnbJQa2J4tI6h7+3LTeh7TE5zW3r25JW1Ni+34qHVS1Ct2E3RbBWmoj0EReMvKaYF1ihm\nmE5YqA1jy64DmJwuZq2rPj6B4dF622j9eWatzXW2ajk8bBWmotNviOA34PICZloDp+gB6AJFm7Xa\nKZ13Xvlnhkfr6NKsZ6U9W92y68CsOso6onrjFLlOJoJfQxny7ac1cNp5odaGuInX0qKd0nnnoYV7\n76ZK6GcxW7UZGwTglksWAmh47LiuRImNX0OnpyVoV1T2Z50NtwuYyctToWzdcIr88KRJHmtGumSJ\naWRADeL526uoEM2so33lxhUYOPfU0qTmEI1fQ5m9XdI0Ubls7oqKbhb35U9dgC03LJ/l1dPbXcXG\na5fNutc0qm/p8BJzlbVvPfJYM9K9g8eZZ/ovrXG7ZdcBbZBVMOPqpZufKE1xJhH8Gspa7ShNE1UZ\nzF1RMM3ibJLsefuzcPv0EnOVsV/95LFmFPZupjludR8ZVvxWmZRFMfVoKNrdKi5pmqjazdyV5MX0\nTETrt+3F2xNTaTcNQPuYewb7a3h6w8fx6uarM8laG/ZupjluTW7Ttse6qCyKxq+hrN4uaWodZdJg\nbNBpil1ERjNLUINULSqmgQTz2hH2bpoKnCzasGMmur23uwoiYPzopPb9jmK6KpNrtAh+A2X0dknT\nRFVWc5eOVUv7cP/u11tsttPMRlNA0gLrtriXNctNwuz3YV5a3ofbn5ZDZw6KogCWSVkUwd9mpKl1\nlEmDCWN4tI6H99S1wtW0CFfWGY6HSwv0SdtiY79XjVsbdGMgigJYFmVRBH+bkabWUSYNJgwbrV0n\n4PPy849a0tEGlxboVW25fdtebNy+v8WDSodNgFhw3EaZSaX1kdd94Fz5CFsJfiJaA+BrACoA7mHm\nzYH9lwP4PoBXm5seYea7bM4VomEzcNLUOsqiwXjo+sfmhe7VCN64GmQUqhXCndcsS/1384qkjdsW\noGFysf0Y2a47+cetKQVykDTMmLqP7UMjr+Ppnx2aOa4+PoGhh4pJ1Bfq1UNEFQDfBHAlgPMB3ExE\n5ysO/TtmXtH8766I5woWtHORc1NiL9ukX18Yfh7rt+1V9o/NC/3ue1PK3/Yqr4VV24pLhQhbblie\nycuvE3hFeBCZPr62XjdxPGdsk6ylZcbUfWz9Qt9j8jhj4/b9ia8ZFRt3zpUAXmbmV5j5GIAHAVxn\n+ftJzhUCtJt7pYfpg6baN/TQPvTf9fisD8HwaF25cOv1j8oFMMjkcXO2zPenjie+VxX+wKO00UUc\nZx2JrCLs42szK4vtZq2w9/RUu7Cgp5p6GdWo5qI4hXySYmPqqQF4w/f3QQAXK477KBE9B6AO4E+Z\neX+Ec0FE6wCsA4CFCxdaNCt7XLHHeRTlXpl1P5g+aEePTbXsmzzOLQUwTqp2aW25vxifsLb76vrS\nlCYgqXunzsRki+n56NqWlUuqiTCTmc2sLM66ky7J2oL585TFiZJSdO4nG9Ja3P1nAAuZ+V0iugrA\nMIDzovwAM28FsBVolF5MqV2xcWlRzCML98owoZ5HP5j8rm2YmJw22t+9/rGx+/r70t83ugF5nFlb\nO9WWw0cncf7/+CH+16cuTD26ulZwwQ8/psjnKGaWqOtOSRWmODWphx7aZ5XRE8hmUT8MG1NPHcA5\nvr/Pbm6bgZnfZuZ3m//eCaBKRKfbnOsqLppV0o4mtlkzyKMfdB+uNKwRulS5YX0Z7Bsdp3RXU1kQ\nPDp5HEMP7Yu8XhP2fFyLQB/sr2H0jivw1RtXWBcLioJqPUj3fBgILRQTZ11tsL+Gk0+y06m7CJks\n6odh07pnAZxHRIvRENo3AfgD/wFEdAaAf2NmJqKVaHxQfg1gPOxcVzFpoUnTrsY1ndhMc6P8to3H\nRx7mpaHVSzD0vX0txVCiWCN6u6t4f+r4rPvxUuXGCbaxDdp6+71JHJtKx9vHW2MwjYXg89XNirzn\n46pLbhbeYrrZz/UX1fDwnrryeYbNYON6RY1b5HJSJQLMi1DBz8xTRHQbgF1ouGTey8z7iejW5v67\nAdwA4E+IaArABICbmJkBKM/N6F5SRfdSEU6YIOrjE1i/bS9u37YXNcsXKqnpxPTCRP1tG6GeR/Tu\nYH8NG7fvT7TI9cnlZ2Lg3FNnCbhVS/vw5EtjWLxhh1LgmfrS9sN2nBvaelroBLmqJnB9fEJrZvI/\nn7K55MZFJ6SffGmskYF11wFl/8YJ3gsbHyb58ZUbVxT+PKyStDHzTmb+LWb+TWb+n81tdzeFPpj5\nG8y8jJmXM/MlzPyM6dwyoJoiq14y729b18osTSdRf9vGNS4vU8FbCT0bnnxpbFZysKHVS/Dwnnps\n19ei0lJ46Zn9eB90VUZQRmuOn7JGVyfFJKS9saGzHpqC96Js9xhavUR5LS8La9FIdk4Nnu+23w4Z\nZnmwEeBZmk6i/raNUFf1Q9rFLgDzi2TjP++Z4Dzb7qbH9if6wNq4f2aBSjCEmZ0YyPz5lAEbIR1V\nkA+tXoJq12wRXu2i0A/rYH8tsudYnkjKBgPBKbJNBGDcKWBw4MVZB4hqlrG1/9qaCpK4fQ6tXoL1\n2/YqX5b58+Zg/rw5xr4PmuB02L50XruTmqDiEGxjWJtrvd2ZuCWWDZU3TVBI61xKj7w/pc/QGlTd\nLZ0OdB5VLiQ5bFuN3zbaMwo2WqDNFDBMy44boRvHLJNW7vSkUcVhGpKp76O4U0Z56Qb7a5g/L3/d\nKNhGU5s71ayjJURIezPYoAullzYiOF637DrQ4nQwOW0O9PNwzaPKT1sK/qxSG/jNHkA826qN6STM\nVq/7qOVlllGRxtqFzrf8rN7ulr73Ik8rRNZCP85Ll/e0XNVG3Uevt7s66/lmoeyUCVshPdhfQ8/c\n1g+6arwmMc0W+T6G0ZamnrQTU6WdaS/MdGIabGGeO0k9OKLek3d8mFuhDauW9uGvd7/esn3Raa0B\nWDb1b7urXTh1/rxEbox5RWFS81o6UxsQ7sbrWsBh3kQR0rbHJvVqc9Wjqi0Ff5oLqKYXKitMgy3L\nbItRhYeN8I1iWnnypTHl9qd/dggrNj2OtyZOVEpSLd4GOalaSWz71sUYpEmFCD/78lXGY8IEiEtZ\nOHUEXVLT9mOPIqR1x3pBXd6HtZ1qUvhpS1NPmrUvdS/Upsf2Z5Yp02QbzNIrKKq5JszbJOoLYrqH\n8YnJWUnabIqd2wTRhDHYX8N8hVkgTdLIm5M07UXWDI/WMfS92c9tfGIyVrSyjig2ddOakf9ddtlc\nk4S20/iHR+s48n5rMWzVAqqNSUP3QqkET1oalmlqrzOrpOEpEPWjYhLUNgFtwWdwSnfVyoPGNgdK\nUHuLS9IYAxtWbHo8kfZrChgy1RPOC5X9HbCLVrYlSpSy/9iwoC5XzTVJIC4gS18YAwMDPDIyEvk8\nnelhQU8Vd15z4qVSHdddrSi/5FGKOACNF+3VzVdHbrstUdruHW9rs9fda4UIx5nR21MFM2ZMLkfe\nn1IKahv3QtV9VCuUiUnF8/qxja4OEnUMJCE4Vv2YnuXwaF3rDuuCu+fiDTu0i/BZvzOAue9Mbfuq\nA1G2thDRHmYesDm2rUw9OtNDz9w5sx7exu32wT266aMuqChrH90oU8+o3k266e80MxiNWY7f5HLk\n2FRLcIuteUf1rCanGfPnVmzdpFvQnRc1ujrIqqV9sdsUlcNH1W6FYc/S9YAh03txVm93ph5JYX1n\napv/uHbymmorU4+NSWJ4tK41J6jO100fARS26BOcenoDUmUWUn3gNm7fb5W4rCsk1/zkNGNBTxU9\nc+dE9pzRPaujx6bxlRtXGD2FugBUArMDb9aj03o9/B/4YF6fH+z75czY8DRvAMYi7VkwMTmN2wM5\noGwWb107Qq62AAAYDUlEQVQOGDItki86rTtTj6SwvjPVCfCPl3bymmorU49uSu6f6pqm7cEpsU2u\n+jjunGkWNjGZfkxC0GYKa5oCe8Sdpid5Vp5QVvWhrVmmu1oJ9QqqVggnz5tjtZCcJWFt9Z6laix4\nZq7e7iqIGgveRWXoHB6t47898lxLUjtdAF5aJqpFG3Zo973WHLvDo3Xcvm2v8hjP1TZsvBZNx5p6\nbFb1TdPeqNGzcaJe0w4uM2kzJk0vSX3TqMeoSPKsDh+d1H44baKrK0RWKZcnp7lwoQ80nqepVKLK\nAwWYLVDHJyZx+Ohk6h5oURjsr2HB/Hkt27M2UdmUnxzsrxnNt0VVv8uKthL8NvZvnaDyQrg9G97n\nvrsvkyyaOkH9ue/ui2U7NA1Ik9kpbn1TP0lMW0melZeXR2frDouuLqLsYFKmmbXPwj8uPWUkLKlg\nUUWFoghKk1LhmTcXbdiB3/z8TiwyvDs25SeHR+s4cqzVG9DL9ZOmi7gLtJWNHwgPdNEFZFx94Zmz\ntusGS70ZPRt3mqwb+N71otoOTUErg/01ZZk7b38Y3vX/+6PP48ixVg35IwtPSWQuiPOsVGaBoK07\nGOEbnB2Y1g9cxbP168wRURO72R6TBv5noFs3Cj5Xk1IRNGmFvTs25Sd17qYnn3TCMaSdArnaSuO3\nQadpPvnSmNX0H0DkabLfG6DLop5g0vTB/gF55zXLEiWKGuyv4T1NoZHdrxzO1NNB9ayimgVU5rhV\nS/tSa2MeeOUjB/trxnxGpr9tzsmCoGlTJfS7qxXccslC6yAplVeeh+rdSWJW9IIAB/truP6i2qwc\nUddf1Kq4lMXzp+00fhtUmuZ6jSalIkqglk47CSNq+mCdzTtKUIsO01Q5a0+H4LOyKZIehi41hB+i\naKUfsyJYPtI2hYDJU0V3ThboXKy92JCo49HklecRfHds3oGwdA/Do3U8vKc+8y5MM+PhPXUMnHuq\nNj7IZc+fjhP8unwhUZNx2QrmsIGvm/pGTR9sGlhJIw8rBrdO0zpIFnVe08idYnp2nqAdOPfU0DxE\nyvNT/GB4bfnS4AUz22w/5MHjTgl49axa2octuw5g/ba9mXr56Pr6OHMsb7C4TglxTcDeuLJxpy1D\nviSPthP8QZuuV3f1F+MT6O2p4q2JSfgj/r18ITeuPKelIHN3tYKTql2xbeRA+MDXuWPmbTs0uZje\nfPE5yqyZOjxNJ47mE+bqGib4bFxwTfEJDOD+3a9j4NxTcf1FtUj3DTSEfrWr8fuWmSX0vwX17MT2\nQ647Lk3NNKy/067ZbKNwGYuqaAgbVzZePWXy/Gkrwa8a0P4XV+eaN3mcZxVk9j/4kZ8fUr78tnbi\nsIFvq8Gl6fsfJEwQfGnwgkgCUOUuGab5DI/WW6pd6QRSXIHm7Q8ztzGATY/tj6252+YSsiGLRei0\nNFObD0ja2S1tZuZeURV/O2wwfVBtPmBpf+SyxErwE9EaAF8DUAFwDzNvDuy/BcCfoTE7fQfAnzDz\nvua+15rbpgFM2QYYxCEsW6SJ+vgENj22vyXARTe1tLETA3YD3zTgGi/Xc5jwLbDWxydw+7a9+Pwj\nz+GkaiVxUE5YVs6obn86wRrUfPy5/HVBPDqBpPPWMQm0KOPDBf99AOjKIFdEGprp8Ggdn/vuvpZn\nrfKwAtIz+6nKK6rwMuimed2w97hMKZxDBT8RVQB8E8AnABwE8CwRbWfmF32HvQrgY8x8mIiuBLAV\nwMW+/auY+VcptltJ0imV97LXxyewfttejPz8UGhRlKi21qgRvqZBPjF5fOaDYJquh7XTlNI3jp1b\nh1/zCWqLptdY9cFQaZq6dnrnuzjlDiPFycMMSTXTsJmTanE1zcV+nYtykMNHJ2e900kWWnVrJuu3\n7cWWXQdmvVNZzczTxEbjXwngZWZ+BQCI6EEA1wGYEfzM/Izv+N0Azk6zkbakWS3Js/X2zK0ofdh7\n5las7aRxB/6WXQcimQ1U2rHNdFzXb7YRrjZ4mk9YxS4VQYGk0+x1i9De+XlV0yoKW3NgUs00bOaU\ntWkjbp2FqN54qrVCb3ZqMkm6KOiD2Pjx1wC84fv7YHObjrUAfuj7mwH8hIj2ENE63UlEtI6IRoho\nZGzMzowSxCZcPwqMRtIwFUeOTRvNI2n488bRUIPn2BRX0fk5pxXh6vllA5jx6Y5CcD3FFARn8tfW\n3Wd3NbtwlgU9VXz1xhWxx6UujUCQKKlAbKKmTZjGZR6mjSQfFpt3StWXf7379ZlxazJJloVURzwR\nrUJD8P+Zb/NlzLwCwJUAPktEv6M6l5m3MvMAMw/09cULsFEN6E9fstCY5ySMqKLPXxc3Tj6eqMFe\nQYIvhY09VycIdMFCUfCSWEW1sfsJrqfoXnx/u/33AQD9dz2O27ftnXV9Lwjn+ouym6C+1dROvXZE\nZeO1y6yOs62e5o0vL27lKzeusM4z5aHr/wpRLtWpVB/wahdhQU915rknSZsed5yWyZRoY+qpAzjH\n9/fZzW2zIKILAdwD4Epm/rW3nZnrzf+/SUSPomE6eipJo014Uy1vqnb/7tfR21PFu+9NzTKbVLsI\nIIQW/jD5sKtIUhc3brCXh0rbsrXn6qaoNvVmdQuzURLkmQiet2ppH+7f/boyxD94H17JP9U9eEE4\n8+Zkp/EfB7D+u3tjeQh92he4FYZtSvI03Dh1pqK8ShLa2NLjuEnHMUP6cdF7R4fNiH8WwHlEtJiI\n5gK4CcB2/wFEtBDAIwD+kJl/6ts+n4g+4P0bwBUAXkir8TqCGvfho5MANVPToqER3LjyHJw878R3\nTzXb765WcPPF57Tu0BBWF9fL86PDVtPornaB0Lgfv5bjvXj+WcPRkGIpJpPUYH94vdlabze+cuOK\nmdmBN7uKknTNf56K4KJwMD8+AcrweUCfg8VjYnLaqtxjEsKE/oKeaouGuqCnioFzT7W+hk0Ssag1\nlXUkNRWlgSoVR1gbr7+oMetUjXW/zIiD6qPicvqGUI2fmaeI6DYAu9Bw57yXmfcT0a3N/XcDuAPA\naQC+RY0X2HPb/BCAR5vb5gD4DjP/KJM7gfmLPTnNmD9vDvbeeYVSG5hTqeDGlbWZYC+/FvHInoMt\nOcSBhgA+df481McnZhZCt+w6YKwda9KwbDXiU+fPm8kB7l+E2rLrAEZ+fmhWINrho5OoVgi93Y3g\nNb83wsbt+3Hk2NSMYFRpgKZ6szotW4WuFnLYzCv4QqmElz/QKbgoV4bF3PPP/AD2/+KdWdu8SlyA\nnTZus2BrOyuwrVnr+iKmv41hs50kruCqcpmup2+w8uNn5p0Adga23e3792cAfEZx3isAlidsoxW6\nert+vAGu03yefGlMWVRhXrWiFPxe8qzgA65WCNUuUnrkeCmYgdYBYCuovPtQDa6gCQQ48dHbeO2y\nWcerPk5Bk5TJ48dWyzPVQgb0PvOqGrlh7rXB/tCZoYLteG/yeGoeTFF5+meHlNujeKGklY/GZWGV\nhDDzq0npqgW8ejzzr6mGs+vpG9omctfmi+0N8KgBLDr3saOTx5WZAr2ShDqBNs2M27ftxabH9s/S\nFMISawXvQ6f9qvBiE2xMzf5+SMOea6qFrOtzApQfYZ3wYkAZUBR2v9UKzZRYVJ1fNLazQL+m3ttT\nxZH3p1p8zNPIR+MKUSPZw975tCtsuZ6+oW3SMod1qH+A6+yhXURKO5zJNq0z6YwfnQz1igkW1g7a\nJXu7q6hW9Pb5qIPIVqT57zcNe67pJYha4MLkshtHaM+fO2fGJHDcMaEP2C0YegvY/jWt8YkT1bbW\nb9uLLww/H/osXRdWHnG85sLGmU3q5ii4XrilbQS/qUODA3xo9ZKWBU/gRJrh4ACK8/A9LSTMfzu4\nuOZftNp75xXYcsNy7YtqqlCVhOD9hi2khWF6CaK+cH7hlQb+NQyTm2ha14uCreDZ9Nh+4wK2F4zo\nJS7TPUvXhZVHnEXqsHGW9oJ12h+StGkbU09kk4RGOvqntrMrB9mHz/sXPYFwE4JJozItounuOU5W\nSY+ealcq03p/353SnLn4hVOwj/wJ2k4KCajy+sSmGHwYfqEWZgrRVb8K0kVApYtC3WBNqBYMAbWJ\nwyZ9AQOhJpu8c83ETTwYZ2ZiswZiu2CddaqWPGgbwR+lo8Nc/FQLhWFWAF1hCe//Jtt9XI3KdM/e\nQlRUJqc5UWlJoHWRcHxicibARpdQ7v2pE4vnh49O4vam19HGa1uFn0caXjvBZHmAfgzZ5IjxBDYQ\nf83g04Ec/B6mxVcbwkw2eQqrJAvJcXMNpeGJFKXdLns+ETto1xwYGOCRkZHMfj9MU1zQU8XbE1OR\nX1r/Kr8q18cP9v2yZU0gq8CXMC+nxjSUZ2X99N9HnAUtD12VLN3v6o732qnrHxtPrjBei1AMxOZ6\n/nuMOiNRFV7xo+snG8+lYNuKJuoY8aMLzsojliBJu7OGiPbYZj9uG40/CiZNsVohvPtedKEPnPj6\nB33p6+MTeHhPfSZsPw+NKqi99fZUwdywaXvX1ZWbjOPXrTs/yXagYXq7XZEB0X+PwTz+tkS123vX\nM5l8/PcSZUaicg20jUlgQOs+7McV+zJgzggbNuMs0oxSlgXwMDpK47cJye6pdil99j1sbP2mNA8m\n398kxBHQYdqLSrPytEvTfZg0eNV5puP9hGn/utz+lS7CdOChJdEQ++963Bh74A+uC5sh6Nph6nvV\nNb16BLp+XNBTxegdV5hvLEfizvKKpl00/srGjRszbk50tm7dunHdOm0iz1h4L9Kho8eMx5m0pu5q\nxWrBznTEO+9N4W9/OoazF3Rj6ZkfDP0tE8Ojdaz9qxF88QcvYtf+f8Xb701FusZp8+fib386hqnj\nsxdd77jmfCw984NY+1cj2v4yXUP1u6bzTMf7mTrOeL7+FtZetrhl39IzP4i1ly3G7f/ht7DotPl4\nvv4W3n1vCrXebmy6dhnWfPiMWduuuvAM3P+Pr+NLP3gRD40cxGnz5+Klf30Ha/9qZNY2Vf+d8cGT\n8ON/+Tflus87703NnDvYX8PZC7pnXfe6FWfh1+8em/n7jmvOVwo4Xd8HfRK85zXYX8PayxZj0Wnz\nlc/0i4MfTjze0sT0zE3PuWjC3pki2bRp0y83bty41ebYjtH4bbVKE1+9cYVVEiebxG5JNYSo9mbT\n7+hmCjY2at01wmZXwfOGR+tWi6cExCrSHWxbsO9UqSM8DylVGo+w9trMjFTt8p6Fqd9rvd2FlelM\nk+HRutZslsZzzooo/ZvnsxAbv4I0bHA2Nl5PWAQLtwexsWWasIlUtrlnk+eBjY1ad40wl0tdlaaw\nD4bKcyPqy6XqO116DX8KDF3BDZVSoTtHh+1CNZE6ormMeDlybDx04ghQXTGVpEI4itunqykw2iaA\nKwwbl8nuakWbx9tbCBzsr2mP8fLXfGnwAqsgI9sc/SpshLpt1Kcug6BNAFrYNaIGBXkBRqriJboM\niFGjOKMoAcHPgSpQKOz3vAVqU4ZG2yRhYRP0JLUgisAm0CnOPZmKqeTVL2llQ82CjhH8yuINzayV\n/ki9jdcuCx2IumP+4veXz/LfD9PMkgyCMIFrE3gT9kIFo2RV9uWwa8SNYLSNpIzzciWNRA0K+t4e\nuypZJmGTlldIlP5wIW3wYH8N119Um0nL7RXH8T/nOM/Y5kOatRB22QOoY0w9UV3AwiL8bH+rFmIu\niTsIVFGWUe3KNkm5/NPaONPtJK53NlNqXd+a+lXVdyobv86LJlgf4N33WtNN69AlPbN1/VwQ8pGx\nFTaumCG8+gremphXHGfg3FMT5RCyfa+yFMJJi9pnSccIfsDeNmdznO1vhWXczCJq1xPQ67ftNQra\nqC9U3EjErCIYvzCsj1o19auu74LbVi3ta1mrUdUHCPOfD6LqX9vMrIePTuLSzU9on6mtsNn0WGtW\n2SIycdooH3EEqO2HNEshnHcKjCh0lOAPI4sVeFOQUdJBoBKoUTQ53ctximWB7yIZHq3jfk0+Iq9O\nggndxyi4beDcU41jIo7GqBI2wY/RKd1VTE4fx5FjrR+C4DMNpmQOBnOpbOY6b6S8zRA2ykccAWrz\nIc1aCLucr0cEf5Msp75Bj5UsB4FOg1JFvw6tXoKhh/a1aKxHjk0lzteTNVt2HdC6PDLSM1eEzVZM\nedyjCiz/OAkTWn77tP/YYMU11TjLcv0jKjbafBwBqjonLa+eKLiar6fjBb/JfTDtqW8eg8Cksanc\nEVW+6JPTHOu+8/RZDquYlBcm4R5X47P18PHKbaoKAXllRnXnme4nT2w/jnHeHVeFrgt0tOCPUq6x\nLITZNoMfM111saj3nfdioe4+bcw8aRIm3OMIH9u+P6sZyBX1N3R919tdzV1QumwOaWc6WvBHKddY\nFmxsmzaJxKLed95l+3ReTbdcsrAQ4ZXmNW0WJj2t2DYAyo9Oy9547bL4jU6AaOb5Y+XHT0RriOgA\nEb1MRBsU+4mIvt7c/xwRfcT23CKJUq6xLAR971UEi4+kUSkob59llZ//V25coU1pXCZsY04G+2ux\nnp9tjITQvoRq/ERUAfBNAJ8AcBDAs0S0nZlf9B12JYDzmv9dDOAvAVxseW5hmDSrrLJo5oFpkTAo\nFNKaahfhs9yummKUZxL3+bVr3wl2hCZpI6J/D2AjM69u/v15AGDmL/uO+b8A/oaZH2j+fQDA5QAW\nhZ2rIutCLB5FFnTIi7wWXDuhLwXBZdJO0lYD8Ibv74NoaPVhx9QszwUAENE6AOsAYOHChRbNSk4n\nLCzlpdl1Ql8KQrvgzOIuM28FsBVoaPx5XVemvOkhfSkI5cBG8NcBnOP7++zmNptjqhbnCoIgCDli\n49XzLIDziGgxEc0FcBOA7YFjtgP4j03vnksAvMXMv7Q8VxAEQciRUI2fmaeI6DYAuwBUANzLzPuJ\n6Nbm/rsB7ARwFYCXARwF8MemczO5E0EQBMGKjim9KAiC0M5E8erpmEIsgiAIQgMnNX4iGgPw8wQ/\ncTqAX6XUnCyRdqZPWdoq7UyfsrQ1q3aey8x9Ngc6KfiTQkQjtlOeIpF2pk9Z2irtTJ+ytNWFdoqp\nRxAEocMQwS8IgtBhtKvg31p0AyyRdqZPWdoq7UyfsrS18Ha2pY1fEARB0NOuGr8gCIKgQQS/IAhC\nh1E6wU9EFSIaJaIfKPY5UwkspJ23NNv3PBE9Q0TLfftea27fS0S5hC+HtPVyInqr2Z69RHSHb59L\nfTrka+MLRDRNRKc29+Xap2HXc2WcWrTTiXFq0U6XxmhYW90Yp8xcqv8A/BcA3wHwA8W+qwD8EI3y\nq5cA+Mfm9gqAnwH4dwDmAtgH4PwC2/lRAAua/77Sa2fz79cAnO5Qn16u2e5UnwaOuwbAE0X1adj1\nXBmnFu10YpxatNOlMWrdL0WO01Jp/ER0NoCrAdyjOeQ6AN/mBrsB9BLRmQBWAniZmV9h5mMAHmwe\nW0g7mfkZZj7c/HM3GumqC8GiT3U41acBbgbwQFZtSQEnxmkYLo3TmDjVnwoKG6elEvwAvgrgvwI4\nrtkfpRJYlhVDwtrpZy0a2p8HA/gJEe2hRlWyrLFp60ebU/4fEtGy5jYn+5SIegCsAfCwb3PefRp2\nPVfGaZR+KXKc2lzLhTEKWPZL0ePUmQpcYRDRJwG8ycx7iOjyotujI0o7iWgVGi/UZb7NlzFznYh+\nA8CPieglZn6qwLb+M4CFzPwuEV0FYBjAeVm0R0fEZ38NgKeZ+ZBvW259WtD14mLVzqLHqcW1Ch+j\nPmz7pdBxWiaN/1IA1xLRa2hM2T5ORH8dOEZXCcymilie7QQRXYiG2eI6Zv61t52Z683/vwngUTSm\nq1kR2lZmfpuZ323+eyeAKhGdDgf7tMlNCEyfc+5Tm+u5ME6t+sWFcRp2LUfGqFVbfRQ7TvNYSEj7\nP+gXc67G7EWzf2punwPgFQCLcWKRZ1mB7VyIRtGajwa2zwfwAd+/nwGwpuA+PQMnAv1WAni92b9O\n9Wlz3ykADgGYX1Sf2lzPhXFq2c7Cx6llO50Yo7b94sI4LY2pRweVpBJYoJ13ADgNwLeICACmuJGt\n70MAHm1umwPgO8z8ozzbqWjrDQD+hIimAEwAuIkbo9O1PgWA3wXwODMf8R2Wd58qr+fgOLVppwvj\n1KadroxRm7YCDoxTSdkgCILQYZTJxi8IgiCkgAh+QRCEDkMEvyAIQochgl8QBKHDEMEvCILQYYjg\nFwRB6DBE8AuCIHQY/x+Qtlxr6/UDtQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11c8f3748>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.scatter(example_table['redshift'], example_table['cigale_dustlumin'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More advanced queries\n",
    "\n",
    "We have written some code in hershelhelp_python to perform more advanced queries. Spefically it is possible to use the depth map tables to get areas that satisfy some depth criteria.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "#This function requires pymoc\n",
    "from pymoc import MOC\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def get_depth_coverage(coverage_dict):\n",
    "    \"\"\"Generate a Multi-Order Coverage maps corresponding to depth constraints\n",
    "    This method connects to HELP Virtual Observatory server to get the list of\n",
    "    HEALPix cells corresponding to depth constraints. The constraints are\n",
    "    expressed as a dictionnary associated to a band name the minimum depth\n",
    "    value in this band.\n",
    "    This method returns the Multi-Order Coverage maps of the corresponding\n",
    "    area.\n",
    "    Parameters\n",
    "    ----------\n",
    "    coverage_dict : dict\n",
    "        Dictionnary associating to a band name the required minimum depth\n",
    "        value.\n",
    "    Returns\n",
    "    -------\n",
    "    pymoc.MOC\n",
    "        The Multi-Order Coverage map corresponding to the depth constraints in\n",
    "        HELP data.\n",
    "    \"\"\"\n",
    "\n",
    "    moc_list = []\n",
    "\n",
    "    for band in coverage_dict:\n",
    "\n",
    "        query = \"select top 100000000 hp_idx_O_13 from \" \\\n",
    "                \"depth.main where ferr_{}_mean <= {}\".format(\n",
    "                    str(band), float(coverage_dict[band]))\n",
    "\n",
    "   \n",
    "        vo_result = vo.tablesearch(\n",
    "            \"https://herschel-vos.phys.sussex.ac.uk/__system__/tap/run/tap\", query).table\n",
    "\n",
    "        moc = MOC()\n",
    "        if len(vo_result):\n",
    "            #for group in vo_result.group_by('healpix_order').groups:\n",
    "            #    order = group['healpix_order'][0]\n",
    "            #    cells = group['healpix_npix']\n",
    "            #    moc.add(order, cells)\n",
    "            moc.add(13, vo_result['hp_idx_o_13'])\n",
    "        moc_list.append(moc)\n",
    "\n",
    "    result_moc = moc_list.pop()\n",
    "    for moc in moc_list:\n",
    "        result_moc = result_moc.intersection(moc)\n",
    "\n",
    "    return result_moc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#We are using our internal code to convert mad to flux \n",
    "#you can comment this and the next cell and set flux_lim = 4.5\n",
    "#if you don't have hershcelhelp_internal installed\n",
    "from herschelhelp_internal.utils import mag_to_flux\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Suppose we want all the regions with IRAC 1 depth higher than 24mag. First lets get 24mag in uJy\n",
    "flux_lim = mag_to_flux([24])[0] * 1.e6\n",
    "#Lets assume a 5 sigma cut\n",
    "flux_lim = flux_lim * 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 4.5600542])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#The limiting flux an object must have in order to be in the catalogue\n",
    "flux_lim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_moc = get_depth_coverage({'irac_i1':flux_lim})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Folder data must exist\n",
    "test_moc.write('./data/irac_i1_deeper_than_24_MOC.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}