{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# HATLAS-SGP master catalogue\n",
    "\n",
    "This notebook presents the merge of the various pristine catalogues to produce HELP mater catalogue on GAMA-09."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This notebook was run with herschelhelp_internal version: \n",
      "0246c5d (Thu Jan 25 17:01:47 2018 +0000) [with local modifications]\n",
      "This notebook was executed on: \n",
      "2018-02-21 22:15:12.935653\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",
    "\n",
    "from herschelhelp_internal.masterlist import merge_catalogues, nb_merge_dist_plot, specz_merge\n",
    "from herschelhelp_internal.utils import coords_to_hpidx, ebv, gen_help_id, inMoc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "TMP_DIR = os.environ.get('TMP_DIR', \"./data_tmp\")\n",
    "OUT_DIR = os.environ.get('OUT_DIR', \"./data\")\n",
    "SUFFIX = os.environ.get('SUFFIX', time.strftime(\"_%Y%m%d\"))\n",
    "\n",
    "try:\n",
    "    os.makedirs(OUT_DIR)\n",
    "except FileExistsError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I - Reading the prepared pristine catalogues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#atlas = Table.read(\"{}/ATLAS.fits\".format(TMP_DIR))[\"atlas_id\", \"atlas_ra\", \"atlas_dec\",\n",
    "#                                                   \"atlas_stellarity\",\n",
    "#                                                   \"atlas_flag_cleaned\",\n",
    "#                                                   \"atlas_flag_gaia\"]\n",
    "#kids = Table.read(\"{}/KIDS.fits\".format(TMP_DIR))[\"kids_id\", \"kids_ra\", \"kids_dec\", \n",
    "#                                                  \"kids_stellarity\", \n",
    "#                                                  \"kids_flag_cleaned\",\n",
    "#                                                  \"kids_flag_gaia\"]\n",
    "omegacam = Table.read(\"{}/omegacam_merged_hatlas-sgp.fits\".format(TMP_DIR))[\n",
    "                                                   \"omegacam_intid\", \"omegacam_ra\", \"omegacam_dec\",\n",
    "                                                   \"omegacam_flag_merged\",\n",
    "                                                   \"atlas_id\",\n",
    "                                                   \"atlas_stellarity\",\n",
    "                                                   \"atlas_flag_cleaned\",\n",
    "                                                   \"atlas_flag_gaia\",\n",
    "                                                   \"kids_id\", \n",
    "                                                   \"kids_stellarity\", \n",
    "                                                   \"kids_flag_cleaned\",\n",
    "                                                   \"kids_flag_gaia\"]\n",
    "ps1 = Table.read(\"{}/PS1.fits\".format(TMP_DIR))[\"ps1_id\",\"ps1_ra\",\"ps1_dec\",\"ps1_flag_gaia\"]\n",
    "viking = Table.read(\"{}/VISTA-VIKING.fits\".format(TMP_DIR))[\"viking_id\",\"viking_ra\",\"viking_dec\",\n",
    "                                                           \"viking_stellarity\", \n",
    "                                                            \"viking_flag_cleaned\",\n",
    "                                                            \"viking_flag_gaia\"]\n",
    "des =  Table.read(\"{}/DES.fits\".format(TMP_DIR))[\"des_id\",\"des_ra\",\"des_dec\",\n",
    "                                                 \"des_flag_gaia\",\n",
    "                                                 \"des_stellarity\",\n",
    "                                                 \"des_flag_cleaned\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II - Merging tables\n",
    "\n",
    "We first merge the optical catalogues and then add the infrared ones: CFHTLenS, CFHTLS, DECaLS, HSC, KIDS, PanSTARRS, UKIDSS-LAS, VISTA-VHS, and VISTA-VIKING.\n",
    "\n",
    "At every step, we look at the distribution of the distances separating the sources from one catalogue to the other (within a maximum radius) to determine the best cross-matching radius."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Omegacam (ATLAS and KIDS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue = omegacam\n",
    "master_catalogue['omegacam_ra'].name = 'ra'\n",
    "master_catalogue['omegacam_dec'].name = 'dec'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add PanSTARRS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAF3CAYAAABnvQURAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt0nHd95/HPdy7SSLYky7Z8iWXHSXDIhZCYmFzIQtOW\ntEmgpWxpCSnXdpuFhRZ6aHtadk/pdrtb9uwpW2go2SyklC5QWmBpSkMpS9MGAnHiJI5J4qQ4cWLL\nsS1fdNdIc/vuH/OMPJZlayQ9M88zo/frnDmaeebRzDcTDvrk+/vN9zF3FwAAAJYuEXUBAAAArYJg\nBQAAEBKCFQAAQEgIVgAAACEhWAEAAISEYAUAABASghUAAEBICFYAAAAhIVgBAACEhGAFAAAQklRU\nb7x27VrfunVrVG8PAABQs0cfffS4u/fNd15kwWrr1q3atWtXVG8PAABQMzN7sZbzWAoEAAAICcEK\nAAAgJAQrAACAkBCsAAAAQkKwAgAACMm8wcrMMmb2sJk9YWZPmdl/nuMcM7NPmtk+M9tjZq+qT7kA\nAADxVcu4hWlJP+Hu42aWlvQ9M/umuz9Udc4tkrYFt2slfTr4CQAAsGzM27HysvHgYTq4+azT3iTp\n88G5D0laZWYbwy0VAAAg3mraY2VmSTPbLWlQ0rfdfeesUzZJOlj1eCA4BgAAsGzUFKzcvejuV0nq\nl3SNmb1iMW9mZneY2S4z23Xs2LHFvAQAAEBsLehbge4+LOl+STfPeuqQpM1Vj/uDY7N//2533+Hu\nO/r65r3cDgAAQFOp5VuBfWa2KrjfIekmSc/MOu1eSe8Mvh14naQRdz8cerUAAAAxVsu3AjdK+gsz\nS6ocxP7a3b9hZu+VJHe/S9J9km6VtE/SpKT31KleAACA2Jo3WLn7Hknb5zh+V9V9l/T+cEtrfl/c\neWDec26/dksDKgEAAI3A5HUAAICQEKwAAABCQrACAAAICcEKAAAgJAQrAACAkBCsAAAAQkKwAgAA\nCAnBCgAAICQEKwAAgJAQrAAAAEJCsAIAAAgJwQoAACAkBCsAAICQEKwAAABCQrACAAAICcEKAAAg\nJAQrAACAkBCsAAAAQkKwAgAACAnBCgAAICQEKwAAgJAQrAAAAEJCsAIAAAgJwQoAACAkBCsAAICQ\nEKwAAABCQrACAAAICcEKAAAgJAQrAACAkBCsAAAAQkKwAgAACAnBCgAAICQEKwAAgJAQrAAAAEJC\nsAIAAAgJwQoAACAkBCsAAICQEKwAAABCQrACAAAICcEKAAAgJAQrAACAkBCsAAAAQkKwAgAACMm8\nwcrMNpvZ/Wb2tJk9ZWYfnOOcG81sxMx2B7ffq0+5AAAA8ZWq4ZyCpA+7+2Nm1iXpUTP7trs/Peu8\n77r7G8MvEQAAoDnM27Fy98Pu/lhwf0zSXkmb6l0YAABAs1nQHisz2yppu6Sdczz9GjPbY2bfNLPL\nQ6gNAACgqdSyFChJMrOVkr4q6UPuPjrr6cckbXH3cTO7VdLXJW2b4zXukHSHJG3ZsmXRRQMAAMRR\nTR0rM0urHKq+4O5fm/28u4+6+3hw/z5JaTNbO8d5d7v7Dnff0dfXt8TSAQAA4qWWbwWapM9K2uvu\nHz/LORuC82Rm1wSveyLMQgEAAOKulqXAGyS9Q9IPzWx3cOwjkrZIkrvfJektkt5nZgVJWUm3ubvX\noV4AAIDYmjdYufv3JNk859wp6c6wigIAAGhGTF4HAAAICcEKAAAgJAQrAACAkBCsAAAAQkKwAgAA\nCAnBqkGGJ3P68iMHNJ0vRl0KAACoE4JVgzx5aERPDIzo+eMTUZcCAADqhGDVIAPDWUnSwZOTEVcC\nAADqhWDVIAND2dN+AgCA1kOwaoDJXEEnJ3JKJkwHhyZV4mo/AAC0JIJVAxwKulRXbOrRdKGk42PT\nEVcEAADqgWDVAJX9VdddsLr8mOVAAABaEsGqAQZOTmrtyjb1r+5Ueyqhg0NsYAcAoBURrBpgYDir\n/t5OJcy0qbeDjhUAAC2KYFVnI9m8xqYK2rSqQ5K0ubdTh0eyyhdLEVcGAADCRrCqs0PBst/m3o6Z\nnyWXXhqmawUAQKshWNXZwFBWCZM2Bh2r/tWdM8cBAEBrIVjV2cBwVuu7M0onyx91dyatno40G9gB\nAGhBBKs6cncdGsqqP1gGrOhnAzsAAC2JYFVHJydyyuaL6l/Vedrxzb2dOjmR0/h0IaLKAABAPRCs\n6qjSldo0u2O1uvz4EMuBAAC0FIJVHQ0MTSqVMK3vzpx2fNOqDpmkgywHAgDQUghWdTQwnNV5qzqU\nTNhpx9tTSa3vzujgSTpWAAC0EoJVnRSKJb00nD1jGbCisoHd3RtcGQAAqBeCVZ3sOzaufNHVv2ru\nYLW5t1PZfFEvnKBrBQBAqyBY1cmegyOSpP7ezjmfr2xgf+LgcMNqAgAA9UWwqpMnBoaVSSe0ZmXb\nnM+v786oLZnQboIVAAAtg2BVJ3sGRrRpVYcSZnM+nzDTeas69DjBCgCAlkGwqoPpQlHPHBk96zJg\nxebVHdr70qimC8UGVQYAAOqJYFUHew+PKV90bTrLxvWK/t5O5Yol7T081qDKAABAPRGs6uDJQ5WN\n6+cOVuu72yVJLxyfqHtNAACg/ghWdTA4OiUzqbsjfc7zejLl54+OTjWiLAAAUGcEqzoYmsyrpyN9\n1o3rFe3ppFa0JXV0dLpBlQEAgHoiWNXB0GROqzvnHrMw2/ruDB0rAABaBMGqDoYmc1rVee5lwAqC\nFQAArYNgVQdDE3n11tyxatfRMYIVAACtgGBVB8OTOa1a0FLgNBdjBgCgBRCs6mBoMq/eGpcC13Vn\nlCuUNDyZr3NVAACg3ghWIZvKF5XNF9W7oraO1YbujCSxHAgAQAsgWIWs0nmqffN6eUgoIxcAAGh+\nBKuQDU3mJGkBm9eDjhXfDAQAoOkRrEJWCVa1dqz6usodq0GCFQAATY9gFbLKUmCtHatMOqlVnWkd\nIVgBAND0CFYhW+hSoFTewM4eKwAAmt+8wcrMNpvZ/Wb2tJk9ZWYfnOMcM7NPmtk+M9tjZq+qT7nx\nNzSxsKVAqTxygaVAAACaXy0dq4KkD7v7ZZKuk/R+M7ts1jm3SNoW3O6Q9OlQq2wiQ5N5dbYllUkn\na/6d9V3tdKwAAGgB8wYrdz/s7o8F98ck7ZW0adZpb5L0eS97SNIqM9sYerVNYGgyt6BlQKn8zcBj\n49Mqlpi+DgBAM1vQHisz2yppu6Sds57aJOlg1eMBnRm+loXhyfyClgElaX1PRsWS68Q4XSsAAJpZ\nzcHKzFZK+qqkD7n76GLezMzuMLNdZrbr2LFji3mJ2FtUx6qLIaEAALSCmoKVmaVVDlVfcPevzXHK\nIUmbqx73B8dO4+53u/sOd9/R19e3mHpjb1EdK4aEAgDQEmr5VqBJ+qykve7+8bOcdq+kdwbfDrxO\n0oi7Hw6xzqax2D1WEtcLBACg2aVqOOcGSe+Q9EMz2x0c+4ikLZLk7ndJuk/SrZL2SZqU9J7wS42/\nYsk1ks2rd4Edq7Ur25Qw6egIwQoAgGY2b7By9+9JsnnOcUnvD6uoZjWazctdWrXAjlUqmdDalYxc\nAACg2TF5PUQzU9dXLKxjJZWXA1kKBACguRGsQjQUXCdwoR0rSVrfTccKAIBmV8seK8zhizsPnHFs\n7+HyFIpH9p/U4eGFdZ/WdWf02IHhUGoDAADRoGMVoslcUZLU2bbwvLqhO6OTEzlNF4phlwUAABqE\nYBWiyVxBktTZVvt1AivWd5eHhB4bYzkQAIBmRbAK0WSuqIRJ7amFf6zrZoaEEqwAAGhWBKsQTeaK\n6mhLqTxTdWHWd5WD1SDT1wEAaFoEqxBN5gqLWgaUTi0FHiFYAQDQtAhWIZrMFRcdrFavaFM6aSwF\nAgDQxAhWIcrmiov6RqAkmZnWdWVYCgQAoIkRrEK0lKVAKRgSyvR1AACaFsEqJO6+pKVAqXxZmyNc\niBkAgKZFsApJvugqlFwrFrkUKJWD1SB7rAAAaFoEq5BMLGE4aMX67ozGpguamC6EVRYAAGggglVI\nTl3OZml7rCRpkOnrAAA0JYJVSCqXs+lY4lKgJB3lm4EAADQlglVIwuxYEawAAGhOBKuQhBOs6FgB\nANDMCFYhmZzZvL74pcCV7Sl1tiWZvg4AQJMiWIVkMldUeyqhZGLhF2CuMDOt787QsQIAoEkRrEKS\nXeJw0Ip1Xe3MsgIAoEkRrEJSvpzN4pcBK9Z3Z3SEjhUAAE2JYBWSpV7OpmJDT3kp0N1DqAoAADQS\nwSokk7miVrQvvWO1rqtd04WSRrNMXwcAoNkQrEIyMV1QRwgdq5mRC2MsBwIA0GyW3mKBiiXXdKG0\nqKXAL+48cNrj/ccnJElf2nlA29Z3zRy//dotSysSAADUHR2rEIQxw6qipyMtSRqdyi/5tQAAQGMR\nrEIQxtT1iq5MOZyNTbHHCgCAZkOwCkGYwSqdTCiTTmiUYAUAQNMhWIUgG+JSoCR1ZdIaYykQAICm\nQ7AKQZgdK0nqzqRYCgQAoAkRrEIQdrDqyqTZvA4AQBMiWIVgMldQMmFqS4bzcVY6VkxfBwCguRCs\nQlC5nI2ZhfJ6XZm0iiVXNuiEAQCA5kCwCsFErqgVIW1cl06NXBidZp8VAADNhGAVgslcOJezqejK\nlIeE8s1AAACaC8EqBJWlwLB0V4aEciFmAACaCsEqBOVgFeZSIB0rAACaEcFqidxd2Vwh1I5VW4rp\n6wAANCOC1RJNF0oqeXgzrCq62pm+DgBAsyFYLdGp4aDhLQVK5W8G0rECAKC5EKyWaHLmOoHhdqy6\nO+hYAQDQbAhWSxT25WwqutqZvg4AQLOZN1iZ2T1mNmhmT57l+RvNbMTMdge33wu/zPg61bEKeSmw\nI61CyTWVL4X6ugAAoH5qSQOfk3SnpM+f45zvuvsbQ6moyVQ6VivC7lhVpq9P5UMdPgoAAOpn3o6V\nuz8g6WQDamlKE9NFmaRM2HusZmZZsYEdAIBmEdYeq9eY2R4z+6aZXR7SazaFbL6gTDqpREgXYK6o\ndKzYwA4AQPMIY2PQY5K2uPu4md0q6euSts11opndIekOSdqyZUsIbx29bK5Yl6W6U0uBdKwAAGgW\nS+5Yufuou48H9++TlDaztWc592533+HuO/r6+pb61rGQzRfVkQ4/WLWnkmpPJehYAQDQRJYcrMxs\ng1l5HczMrgle88RSX7dZ1KtjJZWvGUjHCgCA5jHvUqCZfUnSjZLWmtmApI9KSkuSu98l6S2S3mdm\nBUlZSbf5Mhq+lM0XtaqzrS6v3ZVJ0bECAKCJzBus3P1t8zx/p8rjGJalbK4+S4GS1J1J6eBQti6v\nDQAAwsfk9SVw9/IeqzouBY5N5Zm+DgBAkyBYLUGuWFLJVdeOVb7I9HUAAJoFwWoJssHU9XoFq66Z\nIaHsswIAoBkQrJYgmy8Hq7CnrlcwywoAgOZCsFqCSrCq31IgHSsAAJoJwWoJpipLgXXuWHG9QAAA\nmgPBagnq3bFqTyfVxvR1AACaBsFqCeq9eV2SutpT7LECAKBJEKyWIJsvyiS1p+v3MXZ3pOlYAQDQ\nJAhWS5DNF5VJJ5UoXyqxLsqXtaFjBQBAMyBYLUE9L8Bc0Z1Ja5Tp6wAANAWC1RJk8/W7TmBFVzB9\nfXyarhUAAHFHsFqCel6AuaIyff3o6HRd3wcAACwdwWoJsvlS3aauV1RmWQ2OTdX1fQAAwNIRrJag\nEUuBlenrg3SsAACIPYLVIrm7phqyFEjHCgCAZkGwWqR80VV0r/u3AttTCbUlE+yxAgCgCRCsFqly\nOZvOOneszExdmZQGxwhWAADEHcFqkSqXs6n35nWp/M3AwVGWAgEAiDuC1SLV+wLM1ehYAQDQHAhW\ni5TNlQd21nuPlSR1Z1J0rAAAaAIEq0VqbMcqrYlckenrAADEHMFqkSp7rBq1FCiJrhUAADFHsFqk\nbL4ok9Serv9H2N3BZW0AAGgGBKtFyuaLyqSTSpjV/b262hkSCgBAMyBYLVI2V2zIxnXpVMeKy9oA\nABBvBKtFasR1AivaUwl1pJM6yh4rAABijWC1SNkGXCewwsy0qbdDh4azDXk/AACwOASrRcrmSw2Z\nul7R39uhgSGCFQAAcUawWqRGLgVKlWA12bD3AwAAC0ewWgR311QDlwIlqb+3U0OTeYaEAgAQYwSr\nRcjmiyq6N+xbgVK5YyWJrhUAADFGsFqEkWxeUmOmrlf093ZKkgZOss8KAIC4IlgtwkywomMFAACq\nEKwWYWSy8R2rNSva1JFO8s1AAABijGC1CFEsBZoZIxcAAIg5gtUiRLEUKAUjF4ZZCgQAIK4IVosQ\nRcdKKm9gP8jmdQAAYotgtQij2bxMUnu6sR9ff2+HRrJ5jU7lG/q+AACgNgSrRRjJ5tWeTihh1tD3\nrYxcOMQ+KwAAYolgtQgj2XzDlwElafPqysgFghUAAHFEsFqEkWy+4RvXpaohocyyAgAglghWizCS\nzasznWr4+/Z2ptXZxiwrAADiimC1CMPZvDIRdKwqs6wOnqRjBQBAHM0brMzsHjMbNLMnz/K8mdkn\nzWyfme0xs1eFX2a8jEa0x0oqLwfSsQIAIJ5q6Vh9TtLN53j+Fknbgtsdkj699LLiy90j27wuBUNC\n2WMFAEAszRus3P0BSSfPccqbJH3eyx6StMrMNoZVYNxk80Xlix7J5nVJ2tzbqdGpwsyQUgAAEB9h\n7LHaJOlg1eOB4FhLimrqekV/b3nkArOsAACIn4ZuXjezO8xsl5ntOnbsWCPfOjRRXSewgpELAADE\nVxjB6pCkzVWP+4NjZ3D3u919h7vv6OvrC+GtG29kMh4dKzawAwAQP2EEq3slvTP4duB1kkbc/XAI\nrxtLUS8FrupMa0VbUgfpWAEAEDvzTrk0sy9JulHSWjMbkPRRSWlJcve7JN0n6VZJ+yRNSnpPvYqN\ng6iXAsuzrBi5AABAHM0brNz9bfM875LeH1pFMRd1x0qqjFwgWAEAEDdMXl+g0WxeZlJ7OrqPbvPq\nTjavAwAQQwSrBRrJ5tXVnlLCLLIa+ns7NMYsKwAAYodgtUAj2bx6OtOR1nDqm4F0rQAAiBOC1QKN\nZPPq6Yg6WJVnWR08yT4rAADihGC1QPEIVnSsAACII4LVAsUhWPV0pLWyPcU3AwEAiBmC1QKNZAuR\nB6vyLCtGLgAAEDcEqwVwd41m8+qOOFhJCoaEshQIAECcEKwWYCpfUq5YirxjJZX3WR0ayqo8nxUA\nAMQBwWoBKnOj4hKsxqYLGs0Woi4FAAAECFYLEK9gFYxcYDkQAIDYIFgtQLyCFSMXAACIG4LVAgxP\n5iTFI1htDjpWfDMQAID4IFgtQJw6Vj2daXVlmGUFAECcEKwWoBKsVnW0RVxJ2ebeTu0/PhF1GQAA\nIJCKuoBmMprNy0zqyjT+Y/vizgNnHGtPJfT4gaGZ526/dkujywIAAFXoWC3ASDavrvaUEgmLuhRJ\n0oaejEanCpqcZuQCAABxQLBagJFsXj2d0e+vqtjQnZEkHRmdirgSAAAgEawWZGgyH5v9VVK5YyUR\nrAAAiAuC1QIMjk1rXVd71GXM6MqktaI9pSMjBCsAAOKAYLUAx8amtK47PsFKkjZ2Z3SYYAUAQCwQ\nrGqUL5Z0YiKndV2ZqEs5zYaejI6OTqnExZgBAIgcwapGx8en5a7Ydaw2dGdUKLlOjOeiLgUAgGWP\nYFWjwdFpSYplx0piAzsAAHFAsKrR4Fg5WK2PWceqr6tdCZOOjHBpGwAAokawqtHRoCMUt45VOpnQ\n2pXtfDMQAIAYIFjVaHBsWmbS2pXxmWNVsaEno8MsBQIAEDmCVY0GR6e0ZkW7Usn4fWQbuzMansxr\ndCofdSkAACxr8UsJMRW34aDVKhvYnz0yFnElAAAsbwSrGg3GcDhoxYaeDknSM4dHI64EAIDljWBV\no6Oj01ofs43rFd2ZlDrSSe2lYwUAQKQIVjUollwnxqdj27EyM23oydCxAgAgYgSrGpwYn1bJFds9\nVlJ5AvuzR8ZUKnFpGwAAokKwqsHRytT17nguBUrlDewTuaIODk1GXQoAAMsWwaoGg2OV4aDx7lhJ\n0t7D7LMCACAqBKsaVC5nE+eO1frujMykZ46wzwoAgKgQrGpQuZxN38r4dqzaUgldsGaFnqFjBQBA\nZAhWNRgcm9bqFW1qS8X747pkYxcdKwAAIhTvpBATg6Pxnbpe7ZIN3Xrx5KQmpgtRlwIAwLJEsKpB\neep6fPdXVVyyoUvu0rNHWQ4EACAKBKsaNEvH6tKN3ZLEPisAACJCsJpHqeQ6Nt4cwWrTqg6tbE+x\nzwoAgIgQrOZxYiKnYsm1vgmWAhMJ08s3dOnplwhWAABEoaZgZWY3m9mzZrbPzH5njudvNLMRM9sd\n3H4v/FKj0QzDQavtOL9XTwwMs4EdAIAIzBuszCwp6VOSbpF0maS3mdllc5z6XXe/Krj9Qch1RqYZ\nhoNWe93FfcoXXTv3n4i6FAAAlp1aOlbXSNrn7s+7e07SX0l6U33Lio/B0ebqWF19fq8y6YQe+Nfj\nUZcCAMCyU0uw2iTpYNXjgeDYbK8xsz1m9k0zuzyU6mJgMLgAc1+TBKtMOqlrL1ijB350LOpSAABY\ndsLavP6YpC3u/kpJfyrp63OdZGZ3mNkuM9t17Fhz/OEfHJvWqs60Mulk1KXU7LXb1ur5YxMaGJqM\nuhQAAJaVWoLVIUmbqx73B8dmuPuou48H9++TlDaztbNfyN3vdvcd7r6jr69vCWU3ztHRqaZZBqz4\nsYvLn+33fsRyIAAAjVRLsHpE0jYzu8DM2iTdJune6hPMbIOZWXD/muB1W2L39ODYtNZ1NcfG9YqX\nrVupDd0ZlgMBAGiw1HwnuHvBzD4g6VuSkpLucfenzOy9wfN3SXqLpPeZWUFSVtJt7u51rLthjo1N\n68K+FVGXsSBmptduW6tvPXVExZIrmbCoSwIAYFmYN1hJM8t79806dlfV/Tsl3RluadFz9/J1Apus\nYyVJr724T3/z6ID2DAxr+5beqMsBAGBZYPL6OQxN5pUvetPtsZKkf/OytTKTvss+KwAAGoZgdQ6V\nqevNcDmb2VavaNMVm3r0wL+yzwoAgEYhWJ3D0dHK1PXm61hJ5bELjx8c1uhUPupSAABYFghW59Bs\nU9dne922PhVLrh881xJf0AQAIPYIVucwc53AJty8Lknbt/RqRVuS5UAAABqEYHUOg6NT6sqk1NHW\nPFPXq7WlErr+ojVsYAcAoEEIVudQHg7anMuAFa/d1qcDJyf14omJqEsBAKDlEazOYXBsuim/EVjt\ndcHlbR6gawUAQN3VNCB0uTo6OqUd5zfPcM0v7jxwxjF3V29nWl/ceUBJM91+7ZYIKgMAYHmgY3UW\n5anr01rX5B0rM9O2dV16bnBcU/li1OUAANDSCFZnMZotKFcoNf0eK0nasbVXuWJJjx0YiroUAABa\nGsHqLI4GU9ebvWMlSf29ndrc26EfPHdCpVJLXBsbAIBYIlidxWBl6noLdKwk6fqL1ujERE4P/IiZ\nVgAA1AvB6iya+TqBc3nFph51taf0F99/IepSAABoWQSrszjaYh2rVCKhV1+wWvc/e0z7jzPTCgCA\neiBYncXg2JRWtCW1or11JlJcc8FqpRKmz//ghahLAQCgJRGszqIVhoPO1p1J6w2v3Kiv7BrQxHQh\n6nIAAGg5BKuzGBydUl+LLANWe9drtmpsuqCvPTYQdSkAALQcgtUcSiXXjwbHdf6azqhLCd32zav0\nyv4efe77L8id0QsAAISJYDWH54+Pa3gyrx3nr466lNCZmd51/VY9d2xC39vH9QMBAAgTwWoOj75Y\nnlD+qia6TuBCvPHKjVq7so3RCwAAhIxgNYddLwyptzOti/pWRF1KXbSnkrr9mi36zjODMyESAAAs\nHcFqDo++OKSrz++VmUVdSt3c8WMX6byeDv32V57g4swAAISEYDXLifFpPX98Qle34P6qaivbU/rY\nz1+h545N6E/+34+iLgcAgJZAsJqlsjS2Y2tr7q+q9tptfbrt1Zt19wPPaffB4ajLAQCg6RGsZnn0\nxSG1JRO6YlNP1KU0xEfecKnWd2f0W3/zhKYLLAkCALAUBKtZdr04pFds6lYmnYy6lIbozqT1R//2\nCv1ocFyf/A5LggAALEXrXAgvBFP5on44MKJ337A16lLq5os7D8x5/Ootvfr0Pz8nuem3bn55g6sC\nAKA10LGq8uShEeWKJV3dovOrzuXWKzZqZXtKX3nsoLI5lgQBAFgMglWVXcHG9eUYrDraknrz9n4N\njk7rlz/3iCZzXKQZAICFIlhV2fXCkC5Yu0JrV7bexZdr8fINXfqFHf3auf+E3v3nj2h8mnAFAMBC\nEKwC7q7HDgwty25Vtas29+oTt23Xoy8O6d33PKyxqXzUJQEA0DQIVoHnj0/o5EROO5Z5sJKkn7ny\nPN35tu3afXBY77znYY0SrgAAqAnBKvDoC8tnMGgtbrlioz71S6/Sk4dG9PbP7NTg6FTUJQEAEHsE\nq8CuF09qVWdaF65dGXUpsfHTl2/QXW+/Ws8eGdPrP/4v+utdB+XuUZcFAEBsEawCu14c0tVbepVI\ntO6FlxfjJy9dr/s++Fq9fEOXfvsre/TOex7WwZOTUZcFAEAsMSBU0smJnJ4/NqGff1V/1KXEwlxD\nRN901SZt7OnQPzx1RD/5x/+in7p8vT5x23YlCaIAAMwgWEl6rHLhZTaun1XCTNdduEaXbOjS13cf\n0jf2HNYTA8N692su0C/u6FdXJh11iQAARI6lQJWXAdNJ05WbV0VdSuyt6mzTu67fqrdfu0UbujP6\nL994Wtf/0T/pD/7uaZYIAQDLHh0rSY++eFKXn9ezbC68vFRmpsvO69EfvvkK7RkY1me/t1+f/8EL\n+tz39+v6i9bopkvX6/WXrVd/b2fUpQJAy3F35YuufLEU3FyFUkmFoitXLP/MF0sqlsrH80VXoeqc\nQqmkQqlyzFWselwsuUpevhVLUsldCTOlk6Z0MhHcTD0dafX3dmpTb4d6OlixqLbsg9X9zwxq14tD\n+g83XhQk1lEgAAAODUlEQVR1KU2nshfr2gvW6JIN3Xp4/wk9+dKoHtx3Qr//d09rY09Gv3B1v157\ncZ+u2ERwBRA/7uUwkiuUA0iuUL6fKxY1XblfKCkXhJhcoTrQlJQruorFIJiUXIXgfjF4XCr5zONK\n2MlXhZt84VQ4ylW/bnU9leMz9+P17exMOqHezjat787o4vUrdfG6LnW2zx0vbr92S4Ora7xlHaz+\n9eiYfu1Lj+uyjd16/4+/LOpymlpPR1o3XbZBN122QcfHp7X38Kj2Hh7Vnffv0yf/aZ9SCdOlG7u1\nfcsqbd+ySlds6tH5a1YonWQ1Gogzd1fJpeJMWCjNhIZK96MSFqqfn/34jJt71WueHkCqA0q5i1JS\nPjhW3Z05/X4QVIIQUii68qcFl1nhpc4hxSQlEqaElfeoJsyUSJiSVjlevqUSpuSsWypR7g5l0pXn\nEzPHq89JWNXvBfcTweOESUmrflx1vPI4eD5xWk3les0kU/mxSzP/Xir/ziZzBQ1N5jU8mZv5+aPB\nce0+OCyTtHl1p16+oUuXbuzWhu5MXT7juLKo5hLt2LHDd+3aFcl7S+VvAv7cpx7UZK6oez9wg85b\n1bGg35/rm3M40+R0QS+enNTBk5M6MDSpgaGscoWSJClh0pqV7bp6S69etm6lLuxboU2rOrSpt0Mb\nujNKEboQc9VhoOiuYrH8s1AqqVTSacfmDhalmeWX6lBRWZopVoWX6qWb6tBRWeoph4nSrDBRWRo6\n834lbFR3UM4WmOKgHAJ0RphInBYqdFrASFU9d/otcUZIqfxMJRJKJiv3zzw/lZzjNWcCyqkaE7b8\nvjFdctehoayePTqmZ4+M6dBwVpK0ubdD116wRlf09+hdr9kabZFLYGaPuvuOec9bjsEqVyjpHZ/d\nqccPDuvLd1yn7VsW/m1AgtXilNx1dHRKR0amNDg2rWNj05oqFPXiiUkVS6f+t5hMmDZ0Z3TeqozW\nrmzX2pXtWrOyLbjfpu6OtLozafUEP1dmUox+iFCpVP6jXZz1R3n2H+kzwkNVqDi1/2OOJZVZyyuV\n96mEk2JJMz9LXnm/6uBz6vfzxVOPK52NM95/zvcrndZpKZRccZiXa9KpToVVBYqzhopZ3Y2qDsbM\na8zuYMwcLz+2me6HKZmYoyNTHTROe3yqG5KwMzskM3XMClHLMaS0grGpvPYMjGjn/pM6Pj6tjnRS\nb7tmi26/dotetq75hnHXGqxqWgo0s5slfUJSUtJn3P1js5634PlbJU1Kere7P7bgqhvA3fXRe5/S\nzv0n9T/feuWiQhUWL2GmjT0d2thzeoewUCzNtJOHJ/MaypZ/Hh/P6YUTk8oXSxqePPc1CzvSSXW2\nJdXZntSKtpQ62pLKpJLKpBPKpJPBLaH2VFJtqYTakonyz9SpDZmpREKpZPV/wSZm/nBU/hCZTv/j\nYFL5r1s1L7fPK394XeU/wuVjLp85Jzhe9VwpeLLkCjaQls+p/FEvVX76qeAw89NPLZ2cEQxKp4LP\nnMdndUjmWmYpnKXrUYowYFT+0FcvuZz6w336skvitD/eswJJ8O84kzIl2pKndyHO9nqJMwPCqXrm\nCBKJuc9LJk6vfyZUzP6dmfBE6ED8dWXSuuFla/Wai9Zo//EJ7dx/Un/50Au658H9emV/j372yvP0\nM1eep/UttlQ4b7Ays6SkT0m6SdKApEfM7F53f7rqtFskbQtu10r6dPAzVool158/uF9feviA3nfj\nRXrzdgaCxkUqmVBfV7v6utrPek6hVNLkdFHj0wVN5YuayheVzZeCn0XlCyVNF09tNh3LFnSymFN+\n5lsxpdOXPIpBwGlRycof9Ko/zMnZQaGqQ1EJJ3OFjc621Bkdikr4rOzhSM5+/eCYzepqJBNBd2KO\n4FJ5vdP2h5wRSE4FC1O5bgDxZWa6sG+lLuxbqZsuW6+vP35I9z7xkv7w7/fqv963V9dfuEZvfOV5\nuvr8Xl3Ut6Lpt4HU0rG6RtI+d39ekszsryS9SVJ1sHqTpM97eV3xITNbZWYb3f1w6BUvwODolB4/\nOKzHDwxr98Eh7RkY0WSuqNdfuk6/9VMvj7I0LEIqkVB3R0LdIX61d67NtKWSz3SKSjrVRZrpJgUd\np5ku01niWdDLmmlmVf7+W/DgzON26n4QdCy4n1BVh6x6c6nNCkbLeH8HgPjr62rXr77uQv3q6y7U\nc8fGde/ul/S3uw/pI//3h5Kk9lRCl2zs1uXndeuyjd3q62qf2fLR01ne/tGWLO95i+sl6GoJVpsk\nHax6PKAzu1FznbNJUmTB6tEXh/Tzn/6+JCmVMF12Xrd+4ep+bd/Sq5tfsSG2/0LQWJU9JwCAxrqo\nb6V+46aL9aHXb9O+wXH98NCInnppVE+9NKJvPPFSTXuZKwHr9mu26Pd/9vIGVD2/ho5bMLM7JN0R\nPBw3s2cb9d7PSfq7cF9yraTj4b4kzoHPu/H4zBuLz7vx+Mwb7Jfq9Lr/ObjV2fm1nFRLsDokaXPV\n4/7g2ELPkbvfLenuWgqLOzPbVcu3AxAOPu/G4zNvLD7vxuMzRz3UskPsEUnbzOwCM2uTdJuke2ed\nc6+kd1rZdZJGot5fBQAA0GjzdqzcvWBmH5D0LZXHLdzj7k+Z2XuD5++SdJ/Koxb2qTxu4T31KxkA\nACCeatpj5e73qRyeqo/dVXXfJb0/3NJiryWWNJsIn3fj8Zk3Fp934/GZI3SRTV4HAABoNc09hQsA\nACBGCFYLZGY3m9mzZrbPzH4n6npanZndY2aDZvZk1LUsB2a22czuN7OnzewpM/tg1DW1OjPLmNnD\nZvZE8Jk34FvjMLOkmT1uZt+Iuha0FoLVAlRd3ucWSZdJepuZXRZtVS3vc5JujrqIZaQg6cPufpmk\n6yS9n/+N1920pJ9w9yslXSXp5uDb1aivD0raG3URaD0Eq4WZubyPu+ckVS7vgzpx9wcknYy6juXC\n3Q9XLqDu7mMq/+HZFG1Vrc3LxoOH6eDG5tc6MrN+SW+Q9Jmoa0HrIVgtzNku3QO0HDPbKmm7pJ3R\nVtL6gmWp3ZIGJX3b3fnM6+tPJP22pFLUhaD1EKwAnMHMVkr6qqQPufto1PW0OncvuvtVKl+14hoz\ne0XUNbUqM3ujpEF3fzTqWtCaCFYLU9Ole4BmZmZplUPVF9z9a1HXs5y4+7Ck+8W+wnq6QdLPmtkL\nKm/n+Akz+z/RloRWQrBamFou7wM0LTMzSZ+VtNfdPx51PcuBmfWZ2argfoekmyQ9E21Vrcvdf9fd\n+919q8r/H/5P7v72iMtCCyFYLYC7FyRVLu+zV9Jfu/tT0VbV2szsS5J+IOnlZjZgZr8SdU0t7gZJ\n71D5v+J3B7dboy6qxW2UdL+Z7VH5P96+7e6MAACaFJPXAQAAQkLHCgAAICQEKwAAgJAQrAAAAEJC\nsAIAAAgJwQoAACAkBCsAAICQEKwAzMvMisFMq6fM7Akz+7CZJYLndpjZJ8/xu1vN7PbGVXvGe2eD\n6/DFgpm91cz2mRmzqoAWRLACUIusu1/l7perPBn8FkkflSR33+Xuv36O390qKZJgFXguuA5fzcws\nWa9i3P3Lkv5dvV4fQLQIVgAWxN0HJd0h6QNWdmOl+2JmP1Y1sf1xM+uS9DFJrw2O/UbQRfqumT0W\n3F4T/O6NZvbPZvYVM3vGzL4QXGJHZvZqM/t+0C172My6zCxpZv/DzB4xsz1m9u9rqd/Mvm5mjwbd\ntzuqjo+b2R+b2ROSrj/Le14e3N8dvOe24HffXnX8f1WCmZndHPwzPmFm3wnxXwOAmEpFXQCA5uPu\nzwfhYd2sp35T0vvd/UEzWylpStLvSPpNd3+jJJlZp6Sb3H0qCCZfkrQj+P3tki6X9JKkByXdYGYP\nS/qypLe6+yNm1i0pK+lXJI24+6vNrF3Sg2b2j+6+f57yf9ndTwbX5XvEzL7q7ickrZC0090/HFwL\n9Jk53vO9kj7h7l8Izkma2aWS3irpBnfPm9mfSfolM/umpP8t6XXuvt/MVi/4gwbQdAhWAML0oKSP\nm9kXJH3N3QeCplO1tKQ7zewqSUVJF1c997C7D0hSsC9qq6QRSYfd/RFJcvfR4PmfkvRKM3tL8Ls9\nkrZJmi9Y/bqZvTm4vzn4nRNBLV8Njr/8LO/5A0n/0cz6g3++H5nZT0q6WuWQJkkdkgYlXSfpgUrQ\nc/eT89QFoAUQrAAsmJldqHIQGZR0aeW4u3/MzP5e0q0qd5B+eo5f/w1JRyVdqfJ2hKmq56ar7hd1\n7v+PMkm/5u7fWkDdN0p6vaTr3X3SzP5ZUiZ4esrdi+f6fXf/opntlPQGSfcFy48m6S/c/XdnvdfP\n1FoXgNbBHisAC2JmfZLuknSnz7qKu5ld5O4/dPf/LukRSZdIGpPUVXVaj8rdoJKkd0iab6P4s5I2\nmtmrg/foMrOUpG9Jep+ZpYPjF5vZinleq0fSUBCqLlG5q1TzewaB8nl3/6Skv5X0SknfkfQWM1sX\nnLvazM6X9JCk15nZBZXj89QGoAXQsQJQi45gaS4tqSDpLyV9fI7zPmRmPy6pJOkpSd8M7heDTeGf\nk/Rnkr5qZu+U9A+SJs71xu6eM7O3SvrTYF9UVuWu02dUXip8LNjkfkzSz83zz/EPkt5rZntVDk8P\nLfA9f1HSO8wsL+mIpP8W7Nf6T5L+0cojKPIq7zN7KNgc/7Xg+KDK36gE0MJs1n9wAkDLMLOtkr7h\n7q+IuJTTBEuSMxv6AbQOlgIBtLKipB6L2YBQlbt2Q1HXAiB8dKwAAABCQscKAAAgJAQrAACAkBCs\nAAAAQkKwAgAACAnBCgAAICT/H7PBJQSXCwLlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f34f563fc50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nb_merge_dist_plot(\n",
    "    SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n",
    "    SkyCoord(ps1['ps1_ra'], ps1['ps1_dec'])\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Given the graph above, we use 0.8 arc-second radius\n",
    "master_catalogue = merge_catalogues(master_catalogue, ps1, \"ps1_ra\", \"ps1_dec\", radius=0.8*u.arcsec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add VIKING"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAF3CAYAAACfXf7mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUXWd55/vfc8aaq1RSaZYsG8+WsTHCJthgGwKRmRwC\nCTZT31xYCgk0STqdNKzkQid09+1e3GalGR2HdjukPaQTMHiBjYEkYMADHjC2ZWwQwoPmoVRznfm5\nf5y9S8dSlepU1T5jfT9rnVXn7H2Gt44t+6fnfffzmrsLAAAAtRVr9AAAAACWA0IXAABAHRC6AAAA\n6oDQBQAAUAeELgAAgDogdAEAANQBoQsAAKAOCF0AAAB1QOgCAACoA0IXAABAHSQaPYDZrFq1yrds\n2dLoYQAAAMzrkUceOeLuQ/M9rylD15YtW/Twww83ehgAAADzMrPnqnke04sAAAB1QOgCAACoA0IX\nAABAHRC6AAAA6oDQBQAAUAfzXr1oZjdJerOkQ+6+dZbzfyrp3RXvd56kIXcfNrNnJY1LKkoquPu2\nqAYOAADQSqqpdN0saftcJ939U+5+sbtfLOljkr7v7sMVT7k6OE/gAgAAy9a8ocvd75U0PN/zAtdL\num1JIwIAAGhDka3pMrMulStiX6k47JK+a2aPmNmOqD4LAACg1UTZkf4tkn50wtTiFe6+18xWS/qO\nmT0dVM5OEoSyHZK0efPmCIcFAADQeFFevXidTphadPe9wc9Dku6QdOlcL3b3G919m7tvGxqad/si\nAACAlhJJ6DKzfklXSvp6xbFuM+sN70t6g6Qno/g8AACAVlNNy4jbJF0laZWZ7ZH0CUlJSXL3G4Kn\nvU3St919suKlayTdYWbh59zq7t+KbugAAACtY97Q5e7XV/Gcm1VuLVF5bLekixY7sFZz64PPV/W8\nd13GejUAAJYjOtIDAADUAaELAACgDghdAAAAdUDoAgAAqANCFwAAQB0QugAAAOqA0AUAAFAHhC4A\nAIA6IHQBAADUAaELAACgDghdAAAAdUDoAgAAqANCFwAAQB0QugAAAOqA0AUAAFAHhC4AAIA6SDR6\nAMvNrQ8+P+9z3nXZ5jqMBAAA1BOVLgAAgDogdAEAANQBoQsAAKAOCF0AAAB1QOgCAACoA0IXAABA\nHRC6AAAA6oDQBQAAUAeErjpyd/1w1xGNTucbPRQAAFBnhK46GpnK664n9uvh54YbPRQAAFBnhK46\nGp7KlX9O5Bo8EgAAUG+ErjoaniyHrSMT2QaPBAAA1Buhq47C0HV0kkoXAADLDaGrjsLQNZUrajpX\nbPBoAABAPRG66ujYVE4xK98/OskUIwAAywmhq46GJ3PaPNglSTrCYnoAAJYVQledZPJFTeWKOnN1\nj0zSURbTAwCwrBC66iRcz7W6t0P9nUkW0wMAsMwQuuokDF2D3SkN9qSodAEAsMzMG7rM7CYzO2Rm\nT85x/iozGzWzx4LbxyvObTezZ8xsl5l9NMqBt5rK0LWqO82aLgAAlplqKl03S9o+z3N+4O4XB7e/\nkiQzi0v6vKRrJJ0v6XozO38pg21lw1M5dSbj6kjGtbInpel8UVO5QqOHBQAA6mTe0OXu90pazGaB\nl0ra5e673T0n6XZJ1y7ifdrCscmcBrtTkqRVPWlJ0lGqXQAALBtRrel6lZk9bmZ3m9kFwbENkl6o\neM6e4NiyNFwRusKf9OoCAGD5iCJ0PSpps7u/VNJnJX1tMW9iZjvM7GEze/jw4cMRDKt5lNx1bOrF\noavcNoJKFwAAy8WSQ5e7j7n7RHD/LklJM1slaa+kTRVP3Rgcm+t9bnT3be6+bWhoaKnDaiqj03mV\nXBrsKoeuZDym/i7aRgAAsJwsOXSZ2Vozs+D+pcF7HpX0kKSzzOx0M0tJuk7SnUv9vFYUXrm4Iqh0\nSdLK7pSO0DYCAIBlIzHfE8zsNklXSVplZnskfUJSUpLc/QZJ75D0+2ZWkDQt6Tp3d0kFM/uwpHsk\nxSXd5O47a/JbNLljFe0iQit70npiz2ijhgQAAOps3tDl7tfPc/5zkj43x7m7JN21uKG1j+HJ8kbX\n/Z3JmWOruo+3jehKzfuPAQAAtDg60tfB8FROA10pxWM2c2wlbSMAAFhWCF11UNkuIrQyeMy6LgAA\nlgdCVx0MT+ZmrlwMzbSN4ApGAACWBUJXjWXyRU3liidVuhLxmAa6kmx8DQDAMkHoqrFjUye3iwit\n7E5T6QIAYJkgdNXY8CztIkIre1IspAcAYJkgdNXYTOjqmi10pcttI7KFeg8LAADUGaGrxoYnc+pM\nxtWZip90buYKRqYYAQBoe4SuGputXURoZU/5OIvpAQBof4SuGhuezM26iF4qTznSNgIAgOWB0FVD\nJXeNTOVnXc8lHW8bQYNUAADaH6Grhsam8yq6z6zdms3KnjRXMAIAsAwQumoovHJxrulFqbyY/uhk\nVu5er2EBAIAGIHTV0Kl6dIVW9qSVyZc0lSvWa1gAAKABCF01NDyVU8yk/s7knM9Z1c0VjAAALAeE\nrhoansxpoCuleMzmfM7KnrQkrmAEAKDdEbpq6Nhkbs4rF0MrupMySUdYTA8AQFsjdNXQqXp0hRKx\nmLrTCU1k83UaFQAAaARCV41k80VN5oqnXEQf6krFWUgPAECbI3TVyLHpcuVqRdfci+hDhC4AANof\noatGpoMQ1ZVKzPvcrlRCU7lCrYcEAAAaiNBVI9lCOXR1JOf/iql0AQDQ/ghdNZItlCRJqUT1oYuu\n9AAAtC9CV41k8+XQ1ZGIz/vcrlRCxZIrVyzVelgAAKBBCF01Ek4vpqusdEliihEAgDZG6KqRbKEk\nk5QkdAEAABG6aiabLyqViClmc28BFOoMrnDkCkYAANoXoatGsoVSVVOLEpUuAACWA0JXjWQKJaWT\n8y+ilwhdAAAsB4SuGskViguodDG9CABAuyN01UgmX6qqXYQkxWOmdCJGpQsAgDZG6KqRbKFYVWPU\nUFcqPrN1EAAAaD+ErhrJFkpVbQEUYv9FAADaG6GrRrL5klJVTi9K7L8IAEC7I3TVgLsrWyiqY4HT\ni4QuAADaF6GrBgolV8lVdcsIielFAADa3byhy8xuMrNDZvbkHOffbWaPm9kTZnafmV1Uce7Z4Phj\nZvZwlANvZtlCeePqaltGSFJnKq5MvqRiyWs1LAAA0EDVpIKbJW0/xflfSbrS3S+U9ElJN55w/mp3\nv9jdty1uiK0nm69+s+tQ2CB1Os8UIwAA7Sgx3xPc/V4z23KK8/dVPHxA0salD6u1ZWYqXQubXpRo\nkAoAQLuKek3X+yXdXfHYJX3XzB4xsx0Rf1bTyhaCSteCWkYElS4W0wMA0JbmrXRVy8yuVjl0XVFx\n+Ap332tmqyV9x8yedvd753j9Dkk7JGnz5s1RDashsvmFr+li/0UAANpbJJUuM3uppC9Jutbdj4bH\n3X1v8POQpDskXTrXe7j7je6+zd23DQ0NRTGshgkX0le7DZDE9CIAAO1uyaHLzDZL+qqk97r7zyuO\nd5tZb3hf0hskzXoFZLsJpxdTi5hepNIFAEB7mnd60cxuk3SVpFVmtkfSJyQlJcndb5D0cUkrJX3B\nzCSpEFypuEbSHcGxhKRb3f1bNfgdmk44vbiQSlc6EVPMCF0AALSraq5evH6e8x+Q9IFZju+WdNHJ\nr2h/2UJRJikZt6pfY2Y0SAUAoI3Rkb4GMoWS0smYgipf1dgKCACA9kXoqoFcvrSgHl0hQhcAAO2L\n0FUDmUJxQe0iQkwvAgDQvghdNZArlBYZuqh0AQDQrghdNZDJF5VOLn560Z1NrwEAaDeErhrILrrS\nlVCx5FS7AABoQ4SuGsgWSgvq0RUKG6Qem8pFPSQAANBghK4ayBaKC+pGHwpD18hUPuohAQCABiN0\nRczdlc2X1LGI6cXOYP9FKl0AALQfQlfE8kWXS4vu0yVJx6h0AQDQdghdEQs3u04vaXqRShcAAO2G\n0BWxcLPrxVy92BlWuiapdAEA0G4IXRHLFsLQtfDpxUQspnQixpouAADaEKErYpklTC9K5SnG0Wkq\nXQAAtBtCV8SOTy8uvNIllRukUukCAKD9ELoiFi6kX0zLCKlc6eLqRQAA2g+hK2Lhmq7UIkNXZyrO\n1YsAALQhQlfEwtDVsYgNr6VgenGS0AUAQLshdEUsmy8qZlIiZot6fVcqrrFMQYViKeKRAQCARiJ0\nRSxbKCmdiMts8aFLElcwAgDQZghdEcsWiotuFyGVpxcltgICAKDdELoilsmXFtWNPsRWQAAAtCdC\nV8RywfTiYrHpNQAA7YnQFbFMobjESlc4vUilCwCAdkLoili2UFJ6ke0iJKYXAQBoV4SuiGXzS6t0\npRMxJWLG9CIAAG2G0BWxbKG06C2AJMnMNNCVotIFAECbIXRFqOReXki/hOlFSVrRldSxSSpdAAC0\nE0JXhPKFklxa0vSiJK3oSrGQHgCANkPoilC47+JSWkZI0kBXUiOs6QIAoK0QuiKUKRQlUekCAAAn\nI3RFKBdWupawDZAkDXSXK13uHsWwAABAEyB0RSiTj2h6sTOlXLGkqVwximEBAIAmQOiKUC6y6cWk\nJLrSAwDQTghdEcoE04sdS2wZMdCVkiQW0wMA0EYIXRHK5suVrlRElS5CFwAA7WPedGBmN5nZITN7\nco7zZmafMbNdZva4mV1ScW67mT0TnPtolANvRmHLiKV0pJekFd3lShfTiwAAtI9q0sHNkraf4vw1\nks4KbjskfVGSzCwu6fPB+fMlXW9m5y9lsM0uWygpbqZEfIlXL85UughdAAC0i3nTgbvfK2n4FE+5\nVtKXvewBSQNmtk7SpZJ2uftud89Juj14btvKFopLbhchla9elMSm1wAAtJEo1nRtkPRCxeM9wbG5\njretbL605CsXpfKasJ50gulFAADaSNMspDezHWb2sJk9fPjw4UYPZ1GyhdKSe3SF2AoIAID2EkXo\n2itpU8XjjcGxuY7Pyt1vdPdt7r5taGgogmHVXyai6UVJGuxOaXiSShcAAO0iioRwp6T3BVcxvlLS\nqLvvl/SQpLPM7HQzS0m6Lnhu24pqelEidAEA0G4S8z3BzG6TdJWkVWa2R9InJCUlyd1vkHSXpDdK\n2iVpStLvBucKZvZhSfdIiku6yd131uB3aBrZQkmDQbuHpRrsSukXBycieS8AANB484Yud79+nvMu\n6UNznLtL5VC2LGQLxUgrXUcns5G8FwAAaLymWUjfDrKF0pK3AAoN9qSUyZc0zabXAAC0BUJXREru\nyhVKS94CKLQymKak2gUAQHsgdEUkF9EWQKHB7rQksZgeAIA2QeiKSLjvYlR9ugZnKl2ELgAA2gGh\nKyKZfHntVVR9usLpxeEJQhcAAO2A0BWRXMSVrhVh6KLSBQBAWyB0RSRTCCpdEa3p6utIKBk3DbP/\nIgAAbYHQFZFsPqh0RTS9aGZa0ZViehEAgDZB6IpI1AvppbBBKqELAIB2QOiKSDaYXoyqZYQkrexJ\naZg+XQAAtAVCV0TCSlcqoulFqdyri4X0AAC0B0JXRLL5ohIxUyIWYejqSjK9CABAmyB0RSRbKEV2\n5WJosDut8Uxhph0FAABoXYSuiGQLJaUj2uw6NNhT7tU1QtsIAABaHqErItl8MfJK10q2AgIAoG0Q\nuiKSqcn0Il3pAQBoF4SuiOQKpUh7dElUugAAaCeErohk8sXIutGHZipdE/TqAgCg1RG6IpKtQaVr\noCslM6YXAQBoB4SuiGQLxUi70UtSPGYa6Eyy6TUAAG2A0BWBQrGkfNEj7UYfGuxOUekCAKANELoi\nMJkL912MdnpRklZ2p3V0gtAFAECrI3RFYCJbkKTIW0ZIVLoAAGgXhK4ITIahK+KO9FK5Kz2hCwCA\n1kfoisB4pnaVrpXdKR2byqlU8sjfGwAA1A+hKwK1nF5c0ZVSyaWR6Xzk7w0AAOqH0BWBiZlKVw0W\n0vewFRAAAO2A0BWBiWy5CtVRo5YREqELAIBWR+iKQLimq6MWC+lnQhdbAQEA0MoIXREYyxRkklI1\nWUiflsSm1wAAtDpCVwTGM3mlEjHFzCJ/7xXdSUnSMA1SAQBoaYSuCExkCjWZWpTKi/N70wkqXQAA\ntDhCVwTGM4WatIsIraArPQAALY/QFYHxbL5mlS6pvJj+2BShCwCAVkboisB4plCTdhGhld0pNr0G\nAKDFEboiMJEp1KQxaohNrwEAaH1VhS4z225mz5jZLjP76Czn/9TMHgtuT5pZ0cwGg3PPmtkTwbmH\no/4FmsFYDRfSS8c3vXZn/0UAAFpVYr4nmFlc0uclvV7SHkkPmdmd7v5U+Bx3/5SkTwXPf4ukP3b3\n4Yq3udrdj0Q68iYynsnXfHoxVyxpIltQb0eyZp8DAABqp5qkcKmkXe6+291zkm6XdO0pnn+9pNui\nGFwryBVKyhZKNV5IX26QyhQjAACtq5rQtUHSCxWP9wTHTmJmXZK2S/pKxWGX9F0ze8TMdix2oM1q\nIhtudl27Stdg2CCV0AUAQMuad3pxgd4i6UcnTC1e4e57zWy1pO+Y2dPufu+JLwwC2Q5J2rx5c8TD\nqp3xTLjZNZUuAAAwt2rKM3slbap4vDE4NpvrdMLUorvvDX4eknSHytOVJ3H3G919m7tvGxoaqmJY\nzWFms+saXr24Mtj0mq70AAC0rmpC10OSzjKz080spXKwuvPEJ5lZv6QrJX294li3mfWG9yW9QdKT\nUQy8WYwFla50DRfSDwahi0oXAACta97pRXcvmNmHJd0jKS7pJnffaWYfDM7fEDz1bZK+7e6TFS9f\nI+kOK28EnZB0q7t/K8pfoNEmwkpXDacXu1JxpRMxQhcAAC2sqjVd7n6XpLtOOHbDCY9vlnTzCcd2\nS7poSSNscsenF2tX6TIzutIDANDi6Ei/RPVYSC+FDVKzNf0MAABQO4SuJQorXbVc0yVJK7pSGp7K\n1/QzAABA7RC6lmgiW1A6EVMiVtuvcmU3lS4AAFoZoWuJxjL12ZpnsDutYdZ0AQDQsghdSzSeyauv\nI+oesydb2ZPSZK6oTL5Y888CAADRI3Qt0XimoJ46hC56dQEA0Npqnxba3ES2oN6IQ9etDz5/0rGn\n9o1Kkm558HltGOjUuy5rna2SAAAAla4lG8/k1Zuu/ZqurlQ52E0FG2wDAIDWQuhaovFM9JWu2XSn\ny58xmSN0AQDQighdS1SvNV3d6XLz1cksC+kBAGhFhK4lKJU8WNNV++nFjmRcMZMmmV4EAKAlEbqW\nYCKY6qtHy4iYmbpSCaYXAQBoUYSuJQi3AKrHmi6pPMXI9CIAAK2J0LUE4WbXPXW4elGSulMJphcB\nAGhRhK4lmKhzpau3I6FxQhcAAC2J0LUE9Z5eHOhKaXQqr5J7XT4PAABEh9C1BGPB9GI9rl6UpP7O\npIruMxU2AADQOghdS1D/Slc53I1M5+vyeQAAIDqEriWYyNZ/elGSRqbY9BoAgFZD6FqC8Uxe8Zip\nMxmvy+cNdJYrXaNUugAAaDmEriUI9100s7p8Xkcyro5kTMemCF0AALQaQtcSjGcK6knXZ2oxNNCZ\n0ijTiwAAtBxC1xKUK131uXIx1N+ZZCE9AAAtiNC1BOOZfN0W0YcGupIaYXoRAICWQ+hagvFMQb31\nnl7sSmk6X2Q7IAAAWgyhawnGsw2odAVXMO4bma7r5wIAgKUhdC3BRAPWdIUNUvcSugAAaCmErkVy\n95mWEfUUNkjdN5Kp6+cCAIClIXQtUiZfUqHk6qlz6OrtSChmTC8CANBqCF2LNF7nza5DMTP1dSaZ\nXgQAoMUQuhZpPLh6sK/OlS6pvJie0AUAQGshdC3SeKa+m11XGuhKMb0IAECLIXQtUji92JOu7/Si\nVK50HRjNqFjyun82AABYHELXIjWy0tXflVSh5Do0zhWMAAC0CkLXIk00cnqxM2wbwRQjAACtgtC1\nSGMNunpRqmyQSqULAIBWUVXoMrPtZvaMme0ys4/Ocv4qMxs1s8eC28erfW2rCqcXe+q896LEVkAA\nALSieRODmcUlfV7S6yXtkfSQmd3p7k+d8NQfuPubF/naljOeKag7FVc8ZnX/7HQyrv7OpPYeI3QB\nANAqqql0XSppl7vvdvecpNslXVvl+y/ltU1tIptvyNRiaP1AJ5UuAABaSDWha4OkFyoe7wmOnehV\nZva4md1tZhcs8LUtpxH7LlbaMNBJg1QAAFpIVAvpH5W02d1fKumzkr620Dcwsx1m9rCZPXz48OGI\nhlU745lC3fddrLRhoINKFwAALaSa0LVX0qaKxxuDYzPcfczdJ4L7d0lKmtmqal5b8R43uvs2d982\nNDS0gF+hMcYzjZ9eHMsUZpq0AgCA5lZN6HpI0llmdrqZpSRdJ+nOyieY2Vozs+D+pcH7Hq3mta1q\nPNvY6cX1A52SpH20jQAAoCXMmxrcvWBmH5Z0j6S4pJvcfaeZfTA4f4Okd0j6fTMrSJqWdJ27u6RZ\nX1uj36WuxjOFhmx2HdqwIgxd0zpnbW/DxgEAAKpTVWoIpgzvOuHYDRX3Pyfpc9W+th2MZ/IN6dEV\n2hBUulhMDwBAa6Aj/SLkiyVl8qWGruka6kkrGTcW0wMA0CIIXYvQyH0XQ7GYaW1/B5UuAABaBKFr\nEcZnQlfjKl2StL6fBqkAALQKQtcihJtdN3JNl1ReTM/ViwAAtAZC1yKEla5GXr0olRfTHxjLqFAs\nNXQcAABgfoSuRZjINsn04kCniiXXwfFsQ8cBAADmR+hahLALfCMX0kuVDVJZ1wUAQLMjdC1COL3Y\nyL0XpeO9ughdAAA0P0LXIjRPpatDEg1SAQBoBYSuRRjPFpRKxJROxBs6jq5UQiu6klS6AABoAYSu\nRRjPFNTb4HYRofUDndp7jNAFAECzI3Qtwnim0PCpxdCGAXp1AQDQCghdizCeyTe8XURo44ouPT88\npWLJGz0UAABwCoSuRZhookrXeet6NZ0v6tmjk40eCgAAOAVC1yKMZwoN3wIotHVDvyRp576xBo8E\nAACcCqFrEZppevHM1T1KJWLauXe00UMBAACnQOhahGZaSJ+Mx3Tu2l49uY/QBQBAMyN0LVCp5JrI\nFRq+2XWlC9b3aee+MbmzmB4AgGZF6FqgyVxB7o3fAqjSBev7NTKVpzM9AABNjNC1QOG+i82ypksq\nV7okFtMDANDMCF0LdDx0NU+l67x1fYrHjMX0AAA0MULXAk1kw82um6fS1ZGM6yVD3VS6AABoYoSu\nBRoLKl3N0qcrtHV9P1cwAgDQxAhdCxROLzbT1YuSdP76Ph0cy+rweLbRQwEAALMgdC3QeKb5phel\nys70VLsAAGhGhK4FGpkqh67+zuYKXedzBSMAAE2N0LVA+0amtaIrqc5UvNFDeZG+jqROW9lFpQsA\ngCZF6FqgfSPTWtff2ehhzOqC9X16ci+VLgAAmhGha4H2j2a0fqBZQ1e/nh+e0uh0vtFDAQAAJyB0\nLdDekWltGOho9DBmFXamf4p1XQAANB1C1wKMZ/IazxS0rokrXRJXMAIA0IwIXQuwfzQjSU07vTjU\nm9aavjRXMAIA0IQIXQuwb2RakrS+vzmnF6VyZ3oqXQAANB9C1wLsG2nuSpdUXte169CEpnPFRg8F\nAABUaK69bJrc/tFpxUxa3Ztu9FB064PPz3p8eDKnkkuf+edfaNNgl9512eY6jwwAAMyGStcC7B2Z\n1tq+DiXizfu1hYv89wZToQAAoDlUlR7MbLuZPWNmu8zso7Ocf7eZPW5mT5jZfWZ2UcW5Z4Pjj5nZ\nw1EOvt72jUw39dSiJA10JtWZjM+sPwMAAM1h3ulFM4tL+ryk10vaI+khM7vT3Z+qeNqvJF3p7sfM\n7BpJN0q6rOL81e5+JMJxN8T+0YxeunGg0cM4JTPT+oGOmSstAQBAc6im0nWppF3uvtvdc5Jul3Rt\n5RPc/T53PxY8fEDSxmiH2Xilkmv/SEbrm7QxaqX1A506MJZRoVhq9FAAAECgmtC1QdILFY/3BMfm\n8n5Jd1c8dknfNbNHzGzHwofYHI5O5pQrlrS+SfddrHTGqm4VS65dhycaPRQAABCI9OpFM7ta5dB1\nRcXhK9x9r5mtlvQdM3va3e+d5bU7JO2QpM2bm++Ku5keXU2+pkuSXrK6Rx3JmJ7cS78uAACaRTWV\nrr2SNlU83hgcexEze6mkL0m61t2PhsfdfW/w85CkO1SerjyJu9/o7tvcfdvQ0FD1v0Gd7B8th651\nTdwYNZSIxXTe2j49tX9MuQJTjAAANINqQtdDks4ys9PNLCXpOkl3Vj7BzDZL+qqk97r7zyuOd5tZ\nb3hf0hskPRnV4Otpb9AYdUMLVLok6cIN/crkS/rRrpa/fgEAgLYw7/SiuxfM7MOS7pEUl3STu+80\nsw8G52+Q9HFJKyV9wcwkqeDu2yStkXRHcCwh6VZ3/1ZNfpMa2z8yrc5kXANdyUYPpSpnru5ROhHT\nN5/Yr6vPXd3o4QAAsOxVtabL3e+SdNcJx26ouP8BSR+Y5XW7JV104vFWtG90WusGOhQEyKaXiMd0\n/ro+fXvnAeXedqFSieZt6AoAaB3urrFMQRPZgkolV8ldJZdK7upKxbVujgvO5tpJpVK776LCNkBV\n2juSaZmpxdDWDf36yQsjuu+XR3TVOVS7AACzK5Vcx6ZyOjSe1eHwNnH8/pGJrIYncxqezOnYVE75\nos/5Xuv6O7R1Q78uXN+vVU2wbV4zIXRVaf/ItM45p/kW+J/KWat71JNO6K4n9hO6AGAZyuSLLwpQ\nx0NVZubxobFyqCqUTg5SybiptyOpnnRC3am4Ng926bx1fepOxdWRjCtmJjOVbzKNZ/J6ct+YvvPU\nQX3nqYNa19+hl24c0OVnrlQixowLoasKuUJJhyeyLdEuolIiHtOvn7da337qoP5zsaRkE+8ZCQCo\nTliVOjxRDkyVYepQEKbCkDWeKZz0ejNpZXc6CFQJrR/o1Dlre9XbkZgJWH0dCfV0JJROxBc8vivO\nGtLodF5P7h3VE3tHdc/OAzoyntVvXbKhZZbo1AqhqwoHxzJyV0s0Rj3RGy9cp689tk/3/fKorjy7\ntSp1ALBLMSVbAAAXF0lEQVRcuLvGswUdncjpyERWR4LgdOSEab5nj05pPJPXLEUppeIx9XYk1JNO\nqLcjoa3r+1/0uLcjqZ6OhLpTCcVjtQ0//Z1JXX7mKl1+5ip992cH9S9PH9KqnpSuXOazLoSuKuxt\nocaoJ3rN2UPqSSd09xP7CV0AUEfurpGpvI5MHF8jdSQIVUdfdL9ctZqtr6JJ6g5CU086oZcM9QQB\nKgxTySVVperhdeeu1pGJrO556qBW9qS1dUN/o4fUMISuKoSNUVth38UTdSTjet15q3XPzgP65G9u\nZYoRAJagGEzthWEpDFRhgJq5jed0dDI764LzmKm8RipdDk6re9M6Y1W3eoIgFR7v7Uioqw5VqVoz\nM739ko0amcrrHx95QQNdSW1c0dXoYTUEoasK+4LGqHNdBtvs3njhOn39sX16YPdRvfosql0AEHJ3\nTeWKGp7MzVyhd3Qip6OTuaAaldXRydxMsBqezM46tRc3mwlNPemENgTrpHrSiRcd700n1JEqL0Bf\nTpLxmN7zytP0xe/t0t/f/5x+/6qXaKAr1ehh1R2hqwr7RqY12J1SZ6o5S7fzufLsIXWn4rrrif2E\nLgBtbzp3/Iq9o0FoOv6z3Pbg6GROxyZzGp7KzbldWjJuM2GpO53QlpVd2rqhb+ZYGKZ600l1JGPL\nfpH4fHrSCb3v17bohu//Ul++/zn93pVnNO2UaK0Quqqwb2S6JfZcnEt5inGN7tl5UJ+8tqQEU4wA\nWky+WCqvfRrP6vBE5ngvqRMWmu8bzcwZotKJmLqD1gfdQTXqrDU96k4l1J2OBz+PT+/RVDp6a/o6\n9K5LN+t/3fesfviLI3rdeWsaPaS6InRVYf9oRpsGW3v++S0XrdedP92nrz22T+94+cZGDwcAVCq5\nRqfzLwpNJ4ao8PHwZG7W9+hIxtSbLl+V19uR1LaBTvWmE+oJWh+E1aiuVJw1rU3irDW9Om9tr+4P\nlrwsp3BL6KrC3pFpXXb6YKOHsSSvO3e1Lto0oE/d87Su2bpW3Wn+0QOInns5SJUXmL94cXm4LioM\nU3MtNE/EbObqvJ6OpM5c3RMEqURFwCqfJ0i1pivPHtIN9+7WQ88O6/IzVzV6OHXD/3nnMZ7JazxT\naMl2EZViMdPH33y+3v7F+/Q33/+l/t0bzmn0kAC0iHBq73iAqlhkPpHTkYpF58OTs28RE16xF1af\nKhtyhq0Pwiv20gnWR7W7zSu7tWVlt36464guO2Nw2XSrJ3TNY/9o+crFVg9dkvTy01borRet19/c\nu1vvvHRzy+0lCSA6YUUq3Abm8ERmprv54ZmqVPnnsan8rO+RiB2/Yq87ldDGgS6ds6bvRVfrVU7v\nLbcr9nBqV549pL+7/1k9/sKoLjltRaOHUxeErnkcb4zaugvpK/2Ha87VPTsP6L/d/bQ+c/3LGj0c\nABErlVzDUzkdGitvCRNuD3NwLDNz7OBYds5mnJ3JuDpT8Zmq09lreiuu0ju+yDxcaE5FCot19poe\nrevv0Pd/cVgXbx5YFqGc0DWP/SPtU+mSpA0DndrxmjP02X/ZpX/zqi16+TL52wXQykol11gm/+IG\nnOPH+0eFe+89d3RSE9nCrH2kOpIx9XYk1duR0FBvWmcMdc887u1IqC9YK8XUHurFzPSas4f0Dw+9\noKf3j+n89e3fqZ7QNY99I9OKx0yre1uz0nXrg8+fdGywO6XejoT+6Paf6PeufIne88rTGjAyYHnL\n5I/3kgoDVLg1THi/3LAzp2NTORVnSVIxK/95Hurt0Ore9PG99zqObw8TBisWnKMZbV3fr+90H9T3\nfn5Y563ra/Rwao7QNY99o9Na29fR8tswVEon4vqN89fqnx7do8f3jEgidAFRyRdLOjSe1YHRaR0Y\nLU/rHb+Vp/cOj2c1linM+vp0IjbTjHOmIef6vuBxXD3p460QWCeFVhePmV591ip9/bF9+tWRyUYP\np+YIXfNo9caoc7l484Du331U9+w8qP/41oK6UvyrAJyKu2ssU9ChsYwOjGV0YLQcpMr3szP3j0xk\n5ScUpeIxm6k69XUkdP76/vK0XrBuKlwn1U0LBCxDl2xeoX/+2SF9/+eH9ReNHkyN8X/aeewfzeii\njQONHkbkYmZ604Xr9Lc/2K0P3fKo/ua925ZVgzogVCy5jk5mj1+5F0z5HRorL0Iv38qL0LNzLDzv\n60yovzOp0wa7dOGGfvV3JtXXkVRfZ0J9HUl1peKskwLmkIzHdPmZq3TPzgN6cu+otm5o37VdhK5T\nKJVc+0cyumZreyyiP9GWVd269uIN+tpje/WR236iz73rZWwRhLbh7hrPFnRgNKP9oxkdHC1XovaP\nZmYC1cGgMjXnwvN0eT3UYFdKWwa7g0XnSfV1litWfZ1JKlNABC47fVDfe+aQbvj+L/W5d13S6OHU\nDKHrFI5MZpUrltqmXcRsLj19UBdu7Ncnv/GU/uQff6pP/87FbbV+De1rOlfUvtFp7RuZ1v6RjPaO\nTGv/6LT2ByFr/8i0JnPFk17XnYqrr7McpjYNdumC9X0VV/ElZzqfE6aA+ulIxnXZ6YO664n9emF4\nquW33psLoesUZtpF9LdnpSv0/itOVyZf1KfueUadybj+y9suVIzghQaazhW1f3R6Zu3U/tGM9o1M\n68BoRvtGM9o/Oq2RWRp29qYT6u9Kqr8zqYs2Dai/M1kx1VeuTlHNBZrTq16ySvfvPqov/WC3/vLa\nrY0eTk0Quk5hX9AYdV0bV7pCH7r6TGXyRX32X3apIxnXJ95yPmtQUBPjmfxMeDoQVKYOBtN+YcAa\nnT45UHUm4xroKgeos9f0aiAIVP1dSQ10ptTXmVg2W4kA7aivM6lrL96g//PwHv3Rr5+tFd2pRg8p\ncoSuU9gXbAG0XLbL+XevP1vTuaK+9MNfad/ItP7q2q1a24ZXbqJ2coWSDo6Vp/r2Bbe9I+Uq1f7R\n8jTgePbkVgk96URQkUro3LW9xytUQZWqvzPJhR7AMrDjNWfonx7Zo79/4Dl95HVnNXo4kSN0ncIv\nD0+oKxVXf2ey0UOpCzPTn7/pPK3sSeuvv/tz/fqnv68/236O3n3ZaazzgqRylWrvyLT2Hpsu/xyZ\n1r6RjPYem9K+kYwOjmdOapfQnYoHU34pXbChv1yh6kqqv6P8s7eDChWAsrPX9Oq1567W3933rHa8\n5gx1JOONHlKkCF1zyOSL+sZP9+l1561p+2m2E7vW93cm9eGrz9TXH9unj399p/723t268X3blkW3\n4OWsVHIdnsgG1amgSjUTrsrB6sSGnvGYzYSoDQOdOn99nwY6kxroSs0cZ0E6gIXY8ZozdN2ND+if\nHtnTdjumELrm8M3H92ssU9D1l25q9FAaYmVPWr97+RY99sKIvvnEfr3lsz/U2162Qe98xSa9/LQV\nbR9E2022UJzZ9PjAaFYHxjIzTT7DK/8OjmVUOKF3QjoR04qulAa6kkGgKt8fCI71pBN0RAcQqctO\nH9RFmwb0tz/Yresv3dxWMy2Erjnc+uPndcaqbv3aGSsbPZSGMTO9bPMKnbOmV88OT+qrj+7VPz6y\nRy8Z6tbvbNuk37pko4Z6040e5rJVLLmOTZU3PD4ykX3Rz3AD5PDnbAvTEzGbWTM11JvWmat71N+Z\nnKlQrehKtV1pH0DzMzP93mvO0B/c8qi+vfOArrlwXaOHFBlC1yyeOTCuR547pj9/43lUdCR1pRP6\nT795oT52zXn65uP79Q8Pv6D/9+6n9al7ntFlZwzqFVsGdenpg3rZphXqTPE/6aXIFUoanszp6GRW\nRydyOjKRDW45HQk3Rw6OH52jqWcybupJH9/o+Lx1veoJmnyGXdL7O5LqpEs6gCb1Gxes1Wkru3TD\nvbu1fevatvlvFaFrFrf9+Hml4jG9/eUbGz2UplG57uvtl2zUq89cpUefP6ZfHJrQfbuOyiXFzXTR\npn5dtGlAZwz16CWrunXGUI/W9KXb5g9MtQrFkiazRY1O50+6HZvKaWQqp2NTeR2bzOnYVC4IWjmN\nz7EJciJmM3vz9XYkdFrQ1LMn2LMvDFg96YTSidiy+74BtJd4zPSBV5+h/+drT+rHvxrWZW0y60To\nOsF0rqivPLpH27eu1WAb9giJyuq+Dm3fuk7bVf7Onh+e1K+OTGkyV9DtP35B0/njncC7U3FtXtmt\nNX1prent0Jq+tIb6OjTUkw7aAgQVmI6kejoSdZm/zxdLyhZKyuSLyhZKyuaLyuRLyhSK5WP58rlM\noajpXEnT+fLxTL6o6VxRk7mipnIFTQU/J7NFjWfymsgWNJ4pHz+VZNzUlUqoKxVXV9AhfV1/p7rT\nCXWn4+pOlQNUD0EKwDL12y/fqP/x3Z/rE3fu1B1/cHlbzKQQuk7wzSf2azxT0Lsu29zoobSMzlRc\n56zt0zlry1c3/tbLXGOZwszGwUfGszo2ldMvDk7okeeOaTJbmHVaLJSKx5ROxpROxJVOlO8nYqaY\nmeKx8q1y8bZLkrtcUsldhaKrWCrfCiVXvliaCVm5Qvn+qT7/VEzlzVlTifItnYgpFTzuTCW0oiul\ndCKmjmRc6WRcneEtdfxnVyrOFX0AMI+OZFz/329fpN+9+SH9+dee0H//7Yta/i+fhK4T3Prgczpj\nqFuXnT7Y6KG0LDObaW555uqek86X3DWRLWgiUzhePQorS/niTFAqFIOfJVfJXSUvLx7PF0sn9YIK\n/xzGrBzIEglTLAxqZorHy2EtETv+MxmPKRGPKRkzJeKxmWPJeHgufFwOVol4+XWt/oceAFrFVees\n1h++7iz99Xd/oZeftkLvvqy1W0gQuio8fWBMjz4/or94EwvoaylmNjOdCADAqXzktWfpJ8+P6C/v\nfEpb15fXDbcq5jgq3Pbg80olYnr7JSygBwCgGcRipr9+58Ua6k3rD255VMOTuUYPadEIXYHpXFFf\n/clevXHr2rbcZBMAgFa1ojulL77nEh0ez+oPb/+JiotdmNtgVYUuM9tuZs+Y2S4z++gs583MPhOc\nf9zMLqn2tc3ijp/sDRbQt/Z8MQAA7eilGwf0H996gX7wiyP603/8qV4Ynmr0kBZs3jVdZhaX9HlJ\nr5e0R9JDZnanuz9V8bRrJJ0V3C6T9EVJl1X52oYpFEv67s8O6sv3P6f7fnlU563r0yu2rGj0sAAA\nwCyuv3STnhue1P/8wa/0tcf26jcuWKv3X3F6y2xPV81C+ksl7XL33ZJkZrdLulZSZXC6VtKX3d0l\nPWBmA2a2TtKWKl5bd0cmsrr9x8/rlgef1/7RjDYMdOrPtp+j61+xuSX+oQEAsByZmT52zXn6v161\nRV++/znd+uDzuvvJA7poY7/eevEGrelLa2V3WkO9Ka3sLveCjDXR3o3VhK4Nkl6oeLxH5WrWfM/Z\nUOVr68rd9VtfuE/PD0/pijNX6S/feoFed96attpQEwCAdrauv1P/Yfu5+revPVNfeWSP/tePntUn\nv3FyPacjGdPP/mp70xRUmqZlhJntkLQjeDhhZs/U+jOfk3RLdG+3StKR6N4OVeJ7rz++8/rjO28M\nvvc6e3cN3jP2n2rwpierakF4NaFrr6RNFY83BseqeU6yitdKktz9Rkk3VjGepmRmD7v7tkaPY7nh\ne68/vvP64ztvDL53RK2aqxcfknSWmZ1uZilJ10m684Tn3CnpfcFVjK+UNOru+6t8LQAAQNubt9Ll\n7gUz+7CkeyTFJd3k7jvN7IPB+Rsk3SXpjZJ2SZqS9Lunem1NfhMAAIAmVtWaLne/S+VgVXnshor7\nLulD1b62TbXs1GiL43uvP77z+uM7bwy+d0TK/MSdgwEAABA5tgECAACoA0JXBFplq6N2YmY3mdkh\nM3uy0WNZDsxsk5n9q5k9ZWY7zewPGz2m5cDMOszsx2b20+B7/8tGj2m5MLO4mf3EzL7R6LGgfRC6\nlqhiq6NrJJ0v6XozO7+xo1oWbpa0vdGDWEYKkv7E3c+X9EpJH+Lf87rISnqtu18k6WJJ24MrxFF7\nfyjpZ40eBNoLoWvpZrZJcvecpHCrI9SQu98rabjR41gu3H2/uz8a3B9X+X9GGxo7qvbnZRPBw2Rw\nYyFujZnZRklvkvSlRo8F7YXQtXRzbYEEtCUz2yLpZZIebOxIlodgmusxSYckfcfd+d5r768l/Zmk\nUqMHgvZC6AJQNTPrkfQVSX/k7mONHs9y4O5Fd79Y5R09LjWzrY0eUzszszdLOuTujzR6LGg/hK6l\nq2abJKDlmVlS5cB1i7t/tdHjWW7cfUTSv4q1jLV2uaS3mtmzKi8Xea2Z/e/GDgntgtC1dGx1hLZn\nZibpf0r6mbt/utHjWS7MbMjMBoL7nZJeL+npxo6qvbn7x9x9o7tvUfm/5//i7u9p8LDQJghdS+Tu\nBUnhVkc/k/R/2Oqo9szsNkn3SzrHzPaY2fsbPaY2d7mk96r8t/7HgtsbGz2oZWCdpH81s8dV/gve\nd9ydFgZAi6IjPQAAQB1Q6QIAAKgDQhcAAEAdELoAAADqgNAFAABQB4QuAACAOiB0AQAA1AGhC8CS\nmFkx6Nu108x+amZ/Ymax4Nw2M/vMKV67xczeVb/RnvTZ08G+hk3BzN5pZrvMjF5cQBsidAFYqml3\nv9jdL1C5Y/o1kj4hSe7+sLt/5BSv3SKpIaEr8MtgX8OqmVm8VoNx93+Q9IFavT+AxiJ0AYiMux+S\ntEPSh63sqrBqY2ZXVnSz/4mZ9Ur6r5JeHRz746D69AMzezS4vSp47VVm9j0z+ycze9rMbgm2JpKZ\nvcLM7guqbD82s14zi5vZp8zsITN73Mx+r5rxm9nXzOyRoGq3o+L4hJn9dzP7qaRfm+MzLwjuPxZ8\n5lnBa99TcfxvwtBmZtuD3/GnZvbPEf5jANCkEo0eAID24u67g2Cx+oRT/17Sh9z9R2bWIykj6aOS\n/r27v1mSzKxL0uvdPROEltskbQte/zJJF0jaJ+lHki43sx9L+gdJ73T3h8ysT9K0pPdLGnX3V5hZ\nWtKPzOzb7v6reYb/f7v7cLDP4UNm9hV3PyqpW9KD7v4nwR6rT8/ymR+U9D/c/ZbgOXEzO0/SOyVd\n7u55M/uCpHeb2d2S/lbSa9z9V2Y2uOAvGkDLIXQBqJcfSfq0md0i6avuvicoVlVKSvqcmV0sqSjp\n7IpzP3b3PZIUrMPaImlU0n53f0iS3H0sOP8GSS81s3cEr+2XdJak+ULXR8zsbcH9TcFrjgZj+Upw\n/Jw5PvN+SX9uZhuD3+8XZvY6SS9XOcBJUqekQ5JeKeneMAS6+/A84wLQBghdACJlZmeoHFIOSTov\nPO7u/9XMvinpjSpXnn5jlpf/saSDki5SeflDpuJctuJ+Uaf+75dJ+rfufs8Cxn2VpF+X9GvuPmVm\n35PUEZzOuHvxVK9391vN7EFJb5J0VzClaZL+zt0/dsJnvaXacQFoH6zpAhAZMxuSdIOkz7m7n3Du\nJe7+hLv/N0kPSTpX0rik3oqn9atcRSpJeq+k+RatPyNpnZm9IviMXjNLSLpH0u+bWTI4fraZdc/z\nXv2SjgWB61yVq1FVf2YQNne7+2ckfV3SSyX9s6R3mNnq4LmDZnaapAckvcbMTg+PzzM2AG2ASheA\npeoMpvuSkgqS/l7Sp2d53h+Z2dWSSpJ2Sro7uF8MFqjfLOkLkr5iZu+T9C1Jk6f6YHfPmdk7JX02\nWIc1rXK16ksqTz8+Giy4PyzpN+f5Pb4l6YNm9jOVg9UDC/zM35H0XjPLSzog6b8E68P+QtK3rdxG\nI6/yurYHgoX6Xw2OH1L5yk8AbcxO+MsoACwLZrZF0jfcfWuDh/IiwTTnzMUFANoH04sAlquipH5r\nsuaoKlf7jjV6LACiR6ULAACgDqh0AQAA1AGhCwAAoA4IXQAAAHVA6AIAAKgDQhcAAEAd/P+Tu+OF\nKDN9dwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f34e3812b00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nb_merge_dist_plot(\n",
    "    SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n",
    "    SkyCoord(viking['viking_ra'], viking['viking_dec'])\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Given the graph above, we use 1 arc-second radius\n",
    "master_catalogue = merge_catalogues(master_catalogue, viking, \"viking_ra\", \"viking_dec\", radius=1.*u.arcsec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add DES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAF3CAYAAABnvQURAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0ZGd55/HfU7tUWrol9WL3buPdMTZuryzxEAw2EEyC\nGTCLT4DEgcCAEzJJmOTATJKTkMkMk4AJHiAOycTskOAQOw4hLLYB023HC3bb0N7b7nZ3S2otJdX+\nzh/3XnW1WktJulW3qvT9nKOjUlWp6u3CuH9+3uc+rznnBAAAgJWLRb0AAACATkGwAgAACAnBCgAA\nICQEKwAAgJAQrAAAAEJCsAIAAAgJwQoAACAkBCsAAICQEKwAAABCQrACAAAISSKqNx4aGnLbt2+P\n6u0BAADqds899xx2zq1b7HmRBavt27dr9+7dUb09AABA3czsqXqex1YgAABASAhWAAAAISFYAQAA\nhIRgBQAAEBKCFQAAQEgIVgAAACEhWAEAAISEYAUAABASghUAAEBICFYAAAAhIVgBAACEhGAFAAAQ\nEoIVAABASBJRL6BTff7upxd9zlsu2tqElQAAgGahYgUAABASghUAAEBICFYAAAAhIVgBAACEhGAF\nAAAQEoIVAABASAhWAAAAISFYAQAAhIRgBQAAEBKCFQAAQEgIVgAAACEhWAEAAISEYAUAABASghUA\nAEBICFYAAAAhIVgBAACEhGAFAAAQEoIVAABASAhWAAAAISFYAQAAhIRgBQAAEBKCFQAAQEgIVgAA\nACEhWAEAAISEYAUAABASghUAAEBICFYAAAAhIVgBAACEhGAFAAAQEoIVAABASAhWAAAAISFYAQAA\nhGTRYGVmW8zsO2b2sJk9ZGYfmOM5ZmYfN7O9ZvaAmb2oMcsFAABoXYk6nlOW9EHn3L1m1ivpHjP7\nlnPu4ZrnXCnpFP/rIkmf8r8DAACsGotWrJxz+51z9/q3JyTtkbRp1tOukvR3zvMjSWvM7ITQVwsA\nANDCltRjZWbbJZ0n6e5ZD22S9EzNz/t0fPgCAADoaHUHKzPrkfQ1Sdc758aX82Zmdp2Z7Taz3YcO\nHVrOSwAAALSsuoKVmSXlhaqbnXNfn+Mpz0raUvPzZv++YzjnPu2c2+mc27lu3brlrBcAAKBl1XNV\noEn6a0l7nHMfm+dpt0i61r868GJJY865/SGuEwAAoOXVc1XgiyW9XdKDZnaff99/k7RVkpxzN0q6\nVdKrJe2VNCXpHeEvFQAAoLUtGqycc3dKskWe4yS9N6xFAQAAtCMmrwMAAISEYAUAABASghUAAEBI\nCFYAAAAhIVgBAACEhGAFAAAQEoIVAABASAhWAAAAISFYAQAAhIRgBQAAEBKCFQAAQEgIVgAAACEh\nWAEAAISEYAUAABASglUTOeeiXgIAAGigRNQL6HTOOT1+OKfv/fSQnhrO6bcuP039XcmolwUAABqA\nYNUgVef0yP4Jfe+nB/XM6LSScVOp4nRwPE+wAgCgQxGsGuDp4Sl9/Ns/08GJgtZ2J/W6F56ok9Zl\n9Rf/9jONTZeiXh4AAGgQglUDfOO+Z3VwoqA3nr9Z52xeo3jMVK5WZRLBCgCADkawaoDhXFHpREzn\nbV07c18iFlNPOkGwAgCgg3FVYAOM5IrKpo/PrH1dSYIVAAAdjGDVAMO5grKp+HH39xOsAADoaASr\nBhienLti1d+V1HieYAUAQKciWDXASK6onnmCVb5UVaFUiWBVAACg0QhWIXPOzdtjFcyvYjsQAIDO\nRLAK2fh0WeWqm7PHqi8IVmwHAgDQkQhWIRvOFSRp4YrVFMEKAIBORLAK2UiuKGnuYNXX5d1HxQoA\ngM5EsArZ8ALBKhgSOk6PFQAAHYlgFbLhST9YzdFjJTHLCgCATkawCtnIAj1WEsEKAIBORrAK2bA/\nwyoZn/uj5VgbAAA6F8EqZCO5ogayqXkfnxkSWmZIKAAAnYZgFbLhycWDlcSQUAAAOhHBKmTDuaIG\n6whW49PlZi0JAAA0CcEqZCO5ggZ76qlYFZu1JAAA0CQEqxAF5wQOZNPzPqcv4w8JZSsQAICOQ7AK\n0Xi+rFLFLbgVmIjHlE0nNMZWIAAAHYdgFaLgOJuFmtclqb8rwVYgAAAdiGAVomA46EI9VpLU35Wi\neR0AgA5EsApRcJzN4AI9VlJQsaLHCgCATkOwCtHMVuBiFatMUtOliqaKVK0AAOgkBKsQDeeCitUi\nwarbG7mwfyzf8DUBAIDmIViFaHiyqGwqrkwyvuDz+vxZVgcIVgAAdBSCVYhGcoVFtwElbytQkp47\nMt3oJQEAgCYiWIVoeJHhoAEqVgAAdCaCVYiGJxc+JzCQjMeUTcW1f5xgBQBAJyFYhWhkkQOYa/V3\nJbWfrUAAADoKwSokM+cE1tFjJfnBiq1AAAA6SiLqBbSjz9/99HH35UsVFStVPT08Nefjs/V1JfXo\n8xONWB4AAIgIFauQ5AresM9sqr6s2t+V1JGpkqaLlUYuCwAANBHBKiQzwSpdf7CSpP1j9FkBANAp\nCFYhyfmVp2x64eGggX5GLgAA0HEIViFZbsXqOYIVAAAdg2AVkskl9lgdHRLKViAAAJ2CYBWSXKGs\nVDymVKK+jzQZj2kgm2LkAgAAHYRgFZJcsVJ3f1VgY1+GYAUAQAchWIUkVyjX3V8VOHENwQoAgE5C\nsApJrlCuu78qsLE/Q48VAAAdZNFgZWY3mdlBM/vJPI9fZmZjZnaf//Xh8JfZ+rytwKUFqxP6uzTK\nkFAAADpGPRWrz0m6YpHn3OGcO9f/+sOVL6u9OOf8rcCl9Vid0J+RJB0YZzsQAIBOsGiwcs59X9JI\nE9bStorlqspVt6ytQEnaf4TtQAAAOkFYPVaXmtkDZnabmZ0V0mu2jcklDgcNrO/1gtWhyULoawIA\nAM23tCQwt3slbXXOTZrZqyX9o6RT5nqimV0n6TpJ2rp1awhv3RqC42x6lrgVOJhNSZKGJ4uhrwkA\nADTfiitWzrlx59ykf/tWSUkzG5rnuZ92zu10zu1ct27dSt+6ZSz1OJtAf1dS8ZhpOEfFCgCATrDi\nYGVmG83M/NsX+q85vNLXbSe5JR5nE4jFTGu7UxrJUbECAKATLJoEzOwLki6TNGRm+yR9RFJSkpxz\nN0q6WtJ7zKwsaVrSm51zrmErbkHLrVhJ0lBPiq1AAAA6xKJJwDl3zSKP3yDphtBW1IZyxYqScav7\nnMBaA9mUhqlYAQDQEZi8HoLlHGcTGMiyFQgAQKcgWIUgV1z6cTaBoZ60hhm3AABARyBYhSBXqCx5\n6npgIJvSeL6sYrka8qoAAECzEaxCMFkoq2eZW4GDPd4sq9EptgMBAGh3BKsVmjkncJlbgcGQ0MNs\nBwIA0PYIVitUrPjnBC67eT0tSTSwAwDQAQhWK5QreMfZLLfHKtgKJFgBAND+CFYrtJLhoFLtViDB\nCgCAdkewWqHlHmcT6Mt45wWOcF4gAABtj2C1QrniyipWsZgxJBQAgA5BsFqhlfZYSd52IFuBAAC0\nP4LVCk0WykrETKn48j/KwR4qVgAAdAKC1QpNlyrqTsVlZst+jYEsx9oAANAJCFYrlC9VlE4ufxtQ\n8rYCh6lYAQDQ9ghWK1QoVZVJrOxjHMymNMF5gQAAtD2C1QrlyxV1pVZWsRpgSCgAAB2BYLVC+VJF\n6cRKtwK9Y22GmWUFAEBbI1itUL5UVWalPVZ+xWqYkQsAALQ1gtUK5UsVZZIr+xgHsmwFAgDQCQhW\nK1CuVFWuuhVXrIZmtgIJVgAAtDOC1Qrk/av4VnpVYF9XQomYMcsKAIA2R7BagXzJO85mpRUrM9Na\nzgsEAKDtEaxWIKxgJTEkFACATkCwWoF8yd8KDCNY9aTYCgQAoM0RrFbgaMVq5R/jYDbNViAAAG2O\nYLUChbIfrFY4IFTyRi4wxwoAgPZGsFqB6TC3ArMpTRTKM2ENAAC0H4LVCgRbgekwtgJ7vFlWo7nS\nil8LAABEg2C1AoVSRalETDGzFb9WMH39MA3sAAC0LYLVCuRLVXWFsA0oSUM9HGsDAEC7I1itQL5c\nUXqFU9cDQcVqOEfFCgCAdkWwWgHvAOZwKlaDwXmBXBkIAEDbIlitQL5UDWWGlXT0vEC2AgEAaF8E\nqxUIs2JlZsyyAgCgzRGsViBfqoQyHDQwwHmBAAC0NYLVCuTL4W0FStJQT1ojNK8DANC2CFbLVKpU\nVam60LYCJSpWAAC0O4LVMh09gDm8YDXYk9IIPVYAALQtgtUyFWbOCQzvI+S8QAAA2hvBapnyfvgJ\nt3ndm2XFyAUAANoTwWqZpmcOYA53K1BiSCgAAO2KYLVM+QZtBUqigR0AgDZFsFqmQkOa14OtQEYu\nAADQjghWyxRcFdgV8rgFia1AAADaFcFqmfLlqkxSKhHeR9iXSSgZN7YCAQBoU4moF9Cu8qWKUomY\nYmbLfo3P3/30cfd1JePa9cTIzGNvuWjrsl8fAAA0FxWrZcqXqqH2VwWy6YQmC+XQXxcAADQewWqZ\n8qVKqFcEBrKphHIEKwAA2hLBapnypUqow0ED2XRcuSKT1wEAaEcEq2XKlysN2wqkYgUAQHsiWC2T\n12MV/sfXk06oUK6qVKmG/toAAKCxCFbL5PVYNaBilfIu1KRqBQBA+yFYLYNzToWGXRXovSZ9VgAA\ntB+C1TKUq04V55QJcThoIJumYgUAQLsiWC3DtH+cTbpBzesSwQoAgHZEsFqGfAMOYA70+MGKIaEA\nALQfgtUyFEreFXuNuCownYgpGTdN5AlWAAC0G4LVMgQVq64GVKzMTL2ZpCbypdBfGwAANBbBahny\nZa9i1YgeK0nqTSc0TsUKAIC2Q7BahpkeqwZcFShJvZkEW4EAALShRZOBmd1kZgfN7CfzPG5m9nEz\n22tmD5jZi8JfZmtpZPO6JPV2sRUIAEA7qqfk8jlJVyzw+JWSTvG/rpP0qZUvq7XlSxWZpFSDKlZ9\n/rE2xTLH2gAA0E4WTQbOue9LGlngKVdJ+jvn+ZGkNWZ2QlgLbEX5UlXpZEwxs4a8fm8mKUlUrQAA\naDNhlFw2SXqm5ud9/n0dK1+qKJNozDag5PVYSaLPCgCANtPU5nUzu87MdpvZ7kOHDjXzrUOVLzfm\nnMDATMWKIaEAALSVMILVs5K21Py82b/vOM65Tzvndjrndq5bty6Et45GvlRpyHDQQFCxGp9mKxAA\ngHYSRjq4RdK1/tWBF0sac87tD+F1W1ahVGloxao7FVfcmL4OAEC7SSz2BDP7gqTLJA2Z2T5JH5GU\nlCTn3I2SbpX0akl7JU1JekejFtsq8uWq1jcwWHnT1xM0rwMA0GYWDVbOuWsWedxJem9oK2oD08WK\n0g0atRDozSTosQIAoM0weX2JnHMqlBu7FSiJ8wIBAGhDBKslmi5VVHWNm7oe4FgbAADaD8FqiYKw\n08irAiUvWE0VKyqUKw19HwAAEB6C1RIF23PN2AqUpEMThYa+DwAACA/BaonGg4pVAyevS1KfP8vq\nIMEKAIC2QbBaouZtBXoVq4PjBCsAANoFwWqJgmnozWhel6RDE/mGvg8AAAgPwWqJjlasGhussumE\nTGwFAgDQTghWSzTTvN7gAaExM/VkEmwFAgDQRghWSzSRL8skpRocrCRvO/B5tgIBAGgbBKslmsiX\nlEnGZWYNf6/edJKKFQAAbYRgtUQT+XLDrwgM9HUl6LECAKCNEKyWaDxfbnjjeqA3k9RwrqBypdqU\n9wMAACtDsFqi8XxJ6QYPBw30ZhJyThrOFZvyfgAAYGUIVkvUzK3A3jRDQgEAaCcEqyUKmtebIRgS\n+vw4VwYCANAOCFZLNNHUHivOCwQAoJ0QrJbAOafJQhO3AoPzApllBQBAWyBYLcFUsaJK1SnTpOb1\neMw0mE1RsQIAoE0QrJagWecE1lrXm6Z5HQCANkGwWoLx4JzAJm0FStL6vowOsRUIAEBbIFgtwcwB\nzE2sWK3vTbMVCABAmyBYLcF4sBXYhAOYA+t70zo0UVC16pr2ngAAYHkIVksQRY/V+t60ylWnkSmm\nrwMA0OoIVksQxVbghr6MJKavAwDQDghWSxBJxaovLYlZVgAAtAOC1RJM5EuKx0zJuDXtPdf3+hUr\nGtgBAGh5BKslGJ8uqzeTkFnzgtW6Xq9idYhgBQBAyyNYLcFEvjRzfl+zZJJx9WUSHMQMAEAbIFgt\nwXi+rD7//L5mWt+XoXkdAIA2QLBagpFcUQPZVNPf1xsSSsUKAIBWR7BagtGpotZ2Nz9YbejL0LwO\nAEAbIFgtQbQVq4KcY/o6AACtjGBVp2K5qol8OZJgta43rWK5qvHpctPfGwAA1I9gVacj/pEykVSs\ngunr9FkBANDSCFZ1GokyWPmzrJ7nykAAAFoawapOI5NesIqieT0IVlSsAABobQSrOgUVq8GeKLcC\nqVgBANDKCFZ1GslFV7HqSSfUk07owBgVKwAAWhnBqk5BsFrT3fzJ65K0daBbTw3nInlvAABQH4JV\nnUZzRfV3JZWMR/ORbR/q1pPDU5G8NwAAqA/Bqk7DEQ0HDWwfzOqZkSmVK9XI1gAAABZGsKrT6FT0\nwapcdXr2yHRkawAAAAsjWNVpJFeKpHE9sH0oK0lsBwIA0MIIVnUayRU0kI2mcV3yeqwk6cnDNLAD\nANCqCFZ1cM5pNFfSQDYd2RrW9aSVTcX1JFcGAgDQsghWdZgslFWsVCOtWJmZtg1mqVgBANDCCFZ1\nGM2VJEUzHLQWIxcAAGhtBKs6RHmcTS1GLgAA0NoIVnUYyXln9EVesWLkAgAALY1gVYcRfytwMMLm\ndYmRCwAAtDqCVR1GgwOYI2xelxi5AABAqyNY1WE4V1QybupJJyJdByMXAABobQSrOoz65wSaWaTr\nYOQCAACtjWBVh+FcMfLG9QAjFwAAaF0EqzpEfQBzLUYuAADQughWdQi2AlsBIxcAAGhdBKs6DLdS\nsGLkAgAALYtgtYhypaqx6VLrBKtBRi4AANCq6gpWZnaFmT1qZnvN7PfmePwyMxszs/v8rw+Hv9Ro\nHJn2hoO2SrBa1+uNXHiCYAUAQMtZdDCTmcUlfVLS5ZL2SdplZrc45x6e9dQ7nHOvbcAaIzUSDAdt\nkasCg5ELTzHLCgCAllNPxepCSXudc48754qSvijpqsYuq3UEwWqwRSpWEiMXAABoVfUEq02Snqn5\neZ9/32yXmtkDZnabmZ0VyupawEzFqpWCFSMXAABoSWE1r98raatz7hxJn5D0j3M9ycyuM7PdZrb7\n0KFDIb11YwXBqlV6rCRGLgAA0KrqCVbPStpS8/Nm/74Zzrlx59ykf/tWSUkzG5r9Qs65Tzvndjrn\ndq5bt24Fy26e0RbrsZIYuQAAQKuqJ1jtknSKme0ws5SkN0u6pfYJZrbR/IP0zOxC/3WHw15sFIZz\nRfWmE0olWmcyBSMXAABoTYteFeicK5vZ+yTdLiku6Sbn3ENm9m7/8RslXS3pPWZWljQt6c3OOdfA\ndTfN6FRRAz2tU62SGLkAAECrWjRYSTPbe7fOuu/Gmts3SLoh3KW1hpEWOoA5wMgFAABaU+vsb7Wo\nkRY6zqYWIxcAAGg9BKtFtNIBzLUYuQAAQOupaytwtXLORX4A8+fvfnrO+w+M5VWuOn3qu49psCet\nt1y0tckrAwAAs1GxWsB0qaJCudpyPVaSNNiTluRdtQgAAFoDwWoBrXicTWDQv1JxeLIQ8UoAAECA\nYLWAVjzOJhDM1jpMxQoAgJZBsFpAKx5nEzAzretJ6/nxfNRLAQAAPoLVAkanWjdYSdKWgW7tG5lW\npdoRs1gBAGh7BKsFDE/6waoFm9cladtgt4qVqg6MUbUCAKAVEKwWMDpVVDxm6utqzakU2weDw5iZ\nwA4AQCsgWC0gOM7GP1+65fR3JbWmK6mnRpjADgBAKyBYLcA7ziYZ9TIWtG2wW08N59QhZ14DANDW\nCFYLGM2VWrZxPbBtMKuJfFn7RqejXgoAAKsewWoBw7lCGwSrbknSridHIl4JAAAgWC1gdKr1K1Yb\n+jJKJ2La/dRo1EsBAGDVI1jNo1J1OjJVbNlRC4GYmbYOdOueJwlWAABEjWA1j7HpkqquNY+zmW3b\nYFaPPj+hsalS1EsBAGBVI1jNo5WPs5kt6LO692mqVgAARIlgNY92ClZb1nYrETMa2AEAiBjBah5B\nsFrb4j1WkpRKxHTWiX00sAMAEDGC1TyCA5gHe1o/WEnSzu0Duv+ZIyqWq1EvBQCAVYtgNY92qlhJ\n0s5ta1UoV/WT58aiXgoAAKsWwWoeI7misqm4Msl41Eupy/nb10qSdtNnBQBAZAhW8xjNFdti1EJg\nfW9G2wa7tZt5VgAARIZgNY/hXLEtrgistXPbgO55apQDmQEAiAjBah4juWLb9FcFdm5fq+FcUU8c\nzkW9FAAAViWC1RyqVacnDue03R+82S4uCPqsGLsAAEAkCFZzePbItCYLZZ22sS/qpSzJSUM9WtOd\npIEdAICIEKzmsGf/uCTp9BN6I17J0sRipgu3D+j7Pz2scoV5VgAANBvBag6PHpiQJJ22ob2ClST9\n8os26cB4Xt/76aGolwIAwKpDsJrDIwcmtG2wW9l0IuqlLNkvnLFBQz1pfeHHz0S9FAAAVh2C1Rz2\nHBhvy2qVJCXjMb1x52b9+yPP68BYPurlAACwqhCsZsmXKnrycE6nn9Bejeu13nzBFlWd9JXdVK0A\nAGgmgtUsP3t+UlUnnbGxPStWkrRtMKsXv2BQX9z1jKpVhoUCANAsBKtZ9hzwrgg8rY2DlSRdc+FW\nPXtkWnfsPRz1UgAAWDUIVrM8sn9CmWRM2wazUS9lRS4/c4MGsil94e6no14KAACrBsFqlkef9xrX\n4zGLeikrkk7EdfX5m/Vve57XwQma2AEAaAaCVQ3nnPbsn9DpbTZxfT5vumCLylWnr96zL+qlAACw\nKrTfoKYGOjRZ0Eiu2Jb9VZ+fZ8tvx1BWn73jCfVlknrbxduavCoAQKuoVp0mi2VN5suayJc1kS9p\nPF/SRL6s8emSxvNlTRa8xycL3nMmCyVVq1IsJsXMFI+ZYmbavLZLLz1lSJecPKT+ruTMe8z3d1Gt\nt1y0tZF/zMgRrGo8st+buN5uR9ks5ILtA/ry7mf0+KFc1EsBAKxAuVLVRL6ssemSxqa9UDQ+XZ65\nPTZd0vh0qSY0ed8n8kFIKi/6HjGTMsm40omYMsm4UomYYmZyztvVqTqnqpN++Piwbr77aZmkzWu7\n9IL1vfq5zf3a2Jdp/AfR4ghWNYKjbDplK1CSzjqxT13JuHZxMDMARCqoGI1PHxuIgmrR+KzANPNY\nTTVpIUEo6krGvXCUjCmTiOuE/i7tGIopnYgr49+XTsZmnptOxmZ+JxEzmS3eY1ypOj09MqW9Bye1\n9+CEvvvoQX3/p4f0+vM26fxta8P6yNoSwarGngPjWt+b1kA2FfVSQpOMx7Rz21rdufewdj85op3b\nB6JeEgC0rXyp4gedksamyzO35wpHM5UlP0RN5EtabLRgEHy6UvGZ4LPRD0aZlPdzEIIySe95wX3J\neH2hKAzxmGnHUFY7hrK6/MwNmiyU9aVdT+tr9+7T8+N5XXH2RsWatJZWQ7Cq8cj+ibaeuD6fl5++\nXj95bkwf/Mr9uvX9L23LMxABIAyVqtNk3q8G1VSGgm20ICAFoWj2V7FcXfD1k3FTJhE/JgQNZFM6\ncU1m3lBUWzlq1zDSk07oVy7doX9+cL/u3HtYByfyevMFW5VJxqNeWtPxN6yvXKlq78FJvfSUoaiX\nErp0Mq6rz9+iz975uP70tj3649f/XNRLAoAlc85pqljxmq1nqkReOAp6imp/npijx2ix7TSTZrbG\naqtG2wa6j71vVijK+L+TiK/ei+3jMdPrXniiNvSl9U/3P6dPfe8xXXvxNg32pKNeWlMRrHxPHM6p\nWKl2VON6rR1DWf3qS3boM3c8ocvP3KifP3Vd1EsCsMrM7jGqvSJtJhjN6jE6NkSVVVlkLy0RM78i\n5PUQBf1Em9Z0KZOMKX1M1ehon1FtX1K7Vo1axUU7BrWuJ62b735aN37vMb3v5accc+VgpyNY+fb4\njeunbei8rcDAB195mr776CH9zlfv1+3Xv0xrujunlwxAYznnNF2qzLrizLt97OX7x1eNgmA0kS9r\nsdNLU/HY0VDkB541XUlt7MscWyFKxZVJHH1e8DvJVVwxaiUnrevRdS87SZ/67mP68u5n9K6X7Fg1\ngZVg5Xv0wLgSMdPJ69v7KJuFZJJx/Z83navXf/IuffgbD+nj15wX9ZIANFGhXJnpJZrdWD3TY1R7\nRdoxFaWyyotUi4JttJmeIf8qtHU9aW3xt9K8KlJsZjstk4gfE6Ta/dQLHLWhL6PXvfBEffXeffru\nowf18tM3RL2kpiBY+R7ZP6GT1/UonejsRruzN/Xr/b9wij72rZ/qlWdt0GvPOTHqJQGoU9BjNHPF\n2dSxjdbjswLS7MbrfGnhxutUPKaulDfDKAhBQbXomMpQTRhK1wSlYOYREDhv6xrtPTSpb+85qB1D\nPdox1LnFiwDByvfIgYlVM3vjNy47Wd9+5KD+29cf1An9GZ2/jREMQLMEW2pB2DkyVTouGM13/3i+\npFJl4apROhE75hL8rlRcm9d265T18ePu76rZVuuiWoQGMDNd9cIT9czIlL6062m9/+WnRL2khiNY\nSRqbLunZI9N668WdPWY/kIjHdMM15+nam36saz5zt/786nN01bmbol4W0Dacc8oVKzNVo7mC0VwV\no+DxhcKRScfOJ0rFlU0nNNSTnjMMHb1SzasaUTFCq0kn43rzhVt14/ce01fv3ad3vXRH0+ZtRYFg\nJemnz3uN62d00MT1xWwZ6NbX33Opfv3v79EHvnifHj+U0/WvOKWj/2EHZitXqjoyUx0q6siUd/vI\ndEljU0WvclRTQaoNTgv1Gx0XjpJxdacSGsym/duzwlEqPnM/22noRJvWdOnKszfqmw/s19/c9aTe\n+ZIdUS+pYQhWkh7ZPy6ps84InMtch2O+9pwTVK44/eW3f6bvPHpQX/71S1blQDe0t2B7bXSqpNGc\nF4hGp4pn339YAAAOqElEQVQanSrpSK6okSkvNI1OFTWa8+4/OJFfsOdo3spRb/r47bRZ1SPCEXC8\nS04a1GMHJ/Wnt+3R+dvW6oVb1kS9pIYgWMkbtdCXSazKwyMTsZje8KJNWteT0u0PP69rPvMj/dFV\nZ+vsTf1RLw2r2FSxrOHJokb8UDSa826PThU1kiv54civME17QWmhidjpREzdKa9q1J2Ka213UpvW\ndPn3xdXl39+dOlpdYp4REC4z0xvO36y/uetJ/cbN9+qf3/+Sjhz7Q7CSd/jy6Sf0rdptMDPTz5+2\nXoM9af3TA8/ptZ+4U686a4Ouf8WpOqMDj/hB85UrVY1MFXV4oqjDkwUN5woanizq8GRRw5MFDee8\n74f9MDVdqsz5OjGTulIJZf0Q1J1KaPPabp264Who8oLSsT8nYsw2AlpBdyqhT771RXrjjT/Qb335\nfn322p2KddhFE6s+WD0zMqWfPDumt1y0OhrXF3L2pn797pWn66Y7n9BNdz6h2x+6Q1eevVHXv+JU\nnbaxs7dJsXSlStUPRwUdmizo8IQXjA5NFHR48ujXvtFpTRcrcw6GjMdMPemEsum4sqmE1vemtWMo\nq2zaC0/ZdGLmdncqoUwytmr/AwjoFOduWaM/eM2Z+sgtD+nG7z+m37jsBVEvKVSrOlg55/Shrz+o\nZDymX3vpSVEvpyX0dyX1m5efqne+eIc+e+fj+pu7ntRtPzmgczb365VnbtDlZ27UqRt6+MutQ+VL\nFT8QBRUk7/bM95rQNDpVmvM1UvGYejIJ9aS9r7NP7D/m55mvTELpBEEJWI2uvWSbdj05ov91+6M6\nb8taXXLyYNRLCo05t9gBA42xc+dOt3v37kjeO/ClXU/rd7/2oP749WfrbRdvq/v35moC71RThbJ2\nPTWqh58b0zOj05KkgWxKv3TeJl2wfa3OOrFfm9d28ZdjiyqUKxrJFWf6lYItuBG/Z+nw5NH7hicL\nyhXn3oJLJWLHhaLe2u/phHoySfWkE0ol2HYDML9gh2iyUNbrbrhT49Nl3fqBl2h9b2v3OZvZPc65\nnYs+b7UGqwNjeV3+se/pzBP79IVfu3hJe7yrKVjVGp8uac+Bce3ZP64nD0+pWPGahfu7kjrrxD6d\nvalfJw1ltWWgW5vXdumE/i7+kg2Jc06ThaMTtoPZSUemSxrJFXXEb+o+MlXUcO5ocJoslOd8vZhJ\n2ZS3zRZsxfWka38+GqKyhCUAIaptvXn0wISu+uSdOnfLGv3tOy9s6dNP6g1Wq3Ir0Dmn3/+HB1Wq\nVvVnbzin4xrnGqWvK6mLdgzqoh2DKlWqOjCW13Nj03ruSF5PDU/px0+MHDPbx0za2JfRxv6MBrNp\nretNaTCb1lBPSmuzKfV1JdXflVRfxvve20FbQ8455UtVTRXLmi5VNF2saKpYUa5Y1nSxolyxoqlC\nWbliRZP5snLFsiYLZeUKZe9A20K55oy2kiYLZS10TFsybsoGzdrphNZ2J7V5bZffn+TdH4SkbDrO\nIEkALeG0jb36k1/6Of3Wl+/XG2/8oW645kXaOtgd9bJWZFUGq1vuf07ffuSg/uA1Z2j7Kji3qBGS\n8Zi2DHRry8DR/wNUqk7j+dLMnKBgZtDEdFn7j+RngsNCNdJ4zGaalrv972n/HLJ0IqZ0wvuejMeU\niJuS8ZiScVMiHlPcTDGTYjHzbsdMQXYweTfMJOckJ+d9d973inOqVp0qzqlc9W6Xq06lSlXlilMx\n+F6uqlipet/LVRUqVRVKFRXKVeVLFU2XKsr7Py+lGJyMm1KJ4M8YO3p4bW9aWwa6lEkcnZWUSR69\n8q0r6X1GyTgVJQDt6ZdftFnZdEL/9Sv36zWfuEN/fvU5uuLsE6Je1rLVFazM7ApJfykpLumzzrmP\nznrc/MdfLWlK0q845+4Nea2hODxZ0H+/5SGdu2WN3vHizp38GoV4zLS2O6W1C8wlqfqHyE4Vy8oX\nK5ouVWfCSL5UmQkrM8GlVNVkoaxyxalcqapU9b5Xqk4VJ1WqVVWrXqirOrdgaKtHzKSYeYEsZqZ4\nzJSI2UxYi8dMibgpEYvNPJaIx9SbSWigO1UT9mJKJbzQl5q57YWmlP+VTsRnHuOMNgCr2avO2qgz\nT+jT+z5/r9799/fqVy7drg+9+vSW3hqcz6LByszikj4p6XJJ+yTtMrNbnHMP1zztSkmn+F8XSfqU\n/72lPDMypT/85sPKFSr6n1efw19mEYiZzfTuNErVr0JVZ5WMan808yZry7xqVhCkAADR2DLQra+8\n+1J99LZHdNNdT+iep0b1rpfs0CUnD2pDGw3wrudvtwsl7XXOPS5JZvZFSVdJqg1WV0n6O+d1wv/I\nzNaY2QnOuf2hr3gJDozl9cPHD+sHe4f1w8eHtc+/qu13rzhdp25gLlOnivmpKS6CEgC0k1Qipg//\n4pm66KQB/f4/PKjrv3SfJOmkdVldctKgLjl5UJvWdM1cnRz0kbZSoaSeYLVJ0jM1P+/T8dWouZ6z\nSVJkweqep0b1hk/9QJJ31drFJw3oV1+yQ5e+YIhQBQBAC3vVWRt1+Rkb9PD+cf3wMa848o37ntPN\n81yVf+0l2/SHV53d5FXOranN62Z2naTr/B8nzezRZr33A5I+Hd7LDUk6HN7LoQ585s3HZ958fObN\nx2feZG9twGv+kf/VYHUNvKwnWD0raUvNz5v9+5b6HDnnPq1Q8000zGx3PbMsEB4+8+bjM28+PvPm\n4zNH2Oq5RnuXpFPMbIeZpSS9WdIts55zi6RrzXOxpLGo+6sAAACabdGKlXOubGbvk3S7vHELNznn\nHjKzd/uP3yjpVnmjFvbKG7fwjsYtGQAAoDXV1WPlnLtVXniqve/GmttO0nvDXVpLa/vtzDbEZ958\nfObNx2fefHzmCFVkZwUCAAB0Gs7BAAAACAnBagnM7Aoze9TM9prZ70W9ntXAzG4ys4Nm9pOo17Ia\nmNkWM/uOmT1sZg+Z2QeiXlOnM7OMmf3YzO73P/P/EfWaVgszi5vZf5jZN6NeCzoHwapONUf7XCnp\nTEnXmNmZ0a5qVficpCuiXsQqUpb0QefcmZIulvRe/jlvuIKklzvnXijpXElX+FdXo/E+IGlP1ItA\nZyFY1W/maB/nXFFScLQPGsg5931JI1GvY7Vwzu0PDlB3zk3I+0tnU7Sr6mzOM+n/mPS/aH5tMDPb\nLOk1kj4b9VrQWQhW9Zvv2B6gI5nZdknnSbo72pV0Pn9L6j5JByV9yznHZ954fyHpdyRVo14IOgvB\nCsBxzKxH0tckXe+cG496PZ3OOVdxzp0r79SKC82sNQ4961Bm9lpJB51z90S9FnQeglX96jq2B2h3\nZpaUF6puds59Per1rCbOuSOSviP6ChvtxZJeZ2ZPymvreLmZ/X20S0KnIFjVr56jfYC2ZmYm6a8l\n7XHOfSzq9awGZrbOzNb4t7skXS7pkWhX1dmccx9yzm12zm2X9+/yf3fOvS3iZaFDEKzq5JwrSwqO\n9tkj6cvOuYeiXVXnM7MvSPqhpNPMbJ+ZvSvqNXW4F0t6u7z/gr/P/3p11IvqcCdI+o6ZPSDvP+C+\n5Zzj8n+gTTF5HQAAICRUrAAAAEJCsAIAAAgJwQoAACAkBCsAAICQEKwAAABCQrACAAAICcEKwKLM\nrOLPtHrIzO43sw+aWcx/bKeZfXyB391uZm9p3mqPe+9p/xy+lmBmbzKzvWbGrCqgAxGsANRj2jl3\nrnPuLHmTwa+U9BFJcs7tds69f4Hf3S4pkmDle8w/h69uZhZv1GKcc1+S9KuNen0A0SJYAVgS59xB\nSddJep95LguqL2b28zUT2//DzHolfVTSS/37ftOvIt1hZvf6X5f6v3uZmX3XzL5qZo+Y2c3+ETsy\nswvM7Ad+tezHZtZrZnEz+3Mz22VmD5jZr9ezfjP7RzO7x6++XVdz/6SZ/W8zu1/SJfO851n+7fv8\n9zzF/9231dz/f4NgZmZX+H/G+83s2yH+zwCgRSWiXgCA9uOce9wPD+tnPfTbkt7rnLvLzHok5SX9\nnqTfds69VpLMrFvS5c65vB9MviBpp//750k6S9Jzku6S9GIz+7GkL0l6k3Nul5n1SZqW9C5JY865\nC8wsLekuM/tX59wTiyz/nc65Ef9cvl1m9jXn3LCkrKS7nXMf9M8DfWSO93y3pL90zt3sPyduZmdI\nepOkFzvnSmb2V5Leama3SfqMpJc5554ws4Elf9AA2g7BCkCY7pL0MTO7WdLXnXP7/KJTraSkG8zs\nXEkVSafWPPZj59w+SfL7orZLGpO03zm3S5Kcc+P+46+UdI6ZXe3/br+kUyQtFqzeb2a/5N/e4v/O\nsL+Wr/n3nzbPe/5Q0u+b2Wb/z/czM/sFSefLC2mS1CXpoKSLJX0/CHrOuZFF1gWgAxCsACyZmZ0k\nL4gclHRGcL9z7qNm9s+SXi2vgvSqOX79NyU9L+mF8toR8jWPFWpuV7Twv6NM0n9xzt2+hHVfJukV\nki5xzk2Z2XclZfyH8865ykK/75z7vJndLek1km71tx9N0t865z40671+sd51Aegc9FgBWBIzWyfp\nRkk3uFmnuJvZyc65B51zfyZpl6TTJU1I6q15Wr+8alBV0tslLdYo/qikE8zsAv89es0sIel2Se8x\ns6R//6lmll3ktfoljfqh6nR5VaW639MPlI875z4u6RuSzpH0bUlXm9l6/7kDZrZN0o8kvczMdgT3\nL7I2AB2AihWAenT5W3NJSWVJ/0/Sx+Z43vVm9p8kVSU9JOk2/3bFbwr/nKS/kvQ1M7tW0r9Iyi30\nxs65opm9SdIn/L6oaXlVp8/K2yq8129yPyTp9Yv8Of5F0rvNbI+88PSjJb7nf5b0djMrSTog6U/8\nfq0/kPSv5o2gKMnrM/uR3xz/df/+g/KuqATQwWzWf3ACQMcws+2SvumcOzvipRzD35KcaegH0DnY\nCgTQySqS+q3FBoTKq9qNRr0WAOGjYgUAABASKlYAAAAhIVgBAACEhGAFAAAQEoIVAABASAhWAAAA\nIfn/UY714KZMTbMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f32afc70a58>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nb_merge_dist_plot(\n",
    "    SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n",
    "    SkyCoord(des['des_ra'], des['des_dec'])\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Given the graph above, we use 0.8 arc-second radius\n",
    "master_catalogue = merge_catalogues(master_catalogue, des, \"des_ra\", \"des_dec\", radius=0.8*u.arcsec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cleaning\n",
    "\n",
    "When we merge the catalogues, astropy masks the non-existent values (e.g. when a row comes only from a catalogue and has no counterparts in the other, the columns from the latest are masked for that row). We indicate to use NaN for masked values for floats columns, False for flag columns and -1 for ID columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for col in master_catalogue.colnames:\n",
    "    if \"m_\" in col or \"merr_\" in col or \"f_\" in col or \"ferr_\" in col or \"stellarity\" in col:\n",
    "        master_catalogue[col].fill_value = np.nan\n",
    "    elif \"flag\" in col:\n",
    "        master_catalogue[col].fill_value = 0\n",
    "    elif \"id\" in col:\n",
    "        master_catalogue[col].fill_value = -1\n",
    "        \n",
    "master_catalogue = master_catalogue.filled()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "&lt;Table length=10&gt;\n",
       "<table id=\"table139858957025744-881306\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>idx</th><th>omegacam_intid</th><th>ra</th><th>dec</th><th>omegacam_flag_merged</th><th>atlas_id</th><th>atlas_stellarity</th><th>atlas_flag_cleaned</th><th>atlas_flag_gaia</th><th>kids_id</th><th>kids_stellarity</th><th>kids_flag_cleaned</th><th>kids_flag_gaia</th><th>flag_merged</th><th>ps1_id</th><th>ps1_flag_gaia</th><th>viking_id</th><th>viking_stellarity</th><th>viking_flag_cleaned</th><th>viking_flag_gaia</th><th>des_id</th><th>des_flag_gaia</th><th>des_stellarity</th><th>des_flag_cleaned</th></tr></thead>\n",
       "<thead><tr><th></th><th></th><th>deg</th><th>deg</th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead>\n",
       "<tr><td>0</td><td>730208058155</td><td>19.8749571434</td><td>-28.0734797175</td><td>False</td><td>730208058155</td><td>0.99999899</td><td>True</td><td>3</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>True</td><td>-1</td><td>0</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>-1</td><td>0</td><td>nan</td><td>False</td></tr>\n",
       "<tr><td>1</td><td>730244166129</td><td>340.813431421</td><td>-31.9482257128</td><td>False</td><td>730244166129</td><td>0.99998093</td><td>True</td><td>3</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>False</td><td>-1</td><td>0</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>-1</td><td>0</td><td>nan</td><td>False</td></tr>\n",
       "<tr><td>2</td><td>730211950982</td><td>341.549624302</td><td>-30.0666580713</td><td>False</td><td>730211950982</td><td>0.99999899</td><td>True</td><td>3</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>False</td><td>-1</td><td>0</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>-1</td><td>0</td><td>nan</td><td>False</td></tr>\n",
       "<tr><td>3</td><td>730236841450</td><td>16.4312524322</td><td>-31.951258374</td><td>False</td><td>730236841450</td><td>0.99999899</td><td>False</td><td>2</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>True</td><td>-1</td><td>0</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>-1</td><td>0</td><td>nan</td><td>False</td></tr>\n",
       "<tr><td>4</td><td>730196943452</td><td>350.851313992</td><td>-31.0398321606</td><td>False</td><td>730196943452</td><td>0.99998093</td><td>True</td><td>3</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>False</td><td>-1</td><td>0</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>-1</td><td>0</td><td>nan</td><td>False</td></tr>\n",
       "<tr><td>5</td><td>730240128354</td><td>339.45716753</td><td>-30.0428359203</td><td>False</td><td>730240128354</td><td>0.99999899</td><td>True</td><td>3</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>False</td><td>-1</td><td>0</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>-1</td><td>0</td><td>nan</td><td>False</td></tr>\n",
       "<tr><td>6</td><td>730149898783</td><td>22.3792388871</td><td>-32.0238211482</td><td>False</td><td>730149898783</td><td>0.99999899</td><td>False</td><td>2</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>False</td><td>-1</td><td>0</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>-1</td><td>0</td><td>nan</td><td>False</td></tr>\n",
       "<tr><td>7</td><td>730184690118</td><td>350.382068649</td><td>-35.9833619382</td><td>False</td><td>730184690118</td><td>0.99999899</td><td>True</td><td>3</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>False</td><td>-1</td><td>0</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>-1</td><td>0</td><td>nan</td><td>False</td></tr>\n",
       "<tr><td>8</td><td>730247028434</td><td>5.80262699957</td><td>-31.0384384284</td><td>False</td><td>730247028434</td><td>0.89999998</td><td>False</td><td>0</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>False</td><td>-1</td><td>0</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>-1</td><td>0</td><td>nan</td><td>False</td></tr>\n",
       "<tr><td>9</td><td>730154434261</td><td>345.195018672</td><td>-34.7782563206</td><td>False</td><td>730154434261</td><td>0.99998093</td><td>False</td><td>3</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>False</td><td>-1</td><td>0</td><td>-1</td><td>nan</td><td>False</td><td>0</td><td>-1</td><td>0</td><td>nan</td><td>False</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",
       "require.config({paths: {\n",
       "    datatables: 'https://cdn.datatables.net/1.10.12/js/jquery.dataTables.min'\n",
       "}});\n",
       "require([\"datatables\"], function(){\n",
       "    console.log(\"$('#table139858957025744-881306').dataTable()\");\n",
       "    $('#table139858957025744-881306').dataTable({\n",
       "        \"order\": [],\n",
       "        \"iDisplayLength\": 50,\n",
       "        \"aLengthMenu\": [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        \"pagingType\": \"full_numbers\"\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "master_catalogue[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III - Merging flags and stellarity\n",
    "\n",
    "Each pristine catalogue contains a flag indicating if the source was associated to a another nearby source that was removed during the cleaning process.  We merge these flags in a single one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "flag_cleaned_columns = [column for column in master_catalogue.colnames\n",
    "                        if 'flag_cleaned' in column]\n",
    "\n",
    "flag_column = np.zeros(len(master_catalogue), dtype=bool)\n",
    "for column in flag_cleaned_columns:\n",
    "    flag_column |= master_catalogue[column]\n",
    "    \n",
    "master_catalogue.add_column(Column(data=flag_column, name=\"flag_cleaned\"))\n",
    "master_catalogue.remove_columns(flag_cleaned_columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each pristine catalogue contains a flag indicating the probability of a source being a Gaia object (0: not a Gaia object, 1: possibly, 2: probably, 3: definitely).  We merge these flags taking the highest value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "flag_gaia_columns = [column for column in master_catalogue.colnames\n",
    "                     if 'flag_gaia' in column]\n",
    "\n",
    "master_catalogue.add_column(Column(\n",
    "    data=np.max([master_catalogue[column] for column in flag_gaia_columns], axis=0),\n",
    "    name=\"flag_gaia\"\n",
    "))\n",
    "master_catalogue.remove_columns(flag_gaia_columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each prisitine catalogue may contain one or several stellarity columns indicating the probability (0 to 1) of each source being a star.  We merge these columns taking the highest value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "atlas_stellarity, kids_stellarity, viking_stellarity, des_stellarity\n"
     ]
    }
   ],
   "source": [
    "stellarity_columns = [column for column in master_catalogue.colnames\n",
    "                      if 'stellarity' in column]\n",
    "\n",
    "print(\", \".join(stellarity_columns))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# We create an masked array with all the stellarities and get the maximum value, as well as its\n",
    "# origin.  Some sources may not have an associated stellarity.\n",
    "stellarity_array = np.array([master_catalogue[column] for column in stellarity_columns])\n",
    "stellarity_array = np.ma.masked_array(stellarity_array, np.isnan(stellarity_array))\n",
    "\n",
    "max_stellarity = np.max(stellarity_array, axis=0)\n",
    "max_stellarity.fill_value = np.nan\n",
    "\n",
    "no_stellarity_mask = max_stellarity.mask\n",
    "\n",
    "master_catalogue.add_column(Column(data=max_stellarity.filled(), name=\"stellarity\"))\n",
    "\n",
    "stellarity_origin = np.full(len(master_catalogue), \"NO_INFORMATION\", dtype=\"S20\")\n",
    "stellarity_origin[~no_stellarity_mask] = np.array(stellarity_columns)[np.argmax(stellarity_array, axis=0)[~no_stellarity_mask]]\n",
    "\n",
    "master_catalogue.add_column(Column(data=stellarity_origin, name=\"stellarity_origin\"))\n",
    "\n",
    "master_catalogue.remove_columns(stellarity_columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IV - Adding E(B-V) column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue.add_column(\n",
    "    ebv(master_catalogue['ra'], master_catalogue['dec'])\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## V - Adding HELP unique identifiers and field columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue.add_column(Column(gen_help_id(master_catalogue['ra'], master_catalogue['dec']),\n",
    "                                   name=\"help_id\"))\n",
    "master_catalogue.add_column(Column(np.full(len(master_catalogue), \"HATLAS-SGP\", dtype='<U18'),\n",
    "                                   name=\"field\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OK!\n"
     ]
    }
   ],
   "source": [
    "# Check that the HELP Ids are unique\n",
    "if len(master_catalogue) != len(np.unique(master_catalogue['help_id'])):\n",
    "    print(\"The HELP IDs are not unique!!!\")\n",
    "else:\n",
    "    print(\"OK!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## V.b - Add Specz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "specz =  Table.read(\"../../dmu23/dmu23_SGP/data/SGP-specz-ver2.0.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "specz.rename_column('RA', 'ra')\n",
    "specz.rename_column('DEC', 'dec')\n",
    "specz.rename_column('OBJID', 'specz_id')\n",
    "specz.rename_column('Z_SPEC', 'z_spec')\n",
    "specz.rename_column('Z_SOURCE', 'z_source')\n",
    "specz.rename_column('Z_QUAL', 'z_qual')\n",
    "specz.rename_column('REL', 'rel')\n",
    "specz.rename_column('AGN', 'agn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAF3CAYAAABnvQURAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8XHd97//Xd2akkWa077u1eLfjJfGS2AFCQkjYEuDS\nAKH0ltLmBpqWPnp7W27v7S38ehfae3+3QBuaG5a2FNIQUi44EAiQBJI4iWM73nfZ2vddo22kmfne\nPyS5inFs2ZrRmeX9fDz8QBqdzLyZONJb53zO92ustYiIiIjI0rmcDiAiIiKSLFSsRERERKJExUpE\nREQkSlSsRERERKJExUpEREQkSlSsRERERKJExUpEREQkSlSsRERERKJExUpEREQkSlSsRERERKLE\n49QLFxUV2draWqdeXkRERGTRDh482G+tLb7acY4Vq9raWg4cOODUy4uIiIgsmjGmZTHH6VKgiIiI\nSJSoWImIiIhEiYqViIiISJSoWImIiIhEiYqViIiISJSoWImIiIhEiYqViIiISJSoWImIiIhEiYqV\niIiISJSoWImIiIhEiYqViIiISJSoWImIiIhEiYqViIiISJR4nA4gyeOxfa2LOu7+nTUxTiIiIuKM\nRZ2xMsbcbYw5Y4xpNMZ89jJfv80YM2KMOTz3579EP6qIiIhIfLvqGStjjBt4GLgTaAf2G2P2WGtP\nXnLoi9ba98YgoySZxZ7ZAp3dEhGRxLKYM1Y7gEZr7QVr7TTwOHBvbGOJiIiIJJ7FFKtKoG3B5+1z\nj11qlzHmqDHmx8aYDVFJJyIiIpJAojW8/jpQY60dM8a8G/g+sOrSg4wxDwAPANTU6BKPiIiIJJfF\nnLHqAKoXfF4199hF1tpRa+3Y3MdPA2nGmKJLn8ha+6i1dpu1dltxcfESYouIiIjEn8UUq/3AKmNM\nnTEmHfgIsGfhAcaYMmOMmft4x9zzDkQ7rIiIiEg8u+qlQGttyBjzEPAM4Aa+Ya09YYx5cO7rjwAf\nAj5ljAkBk8BHrLU2hrlFRERE4s6iZqzmLu89fcljjyz4+G+Bv41uNBEREZHEoi1tRERERKJExUpE\nREQkSlSsRERERKJExUpEREQkSlSsRERERKJExUpEREQkSlSsRERERKJExUpEREQkSlSsRERERKJE\nxUpEREQkSlSsRERERKJExUpEREQkSlSsRERERKJExUpEREQkSlSsRERERKJExUqWXefwJG2DE07H\nEBERiTqP0wEkdYQjll+c6eX5M70A3LO5kh11BQ6nEhERiR4VK1kWg+PTPHGgjdbBCbZU5zExHeL7\nhzsYHJ/mnRtKcRnjdEQREZElU7GSmLLWcqh1mD1HO3EZ+PD2ajZX5RGOWJ460skL5/oYmpjmQzdV\nkebWlWkREUlsKlYSU3uOdLKvaZC6Ij+/dlMVeb50ANwuw71bKijwp/OTE92MTs7w6zevwO/VX0kR\nEUlcOkUgMdM/FmRf0yA7agv45K11F0vVPGMMb11dzEe2V9MxPMnXX2oiYq1DaUVERJZOxUpi5uXz\nA7iN4Y51JVecodpUlccHtlbSPTrFme7AMiYUERGJLhUriYnJ6TCvtwyxqSqX7Iy0qx6/qSqP3Mw0\n9jb2L0M6ERGR2FCxkpg42DLIdDjCrpVFizre7TLcUl/Ihf5xOocnY5xOREQkNlSsJOoi1vLKhQFq\nC31U5mUu+p/bXltAmtvw8vmBGKYTERGJHRUribpTXaMMTcywq2FxZ6vmZaa7ubEmnyPtwwSmZmKU\nTkREJHZUrCTqXj4/QJ4vjfUVOdf8z+5uKCIcsexrGoxBMhERkdhSsZKo6hyepKl/nFvqC69rNfWi\nbC9ry7LZd2GAmXAkBglFRERiR8VKourl8wOku11sW3H9ewDuaihifDrM0fbhKCYTERGJPRUriZrA\n1AxH2oe5cUUemenu636ehmI/ZTkZ7G0cwGrBUBERSSAqVhI1rzUPEo5YdtVf29D6pYwx7GoopHt0\nild0h6CIiCQQFSuJikjE8lrTIGtKsynK9i75+TZX5+H3evj6S01RSCciIrI8VKwkKk52jRKYCrGp\nKjcqz5fmdrGjtoBnT/fSOzoVlecUERGJNRUriYoXz81uRbOyJCtqz7mxcna5hufP9EbtOUVERGJJ\nxUqi4qXGPkpzvIvaF3CxynIyqMjN4NlTKlYiIpIYVKxkyaZmwuxvHmJlcfTOVsHsEPvb15bwUmM/\nwVA4qs8tIiISCypWsmT7mweZDkVYWZId9ee+fW0JE9NhXtNK7CIikgBUrGTJXjrXT5rbUFfkj/pz\n72oowutx6XKgiIgkBI/TAST+Pbav9Ypf33Okk6p8H+me6Pf0zHQ3uxoKee50L3/+vvWY69gmR0RE\nZLmoWMmSjAVDdI1Mcef60pg8/2P7WsnOSKN1cIIvPXuOkuyMyx53/86amLy+iIjItdClQFmS831j\nAFEfXF9obdns7NaZ7kDMXkNERCQaVKxkSRp7x8hIc1GZnxmz18jzpVOa41WxEhGRuKdiJdfNWktj\n7xgNxVm4Yjz7tLYsh+aBcaZmtOyCiIjELxUruW79Y9OMTM5EdbX1N7OmNJuIhXO9YzF/LRERkeul\nYiXXrXEZ5qvmVRf4yExzc6Z7NOavJSIicr1UrOS6NfaOke9LozDLG/PXcrsMq0uzONMdIGJtzF9P\nRETkeqhYyXUJRywX+saW5TLgvLVlOYxPh+kYmly21xQREbkWKlZyXdqHJgjGaBubN7OqNAsDnNbl\nQBERiVMqVnJdGnvHMEBDDLaxeTO+dA81hT4tuyAiInFLxUquS2PfGBV5mfi8y7t4/9qyHDpHphid\nnFnW1xUREVkMFSu5ZsFQmLbBiWWdr5o3/5pN/ePL/toiIiJXo2Il16xjeJKIhdpC37K/dnluBl6P\ni+YBFSsREYk/iypWxpi7jTFnjDGNxpjPXuG47caYkDHmQ9GLKPFm/q68yvzlL1YuY1hR6NMZKxER\niUtXLVbGGDfwMPAuYD3wUWPM+jc57i+Bn0Y7pMSX9qFJ8jLTyFrm+ap5tYV+egNBJoIhR15fRETk\nzSzmjNUOoNFae8FaOw08Dtx7meN+D/gXoDeK+SQOdQxPxnTT5aupLZy9E7F5YMKxDCIiIpezmGJV\nCbQt+Lx97rGLjDGVwAeAv4teNIlHE8EQg+PTVDlwGXBeVX4mHpfRnJWIiMSdaA2vfxH4E2tt5EoH\nGWMeMMYcMMYc6Ovri9JLy3JqH56dr6py8IyVx+2iKt+nYiUiInFnMcWqA6he8HnV3GMLbQMeN8Y0\nAx8CvmKMef+lT2StfdRau81au624uPg6I4uTOuaKVUWuc8UKoLbIR+fwJMFQ2NEcIiIiCy2mWO0H\nVhlj6owx6cBHgD0LD7DW1llra621tcCTwKettd+PelpxXPvQJEVZ6WSmux3NUVfoJ2KhdVBzViIi\nEj+uWqystSHgIeAZ4BTwhLX2hDHmQWPMg7EOKPGlY2jC0fmqeTUFPgzQ3K9iJSIi8WNR98tba58G\nnr7ksUfe5NjfXHosiUejkzOMToWozHP2MiCAN81NRV6m5qxERCSuaOV1WbT2IecH1xeqLfTRNjhB\nKHzFeyZERESWjYqVLFrH8AQuA+UOD67PqyvyE4rYiwP1IiIiTlOxkkVrH5qkJDuDdE98/LVZMb9Q\nqLa3ERGROBEfPyEl7llraR+ajJvLgAB+r4eSbC9NmrMSEZE4oWIlizI0McPkTNjRrWwup7bQT8vA\nBOGIdTqKiIiIipUsTvvQ7LIG8bDUwkK1RX6CoQinukadjiIiIqJiJYvTMTSJ22UozfE6HeUNagtn\ni97+5kGHk4iIiKhYySK1D09SnpuBxxVff2XyfOnk+9J4rUnFSkREnBdfPyUlLkXs7JIG8TS4vlBt\noZ/9zYNYqzkrERFxloqVXFVfIMh0KEJVXnzNV82rLfLTPzbNBS27ICIiDlOxkquaX4Az3u4InLdi\nbs7qYPOQw0lERCTVqVjJVbUPTZLudlGcHV+D6/OKsrzkZqZxqE3FSkREnKViJVfVMTRBRV4mLmOc\njnJZLmPYUp3HodZhp6OIiEiKU7GSK5oORegamYrbwfV5W2vyONMTYCwYcjqKiIikMBUruaKzPQFC\nERu381XzttbkYy0cbdNZKxERcY6KlVzR6e4AABW58V2stlTlAXBIxUpERBykYiVXdLprFI/LUJiV\n7nSUK8r1pdFQ7OdQqwbYRUTEOSpWckWnuwOU5mTE7eD6Qltr8jnUOqyFQkVExDEqVnJFp7tHKcvJ\ncDrGotxYk8/A+DStgxNORxERkRSlYiVvqi8QpH9smrLcxChWW2vm5qy07IKIiDhExUre1Jm5wfVE\nKVarS7Pxpbs1ZyUiIo5RsZI3dbp7FIDSBLkU6HYZNlfl6c5AERFxjIqVvKlTXQFKsr1keT1OR1m0\nrTV5nOwcZWom7HQUERFJQSpW8qZOd4+ytjzH6RjXZGtNPqGI5XjHiNNRREQkBalYyWWFwhHO9Yyx\nrizb6SjXRAPsIiLiJBUruaym/nGmwxHWJFixKsryUlPg43UNsIuIiANUrOSy5reyWVuWWJcCYfas\nlc5YiYiIE1Ss5LJOd89uZdNQ4nc6yjXbWp1H9+gUXSOTTkcREZEUo2Ill3W6K0BDcRZej9vpKNds\na00+oDkrERFZfipWclmnuwOsLU+s+ap568pzSPe4tFCoiIgsOxUr+RUjkzN0DE8m5HwVQLrHxQ2V\nuTpjJSIiy07FSn7F/FY2iXrGCmbnrI51jDAdijgdRUREUoiKlfyK+a1s1ibYUgsLba3JJxiKcKpr\n1OkoIiKSQhJnrxJZNqe7A+RmplGWIHsEAjy2r/UNnw9PTAPw1RcvsKuh6A1fu39nzbLlEhGR1KIz\nVvIrTneNsrYsG2OM01GuW25mGjkZHtqHtOSCiIgsHxUreYNIxHKmO8C6BNsj8FLGGKryfbQOTjgd\nRUREUoiKlbxB+9Ak49PhhJ6vmldT4GNwfJrxYMjpKCIikiJUrOQNTs0Prif4GSuAqoJMANqGdNZK\nRESWh4qVvMHprgDGwOrSLKejLFlVng8DtA1qzkpERJaHipW8wenuUWoL/fjSE/+G0XSPi7LcDNp1\nxkpERJaJipW8wZnuAGtKE3++al51vo+2oQki1jodRUREUoCKlVw0OR2maWA8oVdcv1R1QSZTMxH6\nx4JORxERkRSgYiUXne0JYC0Ju0fg5VTl+wDNWYmIyPJQsZKLLu4RmARLLcwrzvbi9bh0Z6CIiCwL\nFSu56ExPgIw0F9UFPqejRI3LmNk5Ky0UKiIiy0DFSi462xNgZUkWblfibmVzOdUFmfSMTjEdijgd\nRUREkpyKlVx0rmeM1SXJcxlwXnW+j4iFjmHNWYmISGypWAkAI5MzdI9OsSqJllqYV1UwP8Cuy4Ei\nIhJbib8KpFyXx/a1vuHzloFxALpGJn/la4kuy+uhwJ+uAXYREYk5nbESAHpGZ9d5Ks3OcDhJbFTn\nZ+qMlYiIxJyKlQDQE5gi3e0i15fmdJSYqC7wMToVYmRyxukoIiKSxFSsBICe0SlKcry4THLdETiv\nOl9zViIiEnuLKlbGmLuNMWeMMY3GmM9e5uv3GmOOGmMOG2MOGGNujX5UiaXe0WDSXgYEKM/NwO0y\nmrMSEZGYuurwujHGDTwM3Am0A/uNMXustScXHPYssMdaa40xm4AngLWxCCzRNx4MMRYMUZLjdTpK\nzHjcLipyM3TGSkREYmoxZ6x2AI3W2gvW2mngceDehQdYa8estXbuUz9gkYTRE5gCoDQnec9Yweyc\nVcfwJKGwFgoVEZHYWEyxqgTaFnzePvfYGxhjPmCMOQ38CPit6MST5dA7f0dgsherfB8zYcvpuT0R\nRUREoi1qw+vW2v9rrV0LvB/4i8sdY4x5YG4G60BfX1+0XlqWqGd0Cq/HRU5Gci9rNr8H4uG2YYeT\niIhIslpMseoAqhd8XjX32GVZa18A6o0xRZf52qPW2m3W2m3FxcXXHFZiozcQpDQnA5OkdwTOy/el\n4fd6eL1lyOkoIiKSpBZTrPYDq4wxdcaYdOAjwJ6FBxhjVpq5n8rGmBsBLzAQ7bASfdZaekanKE3i\nwfV5xhhWFPh4vVXFSkREYuOqxcpaGwIeAp4BTgFPWGtPGGMeNMY8OHfYvwGOG2MOM3sH4YcXDLNL\nHBsLhpiYDlOSxEstLFRT4KN5YIL+saDTUUREJAktaqjGWvs08PQljz2y4OO/BP4yutFkOfQGUmNw\nfd6Kwtk5q4MtQ9y1oczhNCIikmy08nqK6xmdXWohmdewWqgiL5N0t0tzViIiEhMqVimuZzRIZpqb\nbG9y3xE4L83tYmNlDgdVrEREJAZUrFJc79zgerLfEbjQTSvyOdoxQjAUdjqKiIgkGRWrFGatpScw\nRUmKzFfNu2lFPtOhCMc7Rp2OIiIiSUbFKoWNToWYmolQmp0a81XzblyRD6A5KxERiToVqxTWO5oa\newReqiQ7g5oCn+asREQk6lSsUljP3FILqXYpEGYvBx5oGULLrYmISDSpWKWwntEp/OluslLkjsCF\nblyRT/9YkLbBSaejiIhIElGxSmGzdwSm3tkqgG1zc1YHWwcdTiIiIslExSpFzd4RGEzJy4AAq0uz\nyfJ6NGclIiJRpWKVooYnZ5gORVJi8+XLcbsMW2vyONCsYiUiItGjYpWi5u8ITJXNly/nxpp8zvQE\nCEzNOB1FRESShIpViuoZnd98OTXPWAFsq83HWjjcNux0FBERSRIqVimqZ3SK7AwPvvTUuyNw3pbq\nPIxBlwNFRCRqVKxSVG8gmLJ3BM7LzkhjTWk2r7eqWImISHSoWKWgSMTSG5hKua1sLuemFfkcah0m\nHNFCoSIisnQqVimobWiCmbBN+TNWMDtnNRYMcbYn4HQUERFJAipWKehM92yJULGCm2oKADig9axE\nRCQKVKxS0PzZmRJdCqS6IJOiLC8Hm7UCu4iILJ2KVQo62zNGni8Nb5rb6SiOM8aws66AVy8MakNm\nERFZMhWrFHS2J0BpCi8MeqmbGwrpHp2iZWDC6SgiIpLgVKxSzEw4woW+8ZReGPRSt9QXAvDKhQGH\nk4iISKJTsUoxLQPjTIcjGlxfoKHYT3G2l1dVrEREZIlUrFLM2Z4xAEpUrC4yxnBzfSGvnB/QnJWI\niCyJilWKOdMdwBjdEXipW+oL6Q0EudA/7nQUERFJYCpWKeZsT4AVBT7S3PpXv9DN9bPrWelyoIiI\nLIV+uqaYsz0BVpdmOx0j7tQV+SnN8fLKeRUrERG5fipWKSQYCtM8MMGaMhWrS83PWWk9KxERWQoV\nqxRyoW+ccMSySmesLuuW+kL6x4Kc7xtzOoqIiCQoFasUMr+VzRoVq8u6pWF+PSttbyMiItfH43QA\nWT5newJ4XIa6Ij8HU3jT4cf2tV72cWstuZlpfGd/G25juH9nzTInExGRRKdilULOdI9RV+Qn3aMT\nlZdjjKG+yM/ZnoDmrEQkob3ZL5CX0i+Q0adilULO9gS4oSrX6Rhxrb7Yz6G2YXoDQaejiIjE3GIL\nGKiELZZOXaSIiekQbUMTrC7RfNWV1BVlAXBBA+wiInIddMYqRTT2jmEtrCnLcjpKXCvwp5PnS9MK\n7CKybHTZLrnojFWKmN8jUEstXF19URZN/eNEIpqzEhGRa6MzVinibE+AdI+LFQU+p6PEvfoiP6+3\nDnG2N8Dashyn44hICrDW0jcWpLl/goHxIOPBEGPzf6ZCGGN4+lgXDcV+GkqyaCjOYl15DgX+dKej\nyyVUrFLE2Z4ADcVZeLRH4FXVFfsBeOX8gIqViMRMY2+Al8718+TBdpoGJhgPhgDwuAxZXg9+r4ds\nbxrlOZmEraV5YJz9zYMEQ5GLx71lVRFvW12iu73jiIpVijjbHWBHXYHTMRJCvi+dfF8ar5wf4BO7\n65yOIyIJ6M3mpoYnpjnaPsLhtmG6R6cAyPOlsboki9oiP3WFfgqz0jHGXPaft9YSmArRGwjyeusQ\nz5/p4/XWYe7eUMamqtw3/edk+ahYpYDA1AydI1Os1h6Bi7ayJJuXzw8wHYroN0ERWZLpUISj7cMc\nahumuX8cC1TnZ/LeTeWsK88h37f4y3nGGHIy08jJTGNlSRY76wp46mgn3znQxqtNA7xvUwUVeZmx\n+z8jV6VilQLmB9e11MLirSnNZn/zIAeaB9m1ssjpOCKSgIYmptl3YYD9zUNMzoQpykrn9nUlbKnK\nozDLG5XXWFHo59O3reRgyxA/PdHNV37RyCdvraeuyB+V55drp2KVAi7uEagzVovWUOIn3e3iudO9\nKlYismjWWvY1DfLtfS2c7BzFGFhfnsMtDUXUFvpicqnOZQzbawvYUJHDI788z+OvtfK7t68kJyMt\n6q8lV6drHCngbE+AzDQ3lTo9vGhej5ud9QU8d6bX6SgikgBC4QhPHenk3of38pFHX+VC3zhvXV3M\nH71zDffvXEFdkT/m80++dA8f27mCqVCYf36tlbCWjHGEzlilgNNdAdaUZeNyaajxWty+toTPP3WS\nloFxVhTqtLqI/KqJ6RBP7G/j63ubaBucpK7Iz3/7wEZmQtaR+czSnAw+uLWK7xxo4yfHu3jPpopl\nz5DqdMYqyVlrOd09yrpyXQa8VrevLQHgudM6ayUib9Q1Mslf/eQ0u77wHJ976iQl2Rn8n4/fxLN/\n+DY+tnOFoze9bK7O45aGQvaeH+Bo+7BjOVKVzlglud5AkKGJGa3HdB1WFPppKPbz3OleLbsgIgAc\nah3iG3ub+fGxLiLWcuf6Uh54az03rYiv5WzetbGMjqFJvvd6B6U5GZTmZDgdKWWoWCW5k12jAKzV\n4Pp1uX1tCf/4cgvjwRB+r/5zEUlFk9Nhnj7Wxbf2tXCodZhsr4ff3FXLv91VS3Wc7mbhcbm4f0cN\nf/t8I9/e18pDb1+ppWOWiX5SJLnTXbN3BOqM1fV5+9oSvvpiE3sb+3nnhjKn44jIMrHWcrxjlMf3\nt7LncCeBYIi6Ij/v3VTOTTX5eNPcvHiu3+mYV5STmcZ926r5xt4m9jcPslt3OC8LFaskd7p7lMq8\nTHJ9uu32emyvLSDb6+H5M70qViIpoGtkkh8d7eJ7r3dwsmsUr8fFe24o58Pbq9lRV8A/v9bmdMRr\nsrIkixWFPl4+388tDYW4tDJ7zKlYJbnTXQFdBlyCNLeLt6wu4vnTfVhrtV2ESBLqCwT58fEufnik\ni9eaBwG4oTKXv3j/Ru7ZXEFuZmL/Yrq7oYjHXmvlZOcoGytznY6T9FSsklgwFOZ83xjvWF/idJSE\n9vY1JTx9rJuTXaNsqNA3JZFk0D40wc9O9vDMiW72XRjEAiXZXt6xrpRNVbkUza2M/qOjXc4GjYL1\nFTnk+9LY29ivYrUMVKyS2PnecUIRq/mqJbptzWwxff50r4qVSIKaXXomwE9P9PDTk92c6Jy9sWd1\naRa3rSnmhqo8ypL0zjmXMexeWcQPj3bRNjgRtwP3yULFKomd7p79xqE1rJamONvL5qpcnjvdy0O3\nr3I6jogsUigc4UDLED87OVum2gYnMQZurMnnT9+9ljvXl1FX5Oexfa1OR425m2ry+dnJHl5q7Oej\nO2qcjpPUFlWsjDF3A18C3MDXrLVfuOTrHwP+BDBAAPiUtfZIlLPKNTrVNUq6x0WtVg1fsrevLeFL\nz55jcHyaAv/id6IXkei7UhEKRSKc7x3jeMcoF/rHGJqYId3tYvfKQj5920ruWFdCSXZynpm6Em+a\nmx21Bew938/wxDR5Pn0fi5WrFitjjBt4GLgTaAf2G2P2WGtPLjisCXibtXbIGPMu4FFgZywCy+Kd\n7g6wujQLj1trlyzV7WtL+OLPz/GLM7188MYqp+OIyAKhcITG3jGOdYxwqnuUqZkIXo+LtWXZ3L0x\nl9UlWXjT3FgLPz+ZujspzK7G3s8r5wd41w3lTsdJWos5Y7UDaLTWXgAwxjwO3AtcLFbW2pcXHP8q\noJ88ceBUV4Db1hQ7HSMpbKyYHWZ97rSKlUg8iFhL88A4h1uHOd45wtRMhIw0F+vLc9hYmcvKYv1S\neak8XzobKnJ5rXmQ29eW4E1zOx0pKS2mWFUCCxfuaOfKZ6M+Cfx4KaFk6foCQfrHglpqIUpcLsPt\na4v58fFupmbCZOgbkogjzvUEeOZEN0fahhmenL3Mt6Eih01VeTSU+PG4VKau5NaVRRzrGOFAy5AW\nDI2RqA6vG2PezmyxuvVNvv4A8ABATY2G52LpTPfsiuvry3VHYLS8d1MFTxxo5xdnerl7o06jiyyX\n8WCIHx7t5PH9bRxqHcZlZhe+vGtDGevKc7RVyzWoLvCxokALhsbSYopVB1C94POqucfewBizCfga\n8C5r7cDlnsha+yiz81ds27bNXnNaWbT5OwLX6IxV1OxqKKQoy8v3D3WqWInEmLWWw23DfGd/G08d\n6WR8OszKkiz+83vWEY5YsjMSe9FOJ+1eObtg6JnuAOv0y3fULaZY7QdWGWPqmC1UHwHuX3iAMaYG\n+B7wcWvt2ainlGt2qitASbaXwrlF7mTpPG4X79tczrdfbWVkcibhV2MWiUeT02F+cLiDf3q1hROd\no2SmuXnf5nI+vL2GG2vyMMakxPIIsbSuPIeMNBcnOkdVrGLgqsXKWhsyxjwEPMPscgvfsNaeMMY8\nOPf1R4D/AhQCX5nb8iNkrd0Wu9jyZua/4bx8vp/czDR9A4qy92+p5O/3NvOT4118eLsuZ4tES1P/\nON96tYXvHmhjdCrEmtJs/uv7N3LvlgqdnYoyt8uwujSbM92jRKzV5cAoW9SMlbX2aeDpSx57ZMHH\nvw38dnSjyfUKRyy9gSArS7KcjpJ0NlXlUlfk5/uHOlWsRBbhSr/cWWtp7Bvj5cYBzvQE8LgMd28s\n4zduqWV7bb725oyhdWU5HG0foW1wghVa6zCqtPJ6EuofCxKO2KTdnsFJxhju3VLBl549R9fIJOW5\nmU5HEkk406EIR9qG2Xu+n95AkCyvhzvWlrC9roCcjDQae8do7B1zOmZSW12ajcvMrneoYhVdupUi\nCXWPTAHAqhi6AAAgAElEQVRQlqtiFQvv31KJtbDncKfTUUQSSmBqhp+e7OYvf3Ka/3u4A4/L8KGb\nqvjju9Zwx7pScnTJb9lkprupLfRzqmvU6ShJR2esklD36BRuYyjO1uB6LNQW+dlcncf3D3fy797W\n4HQckbjXMzrFS439HG4bJhKxrCvPYffKImoLfbrc56C15Tk8faxLW3VFmYpVEuoemaI426uF8mLo\n/Vsq+PxTJznbE2B1qZa0ELmUtZZXzg/wDy83cbZnjDS3YduKfHavLKJIdyvHhXVl2Tx9rItTXaNa\nLDSK9JM3CXWPTukyYIy9d1MFbpfh+4d+ZUk3kZQ2E47wg8MdvPdvXuL+r+2jY3iKd6wr5Y/vWsu9\nWypVquJIYZaX4mzvxXUPJTp0xirJTEyHGJmc0eB6jBVne9m9sogfHO7kj965BpdLlzMktY0FQzz+\nWit/v7eZjuFJGor9fOGDNxAMRUjTnn1xa11ZNi819murrihSsUoyGlyPnqutAVaa7eWFs3184cen\n+dP3rFumVCLxpW1wgn94uZkn9rcRCIbYUVfA/3fvBt6+pgSXS4t5xrt15Tm8cK6fsz0BNlXlOR0n\nKahYJZnuURWr5bK+Ioe0w4bD7cNORxFZVt9+tYXmgQn2NvZzqmsUY2BjZS67G4qoLvDRMxrk8f1t\nTseURagu8OFLd3O6W8UqWlSskkz3yBS+dDfZXv2rjTWvx8268hyOtY/oNLqkhLFgiO8f6uDh5xvp\nGpkiM83N21YXs7O+UFs8JSiXMawty+ZUV4BwxOLWWMOS6advkukcmaQ8N0O3MC+T7bUFHG0f4QeH\nO7QSuyStY+0jPPZaCz843MnEdJjy3Azev6WSLdV5pHs0P5Xo1pbl8HrrMC2D49QXaceOpVKxSiLB\nUJiekSC7VxY6HSVl1Bf5Kc/N4GsvNnHftmoVWkka/WNBnjrSyfde7+BYxwgZaS7et6mC+3fWcLJz\nVH/Xk8iqkizcLsPproCKVRSoWCWRs91jhK2lIk/brCwXYwy3riziuwfb+cXZPt6+psTpSCLXbWom\nzM9O9vB/D3Xwy7N9hCOW9eU5fO596/nAjVUXL/ed6go4nFSiyZvmpr5odhX2d99Q7nSchKdilUSO\ndYwAUKlitaxuqMrlhXN9fP3FJhUrSTjBUJiXzvXzo6Nd/OxkD4FgiPLcDH7nLfV8YGsla8q0AG4q\nWFuew1NHOukLBLVrxxKpWCWR452zp+u1NcHy8rhc/NtdtfzVT85wsnOU9RU5TkcSuaLpUIS9jf38\n8GgXPz3ZTWAqRG5mGqvLstlSnUddkR+XMRxsGeJgy5DTcWUZrCvL5qkjcLp7lOLsYqfjJDQVqyRy\nvGOEirxMzT444GM7VvC3zzXytZcu8L/v2+J0HJFfMROeLVM/OtrFMye6GZ0KkZPh4e4NZbxnUzm7\nVxbx3QPtTscUh+T50in0p9PUP85bVqlYLYWKVZKYDkU43RVgZ32B01FSUq4vjfu2VfPtfS38yd1r\nKdXK9xIHwhHLf/vRKY62D3Oic5TJmTAZaS7Wl+ewsTKXlSVZeFwuOoenVKqEuiI/JzpHiViLS7+g\nXzcVqyRxtifAdDiiwXUH/dbuOr75SjP/+HIzf3z3WqfjSAo71TXK9w918IPDnXSPTuH1zJapG+bL\nlLaYkcuoLfJzoGWI3tGgFpleAhWrJHFcg+uOqyn0cdeGMr69r5WHbl+JL13/ecnyGRqf5l9eb+fJ\ng+2c7g7gcRluW1PC29eWsLYsW/v1yVXVFfoBaBoYV7FaAn3nTxLHOkbI9no0uO6w335LHT8+3s13\nD7Tzb3fVOh1Hkthj+1qx1tI0MM7+pkGOd44Sjliq8zO5Z3MFN1Tm4tcODHIN8nxp5Gam0dw/zi31\nWg/xeum/uiRxvGOEDZU5ui7usJtWFLC1Jo+vv9TEr9+8QttDSEwEpmbY29jPvqZB+seCZKS52F5b\nwI7aAp1pkOtmjKGuyM/53jGstboR6jrp3HASmAlHONUd4IbKXKejCPDv3lpP6+AETxzQJrQSXU39\n43xuzwlu/u/P8qNjXWSmufjQjVV89u513LO5QqVKlqy20E8gGGJgfNrpKAlLZ6ySwLmeMaZDETZW\n5jIeDDsdJ+XdtaGM7bX5/M9nzvDuG8q1Oa0sibWWvY0DfGNvE8+f6cXjMrx3UwXluRlU5fucjidJ\nprZo9u9Uc/84RVlaKPR6qFglgfnB9Rsqc3n1wqDDacQYw+fu2cD7/uYlvvjzs/z5+zY4HUkSxGP7\nWi9+HIpEONY+wkuN/XSNTOH3enj7mhJ21BWQk6GyLrFRnOXFn+6mqX+cbbVavud6qFglgWMdI2R5\nPdQW+lWs4sSGilw+uqOGb77Swkd31LC6VNuCyOJMTofZ3zzIy+f7GZ0KUZLt5YNbK9lSnadlEiTm\njDHUFvlpHhh3OkrCUrFKAsc6RthQkYNLg9Jx5d+/cw0/PNrF5586wbc+uVODoHJFbYMT/PBoJwda\nhpgORWgo9vPBG6tYVZKlvzuyrOYXCh2emCbPpzvNr5WKVYILhSOc6hrl129e4XQUuUSBP50/vHM1\nf77nBM+c6Obujdo1Xn7V661DfO3FC/zkeDcAm6vy2L2ySIv9imNq59azah4YZ4uK1TVTsUpwjX1j\nBEMR3REYpz62s4bH9rXyFz88xW1rSshIczsdSeJAKBzhmRM9fGNvEwdbhsjJ8PDAWxvIzUzTzQ7i\nuLLcDLweF039E2ypznc6TsJRsUpwx9pnB9c3qljFJY/bxZ/fs577v7qPR1+4wO/fscrpSOKg4Ylp\nHt/fxjdfbqZzZIqaAh9//r713LetGr/X84bhdRGnuIyhttBPc7/mrK6HilWCO94xgj/dTX2R3+ko\n8iZ2NRTxnhvK+covGnn/lkpqCnWLfKr565+d5ZULAxxqHWImbKkv8vPrO1ewtjwblzH84HCn0xFF\n3qC2yM+ZngBjwRBZWsH/mujdSnCzg+u5GlyPc3/6nnW8cK6PTz92kCcf3KVLgikgFI7w81O9fPOV\nZl4+P4DHZdhclceulYWU52p+SuJbXeG/rmelKyLXRsUqgYXCEU52jXL/Dg2ux7vKvEy++OEtfPIf\nD/Bn3z/OX31ok+70SlIDY0Ee39/Gt19toXNkisq8TO5aX8q22gLt3ScJoyI/kzS3oXlAxepa6b/y\nBHa+b5ypmQg3VOU4HSXlLXY25vfvWMWXnz3H1pp87t9ZE+NUslystRxqG+afXmnhR0e7mA5H2L2y\nkD+/ZwN3rC3hiQPtTkcUuSYel4vqAp/mrK6DilUCO7ZgxXVJDJ+5YxVH2ob53J4TrK/IYUt1ntOR\n5Do9tq+V6VCEo+3DvNo0QOfwFF6PixtX5LOzroDSnAwGxqZVqiRh1RX6ee50L1MzYY0vXAMVqwR2\nvGMEX7qbuqIsp6PIIrldhi99ZAvv/ZuX+NS3DvLD37uVQu3HlXBaBsZ5+lgXB1uGmJwJU5Lt5Z7N\nFWytzsOrH0CSJGqL/Fhm/76vKdOVkcVSsUpgxztGWF+eg1uD6wklz5fOI79+E//m717m9/75EN/8\nrR3aqiQBRCKWX57t45uvNPOLs30YYH1FLjfXF1BX6NfMnCSd6nwfbmNo6lexuhYqVglqOhThWMcI\nH9eK6wlpY2Uu//X9G/kPTx7lj588yl99aJPKVZwanpjmuwfa+adXW2gdnKA428vv376KzDQ3OVrM\nU5JYusdFRV4GLYMTTkdJKCpWCepE5wjBUIRttVoVN1H92rZqekan+F8/Pcv4dIgvf3QrXo8uIzlt\n/kaEzuFJXr0wwJH2YWbCltpCHx/ZXs36ihw8LpVgSQ0rCv28emGAUCTidJSEoWKVoA62DAFwY42K\nVSJ76PZV+L0ePv/USX77Hw/w6Me3kZmucuWUqZkwh1qH2Nc0SOvgBGluw5bqfG6uL9DaU5KSagp8\nvNTYT9fwlNNREoaKVYI62DJEdUEmJTkZTkeRJfrE7jr86R4++72j/MY39vH139xOToYuMS2n1oEJ\nvv1aC0/sb2NoYoairHTefUM5N9Xkq+hKSqspmF0oVJcDF0/FKgFZaznQMsStK4ucjiJRct/2anxe\nN3/w+GHu/+qr/OMnduhuwRgLhsL8/GQv3znQxovn+nAZw53rSqnIy6ShWMPoIgA5mWnk+9JoHdB6\nVoulYpWA2ocm6QsEuXGFLgMmk/duqsCf7uHBbx3kXV96kf/1a5t56+pip2MlndPdo3xnfxvfP9TB\n0MQMFbkZ/P7tq/jojhrKcjO0EbLIJWoKfDT1j2Ot1S8ci6BilYAOtAwCsE3FKuEs5of277ylnu8c\naOM3vvEan9hdy5/cvVaL8y1R5/AkPzraxZ4jnRzrGCHNbXjn+jLu217NrSuLtGSJyBXUFPo50j5C\nx/AkVfnaRP5qVKwS0MGWIbK8HlaXZjsdRWKgIi+Th96+kqb+cf5+bzN7G/v54oe3sr5C68hci/6x\nID8+Nlum9jfP3uyxqSqXP3vvej6wtZICf7rDCUUSw4q5OauDLUMqVougYpWADjQPsbUmT79lJ7E0\nt4vP3bOB29YU8x+ePMr7H97Lp25r4HfeWk+WNvK9LGstJzpHef50L8+d6eVw2zDWwqqSLN6xrpRN\nVbkUzc2t/eR4t8NpRRJHaU4G6W4Xr7cMce+WSqfjxD19h04wgakZzvQEuHtjmdNRZBnctqaEZ/7g\nrfzZD47zpWfP8a1XW/j021fysZ01ujzI7CW+15oGeeX8AL8420vPaBCAzVW5fOaOVbxrYzlryrI1\nNyWyBG6Xoaogk4OtQ05HSQgqVgnmUOvsb+E3ab4q6S0sA7sbiqjJ9/HTk938xQ9P8uVnz3HH2hK2\n1uTz8VtSY/X9UDjC+b5xjrQNs69pkNeaB2gbnAQgI81FQ3EWt64sZnVpFtlzy1UcbBm6uOabiFy/\nFQU+XjjXz3gwhF9nza9I706COdgyhMvAluo8p6PIMqsu8PHJW+tp7B3jpye7+d6hDp4708vAeJBf\n21ZNZV7yLGA5NROmsXeM090BjneMcKxjhBOdI0zNzK7+nO9LY0ddAZsq86gr8lOWm4FLdyuJxExN\ngZ9wpI8j7cPsatBSP1eiYpVgDrYMsaYs5+Jv5JJ6VpZk0VDcwKmuAK9eGOCLPz/Hl549x60ri/jw\n9mruXF8al1vjXO5y3OR0mP6xIP1jQfrGgvSOBukZnWJoYpqInT0mM83NxsocPrqjhk1VudxQmUt9\nURYul9ElPpFlMr9Q6OstQypWV6FilUDCEcuh1iE+eGOV01HEYcYY1lfksL4ih7esKuK7B9t58kAb\nDz12CH+6m5vrC9m9sohbVxWxqiTLsbVnrLX0jQVpG5zgUOsQgxPTDI1PMzA2Tf9YkPHp8MVjXQYK\n/V7KcjPYXJ1HaU4GpdleirK9F89GTU5HeK1piNeadHlPZDllprtZVZKlS+uLoGKVQM50BxifDmu+\nSt6gusDHH965ms/csYqXGvv52clu9jYO8OzpXgCKs73sqCtgdUk29cV+GoqzqCvyL3mrllA4wsD4\nNH2BIL2BKfoCQbpGpugcnqRzeIqO4Uk6hycJhv5181bD/ErO6ayvyKEoy0tRlpfCrHQK/Ona3Fgk\njt20Ip8fH+8mErG4dFf6m1KxSiAH5xYGVbGSy3G7DG9bXczb5lZrbx+a4OXGAfae7+dgyxBPH+vC\n2n89vjw3gzxfOtkZHnIy0sjJ9JA9N5QatpaIhUjEEo5YxqdDjE6GGJ2aYXRyhtGpEEMT0294PgBj\noCTbS0VeJusrcrhzfSmVeZnUFPo42jZCni+NNLfKk0giunFFPo/vb+NC/xgrS7SO4ptRsUoA83Mk\n//J6B9kZHl4426dtBeSqqvJ93Lfdx33bq4HZgfCm/nHO941xoW+cloEJRiZnCEzN0DE8yamuGcaC\nIYwBlzFzf2YLm9/rISfDQ4E/nRWFfnIyPBRmeSnO9nKma5SsjDSyvR6yMzx4LlOcuoanKM7W3oci\niWz+l/qDLUMqVlewqGJljLkb+BLgBr5mrf3CJV9fC/w9cCPwn6y1/yvaQQVaBsapKfCpVMkbXM8A\n9/wluGhYX5EblecRkfhWX+Qnz5fGwZYhPry9xuk4ceuq5+SNMW7gYeBdwHrgo8aY9ZccNgj8PqBC\nFSOjkzMMTcywotDvdBQREUlBxhhuqsnXAPtVLGbYYQfQaK29YK2dBh4H7l14gLW211q7H5iJQUYB\nWgYngH/ds0lERGS53bgin/N94wyNTzsdJW4tplhVAm0LPm+fe0yWUevAOB6XoTwvw+koIiKSoubn\nrA616azVm1nW23OMMQ8YYw4YYw709fUt50snvJbBCaryM3U7uoiIOGZzVR5ul9HlwCtYzE/pDqB6\nwedVc49dM2vto9babdbabcXFxdfzFCkpOBOmc3hS81UiIuKozHQ3Gypy2N+sYvVmFlOs9gOrjDF1\nxph04CPAntjGkoUu9I8TsbNbmYiIiDhpR20Bh9uGmZoJX/3gFHTVYmWtDQEPAc8Ap4AnrLUnjDEP\nGmMeBDDGlBlj2oE/BP6zMabdGJMTy+Cp5FxvgDS30eC6iIg4bmd9IdOhCIfbhp2OEpcWtY6VtfZp\n4OlLHntkwcfdzF4ilBg41zNGfVHWZRdeFBERWU47agswBvZdGOTm+kKn48Qd/aSOc22DEwyMT7Oq\nVJcBRUTEebm+NNaV5fDqhQGno8QlFas498K52bsnNV8lIiLx4ub6Ql5vHSIY0pzVpVSs4tyLZ/vJ\nzUyjOErbj4iIiCzVzvoCgqEIR9tHnI4Sd1Ss4lgoHGHv+X5WlWRpf0AREYkbO2oLAHj1vC4HXkrF\nKo4daR8hMBViVal2ERcRkfiR709nbVk2+5oGnY4Sd1Ss4tgLZ/twGWgo1sKgIiISX26uL+RAyyDT\noYjTUeKKilUce/FcH5uq8vClL2pVDBERkWVzc30BUzMRjnVoPauFVKzi1MjkDIfbhnnrqiKno4iI\niPyKHXWza1i9ekGXAxdSsYpTr5zvJ2LhLau1p6KIiMSfAn86q0uztJ7VJVSs4tQvz/aT7fWwpTrP\n6SgiIiKXdXN9IQdbhpgJa85qnopVHLLW8sLZPm5pKCRN29iIiEic2llXyMR0mGMdWs9qnn5qx6Hm\ngQk6hid1GVBEROLajrrZ9az2ac7qIhWrOPTi3DY2b1ulYiUiIvGrONvLypIs9jVpzmqeilUceuFs\nHysKfdQU+pyOIiIickU76wrY3zRISHNWgIpV3JkORXjl/ABv0TILIiKSAG6uL2R8OsyJzlGno8QF\nFas488LZPsanw9yxttTpKCIiIle1s35uzkqXAwEVq7iz50gn+b40btUZKxERSQAl2RnUF/m1UOgc\nFas4Mh4M8bOTPbxnU7mWWRARkYSxa2Uhr5wfYGom7HQUx+mndxz5+akeJmfC3LO50ukoIiIii/bO\n9WVMzoR56Vy/01Ecp2IVR/Yc7qQiN4NtK/KdjiIiIrJoN9cXkp3h4ZkT3U5HcZyKVZwYGp/ml2f7\neN+WClwu43QcERGRRUv3uLh9bQk/P9WT8ssuqFjFiaePdxGKWO7ZXOF0FBERkWt214YyhiZm2N88\n5HQUR6lYxYk9hztZWZLF+vIcp6OIiIhcs7etLibd4+KnJ1P7cqCKVRzoHJ7kteZB7t1cgTG6DCgi\nIonH7/Xw1lVF/PRED9Zap+M4RsUqDvzwaCfWwvt0GVBERBLYO9eX0TE8mdKrsKtYxYE9RzrZXJ1H\nbZHf6SgiIiLX7Y51JbgMKX13oIqVwxp7xzjeMaqhdRERSXiFWV621xbw0xM9TkdxjIqVw/Yc6cQY\neN+mcqejiIiILNk7N5RxpidAc/+401EcoWLlIGstTx3pZFdDISU5GU7HERERWbJ3ri8FUvdyoIqV\ngw60DNHUP67LgCIikjSqC3xsqMhRsZLl9/DzjRT403U3oIiIJJW7NpRxqG2Y3tEpp6MsOxUrhxzv\nGOEXZ/r45K11+NI9TscRERGJmnduKMVa+Nmp1BtiV7FyyMPPN5Kd4eHjt6xwOoqIiEhUrSnNZkWh\nj2dS8O5AFSsHNPYG+MmJbn5zVy05GWlOxxEREYkqYwx3byzj5cZ+OocnnY6zrFSsHPCV58+T4XHz\nid11TkcRERGJiY/fPHtF5msvNjmcZHmpWC2z1oEJfnCkk4/trKHAn+50HBERkZioyvdxz5YKHt/f\nytD4tNNxlo2mppfZIy+cx20Mv/PWeh7b1+p0HBERkZh58G0NfO/1Dr75Sgufeccqp+MsC52xWkbd\nI1M8eaCdX9tWRakWBBURkSS3ujSbd6wr4R9ebmJiOuR0nGWhYrWMvvriBcLW8uDbGpyOIiIisiw+\ndVsDQxMzfGd/m9NRloWK1TIZGAvy2L5W7t1SQXWBz+k4IiIiy+KmFQXsqC3gqy9cYCYccTpOzKlY\nLQNrLX/2g+PMhCN8+raVTscRERFZVp+6rYHOkSn2HO50OkrMqVgtgycOtPH0sW7+6K41rCzJcjqO\niIjIsrptTTFry7J55JfniUSs03FiSsUqxhp7x/jcnpPsXlnIA2+pdzqOiIjIsjPG8KnbGjjXO8az\np3udjhNTKlYxFAyF+czjh8hIc/G/79uCy2WcjiQiIuKI99xQTlV+Jl/5RSPWJu9ZKxWrGPqfPznD\nic5R/ueHNmt5BRERSWket4tP3dbAodZh/n5vs9NxYkbFKkZ+caaXr73UxG/csoJ3rC91Oo6IiIjj\nPrq9hrs2lPLfnj7FS+f6nY4TEypWMdAbmOKPvnuUNaXZ/Om71zkdR0REJC64XIb//74trCzO4ncf\ne52WgXGnI0WdilWUNfYG+NDfvUJgaoYvf3QrGWlupyOJiIjEjSyvh0d/4yYAfuebBxgLJteK7CpW\nUfTLs3184OGXmZgO888P3MyasmynI4mIiMSdFYV+Hr7/Rhp7x/j3TxxOqiUYVKyiwFrLP+xt4hN/\n/xpVBT5+8NBubqzJdzqWiIhI3Lp1VRH/6T3reeZED3/zXKPTcaLG43SARDcTjvD5p07wrVdbece6\nUr70kS34vXpbRURErua3dtdysnOUv/75WbpHp/iTu9eQ50t3OtaSqAFcp6mZME8ebOcbLzVxoX+c\nB9/WwB/ftUZrVYmIiCySMYb//sGN5PnS+IeXm3nmRDf/8V1r+dBNVRiTmD9PF3Up0BhztzHmjDGm\n0Rjz2ct83Rhjvjz39aPGmBujHzU+9I8F+eufnWXXF57jP3//ONkZHh79+E189l1rVapERESukdfj\n5s/eu56nHrqV2kIf/+HJo9z3f17hTHfA6WjX5apnrIwxbuBh4E6gHdhvjNljrT254LB3Aavm/uwE\n/m7ufxPexHSIY+0jHGob5lDrEM+f6WM6FOEd60r5nbfUsaOuIGFbtYiISLxYX5HDkw/u4smD7fyP\nH5/i3V9+kRtr8tjVUMSuhkK21uST7on/0fDFXArcATRaay8AGGMeB+4FFhare4Fv2tk16l81xuQZ\nY8qttV1RT3ydIhFL2FrCEct0OMLUTJjgTIRgKMzUTIShiWn6x4L0BYL0j03TFwhypjvAmZ4A4bm7\nFWoKfNy3rYpP7K6joVibKYuIiESTy2W4b3s1d64v5esvNfHCuT6+/Nw5vvTsOTLT3GyrzaehOIvi\nbO+//snyUpmXSb4/PmazFlOsKoG2BZ+386tnoy53TCXgWLE62TnK+x/ee7FMXQuvx0VRlpe6Ij+f\nvq2BrTV5bK7KozDLG6O0IiIiMi/fn84f3bWGP7prDSMTM+xrGuDl8wO8emGAw23DBKbeuPbVR3fU\n8D8+eINDad9oWYfXjTEPAA/MfTpmjDmznK9/Lc4CLwPfjt1LFAHJuZ6/s/S+xobe19jQ+xobel9j\n4GNOB7iCL8z9ibEVizloMcWqA6he8HnV3GPXegzW2keBRxcTLNkZYw5Ya7c5nSPZ6H2NDb2vsaH3\nNTb0voqTFjMFth9YZYypM8akAx8B9lxyzB7gN+buDrwZGImn+SoRERGR5XDVM1bW2pAx5iHgGcAN\nfMNae8IY8+Dc1x8BngbeDTQCE8AnYhdZREREJD4tasbKWvs0s+Vp4WOPLPjYAr8b3WhJT5dEY0Pv\na2zofY0Nva+xofdVHGNmO5GIiIiILFX8r7QlIiIikiBUrBxwtS2C5NoZY75hjOk1xhx3OksyMcZU\nG2OeN8acNMacMMZ8xulMycAYk2GMec0Yc2Tuff2805mShTHGbYw5ZIz5odNZJDWpWC2zBVsEvQtY\nD3zUGLPe2VRJ4R+Au50OkYRCwL+31q4HbgZ+V39foyII3G6t3QxsAe6eu6Nalu4zwCmnQ0jqUrFa\nfhe3CLLWTgPzWwTJElhrXwAGnc6RbKy1Xdba1+c+DjD7A6vS2VSJz84am/s0be6PBl6XyBhTBbwH\n+JrTWSR1qVgtvzfb/kckrhljaoGtwD5nkySHuUtWh4Fe4GfWWr2vS/dF4I+BiNNBJHWpWInIVRlj\nsoB/Af7AWjvqdJ5kYK0NW2u3MLtTxQ5jzEanMyUyY8x7gV5r7UGns0hqU7Fafova/kckXhhj0pgt\nVd+21n7P6TzJxlo7DDyPZgSXajdwjzGmmdkRi9uNMd9yNpKkIhWr5beYLYJE4oIxxgBfB05Za/+3\n03mShTGm2BiTN/dxJnAncNrZVInNWvsfrbVV1tpaZr+vPmet/XWHY0kKUrFaZtbaEDC/RdAp4Alr\n7QlnUyU+Y8w/A68Aa4wx7caYTzqdKUnsBj7O7G//h+f+vNvpUEmgHHjeGHOU2V+2fmat1fIAIklA\nK6+LiIiIRInOWImIiIhEiYqViIiISJSoWMn/a+9uQ7Os4jiOf39NKbM1EBKColVoimFGCplkRs9P\nUCANKntRYUYlib4oCnoXRiRkIkW9KMJCSCmofABBipW6siwko9I30sNCxZ62svnrxXUWd0t3b3Db\n9NtIuAoAAAPRSURBVOb3gcHFOde5z//aYPz3v87OiYiIiAZJYhURERHRIEmsIiIiIhokiVVERERE\ngySxioi6JPWVPax2StohabGkk0rfdEnLBxnbLumO/y/a/8zdU87kOy5I6pD0jaTsWxXRhJJYRcRQ\n9NieZnsK1S7hNwBPAtj+2PbCQca2AyOSWBXfljP5hkxSy7EKxvZq4L5j9fkRMbKSWEXEsNjuBuYD\nD6kyp7/6IumKmh3aP5XUCiwFLi9ti0oV6QNJ28vXZWXsHEmbJb0paZekVeVIHSTNkPRhqZZtk9Qq\nqUXSM5K6JH0u6f6hxC/pLUmflOrb/Jr2XyU9K2kHMPMoc04p15+VOSeUsXfVtL/Yn5hJur484w5J\nmxr4Y4iI49SokQ4gIk48tneX5GH8gK4lwIO2OyWdBvQCjwJLbN8MIOlU4BrbvSUxeQOYXsZfDEwB\nvgM6gVmStgGrgQ7bXZJOB3qAe4GDtmdIOhnolLTR9p464d9je385o69L0hrb+4CxwFbbi8s5nruO\nMOcC4Dnbq8o9LZImAx3ALNuHJK0E7pS0DngJmG17j6Rxw/5GR8QJJ4lVRDRSJ7BM0ipgre29pehU\nazSwQtI0oA+YWNO3zfZegLIuqh04CHxvuwvA9s+l/1pgqqS5ZWwbMAGol1gtlHRbuT67jNlXYllT\n2i84ypwfAY9LOqs839eSrgIuoUrSAMYA3cClwPv9iZ7t/XXiiogmkMQqIoZN0nlUiUg3MLm/3fZS\nSe8CN1JVkK47wvBFwI/ARVTLEXpr+v6oue5j8N9RAh62vWEYcc8BrgZm2v5d0mbglNLda7tvsPG2\nX5e0FbgJeK+8fhTwqu3HBsx1y1DjiojmkTVWETEsks4AXgBWeMAp7pLOt/2F7aeBLmAS8AvQWnNb\nG1U16DAwD6i3UPwr4ExJM8ocrZJGARuABySNLu0TJY2t81ltwIGSVE2iqioNec6SUO62vRx4G5gK\nbALmShpf7h0n6RxgCzBb0rn97XVii4gmkIpVRAzFmPJqbjTwF/AasOwI9z0i6UrgMLATWFeu+8qi\n8FeAlcAaSXcD64HfBpvY9p+SOoDny7qoHqqq08tUrwq3l0XuPwG31nmO9cACSV9SJU9bhjnn7cA8\nSYeAH4CnynqtJ4CNqragOES1zmxLWRy/trR3U/1HZUQ0MQ34gzMiomlIagfesX3hCIfyL+WV5D8L\n+iOieeRVYEQ0sz6gTcfZBqFUVbsDIx1LRDReKlYRERERDZKKVURERESDJLGKiIiIaJAkVhEREREN\nksQqIiIiokGSWEVEREQ0yN/oSZQVzOjtdAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f32f94ea710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nb_merge_dist_plot(\n",
    "    SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n",
    "    SkyCoord(specz['ra'] * u.deg, specz['dec'] * u.deg)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue = specz_merge(master_catalogue, specz, radius=1.5 * u.arcsec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VII.a Wavelength domain coverage\n",
    "\n",
    "We add a binary `flag_optnir_obs` indicating that a source was observed in a given wavelength domain:\n",
    "\n",
    "- 1 for observation in optical;\n",
    "- 2 for observation in near-infrared;\n",
    "- 4 for observation in mid-infrared (IRAC).\n",
    "\n",
    "It's an integer binary flag, so a source observed both in optical and near-infrared by not in mid-infrared would have this flag at 1 + 2 = 3.\n",
    "\n",
    "*Note 1: The observation flag is based on the creation of multi-order coverage maps from the catalogues, this may not be accurate, especially on the edges of the coverage.*\n",
    "\n",
    "*Note 2: Being on the observation coverage does not mean having fluxes in that wavelength domain. For sources observed in one domain but having no flux in it, one must take into consideration de different depths in the catalogue we are using.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "atlas_moc = MOC(filename=\"../../dmu0/dmu0_ATLAS/data/ATLAS_HATLAS-SGP_MOC.fits\")\n",
    "kids_moc = MOC(filename=\"../../dmu0/dmu0_KIDS/data/KIDS-DR3_HATLAS-SGP_MOC.fits\")\n",
    "ps1_moc = MOC(filename=\"../../dmu0/dmu0_PanSTARRS1-3SS/data/PanSTARRS1-3SS_HATLAS-SGP_MOC.fits\")\n",
    "viking_moc = MOC(filename=\"../../dmu0/dmu0_VISTA-VIKING/data/VIKING_SGP_MOC.fits\")\n",
    "des_moc = MOC(filename=\"../../dmu0/dmu0_DES/data/DES-DR1_HATLAS-SGP_MOC.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "was_observed_optical = inMoc(\n",
    "    master_catalogue['ra'], master_catalogue['dec'],\n",
    "    atlas_moc + kids_moc + ps1_moc + des_moc) \n",
    "\n",
    "was_observed_nir = inMoc(\n",
    "    master_catalogue['ra'], master_catalogue['dec'],\n",
    "    viking_moc\n",
    ")\n",
    "\n",
    "was_observed_mir = np.zeros(len(master_catalogue), dtype=bool)\n",
    "\n",
    "#was_observed_mir = inMoc(\n",
    "#    master_catalogue['ra'], master_catalogue['dec'],   \n",
    "#)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue.add_column(\n",
    "    Column(\n",
    "        1 * was_observed_optical + 2 * was_observed_nir + 4 * was_observed_mir,\n",
    "        name=\"flag_optnir_obs\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VII.b Wavelength domain detection\n",
    "\n",
    "We add a binary `flag_optnir_det` indicating that a source was detected in a given wavelength domain:\n",
    "\n",
    "- 1 for detection in optical;\n",
    "- 2 for detection in near-infrared;\n",
    "- 4 for detection in mid-infrared (IRAC).\n",
    "\n",
    "It's an integer binary flag, so a source detected both in optical and near-infrared by not in mid-infrared would have this flag at 1 + 2 = 3.\n",
    "\n",
    "*Note 1: We use the total flux columns to know if the source has flux, in some catalogues, we may have aperture flux and no total flux.*\n",
    "\n",
    "To get rid of artefacts (chip edges, star flares, etc.) we consider that a source is detected in one wavelength domain when it has a flux value in **at least two bands**. That means that good sources will be excluded from this flag when they are on the coverage of only one band.\n",
    "\n",
    "This now takes place in the fonal stage when photometry is folded in\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VIII - Cross-identification table\n",
    "\n",
    "We are producing a table associating to each HELP identifier, the identifiers of the sources in the pristine catalogue. This can be used to easily get additional information from them.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['omegacam_intid', 'atlas_id', 'kids_id', 'ps1_id', 'viking_id', 'des_id', 'help_id', 'specz_id']\n"
     ]
    }
   ],
   "source": [
    "id_names = []\n",
    "for col in master_catalogue.colnames:\n",
    "    if '_id' in col:\n",
    "        id_names += [col]\n",
    "    if '_intid' in col:\n",
    "        id_names += [col]\n",
    "        \n",
    "print(id_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue[id_names].write(\n",
    "    \"{}/master_list_cross_ident_sgp{}.fits\".format(OUT_DIR, SUFFIX), overwrite=True)\n",
    "id_names.remove('help_id')\n",
    "\n",
    "master_catalogue.remove_columns(['atlas_id', 'kids_id'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IX - Adding HEALPix index\n",
    "\n",
    "We are adding a column with a HEALPix index at order 13 associated with each source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue.add_column(Column(\n",
    "    data=coords_to_hpidx(master_catalogue['ra'], master_catalogue['dec'], order=13),\n",
    "    name=\"hp_idx\"\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## X - Saving the catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "columns = [\"help_id\", \"field\", \"ra\", \"dec\", \"hp_idx\"]\n",
    "\n",
    "bands = [column[5:] for column in master_catalogue.colnames if 'f_ap' in column]\n",
    "for band in bands:\n",
    "    columns += [\"f_ap_{}\".format(band), \"ferr_ap_{}\".format(band),\n",
    "                \"m_ap_{}\".format(band), \"merr_ap_{}\".format(band),\n",
    "                \"f_{}\".format(band), \"ferr_{}\".format(band),\n",
    "                \"m_{}\".format(band), \"merr_{}\".format(band),\n",
    "                \"flag_{}\".format(band)]    \n",
    "    \n",
    "columns += [\"stellarity\", \"flag_cleaned\", \"flag_merged\", \"flag_gaia\", \"flag_optnir_obs\", \"flag_optnir_det\", \"ebv\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Missing columns: {'specz_id', 'omegacam_intid', 'zspec_association_flag', 'viking_id', 'omegacam_flag_merged', 'zspec_qual', 'zspec', 'stellarity_origin', 'ps1_id', 'des_id'}\n"
     ]
    }
   ],
   "source": [
    "# We check for columns in the master catalogue that we will not save to disk.\n",
    "print(\"Missing columns: {}\".format(set(master_catalogue.colnames) - set(columns)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#master_catalogue[columns].write(\"{}/master_catalogue_hatlas-sgp{}.fits\".format(OUT_DIR, SUFFIX))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## XI - Folding in photometry\n",
    "On HATLAS-SGP there is too much data to load all in to memory at once so we perform the cross matching without photometry columns. Only now do we fold in the photometry data by first cutting the catalogue up in to manageable sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "298\n"
     ]
    }
   ],
   "source": [
    "split_length = 100000 #number of sources to include in every sub catalogue\n",
    "num_files = int(np.ceil(len(master_catalogue)/split_length))\n",
    "print(num_files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "surveys = [\n",
    "    #['atlas',    \"ATLAS.fits\" ,     Table.read(\"{}/ATLAS.fits\".format(TMP_DIR)   ), \"atlas_id\"],    \n",
    "    #['kids',    \"KIDS.fits\" ,   Table.read(\"{}/KIDS.fits\".format(TMP_DIR) ), \"kids_id\"],     \n",
    "    ['omegacam', \"omegacam_merged_hatlas-sgp.fits\" ,      Table.read(\"{}/omegacam_merged_hatlas-sgp.fits\".format(TMP_DIR)    ), \"omegacam_intid\"],       \n",
    "    ['ps1',    \"PS1.fits\" ,         Table.read(\"{}/PS1.fits\".format(TMP_DIR)       ), \"ps1_id\"],\n",
    "    ['viking',    \"VISTA-VIKING.fits\" ,  Table.read(\"{}/VISTA-VIKING.fits\".format(TMP_DIR)), \"viking_id\"],\n",
    "    ['des',    \"DES.fits\" ,  Table.read(\"{}/DES.fits\".format(TMP_DIR)), \"des_id\"]\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Sort catalogue by HELP id so that it is split up in RA strips\n",
    "master_catalogue.sort('help_id')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n=0\n",
    "for sub_file in range(num_files):\n",
    "    # the following used to have a -1 which was wrong as it left out objects\n",
    "    sub_catalogue = master_catalogue[n*split_length:(n+1)*split_length] \n",
    "    #print(n)\n",
    "    for survey in surveys:\n",
    "        #print(survey[0])\n",
    "        sub_catalogue = join(sub_catalogue, \n",
    "                               survey[2], #Table.read(\"{}/{}\".format(TMP_DIR, survey[1])),\n",
    "                               join_type='left',\n",
    "                               metadata_conflicts='silent',\n",
    "                               keys=survey[3]\n",
    "                            )\n",
    "    #print('Finished join')\n",
    "    #sub_catalogue.remove_columns(id_names)\n",
    "    \n",
    "    #We merge all the flag_merge columns\n",
    "    sub_catalogue['flag_merged'].name = 'flag_merged_tmp'\n",
    "    flag_merged_columns = [column for column in sub_catalogue.colnames\n",
    "                            if 'flag_merged' in column]\n",
    "    \n",
    "    flag_merged_column = np.zeros(len(sub_catalogue), dtype=bool)\n",
    "    for column in flag_merged_columns:\n",
    "        flag_merged_column |= sub_catalogue[column]\n",
    "    \n",
    "    sub_catalogue.add_column(Column(data=flag_merged_column, name=\"flag_merged\"))\n",
    "    sub_catalogue.remove_columns(flag_merged_columns)\n",
    "    sub_catalogue['flag_merged'].name = 'flag_merged_tmp'\n",
    "    flag_merged_columns = [column for column in sub_catalogue.colnames\n",
    "                            if 'flag_merged' in column]\n",
    "    \n",
    "    flag_merged_column = np.zeros(len(sub_catalogue), dtype=bool)\n",
    "    for column in flag_merged_columns:\n",
    "        flag_merged_column |= sub_catalogue[column]\n",
    "    \n",
    "    sub_catalogue.add_column(Column(data=flag_merged_column, name=\"flag_merged\"))\n",
    "    sub_catalogue.remove_columns(flag_merged_columns)\n",
    "    \n",
    "    \n",
    "    for col in sub_catalogue.colnames:\n",
    "        if (col.startswith(\"m_\") or col.startswith(\"merr_\") or col.startswith(\"f_\") or col.startswith(\"ferr_\")):\n",
    "            sub_catalogue[col] = sub_catalogue[col].astype(float)\n",
    "            sub_catalogue[col].fill_value = np.nan\n",
    "        if \"stellarity\" in col:\n",
    "            sub_catalogue[col].fill_value = np.nan\n",
    "        elif \"flag\" in col:\n",
    "            sub_catalogue[col].fill_value = 0\n",
    "        elif \"id\" in col:\n",
    "            sub_catalogue[col].fill_value = -1\n",
    "        #Remove residual ra decs from join\n",
    "        if col.endswith('_ra'):\n",
    "            sub_catalogue.remove_column(col)\n",
    "        if col.endswith('_dec'):\n",
    "            sub_catalogue.remove_column(col)\n",
    "        \n",
    "    sub_catalogue = sub_catalogue.filled()\n",
    "    \n",
    "    \n",
    "    #Adding detection flag\n",
    "    \n",
    "    # SpARCS is a catalogue of sources detected in r (with fluxes measured at \n",
    "    # this prior position in the other bands).  Thus, we are only using the r\n",
    "    # CFHT band.\n",
    "    # Check to use catalogue flags from HSC and PanSTARRS.\n",
    "    nb_optical_flux = (\n",
    "        #Omegacam \n",
    "        #1 * ~np.isnan(sub_catalogue['f_omegacam_ul']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_omegacam_u']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_omegacam_g']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_omegacam_r']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_omegacam_i']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_omegacam_z']) +\n",
    "        # PanSTARRS\n",
    "        1 * ~np.isnan(sub_catalogue['f_gpc1_g']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_gpc1_r']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_gpc1_i']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_gpc1_z']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_gpc1_y']) +\n",
    "        #DES\n",
    "        1 * ~np.isnan(sub_catalogue['f_decam_g']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_decam_r']) +    \n",
    "        1 * ~np.isnan(sub_catalogue['f_decam_i']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_decam_z']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_decam_y']) \n",
    "    )\n",
    "\n",
    "    nb_nir_flux = (\n",
    "        1 * ~np.isnan(sub_catalogue['f_vista_z']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_vista_y']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_vista_j']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_vista_h']) +\n",
    "        1 * ~np.isnan(sub_catalogue['f_vista_ks'])\n",
    "    )\n",
    "\n",
    "    nb_mir_flux = np.zeros(len(sub_catalogue), dtype=bool)\n",
    "\n",
    "\n",
    "    has_optical_flux = nb_optical_flux >= 2\n",
    "    has_nir_flux = nb_nir_flux >= 2\n",
    "    has_mir_flux = nb_mir_flux >= 2\n",
    "\n",
    "    sub_catalogue.add_column(\n",
    "    Column(\n",
    "        1 * has_optical_flux + 2 * has_nir_flux + 4 * has_mir_flux,\n",
    "        name=\"flag_optnir_det\")\n",
    "        )\n",
    "    \n",
    "    \n",
    "    sub_catalogue.remove_columns(id_names)\n",
    "    \n",
    "    sub_catalogue.write(\"{}/tiles/sub_catalogue_sgp{}_{}.fits\".format(OUT_DIR, SUFFIX, n), overwrite=True)\n",
    "    n += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to generate final catalogue¶\n",
    "After this notebook has been run there will be a set of sub catalogues in data/tiles/\n",
    "\n",
    "These need to be stacked using stilts:\n",
    "<pre>\n",
    "suffix=20180111\n",
    "\n",
    "ls ./data/tiles/sub_catalogue_sgp_$suffix_*.fits > ./data/tiles/fits_list_$suffix.lis\n",
    "\n",
    "stilts tcat in=@./data/tiles/fits_list_$suffix.lis out=./data/master_catalogue_sgp_$suffix.fits\n",
    "</pre>\n",
    "For many purposes this file may be too large. In order to run checks and diagnostics we typically take a subset using something like:\n",
    "<pre>\n",
    "stilts tpipe cmd='every 10' ./data/master_catalogue_sgp_$suffix.fits omode=out out=./data/master_catalogue_sgp_RANDOM10PCSAMPLE_$suffix.fits\n",
    "</pre>"
   ]
  }
 ],
 "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
}
