{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pylab\n",
    "import pymoc\n",
    "import xidplus\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from astropy.table import Table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook uses all the raw data from the XID+MIPS catalogue, maps, PSF and relevant MOCs to create XID+ prior object and relevant tiling scheme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in MOCs\n",
    "The selection functions required are the main MOC associated with the masterlist. As the prior for XID+ is based on IRAC detected sources coming from two different surveys at different depths (SERVS and SWIRE) I will split the XID+ run into two different runs. Here we use the SERVS depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Sel_func=pymoc.MOC()\n",
    "Sel_func.read('../data/ELAIS_N1/MOCs/holes_ELAIS-N1_irac1_O16_MOC.fits')\n",
    "SERVS_MOC=pymoc.MOC()\n",
    "SERVS_MOC.read('../data/ELAIS_N1/MOCs/DF-SERVS_ELAIS-N1_MOC.fits')\n",
    "Final=Sel_func.intersection(SERVS_MOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in XID+MIPS catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: UnitsWarning: 'degrees' did not parse as fits unit: At col 0, Unit 'degrees' not supported by the FITS standard.  [astropy.units.core]\n",
      "WARNING: UnitsWarning: 'muJy' did not parse as fits unit: At col 0, Unit 'muJy' not supported by the FITS standard. Did you mean MJy, mJy or uJy? [astropy.units.core]\n"
     ]
    }
   ],
   "source": [
    "XID_MIPS=Table.read('../data/ELAIS_N1/MIPS/dmu26_XID+MIPS_ELAIS-N1_SERVS_cat_20170710.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<Table length=10>\n",
       "<table id=\"table4565526232\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>help_id</th><th>RA</th><th>Dec</th><th>F_MIPS_24</th><th>FErr_MIPS_24_u</th><th>FErr_MIPS_24_l</th><th>Bkg_MIPS_24</th><th>Sig_conf_MIPS_24</th><th>Rhat_MIPS_24</th><th>n_eff_MIPS_24</th><th>Pval_res_24</th></tr></thead>\n",
       "<thead><tr><th></th><th>degrees</th><th>degrees</th><th>muJy</th><th>muJy</th><th>muJy</th><th>MJy / sr</th><th>MJy / sr</th><th></th><th></th><th></th></tr></thead>\n",
       "<thead><tr><th>str100</th><th>float64</th><th>float64</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th></tr></thead>\n",
       "<tr><td>HELP_J161006.621+533656.872</td><td>242.527587012</td><td>53.6157979058</td><td>10.4667</td><td>22.6099</td><td>3.04599</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161005.559+533642.824</td><td>242.523161793</td><td>53.6118954771</td><td>12.3774</td><td>25.8053</td><td>3.60663</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161008.641+533615.951</td><td>242.536005913</td><td>53.6044309555</td><td>116.806</td><td>131.279</td><td>103.169</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161002.817+533625.213</td><td>242.511735883</td><td>53.6070036615</td><td>70.9787</td><td>88.4681</td><td>52.4919</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161009.708+533610.215</td><td>242.540451147</td><td>53.6028373869</td><td>132.994</td><td>146.192</td><td>120.556</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161000.783+533454.820</td><td>242.503263565</td><td>53.5818943665</td><td>5.79204</td><td>13.5911</td><td>1.59307</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161002.558+533606.162</td><td>242.510656314</td><td>53.6017117556</td><td>18.1655</td><td>31.734</td><td>6.73982</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161005.420+533611.962</td><td>242.522583128</td><td>53.6033228215</td><td>12.266</td><td>24.3795</td><td>3.96191</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J160959.123+533542.063</td><td>242.496347238</td><td>53.5950173915</td><td>17.7358</td><td>31.0833</td><td>7.03387</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161000.625+533444.663</td><td>242.50260411</td><td>53.5790731835</td><td>52.7274</td><td>68.7187</td><td>37.6159</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=10>\n",
       "                                              help_id                                                ...\n",
       "                                                                                                     ...\n",
       "                                               str100                                                ...\n",
       "---------------------------------------------------------------------------------------------------- ...\n",
       "HELP_J161006.621+533656.872                                                                          ...\n",
       "HELP_J161005.559+533642.824                                                                          ...\n",
       "HELP_J161008.641+533615.951                                                                          ...\n",
       "HELP_J161002.817+533625.213                                                                          ...\n",
       "HELP_J161009.708+533610.215                                                                          ...\n",
       "HELP_J161000.783+533454.820                                                                          ...\n",
       "HELP_J161002.558+533606.162                                                                          ...\n",
       "HELP_J161005.420+533611.962                                                                          ...\n",
       "HELP_J160959.123+533542.063                                                                          ...\n",
       "HELP_J161000.625+533444.663                                                                          ..."
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XID_MIPS[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGoCAYAAAAAZTE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXm8bFlZ3/1de6iqM93bl+Y2TSND8zYupI0NAi9GJCIJ\nGBDzgagYsOkXyauoDNGoCSC80KgMBiOYjwoyRZBWMUAEEqII2gkGwyQ4EHaAFgembprbfc9UVXtY\n7x9r76o6555Tp6pO1aldp37fz6f7rqq99q7nV7tOPfU861lrGeccQgghxEkTzNsAIYQQy4kckBBC\niLkgBySEEGIuyAEJIYSYC3JAQggh5kI0bwMqbrtt89jleOfOrXLhws40zFlIpF/6pb9++s+f3zDz\ntqGunKoIKIrCeZswV6Rf+peZZde/iJwqBySEEGJxqE0KTgiAP/7EFw499sgH3uMELRFCzBpFQEII\nIeaCIiAxF4ZFOkKI5UAOSCwMhzktpeaEWEyUghNCCDEX5ICEEELMBaXgxEzRWI8Q4jAUAQkhhJgL\nckBCCCHmghyQEEKIuaAxIDEV5jnWU732xnqLza1273mVZwtRbxQBCSGEmAtyQEIIIeaCHJAQQoi5\nIAckhBBiLqgIQYyFJpYKIaaFHNAJMe4Xtyq4hBCnHTmgmjLPjdkU5QghTgI5oClzEl/eh73G9z36\n/jN/7UVC2zcIUW/kgJYYRTpCiHkiB3SK+G8f+vyelQCEEKLOyAFNgCIHIYQ4PpoHJIQQYi7IAQkh\nhJgLckBCCCHmghyQEEKIuaAiBLF0zHOSrxCijxzQEFTtJoQQs0MpOCGEEHNBDkgIIcRckAMSQggx\nFzQGJMQAWsBUiJNDEZAQQoi5IAckhBBiLsgBCSGEmAsaA0LzfYQQYh4oAhJCCDEXFAEJMQKqjhNi\n+igCEkIIMRfkgIQQQswFOSAhhBBzQWNAQhwDjQ0JMTlL5YBUbi2EEPVBKTghhBBzQQ5ICCHEXFiq\nFJwQJ4XGhoQ4GkVAQggh5sKpjIBUbCDqyrDPpqIjsWwoAhJCCDEXTmUEJMQionEjsWwsrAM66I91\nY73F5lZ7DtYIMTtGTSlXn385LLEo1N4BaTxHiPGY5t+MnJmYJcY5N28bhBBCLCEqQhBCCDEX5ICE\nEELMBTkgIYQQc0EOSAghxFyQAxJCCDEX5ICEEELMBTkgIYQQc2FmE1GttU8DnlY+bAEPBK5MkuSO\nWb2mEEKIxeFEJqJaa38F+GSSJL8+8xcTQgixEMw8BWetfQhwrZyPEEKIQU5iLbjnAzce1ck554wx\nJ2COEEKcKCN/sWVZ7qIonKUt8+BQ/TN1QNbaywCbJMkfjdL/uOlAY8yxr7HISL/0S3/99I/zw/rC\nhZ0ZWjIfzp/fOPTYrFNw/wh4/4xfQwghxAIyawdkgVtG6Tj4K6EO7UntnfS8urWPYtk1L7v+cfrW\nuX0U45w36WssM7XZjsFNwZC6huAnhfRLv/TXT78Zwxvddttm/QQck/PnNw7Vr4moQggh5oIckBBC\niLlQKwc0Tv72qHxrHfLLs8xZL6P+aeTjF7E97Phh70Md7D6Je74o+sXBaAzoFCH90i/99dOvMaDD\nx4BOYiKqEEIcyic+8Qle+tKXEoYhD3/4w3nWs5615/jm5ib/+l//a3Z2dmg0GvzCL/wC58+f54Mf\n/CC/+Iu/yMrKCo94xCP4sR/7sbnaKcanVik4IcTy8eIXv5hXvvKV3HTTTfz5n/85n/rUp/Ycf+c7\n38nXf/3X89a3vpXHPvaxvOENb6AoCl74whfyy7/8y9x0003ccsstfPSjH52rnWJ8ahMBDYbP82pP\nw/ZJz1sWzYddYxH0z/IaR7Xf8Y538P73v5/t7W0uXLjAM5/5TB7zmMfw4Q9/mFe96lUEQcC97nUv\nXvKSl9But3nBC17A5uYmt956K095ylN4ylOewvXXX8/ll1/OnXfeyQtf+EJ+5md+hjAMcc7xyle+\nkrvf/e68/OUv5+Mf/zgAj3/847nhhht47nOfS6PR4Itf/CK33norL3vZy7j22mt51KMexdVXX801\n11zD8573vJ7GZzzjGezs9Gf0X3PNNbzoRS86UNvm5ibdbpd73/veOOf4tm/7Nj70oQ/xgAc8oNfn\n67/+67nlllswxrC1tUUcx9xxxx2cOXOGe97znhhj+OZv/mY+/vGP85CHPITrr7+et7zlLXte64Yb\nbuDqq6/mlltuwTnHL/3SL3H+/Pmejb/5m7/J7//+7++5V694xSu4xz3ugXOOra0tut0u97rXvQAu\nsVNMRm0c0OAf87zakzLpNZZR82HXWAT9s7zGKO3d3V3e+MY38rWvfY0nPelJfMd3fAcvfOELuemm\nm7j88st59atfzTve8Q6uvfZaHve4x/GYxzyGr3zlK9xwww08+clPBuC7vuu7ePSjH81b3/pWvumb\nvomf+qmf4mMf+xibm5t8+tOf5gtf+AK/8zu/Q5Zl/MAP/AAPe9jDALjqqqt4yUtewtve9jbe9ra3\nceONN/KlL32Jt7/97Zw7d26Pxte+9rUj69/a2mJ9fb33/NraGn/3d3+3p89ll13Gn/zJn/C4xz2O\nO++8k9/8zd/k3LlztNttbrnlFu5973tz88038w3f8A0453jLW95y4Gs96EEP4sYbb+Smm27ita99\nLS94wQt6x6+//nquv/76Q+2t7KzYb6eYjNo4ICHEcB760IcSBAF3vetdOXPmDLfeeiu33XYbP/7j\nPw5Ap9PhW7/1W/n2b/923vzmN/O+972P9fV10jTtXePqq68G4Hu/93t53etexw/90A+xsbHBT/zE\nT3DLLbfw4Ac/GGMMcRxz3XXX8bnPfQ6g90v/7ne/ey9COnfu3CXOB4ZHQLA32nj5y1/O9vZ279j2\n9jYbG3vXDvuVX/kV/uW//Jf8i3/xL0iShOc85zm8613v4hWveAUvfvGLaTQa3O9+9zvQlkG+5Vu+\nBYAHPehBvP/9e1cIOywCuuqqqwBYX18/0k4xPnJAQiwIf/VXfwXAV7/6Vba2trjyyiu58sor+dVf\n/VU2Njb4wAc+wOrqKm984xt54AMfyJOf/GT+9E//lJtvvrl3jSDww77vf//7echDHsKznvUs3vOe\n9/C6172OxzzmMbzjHe/gaU97Gmma8md/9mc84QlPAA4uKa6utZ9hERBcGm3Ecczf/u3fcs973pMP\nfvCDPPOZz9zT/8yZM70v+7vc5S5sbW0B8MEPfpDXv/71xHHMs5/9bL7ne77nyPfvyiuv5OMf/zjX\nXHPNUJv2s76+fqSdYnzkgIRYEG677Tae9rSnsbm5yYte9CLCMOT5z38+z3jGMyiKgvX1dV7xilcA\n8PM///P8l//yXzhz5gxhGNLtdvdc6xu/8Rt57nOfy6/92q+R5znPe97zuPbaa/nwhz/M93//95Om\nKY997GO59tprZ67rxS9+MT/90z9Nnuc8/OEP57rrrgPg6U9/Oq95zWt4znOewwtf+EJ+67d+izRN\n+dmf/VkArrjiCp70pCfRbDb57u/+bu53v/vhnOOpT31qLw03yDvf+U7e9KY3sbq62nufpmGnmBzN\nAzpFSP/p1f+Od7yDv/7rv+Ynf/InD+1zmvWPQqX/pS99Kc9//vP3HHvqU5/KjTfeyH3ve9952KV5\nQIegMmwhxKniB3/wB+dtghgRRUCnCOmXfumvn35FQIqAhBBC1Aw5oDFY9gUIpVn6D2qfNpb9np8k\nckBjcFKTPOuKNEv/Qe3TxrLf85NEDkgIIcRckAMSQggxF2rjgOqwadQ4ed1pn1eH9jhM47w6tKdh\n+6Tn1aE9Drrnk91zcTgqwz5FSL/0S3/99KsMW2XYQgghaoYckBBCiLlQKwc0Ti73qNxsHXLNs8xf\nL6P+cfLx87Z1mu1hxw97H+pg90nc80XRLw5GY0CnCOmXfumvn36NAWkMSAghRM2QAxJCCDEXauOA\n6pCnnbTOfxrn1aE9DtM4r27tSW2f9Lw6tMdB93yy88ThaAzoFCH90i/99dOvMSCNAQkhhKgZckBC\nCCHmghzQGCxj7f+yz32QZuk/qi0mRw5oDJZxn5BxNS+7/tOA7vnyaZ4XckBCCCHmghyQEEKIuSAH\nJIQQYi7IAQkhhJgL0Swvbq19HvDPgAbwq0mSvGGWryeEEGJxmFkEZK19JPCtwMOBbwfuOavXEkII\nsXjMMgL6TuAvgHcCZ4CfnuFrCSGEWDBm6YDuCtwbeDxwNfAua+39kyQ5sIB+WhO7ln2CmPRL/zKz\n6PrPnVslisJ5m3FizNIB3Q58OkmSLpBYa9vAeeDWgzpPY2JXXRcjPCmkX/qlv376x3GKFy7szNCS\n+XD+/Mahx2ZZBfdB4J9aa4219ipgDe+UhBBCiNk5oCRJ3gP8GfBh4N3AM5MkyWf1ekIIIRaLWu0H\nNBhCH9Ue9bnT1q5YRv2TvCenoT1MXx3sO2nNoxyvU9uMkYNbtv2AauWAjnuNuuaATwrpl37pr59+\nOSBtSCeEEKJmyAEJIYSYC7VxQONuAnWS7XHsHadvndtHMWnfRWmPo2OcvnVuH4Xu+fT6Co/GgE4R\n0i/90l8//RoD0hiQEEKImiEHJIQQYi7IAR3BSeXX68Qyah5kGfUvo+ZBll3/vJADOoLBnPJx2ovE\nMmoeZBn1L6PmQZZd/7yQAxJCCDEX5ICEEELMBTkgIYQQc0EOSAghxFyQAxJCCDEX5ICEEELMBTkg\nIYQQc6GWDmjeCxDOg0XTPMvrLYL+abBomhftnov6U0sHNK1JYYs0cWzRNM/yeougfxosmuZFu+ei\n/tTSAQkhhDj9yAEJIYSYC7VyQOPkdSc9vujt4x5fxPaox0fpu4jtUY+P0ndR2pMcH+W8eesRe9GG\ndKcI6Zd+6a+ffjOGJ9KGdEIIIcQJIAckhBBiLtTGAdUhT6sxmMnGGI46vujt4x5fxPZB6J5Pds/F\n4WgM6BQh/dIv/fXTrzEgjQEJIYSoGXJAQggh5kI0Sidr7cOBfwC8CXhYkiT/faZWjcm0Q+/B682i\nXUdmrXmW92ja15PmeuqfNsuouW4cOQZkrf1XwBOAewD/EPgg8IYkSV45TUM0BnR8pF/6pb9++jUG\ndLwxoKcB3wlsJ0lyO/BQ4OnTMU0IIcSyMooDypMk6Q48bgP5jOwRQgixJIzigG621r4SWLPWPgF4\nF/D+2ZolhBDitDOKA/pp4DPAJ4EbgP8K/NQsjRJCCHH6ObQKzlp7r4GH7y3/q7gK+NtZGSWEEOL0\nM6wM+2bAAYMVDNVjB9x3hnYJIYQ45RzqgJIkufokDRFCCLFcDEvBvThJkhdba9+Ej3j2kCTJiZRi\nL+MEsUXQfFLXrqv+abMImnXPxbQZloL7WPnvH096cWvtx4GL5cO/TpLkB8e9xuCHaVrturMImk/q\n2nXVP20WQbPuuZg2w1Jw7y6bVyVJ8rLBY9balx51YWttCzBJkjzyWBYKIYQ4lQxLwb0cuAL4Z9ba\n++0751uA5x9x7euAVWvtH5TnPD9Jkj89pr1CCCFOCcNScG8HHgD8Y3xFXEUG/OwI194BXgm8Hrgf\n8F5rrU2SJDvshGls4rTsG0FJv/QvM4uu/9y5VaIonLcZJ8Yoi5GeTZLkznEvbK1tAkGSJLvl4w8D\n35Mkyd8d1F+LkR4f6Zd+6a+ffi1GevhipKNsx/AEa+0vAufKxwZwSZIc5aafjt/C4cestVcBZ4Av\njfB6QgghloBRHNCLgEcmSfKXY177DcB/tNZ+EF/G/fRh6TchhBDLxSgO6AsTOB/KFbSfMmr/k9wc\naxpzCyY9vujtg/Qto+Zl17+MmifVLw5nlDGgVwFfB/wBfisGAJIkefM0DdEY0PGRfumX/vrp1xjQ\n8caAzgKb+N1QKxwwVQckhBBiuTgyAqqw1p5LkuTCrAxRBHR8pF/6pb9++hUBHWNLbmvtddbaTwOf\ntNZeZa39rLX2m6dq4TGZ9vyhk2xPw95psAiaB9E9Pz6LoHmQZbznp51RNqT7D8ATgduTJPki8KPA\na2Zq1ZhM41fPtNafGrc9DXunwSJoHkT3/PgsguZBlvGen3ZGcUCrSZL87+pBkiTvA5qzM0kIIcQy\nMIoD+pq19jrKLRmstT8AfG2mVgkhhDj1jFIF96PAbwDXWmvvAD4DXD9Tq4QQQpx6RqqCs9beDdgC\nQuCKJEk+O21DVAV3fKRf+qW/fvpVBXe8KrjnAO9NkmQbvx7cu621PzxF+4QQQiwho4wB/TDwCIAk\nSf4GeDDw7FkaJYQQ4vQzigOKgc7A4y5lQYIQQggxKaMUIfxn4APW2reVj/858HuzM+lwprG44KIh\nzdIvzadf/7IyigN6Pn4i6rcDKfDLSZL855ladQjLOBFMmqV/nPaisuz3fFkZxQF9JEmSbwb+06yN\nEUIIsTyMMgb0FWvtI8ottoUQQoipMIoDeghwM7Brrc2ttYW1Np+FMZMu3jfOefNajHDaixdO+7w6\ntMfhtGieVMc41EHbsn/OxcGMvB3DrNFE1OMj/dIv/fXTr4mox9iQzlrbAH4KsPj5Pz8OvLzcclsI\nIYSYiFFScL8CrOMnoGbANcAbZmmUEEKI088oDujBSZI8H0iTJNkB/h/gQdM2pA552nHaR+V6523f\nLNrHPb6o7eMeX+S27vnxNYvDGcUBuTINV+Um78oMVkKY10ZRk7aPmocwb/tm0T7u8UVtH/f4Ird1\nz4+vWRzOKA7oVcAfAne31r4K+CjwSzO1SgghxKln1O0YHgB8B95h3ZwkyZ9P2xBVwR0f6Zd+6a+f\nflXBHW87hhh4DPBP8U7oYdbaWiU4J8231iHXPG3bp3GNumoeZJE1T8P2SVk0zYMs8j0XBzPKUjyv\nB1aAX8c7rBuAa/Hl2LVg0l89dcg1T9v2aVyjrpoHWWTN07B9UhZN8yCLfM/FwYzigB6WJMn9qwfW\n2ncDfzk7k4QQQiwDoxQh/J219pqBx3cDvjAje4QQQiwJo0RAMfBJa+1/x09E/TbgS9baDwAkSfKo\nGdp3KIMDjoe1TxvLqHmQZdS/jJoHWUbNy8QoDuhF+x6/chaGjMsy5l2XUfMgy6h/GTUPsoyalwkt\nRnqKkH7pl/766VcZ9jHKsIUQQohZIAckhBBiLoyyHcPVwOOB+wEF8Fng3UmS/M2MbRvKMg7OLqPm\nQZZR/zJqHmQZNS8Th44BWWvvjl8H7t7AnwB/A6TA1cC3A58HfjJJkr+fhiEaAzo+0i/90l8//RoD\nmmxDupcDNyZJ8qmDDlprrwNeBjz1eOYJIYRYRlQFd4qQfumX/vrpVwR0vC257ws8A78PUO9CSZI8\nfYRzrwA+Bjw6SZJPH9V/0g/QNM6rQ3sati+aZt1zaZ5Ex6TnLZrm084oE1Hfjt8P6H/A6BvRlato\nvxbYHfWcSW/WNM6rQ3satk96Xh3a47Ds93zS8+rQnpRF1i8OZhQHZJIk+ekJrv1K4DXA8yY4Vwgh\nxClnFAf0P621TwR+L0mSYpSLWmufBtyWJMnvW2tHdkDT2D9j2ffgkH7pX2YWXf+5c6tEUThvM06M\nYWXYBT7lVt3RqqMBXJIkh75L5cKlrvzvgcD/Af5ZkiRfHmKLO26utfrwnURe96hc7zxyzSehf5i+\neefaZ6X/qHt61PHTrH/emmelf8qaVYRwCBNVwVlrm0mSdEbs+8fAjxxVhOAmMWQfyz7oJ/3SL/31\n0y8HdLwtuT+073EAfHQKdgkhhFhiDh0DKvf7eWTZHhz7yYB3jfoCSZI8ckLbhBBCnGKOTMFZa1+d\nJMm/mrUhzjl33Dr/Uc6fd5561Paomk+T/nE4Kl+/KJonvefjHK+Dtml+zsftO++2UnATjgFZa78T\n+F9Jktxhrb0B+L+BjyZJ8h+nbaTGgI6P9Eu/9NdPvxzQBGNA1tpXAc8HWtbanwV+APgr4J9ba395\n6lYKIYRYKoYVITwGeFRZOv14fBn1rwFPBB59EsYJIYQ4vQxzQDvAFWX7K8Ba2V7DFyKcOIOR7Ljt\nRWUZNQ+yjPqXUfMgy65/mRi2EsKNwEestb8NfBq42Vr7h8B3Ar9wEsbtZzC/O257UVlGzYMso/5l\n1DzIsutfJg6NgJIkeTfwCOCLQAP4ELAJPG0WRQhCCCGWC+0HdIqQfumX/vrpVxXcMVZCEEIIIWbB\nsJUQvgjc7YBDRy5GKoQQQhzFsCKEhwIfAJ6YJMmnTsgeIYQQS8KwIoQvAD+Br4YTQgghpoqKEE4R\n0i/90l8//SpCmFIRgrX2h49vzvSZxkS040x+m8fEuUXWvOz6p2HvNK4h/fXXfNoZtwruR2ZixTGZ\nxq+e40x+m8fEuUXWvOz6p2HvNK4h/fXXfNoZ1wHJpQshhJgK4zqg18zECiGEEEvHKBvSxcA/Ae7K\nQASUJMmbp2xLrw7hGBs/VRc61nWGtSsmPT7L9rz1z0PzSeg/6p4edfw065+35lnpP+qeHnV8X1tF\nCIcwbB5Qxe8Cdwf+N1C9OQ6YqgOqQ256nFzupMcXvT3seB3s0z3XPZ9G+7jHD+orLmUUB3T/JEnu\nP3NLhBBCLBWjjAF9zlp7r5lbIoQQYqkYthbcH+FTbVcAf2Gt/SQDG9ElSfKoaRtzVL512n0Xpb2M\n+pdR87LrP67mOusXBzMsBffikzKiYpwc6jT6Lkr7IE67/nF0jNN3UdoHcdr1j6PjKOqgR87naA51\nQEmS3Axgrf0PSZI8e/CYtfY3gJtnbJsQQohTzLAU3OuB+wIPsdZeu++cy2ZtmBBCiNPNsBTczwH3\nAV7N3hWxM3xJ9omw7PnYumo+qWvXVf+0WQTNi3DPxWIxdCKqtTYAHg7cAyiALwIfSZKkM21D3BQ+\nQcv+QZR+6Zf++unXRNQJJqJaa78V+A3g88CXy6fvDtzPWvv0JEneP00jhRBCLBfDUnC/DjwuSZLP\nDD5prb0GeCfwD2ZpmBBCiNPNsImo4X7nU3ILaFVsIYQQx2NYBPQea+27gd8GvlQ+dyXwA8B/nbVh\nQgghTjdHFSF8D/BdwFX4qOcLwH9NkuQ/TdsQFSEcH+mXfumvn34VIRxehHDkdgzQ25LhG4BukiSf\nnqJtPeSAjo/0S7/010+/HNDhDujQMSBr7XvLf78RP+/nN4Dfsdb+RfmcEEIIMTHDihDuVv77KuDZ\nSZI8KEmS64BnAW+YuWVjMMYPjLGvN4t2HZm15lneo2lfT5qlfxptcTSjbMdwNkmS91YPyjXiVmdn\n0vhMO+xexoUJF21hxmW854umedn1i6MZ5oCusdb+GtCx1v4QgLX2nLX2p+hXxQkhhBATMawM+1rg\nocAd+BUQwJdgPwx42mzNEkIIcdo5tArOWttKkqQ97ORhfay1IfA6wAIO+JEkSf5yyOV6hXCTLkZY\n5V+nvajjQdU1kx6fZXtW+g9inOOnRX8d7/ms9B9EHe/5rPVP6Z6rCu4QhqXg3mqt/SFr7cb+A9ba\nDWvtM4HfGnL+dwMkSfJw4AXAzw8z8iTz0dPI6056fBHbBzHO8UVvH/f4IrYPQvd8snsuDmdYCu77\ngB8FPmKtvQP4e/xWDPcBLsdv0/B9h52cJMl/tta+p3x4b3wqTwghhACGpOAGsdZeB9wPvyXD55Ik\n+eSoL1DunvpE4HuTJPmDw/o555xKGIUQp5CRv9iyLHdRFM7Slnkw/koI1to/AT4DvBf4gyRJLkz6\n6tbaK4H/BTwgSZLtg/pUDug4udnDctTzzk1Psz2q/tPUPuo+HnV80dvLqP80adYY0IRL8ZRbLzwW\neAywAvwx8N4kST521Itaa58KfF2SJC+z1p4BPol3QLsH9XdTSJoeNUh62pF+6Zf++umXAzrmWnDg\nK96AR+Id0oOTJPm2I/qvAW/Cr6AdAy9PkuT3DusvB3R8pF/6pb9++uWAJtgRtcJaezXwePpjQJ8B\n/v1R55WptieNbqYQQohlYthipHe31v4OvtT6PsBngQS4F/A2a+3vWGu/bhZGDf5gmEW77iyC/lle\nr66aZ8ki6Nc9F9NmWAT0cuDGJEk+ddDBsjLuZcBTp23UKLX1x2nXnUXQP8vr1VXzLFkE/brnYtqM\nPAY0azQGdHykX/qlv376NQZ0vDGg+wLPAO4K/XruJEmePhXrhBBCLCVHOiDg7cAfAv8DOHXeWQgh\nxHwYxQGZJEl+euaWCCGEWCpG2ZDuf1prn2itHaWvEEIIMRLDluIp8Cm3atyn6mgAlyTJVBcsUhHC\n8ZF+6Zf++ulXEcIERQhJkgybI9Q8rlFCCCGWmyPTatbaD+17HAAfnZlFNWMZJ8Ut+6RAaZb+Sdti\nPA6NgKy1H8Cv/Val4yoy4F2zNas+LOOkuGWfFCjN0j9pW4zHkRNRrbWvTpLkX83aEI0BHR/pl37p\nr59+jQFNOBHVWnt/4KvW2tfgFyL9IiNuxyCEEEIMY1ihwY8Bv413PB8BKqfzemvtT07bkFkvTDjt\nvO40+ta5fRST9l2U9jg6xulb5/ZRnEbNJ3HPxeEMK8P+P8ADkyTZ2ff8KvDxJEnuP01D5p2Ccw4W\n/TNT1xTESSH90l9H/UrBHZ6CG1YFl+I3ktvPSnms9jjnDv1A7n/eAYXzjmjY9YQQQkyHYWNAPw/8\nmbX2/cCXyufuDjwK+JlZG3YU5V7rlzx3WF+4NCy+pH857bZ6uuq/v+JF4bUQQhyfQyOgJEluAr4N\nvwjpDrBbth+RJMlvz8KYcXK2o+RbjTG9/8A7j6JwFM7hymininoqH1O5msJB4Rx54QaOmz2vNRhh\nHVWWOe089aTHF7E96vFR+i5ie9jxYZ/9RW5PcnyU8+atR+xl2DygeyRJ8gXgzdZaC3wnsApMdQme\nQSatuR+WGus7CHA4cNV4T98pOSAM+h+WrPAnVJ8fR+nMcJek6Eb9o5j2PINJjy9ie9Tjo/RdxPaw\n46OkmBexPcnxUc6btx6xl2FjQO8GsNZ+L/A+4GrgG4Cby+dqx0i/OA7ps/+jMvhYv2SEEGL6jLId\nw3Pxabe/AbDWvhT4I+A/zdKwSRj8xVHs+SVStUzvcVFGOAZI08K3jcEEBlc4wGGikHywMiEwGNNP\n3ZVXJzB5owyxAAAgAElEQVT0nq+eM5iFr6oTQohZMsoWC9vA3w08vpVLA4ZjM2oedX+BwMF9KvPM\nnpRZb8zG+AsFxvgxoCPsqa4X7HMo1fWq9N3eqrvxQvFFy01P+xqLoH+W16ir5kEW+Z5Pw96T4I8/\n8YUTfb15M8wBnbXWfgq4G/DvAay1DwZ+D1+MMFVGzaP2xmX29XFQFhf0iwYqh1BFQ85BMVBUkGWF\n90VUxQb9/kXhyIuCIi961+pmBa7wffO8oCjDoDQvyqiJ8nxvV1XscJCOSfTXKTc97Wssgv5ZXqOu\nmgdZ5Hs+DXvF9Bm2HcP/Za29C/BQoFU+fSXw+8BrTsC2kakcTF6AT4l5L1UUBUUBJjBkRQGuqnor\nCENDUUCRO0JjKCjTarl3NABB1k+1BWWKrlu+Zhz6/J3zl6UAjHFQFitkDgyGwLg9TlPjSUII4Rk6\nBpQkydfwDqfiHkmS/MpsTZqQXqQB7P+Od4P/9B2C7+8OjKooo6rqRMfeyrn9L+LKSUSmsgHvkPrP\nqphBCCEGGaUIYZAfAX59FoZMgk+99cd5fDoN8rzwjiUIMMan0gZTYXnhnU41x6coHEVREASGMAwo\nckfhCrqZ7xcFAcZAUeQ9J9bpGsIgIAh90UJYVSLgS7qDwJA7fKQVeNtMv0svSlNUJIRYVsZ1QDX7\npqxcTv+L3DhwBowJMEAQBqRZTuF8tVp/rMen7Iq82LMETxD4Pt2soJMWrDRCiqLAGENW+PSdc9CI\nDHlRkDtDIwr8MQoCY3p2BGVlnCscURiUdlbP+xeU8xFCLCvjOqBajf3sr18bLLf21W4Gh3cqrth/\nbq9rr4AgLwqisi6jilDy3EFoCI1/vaogIS+jHh8Zeefnh5l8AUMYVJGZP1aa07NabkcIcRD7K+Ee\n+cB7zMmS2XOkA7LWxsA/Ae4K7FprbwBIkuTNM7ZtKL15PmVVdc+/GEMAZEVBlrlehRy4cnkdSkfh\nU3VZVtApiw6cc+x0cgLjHVg3LUh92FNGSX2H14hDjDHkeUEUBRhjfOrNQBAY4ijwTs9AGAZEobcr\nDAOCwEc+QeBHh4Kg7/CEEGJZGCUC+l38IqT/m37I4YC5OqDKir0TQPtFBVU00n+e0vp+tUI1IbV3\nsfJ5Exg/AXWgH/v6BYHplWFDNbXI9YoVXFm4UDmz8kHZLivonC9UMIqHhBBLyCgO6P7T3vvnuPQi\nETOYhBucOOqjiiI3FK4oVz0o24XDuer58qwypVYUPmSpig3yvOjNLaoKFqoIZ3u3SxgG5FlBmpky\nyvHjTmmW04gjPxZUQBQGpDiCwtCIQvLcYQoII29/ljvC8lyo17iQMWbPeNU47dOANC+H/mXUXAdG\ncUCfs9beK0mSv525NSNiBsKeKrAJTH/Ojh+nMRh88YArnC82oJy34wq/5I5zZIXzaTZ8UUKW5zTj\nkG6as9vJCENDljuy3NFqBGWfclJqacNKM8R1/TWikNKZdblso0EYBKRBQSMKaMShn4+ET9UVmZ+z\nFIYBLne9FF6dChSWfQFGaV4O/cuouQ4MWw37j/Df7VcAf2Gt/SSQVceTJHnU7M07HoOOaqBgrn+c\nvcUAVXuwYMBfp/+8/3fvB7Fyhf7pfljmM2/7V0MwvTlGrkzFDZYlqEBBCLEsDIuAXnxSRkzC3nB4\nn3Mox3oCYyCovtIDiqKgW5TL9uT9deGCoHIWEIWGTpqSpgVZXpDmlecytDuZ3yMoL5fjKRxZXtDu\nZL4QAZ9+c8BqK+biliOOfGpuc6dLaAJWVyKiwBBGAXHYL15oxCEu9CXc1aoLVKXl7HWIQojl4TRX\nxQ1biudmAGvtf0iS5NmDx6y1vwHcPE1DJs27mt7/9/apvrQDDEX57R0EAUVVmkYVhfSLBipH1E29\nkykGyqgD48O/ouin+vLCkeaOtZWwtzbcbreg2QjJsoIiCimcKecVOVqxI00LwmYIzqf8gsARmYDC\n+Y2WDAel3rwDPE6eetp57Un71rk9qY5l1L+Mmmf1Xi0zw1JwrwfuCzzEWnvtvnMum7Yhk+Rdq5vs\n9ta79Y6x7zmAKDCkeX9SUBgYstz3MXjnEgaGqkt1laoYgXKeT1EUmPK6nW5OWFbFGSBN8zKSynBA\nMw4BRycrCMOALPNOMDYOV0CWF2VVHRA4jAuAausH1/v3OO/XtPPak/atc/sopn1e3dpHseyaZ/Ve\nLTPDUnA/B9wHeDVw48DzGb4kuxb4VBVgnJ8s6gJcUI7NOEea9dNzAHEcEEY+4umkOUUBa62Ydjcj\nL5yvnnNhmfIqeqsZFA46uaOb5rTTgjT1TqfZCMnygk7XUdBf/WBrJy0jp5QoDHrpve3YRzutRsja\nSoOwnDPUaoS0mhFhYAjLooXKiUaRnxxrcIdESEIIsXgMS8F9Hvg8cF31nLX2yUmS/NYoFy4nsL4R\n78SawM8lSfKuY9h6CT46GNj/g711Bm7/g/19C3pzcKq+/Tk8g1UIvp9/vk8Q7O0zeLA3hhOY3kKl\nQW91BD8h1TtJ8Kv0DJxcdro06hk+FqRfXUKIRWJYCu6GA55+SelYRlkJ4Xrg9iRJnlpu6/AJYKgD\nGjVnWzkeB34ZHePr2XyCrH9OEIDJHWEAhatGi0pHYhzNOCDN/RlxHNLp+iK/ZhT0djnF+UKDrHA0\nooAo9M+3GoGfz2P8mnBF4UjLcSC/vYP3VsbkvsDAQafraDZCQhOQZTln15sEQUA3y9lYiclzn/5r\nNv1OrCYwvlChTNH5wnIwZVpuWCRUh/z6tNvj5OPnbes028M0HfY+1MHuk7jni6JfHMywFNyPAPcD\n3k1/KGQD+A4YaSWE36W/bbdhoIT7MEbN2VZfvKb3vzJiGIhmwsA/bsQBWe4wRTlSZMr125wPWWJT\nju04g4sCullBoxESu5CWc2zvphgDYegnlGZ5QRiGhIEhzQq6ac6FzbJPWfnQyXLyAqIAwiCg0y3K\nCac52+2Ms2t+ftDXLnaIooD1lZhtkxGFBa1GhMPRjEOiKCCjIC4dmMMR4IspTK/yz+15Tw76Ujot\n7XHy8fO2dZrtYccPex/qYPdJ3PNF0S8OZpgDegTwIuAbgWckSXKbtfbPkiT5wVEunCTJFoC1dgPv\niF5w1DnjjG1UX7R7n6uuM5DQMgZjqoF80zvOnr7eKVUpNVc+b8q0WeH6dd5m4OTq5YPADERivvy7\nv3FD9YfUt9H/158flOUFzvlxnrwoCF0I+EVOD14jzlEVKcDe4z3nvKTjRMuqu0L6F1v/2mqDIBi2\nUTWcP79xQtbMnmFjQDnw/1lrHw68y1r78+wdAjkSa+09gXcCv5okyU1H9T/uL4bAQFGuhD14KR/B\nGEL8KgZUO5UGgAnKtFnQW62A0rnleUFA5FNwuV9VAWMo8oLtdkZgDFHoCwjS3LG926UofGovKyvl\nssz5OUBh0HNwd253CYL+fJ+VRkS20SAwfj7QajNicwdW4pCNtQZp5ndwjaPQL3LqfMn4oL5qG4ij\nfhmfZpY95SH99dQ/jlPc3uke2ee22zaPY86JM8xhDne1QJIkfwJ8J/Ak/KoII2GtvRvwB8C/TZLk\njUf131NMcJw2e2OCKlLyaSvT6zH4XFC2q2uF5d49gV+2uoyiAqIoJAoMcexn7ASBn0QahkFvz6Aw\nDDCBn3wamIBOVq45d8DUnioaajTCXlQUlxVzzkGzGfX+qIJy6wczILD/x7Y3/TbpezcO0zivDu1p\n2D6Nayya/knPq0N7GraL6WBG/cVgrb0X8I+ADyVJ8rkR+r8a+H7g0wNPPzZJkt2D+rsp/XTxYyW9\na15yPMsLP9/G9+g5gbzw2y5kvRWuHUXhVzZIcz+BtHICReHY7WR+PbhyNYQ0d+y0U5zz84TyArIs\nJ80LNnfSXirNBNW23QVR6Mu9ozDgzHqD0EAzjlhfawAQh4a7nGn1HFszDoijoHR85W+HMvoJfM5w\nqZfxqesv4JNC+uup34zhuX73fZ8+UsCirYRw/vzGofoPdUDW2n8M/EdgB/i3wC8Dfwo8GHhekiRv\nm6aR03BAwz6A1fMDOygMzPGpxmv846zcA6hKmVUVblXbL2Ba7azq6KQ53TSjmxZl6s47pc2dlHYn\npd3NKZx3ZnnhV1hodzOyzHFuo0FWVuKtlhFPYAxX3GWFMDBEoWFjtcFKM6IRhwQBtBoRcZkurFZr\n6Dkh+lHQQeNkp5m6fgGdFNJfT/1yQIc7oGFFCP8OeCywDvwR8E1JknzGWnseeB8wVQd0ori9ObED\nP7gDb5ljMPyuvtT3f7mb/mMzUJU2zA5zYHNfDrH/XFWtd0k2z+2/gBBC1J9hDihKkuQvrbUBcGeS\nJJ8BKKvhxt3Ke64MRgLGVZND6f0HlPNIfcVcGIIr+g4lLHcwdRQEJiB05ZI+5fjO2kpEqxHSzXLy\nDAoKOt2CdRcTRwa3Ct28IC/TeXlRsJaG7HZ9cUPh/JI8u+0Uh5+/dCH0BQ5xHJAVBTu7IasrMWFo\n2O1kbKw2fJVc6COfokwlVts5GDP9nPWyz3eQ5uXQX1fNixb5jMIwR/Ln1tqbgDX8nkD/DngT8ETg\nsydh3DTYn4bqrV7g+pNZiyrd5vxio2EQQDAQxfTSdIGvjCsceXluXjhfzRYbVlzE1m5GXgREgXcu\njdiP51zY7NJohaxSbvUdFbSavkLOp+gcRe7opBlBYFhrxWS536soCgOMCXC7GWEIZ9aadLKcOA4J\nCnCB32TPlM6nV06+p3R8Ou/lUe3TjDQvh/5l1DwvhlXB/SDwh8B7gEcCXXza7RuAH565ZVPisC/f\nS6tbytRW1d53XvU46C2RQH+OThVdGV8VV10nCvv9orA/s9uYqozaO69w4C4EptrQzvUKGoqy0KEo\ny8u7WU6eF33H2Jtj5AaWEFJOTghRb4ZFQOG+8umfKf/rYa1tJUnSnollJ4QZqBwzxuzxyK4cKzJV\nuyryjuhVyxkgzXwyLzCGlUZIXhSkmeOy9YbfMygruGyjSTct6HQz2t2CLMt95Z0xrJRzhnY7Od3U\nr8R9cSf1Tgi4uJuy0gw5s9pgtRXRzXzf3U7OSjNirRURBD5lF4VBWfZdlpfLDwkhasowB/RWa+1/\nA347SZI9M5/K1Q1uAP4JPiV3ajn8+9uPGvUeDURBVY2DoYp6fP896UDjx5z6Uf2lRQvVigzO+aV5\nqlPDoEq34SfOlgTlNuM4H4ntZ5yquGWroBNCnDzDHND3AT8KfMRaewfw9/j13O4DXI7fpuH7Zm3g\nPNmfpitcP8HlB/t9wUAUUi4e6p1CFAZAgcsdUWjICYjDAtMIe2m7LDN0s5yzqzGdrCBNM4rCXyAt\n5yplRYEzfp+hze0uRR7hnCPLctZWGhjgaxd9pFVt2RBHAUEBGTlh2F9Drj8Jd3rvjXLiy8Uy3vNl\n1HySjDQR1Vp7HX5h0gL4XJIkn5y2IbOeBzQLqi/3g17SOT/fJ88Lv1tqb/03x9ZuRjfNycpdV7d3\nu3TSnN12Rrub+Uq53PXSdGnmnVIQGM5ftkIUBuUcoYCz6w1azYiVZsRl6w0aUUij0Y+KqjEpB/2t\nvqf6HtQnUlr2Lwvpr6f+acwDWuQKuEnnAfUoHc7Unc5poZoXVD7y/7h96bkDl8nxfQJj9oxD9VdK\n7T82DLxEbwUHV5Zw9yfZOi79/M7yT7IuzkcIsXgs1HyeuuF9g6nWpgaqdarBBAG4ggBDEQZERdEr\n5b5sw5dzVwUK2WpEO819eXaas7nTZXMnJctDsjxnpem3/d7aTbn9zjaFc7TikI21mG6as9qK6WY5\nu+2UjdUGrWZOGEAjisp5QmVaMDC9bcRNbzVtIYSYD0cuRnpSLNoihf3z+0UBVVm1MdXzAUGVLotC\n4sjv8dNshDTiiGYjIo7982vNmGYjIgwDWs24nJBqev26WdErcMjzMiVnfJ+88Cm7OAr8fKLCEQT9\ndhj2t/f2gZmPqo6zcrYWplxc/ZOyyJqnYbuYPrVxQIu8adSwCZ+DC+cYY8rHplf+HZRzgozx7SgK\niENDMw56Wy4451hpRr1zQgNpVvj/0pw0z8lyX5qdZn4CbKebUW1L4Vds6DubqqS8yplP8kc26ftW\nh/s8bduncY26ah5kkTVPw3YxfZSCmzHGQHjAF3wchhRF4Xc9zfr7Da00Qtxag/PnVtntZGy3M279\n2g5xFLK+EnPHVofN7Zx2N+dinALecXW6OautlGYc0WyEbKw2WG0VrDRCWo0IAxSGMhraN14lhBBz\noDYR0GnmqF9RvawYVZTkHVde7K2wq5pp7tezq6iW3gn7IVNZCl5eM6BXFDFY46Bfd0LUn0WugDuK\nWkVAoywCOJg2GrZI4DjXmnX7oBRX9XwYAKEhMAE0Ic8KHH5cZ7UZ4RxcfrZJmuZsdzKgQbMRcOFi\nt5x0CqlztDsZWe787qyNkLQoOLPaJMtC2mnOajOiUY4RRWFQOryBooR9th6UVqzD+zjsntf1/h+3\nPUzTYe9DHew+iXu+KPrFwdTKAY2Tyz0qv1uHXPNRX5Dgx3PCIMQ5R5gH0PDnddMC5zLWWiErzRXa\nnYyVVs6FoEMYBERhyMXtbrl3Eex2clrOn+vHhxxZ5nAbTdYDw27HzztaD8OySKF0fux1ONUf/WFO\nc57tcXL687Z1mu1hxw97H+pg90nc80XRLw5GKbga46vrvCMIAtPbkTWs2pQpNcDRX7jUQdku/Ird\nuaNw1YZ6RX/H2JMWJIQQA9QqAlpmjDHEYekUnF9ItBEFpHEI+Khmp1PO+UlzdnZTNlZidjopO7sZ\nO52UrXZKnjs21mJ2MAQ7hs2dlNVmxJm1BmfXW7S7Oa1mSCMKCUNDHAQ0m/115qq9hXx7Tm+GEGIp\nqE0EVIf5AeOUI0/7vCrtFZT/huWmeWEYEAQBWeEjnzjy67sFYUAYBvgFtQ3GBGS5IwgDAlOO7Dgf\n/cSxr/Pujf/g24ExmLA/P2jvXqvD46Npaa5T+yhmcc/r1D6KZdc8q/dqmalNBFSHPO04Odtpn3fJ\nNQxEoSHN/JhMMw5oOz/es9KMKJyfeLq+GpMXBc084Oxagzu3u37eT+HIA+/Q7tzsYAw045DN7Q7r\nqw0CE9JNcz/2FPRXRQh6K8cN/wOatuY6tI9ilve8Du2jWHbNs3qvlpnaOCDRxxgfjQShIQz8B7kR\nB6w0Q9LMTyQ9sxYTRwFbO12acUia5dy5lRKFAV/b7PidWsso6sJWzlfv7HD3u+Zctt5kt+v3EdpY\njcnygCwvaDUiv2leOe5kcP0BJsoIST/qhBBTpDYpODE6vrzTt/upswNKVKn69J3RoBPZ395z+iXO\nRr/ohBDTRRHQEcx7HkCVSzbOr+8WR766LY5Dzqw1WGtF7LQz2mVUc+FimzgOSNOCnXbq+4YBaZaz\nudWh08lYX41ZbcVsbndYW4lZW4lpdzNWmzGtcisHv4hpNS61XDnted/zeSPNy6d/XsgBHUFd5gFU\nu52GUUCWGzAFa62YLC9oRCHttCArdjiz1iCMQm6/s+0XMU1znHFEgaGTFTQaEZ20wLmUZiP0exYV\njrueXfGLmhZ+TKgoHM5AI+rPE6qc0GD7NFKXez4vpHn59M8LpeAWiF40ZPa2g7KSrRmHfq03oBkH\ne5xEUfixo26akQ/MD8py3+6kfkFTwC9cWv5NFQf8cZ1m5yOEODkUAS0g1S6ncTMkyw1pmkMzZLW5\nwvnLWmztpFzcTtnaTblzq0O3XDX7jq0ut9/Z7q0xd3at3AIiCLj8shYbKzFn15ust8ptIuKAuCwD\n96t2+y3Iq9ScEEIcB0VAC0oVBfkFDkzpHAI/Ryjwc4KqNtCbXOq3Bvf/RZEf73E4mnHoq+CcI4qC\nXlGCKRs+7TYfrUIsM3/8iS/M24SZoQhoQemtJRcZv2G38bloFwZEYcFlGw3CgLIcG3COODK0XEi7\nk5Hmjtvv2GV9LSYyAV/+6haXbTRxDvKiYH0lptWIvEMKfSSUl1tGEPqiCDOw+rYiIiHEuMgBLThh\nEBA2ApxzZUGBr5BrdDJazYizGwW3XdihEYVEUchXvrZDHEeYoODidoeLOylF4dhYjfnaxS6Xn+1y\n+WUt7tzq0ogjLj/botUIWWmGRKFPxRWFI45DKPpR0f7ChNNeqCCEOD5KwZ0S+kUJg7uuDjy35/je\nc8Cn2/q7ppbRVPnfYCGCG3k+kJyPEGI4ioBOEQZfINCIIQpD0qza3nuN7XbG5k7MXc40ubjtixHW\nV0J22zkXNjt+p9WdLhe3u9zyxTs5s9bg8stW+fvbNlltxlx+tsVqM2SlFRMFAa1mQCP2qydUa9dF\nYbDHPRm3t2JPCCEGqaUDmtaksEWaODYVPaaKfYAgJAwLotDhMKw0yxdykGZ+6Z3CQTsowBg6ae7n\n/ziII0Mc5WztdDHGUJQLobYafpvwMAhYWVkpq+kccRRCb2u7/ryIaqfWS1ZoGIjGpnk/q7TfUt3z\nCdrTsHca11gEzWK21NIB1WkxwpNiFnqCoJwnZCAKfWVbIw566741oqrdZaUV0elkFIWjmzp22imr\nLV+EEBrD+krM1k7K2kqEAXZ2u7Sa/jpZlpevF/iiB7+YXa9tzN4xoaIotxR35ZbjBt/uNV1vNe/B\nJYAqnb1/h/TZ27/fZ++E2v5yRIPtk+IkP+fT+MxP+xqLoLkuHKcSrs5betfSAYnJ6O/fY/xOq4RE\nKwGF8xvRbaxE7HQz0rTgrpe1uP3ONrvtlCvvssJtd7S5uN3hy7fv0E1zNne6fPbvL7Kx6hc9jaOA\nK86t0mqGrLZizq41aDZCzqw1WWtF5WZ5hjOrca8kHPyK3q1GSGAMeelw4tD0Vl2IQj9mVc0xKhf8\nxhgoynYY9NepM8YXXhgzOBrlyPP++xCGvVXwcIV/3HN+pTOq3itTeT3AOdNbC5wDVgTvj50d5y4J\nISrkgE4RR6UcqnlDrmzn5e6pxhjyvCDLfVWbw0cvwJ6tGqov8mpOUBj028AlbfBjUs45nDG94oZq\n4VT/Xz9Cwvk0Xq8P9Pow8Pz+yrv+dfuLqg4+b0ywxyZc/z0YDKD8CuDVtQcO9N6/wWOXvrdHMZga\nGrR/WpymX/xiOaiVAxonr3tUrndeueaTymWPetwY4+fsGP/bvtWICLMMnONud1nh4naXRpQShRus\n3LlLFMJuJ/MrJmQF7W5GYAxhFPDXX7iDZiNgfbXJl7+6xWqrwV3OtsrdWw1RFNJshKw1Q6LYp+rC\nMKAR+bGlAEMU+fRcGBjSrCgLGHzRRF4UmNLBFEXl5EKyoiAAoijw24vjaMUheak5CgPSzAF+zbso\n8pvyhaEhz32fOCqdEI5enWDgiyfKSVS9yMs53yeMgj2O1FcT+kVhq/ShCaqVyV3pkH07CPxmgT7a\nKp1mOUbnnC+Z7zveftoSyjQk/Whs/w5NVbJyv+sa15fV4fM8rc953bWJg6mVAxonr3tUrndeueaT\nymWPczwIDMaVX56hj2QacUQnzWk1ItprOV+8bYu7nVvh8jMtPv+li/hII6Ddycjygm5a8MUv30YY\nhkRRTBRFtJoRzUZMHPodVl0B19zzMtIymrr8TLP3C//MWoMsL3qRVZblFAW0mmHvD7WbFr2Vt4Mw\noBkFYGB7N2O9TPN1soJGFPiqu3LFhryATjdnY8UvIbTSDHu7whaFY6Xpzw337TG+2oqh7FP9W5Wc\nR1FQrgxuBvpHpZ6+nYORV1BqDQIDLiudjHc+rUYIzlH07g97Aqy8KMrdavspwF7EVfZ3VFFnGUEO\naPHSDv8yrj4TgxFXHT7P0/ycj3LevPWIvWge0BIxWH3m/+1XqsVR4AsDynYchZQ/2nv9wzDsLWrq\nXLmQqfNf2kUBuXN0s7yMDpyPVoqibBfeSTlHXhR+XMo5v3Nr2acofL+qT7WcUFGU1+XSPnnuI4jC\nOdKswBWOvChTiPuuWZTjTpX9RVEArndssO3tZo/ewbZ3Cv0fQ9V/DDyPq47tTduN+r1UdRu5v77w\nxIIx0wjIWvsw4BVJkjxylq8jjqbMwBEaiCKDCwOasS+xXm1G3GWjQSct2G5n3OP8Glu7KRe3u1zc\n6XLHZpfb7tjhbuda3LHV5u+/fAfdbofd7ZQvbV8k7e5CnpJnKZ/8SJvAOEwQkXZ3KbIOeWcLR+j3\nM2o0yYqCMAiJG00uXvgKFAUra+t0u13Szi6RyWl3dgmjJmtn7sL25gWMCVk/c45ut41z0Gqt0Fi5\nDMKYlZVV8iwjd45mHJGnXQgiVlbXaKycBROystoijmMKByutJmfXW2BMbxyrFYe0mhHgWGnGrJSR\nWZrmEAQ044BGHNCIAhpRyForJop8Mq+bFcRRwGUbPuILy/RfFf0UhcOVablWHBAGQbnbbUAYGlzh\nnVschWVaEAbHyZpx6JdbGriXgz8MALwr9URBf/zMVINiFc5HVkdl6/pFGz7dFxxwQnX9/ZHVUYzb\nX5xeZuaArLX/BngqsD1K/zrkacfJ306j70m3fVpn8IvMPxeFAbkxRAW0GhGdbk4zDtlYbbC5mxGG\nhpVmTDctyLICMDQaDbY6bfI8J087ZGmHIktJ25uYIAQTkHW26GzfQdrdIYqamDDurR8XxS3YvsjW\n175EGK+wu7sNRUa3vU2Rdci6O6yduwfd7pdpb91Oa+Nyul/tkKdtwrhFu91mJTXErXXa7Q7OFRgT\nslXkZGmb5upZdlNDK4uJGi3aGTQaBVEUkRcpDu8AKqKzLbZ2U5zzqbndbo5z0E0LwqAgzQLWiSgK\nSDPHaisut0cvaKcFjTgky30aLQeKHOLQb3lRlOmz2DjS3PlUqPPjRQbIqpVijS+EqAozfBl7ZaF3\nBM65XsXhINWnr3JOznHJve5/Rvu1fod9dvc+5197f99LI+rR/pYOcj51+Bs5iFn1PUnmtZjpKOXf\ns4yAPgf8c+Ato3SuQ552nPztNPrOU0/1JRAEEBSOwvjUW5Y7Vpp+lew8Lzi33sCU7SiAwhVsbXdw\nwMHcWjEAAA6OSURBVPrGOnnWJQwcnV1Ii5wgbOBcTpGluCInaq6SpW2KIgdXYExAnqWknR3CsEHW\nbZN1d4maa+RpB1fkpN0d0vYWxe1/Q2v9PGlnB1cUxCsbZJ0dTHubqLkGDjrtLRqtNX99IIpXyLMu\n7e0LNHGEUUyedSmaa+By0q7BuTUMBcY4mk0ftdxxsc3qSowxhgsX/U6xPkwoCwbwzggHQRywtZv6\nsajQFwNst1NajRBX5H77itCAqyr2TPl+OsLA+dQiBQ1jyHN/tKjSfYEhxPWClrxwZIUjLGslqogq\nCHwZgqu8GFU1YX/l8yr9BxCG5XP+gO9btqt/BseceiUOvev7Uv6DxpMGI6FRPnsHUZe/i1HtOm5f\n4ZmZA0qS5O3W2vuMc840wvJlD+3H1W+MIY4NMf6PptkIObvRIM/92Mlu129Ul6Y5F3dSdnZTLn7j\nVVzc7rK53eXCZpsvX9hlt53S7uZ8+dbbabd36LZ3wDnSzjbd3U12L95GkacUua/A67YvstvewgQh\nYeydQHPtMtL2JkXeJA/bFEXB7ubtRA3v1DpbF8i62zRWLiPPunR2LlDkGbthzPpdvo7W2jnCuElj\n9SxhGBNEDcJ4hdbqGmEY9fSmaYpzjnNnVmnGEWsrMUX5JbrailhrRcShodGI/I6zccjZtcae980B\naV5QuDLVBlzc7nLV+TW6aU67k7O2Ur4mfj5UGHqHhfNzmYIA0txHOYGpvvh9UYUrixCM8SXyhAEh\nfvNBM+gpAFdAFO11PD791496B+dUDcRF/XGm0un1ihlKW/qrWVwaUU36mZv0nHle96RYW2305tAt\nOufPbxzZp7ZVcJNQp7B3HhxXvxu4huGgFE6/T9U/KH9KG2PKbcMpS6kDoCi/xAKqL8Ig8Ok+HARB\nCMYQBH58qJqv49N0hsCEmPJ5YwwE5W/zMCQIwt6XbBBGmMA/X51rgqC8piEoz/ffrab/fBD001mO\nXrvS0iu+CPZGMoPvSS/9VD4OTFWKzQHvYRm1mH60cVi/3jWq6GbAWew50e291p5QhoFIprxj+x/v\n/bhUbsf1mwPF3/v1+5cbTPEezEFjPrP4O63r3/84TnF7pztDS06W227bBIY7olo5oLowj1xzHehV\nvWF8uXZgiMLYR0BZwWrLl27vtjPyYo2d3ZQLm10ecF/Hbjfjwp1tUnuerZ02F7d2iMKIIs+5cOdF\nmo2IokjZ3dkhJybPOuxevJ14ZQMK2Om0WVndwLiCdnuLVrOFKzJ2d7ZZWTuDMYb2zjat1RUMAVmW\nsnH2LFEYs3nxDi4/fwVxHNPupKytrhBHEVs7HS47u06r2aCb5sRRQLMZk+U5a60GK82YvChoxJEf\nw8lyVloxrUZImhZEUUAUBhSuYLUVe0fkwBmIwwBHWThgIMscYWRoNnwEleY+mqy+9H2JtSMIg97k\n3DgKelGMKSsS/RhP9QvYYfDnAQTldutUjqxMpVXVi740PMAv7zBwU6kiLO8cB9fo60dT3nGVn9be\nWFLvMvu+RPcXH7j9jnBf3zqyrH/ndUIO6ADmkWuuA/3BZBj8tRtFAc2GT/Gs5gVn13wJdpY7rrqi\noNPNaXcyOmnO1m7m000BtLsFO+2UNPN7FPk/zoJbv7bLSjMijgK6WcFuJ/NjTKEfgwJYa/U/mmFo\nSDOHwX+pd1MfWTUaYW85n7NrTdrdDPBjWWnmn19pxb1S7LwoOLPW9I6oUY1z+dLsatxrMKIZLFFv\nxCG7ndxPmi0ntUZB4HWmRT+6gV5KLY7CPSksvzZf/z2uIh9fMe5685SqL/Yg6E9c3Y8DnPOFFK48\nf/C1+n24ZP5T9RrBXg8z8NqHf0YGzx+kpj5mKMv6d14nZuqAkiT5PPAts3wNMXsM5bdlb2yh+sUO\nhTO9VQdMZnrbehvjJ4tWhQ3VlynGlzNXYy5R6Eub24UjDH3aLs3LVQUY+AINIc3pXTPNijJN5teM\nC8rnu6l/Pgr986ExmMBXljkMRdEft/Ff2v6LvIoMqnGTasyjKhAwxvTWpKscTeEckfGb9DkHpr9B\nbC8qqRyPdxJ9pwb9L20z+ICBiKEMKy6NNlzvnvS+70x/3GbPGNHA2NKee3qIx6hrtLJs1HkB0Wmi\nCEgcSfmd3JtDBH42vj9WVWKZcqJm/5yi+qLETxL1qT2/akLPu5i982XAP++MH18KAkNWOZUoKKu8\nXNWNIPBOpiplDga+lMuhHX8dA/nA2Iz/pWou+eVene8G2lUxgFuNel/QxcA42Wrr0j+jg8Y9Dnl3\nqartjhpv6Z0x6GDo359LnczBzwtRF+SAxEj0voT3VFhVx0xvbbPB79xwYHwg9HXDPpIwfvKlw49r\n9PqHA88H1Twd16seqxwVrj/+4SOuMr01YGvPNuhFU/triy4Z16jONQxsptfvO9g/HMG5jBdNmEPa\nl1ag9a+/t33Q+IMCGlFn5IDEWByauvEHDz3HOT9X31SpoqAfXVRfuAZw1dYLxkdKYZnqc4Xzc2vw\nkUlRxlZFUQ6+HxDNVBQDUQwMlCn3Db9Ei+s5ntHmfxw0SD8t5ETEaUUOSEyFalD9MILokvjjkDTV\npeMdobl0XkTI3sqrw2ZOHPgao0Qve7oP769xEyEm43TMeBJLySjf+3IOQtQXU5eyQTcFQ5a9Bn/Z\n9S87y37/66rfjPEr6LbbNusn4JicP79xqP5aRkCjLPlxnHbdWQT9s7zepO2DiiMWBd3zemoWs6WW\nDmjZJ4gtgv5ZXq+ummfJIujXPRfTppYOSAghxOlHDkgIIcRcqJUDGievO07ed9a55nnksid9fxax\nPerxUfouSvsgllH/pOM7ddCg8aijURXcKUL6pV/666ffjOGJVAUnhBBCnAByQEIIIeZCbRxQHfK0\n08xN18HWkxqPOOr4orePe3wR2wehez7ZPReHozGgU4T0S7/010+/xoA0BiSEEKJmyAEJIYSYC6fC\nAU073zrvXPu4Nk6DRdA8yKLpn7aN076e9NdT82nnVDigRVujahr2LqPmQRZN/7RtnPb1pL+emk87\np8IBCSGEWDzkgIQQQswFOSAhhBBzQQ5ICCHEXJADEkIIMRfkgIQQQswFOSAhhBBzofYOaBkniC2j\n5kGkWfonbYvFovYOaBkniC2j5kGkWfonbYvFovYOSAghxOlEDkgIIcRcqJUDmjSvO2nfOrePYpzz\n6qBnlvqPOr4o7aNYxns+7feqTtoE2pDuNCH90i/99dNvxvBE2pBOCCGEOAHkgIQQQsyF2jigOuRp\nx2kfleudt32zaB/3+KK2j3t8kdu658fXLA5HY0CnCOmXfumvn36NAR0+BlQbBySEEGK5qE0KTggh\nxHIhBySEEGIuyAEJIYSYC3JAQggh5oIckBBCiLkgBySEEGIuyAEJIYSYC9G8DZgG1toA+FXgOqAD\n/L9Jknx2vladPNbahwGvSJLkkfO25SSx1sbAG4H7AE3g55IkeddcjTpBrLUh8DrAAg74kSRJ/nK+\nVp0s1torgI8Bj06S5NPztkeMxmmJgJ4AtJIk+YfAc4FfnLM9J4619t8Arwda87ZlDlwP3J4kySPg\n/2/v7kKkquMwjn81KnoxyJDsxSyCnhsJLyyid6iFEoySaKMw16jowgu9tFSQKCEjg4wyLyTDLHtF\niDQQTSiiIIIieIIiUOzFqDbLLnK1i3MWxm1n3dmL+e/MPB9YmHPmzH9/s7A853/+h9/hNmBD4Xra\nbQGA7euAlcCTZctpr/oEZCPwT+laojXdEkDXAzsBbH8KzCtbThHfAQtLF1HIm8Cq+vUU4GjBWtrO\n9nvAI/XmbOCPguWU8AzwEnCwdCHRmm4JoHOAwYbtIUldcXlxvGy/Dfxbuo4SbP9l+7CkacBbVLOA\nnmL7qKRXgOeBraXraRdJA8Ah27tK1xKt65YA+hOY1rA91XZPnQX3OkmzgD3Aq7ZfK11PCbYXA1cA\nmySdVbqeNnkQ6JO0F5gLbJE0s2xJMV7dMkv4mOo6+HZJ1wBfFa4n2kjS+cCHwFLbu0vX026SFgEX\n214LHAGO1T9dz/aNw6/rEHrU9k/lKopWdEsAvUt1FvQJ1RrAksL1RHs9BpwLrJI0vBZ0u+1eWZR+\nB9gsaR9wKrCsh757dLA8jiEiIoroljWgiIjoMAmgiIgoIgEUERFFJIAiIqKIBFBERBTRLbdhRweT\ndCnwLfDNiLcW2N4/yvEDwGbgPtvbGvYvA9YDl9n+QdJx21NGjH8cOI2qbcsS2wckXQK8QNXGZmp9\n3FLbv4xR8ylUPeduoLr1f5Pt50Ycsw6YYXtgnH+KiJ6SAIrJ4qDtuS0cfwC4G9jWsG8hzfugnTC+\npLVUbWvuompkuWU4zCStoOotNlZvvSXAecCVwBnA55L22f6iHuMWYAB4v4XvFNFTcgkuOtVHwLzh\nljOSZgOHObEn4Fj2UbWtAZgJnNnw3gZO3lH7a2CN7WO2/wa+B2bVtUyn6kj91DhriehJmQHFZHGh\npC8btrfaXjfG8UeBXcB8qm7Y9wDbgTUn+0V1+/5+qhZOACuArZLWALuBD+qxmqq7rg+Pdy1wNbCo\n3rUReJw6kCJidAmgmCxavQQHVUg8TBVAd1KFUbMAagy404HPqJ4dhe2dki4CbgZuBZ4G7q3HHJOk\nm4DXgftt/y7pIWC/7d31WlVENJEAik62h6rz8xzgV9uDkpodO2rA1ZfLVtleTvVMqZ2SngB+lDTD\n9qFmA0paCLwI9NveW+/uBy6ow246cLak9fX4EdEga0DRsWwPUXXBfhl4Y4LDDAJ3SHqgYd/lwM/A\nb80+JOkqqvDpawgfbPfZnlOH3WpgR8InYnSZAUWn20619rJjIh+2PSRpPvBsPfM5QnWL9oI64JpZ\nSfX/s6Vh1rXa9oTqiOhF6YYdERFFZAYUk5Kk5cDiUd46aHt+m2rop7pD7n8mcMNERIyQGVBERBSR\nmxAiIqKIBFBERBSRAIqIiCISQBERUUQCKCIiivgPRwRXn0A9w6MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11f207588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "skew=(XID_MIPS['FErr_MIPS_24_u']-XID_MIPS['F_MIPS_24'])/(XID_MIPS['F_MIPS_24']-XID_MIPS['FErr_MIPS_24_l'])\n",
    "skew.name='(84th-50th)/(50th-16th) percentile'\n",
    "g=sns.jointplot(x=np.log10(XID_MIPS['F_MIPS_24']),y=skew, kind='hex')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The uncertianties become Gaussian by $\\sim 20 \\mathrm{\\mu Jy}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "good=XID_MIPS['F_MIPS_24']>20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "120282"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "good.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "im100fits='../data/ELAIS_N1/PACS/ELAIS-N1-100um-img_wgls.fits'#PACS 100 map\n",
    "nim100fits='../data/ELAIS_N1/PACS/ELAIS-N1-100um-img_noise.fits'#PACS 100 noise map\n",
    "im160fits='../data/ELAIS_N1/PACS/ELAIS-N1-160um-img_wgls.fits'#PACS 160 map\n",
    "nim160fits='../data/ELAIS_N1/PACS/ELAIS-N1-160um-img_noise.fits'#PACS 100 noise map\n",
    "#output folder\n",
    "output_folder='./'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "from astropy import wcs\n",
    "\n",
    "#-----100-------------\n",
    "hdulist = fits.open(im100fits)\n",
    "im100phdu=hdulist[0].header\n",
    "im100hdu=hdulist[0].header\n",
    "im100=hdulist[0].data\n",
    "w_100 = wcs.WCS(hdulist[0].header)\n",
    "pixsize100=3600.0*np.abs(hdulist[0].header['CDELT1']) #pixel size (in arcseconds)\n",
    "hdulist.close()\n",
    "\n",
    "hdulist = fits.open(nim100fits)\n",
    "nim100=hdulist[0].data\n",
    "hdulist.close()\n",
    "\n",
    "#-----160-------------\n",
    "hdulist = fits.open(im160fits)\n",
    "im160phdu=hdulist[0].header\n",
    "im160hdu=hdulist[0].header\n",
    "\n",
    "im160=hdulist[0].data #convert to mJy\n",
    "w_160 = wcs.WCS(hdulist[0].header)\n",
    "pixsize160=3600.0*np.abs(hdulist[0].header['CDELT1']) #pixel size (in arcseconds)\n",
    "hdulist.close()\n",
    "\n",
    "hdulist = fits.open(nim160fits)\n",
    "nim160=hdulist[0].data\n",
    "hdulist.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in PSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pacs100_psf=fits.open('../../dmu18/dmu18_ELAIS-N1/dmu18_PACS_100_PSF_ELAIS-N1_20170720.fits')\n",
    "pacs160_psf=fits.open('../../dmu18/dmu18_ELAIS-N1/dmu18_PACS_160_PSF_ELAIS-N1_20170720.fits')\n",
    "\n",
    "centre100=np.long((pacs100_psf[1].header['NAXIS1']-1)/2)\n",
    "radius100=15\n",
    "centre160=np.long((pacs160_psf[1].header['NAXIS1']-1)/2)\n",
    "radius160=25\n",
    "\n",
    "pind100=np.arange(0,radius100+1+radius100,1)*3600*np.abs(pacs100_psf[1].header['CDELT1'])/pixsize100 #get 100 scale in terms of pixel scale of map\n",
    "pind160=np.arange(0,radius160+1+radius160,1)*3600*np.abs(pacs160_psf[1].header['CDELT1'])/pixsize160 #get 160 scale in terms of pixel scale of map\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  0.           0.33333333   0.66666667   1.           1.33333333\n",
      "   1.66666667   2.           2.33333333   2.66666667   3.           3.33333333\n",
      "   3.66666667   4.           4.33333333   4.66666667   5.           5.33333333\n",
      "   5.66666667   6.           6.33333333   6.66666667   7.           7.33333333\n",
      "   7.66666667   8.           8.33333333   8.66666667   9.           9.33333333\n",
      "   9.66666667  10.        ]\n"
     ]
    }
   ],
   "source": [
    "print(pind100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x121cbc940>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABGsAAAItCAYAAACO+PfWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3X2MpFl9H/rvdE/39MzO7M4uC5tr2TFxDEf8RRAkgIyX\nlWJQcMTFiaL8gcKNnVgIC/via+fFL4DlyCiKLmD5EhHHG9Y4lizFJtnIWMJYimIw2I51MZGMAsde\n7NiRL2Zh2deZnel56ftH9Qy1M089Z+qprqqnez4faaTu59ennqdOPVV15pzfOefY3t5eAAAAABiH\njXVfAAAAAABfp7MGAAAAYER01gAAAACMiM4aAAAAgBHRWQMAAAAwIsfXfQEAMFbHjh1b2ZaJe3t7\nx1Z1LgCAsVpl+ysZbxtMZg0AAADAiOisAQAAABgR06AAYIZjx0aZFQsAcGSNqf1VStlI8sEkL01y\nMcn31lofmYq/Mcm7k1xO8lCt9cGp2AuSfCbJ62qtXyilfGuSDyfZS/K5JG+vtV6ddW6ZNQAAAAA3\n+64kO7XWVyf5kSTvuxYopWwl+ekkr0/y2iRvLaXcNxX7t0menXqs9yd5Z63125McS/KmvhPrrAGA\nGY4dO7ayfwAArLb9dQttsNck+fUkqbX+bpJXTMVekuSRWuvjtdbdJJ9Kcv9+7L1JfjbJ/zf19y9P\n8on9nz+W5Dv6TqyzBgAAAOBmdyZ5cur3K6WU4zNiTye5q5Ty3Um+Umv9+A2PdazWujf9t30ntmYN\nAMwg4wUAYLVG1v56KsmZqd83aq2XZ8TOJHkiyf+ZZK+U8h1J/lqSf19K+d+TXO3425lk1gAAAADc\n7NNJvjNJSimvSvIHU7HPJ3lRKeWeUsp2JlOgfqfWen+t9bW11geS/Pck/0et9S+SfLaU8sB+2Tck\n+a2+E8usAYAZNjaMaQAArNLI2l8PJ3ldKeW3M1kU+HtKKW9OcrrW+nOllB9K8vFMEmEeqrX+ec9j\n/XCSB/c7dj6f5CN9Jz62t7fXFweA29b29vbKviR3d3dHlfMLALAOq2x/JeNtg8msAYAZRjZnGgDg\nyNP+mhhVfhEAAADA7U5mDQDMYGQHAGC1tL8mZNYAAAAAjIjOGgAAAIARMQ0KAGaQhgsAsFraXxMy\nawAAAABGRGYNAMxgZAcAYLW0vyZk1gAAAACMiMwaAJjByA4AwGppf03IrAEAAAAYEZk1ADDDxoYx\nDQCAVdL+mlALAAAAACMiswYAZjBnGgBgtbS/JmTWAAAAAIyIzBoAmMHIDgDAaml/TcisAQAAABgR\nmTUAMIORHQCA1dL+mpBZAwAAADAiMmsAYAYjOwAAq6X9NSGzBgAAAGBEdNYAAAAAjIhpUAAwgzRc\nAIDV0v6akFkDAAAAMCIyawBgho0NYxoAAKuk/TWhFgAAAABGRGYNAMxgzjQAwGppf03IrAEAAAAY\nEZk1ADCDkR0AgNXS/pqQWQMAAAAwIjJrAGAGIzsAAKul/TUhswYAAABgRGTWAMAMRnYAAFZL+2tC\nZg0AAADAiMisAYAZjOwAAKyW9teEzBoAAACAEZFZAwAzbGwY0wAAWCXtrwm1AAAAADAiOmsAAAAA\nRsQ0KACYwQJ3AACrpf01IbMGAAAAYERk1gDADEZ2AABWS/trQmYNAAAAwIjIrAGAGYzsAACslvbX\nhMwaAAAAgBGRWQMAMxjZAQBYLe2vCZk1AAAAACMiswYAZjCyAwCwWtpfEzJrAAAAAEZEZg0AzLCx\nYUwDAGCVtL8m1AIAAADAiMisAYAZzJkGAFgt7a8JmTUAAAAAIyKzBgBmMLIDALBa2l8TMmsAAAAA\nRkRnDQAAAMCImAYFADNIwwUAWC3trwmdNQAAAAA3KKVsJPlgkpcmuZjke2utj0zF35jk3UkuJ3mo\n1vpgKWUzyYNJSpK9JG+rtX6ulPKyJL+W5I/2i/+bWut/mHXulXTWXLp0aW/mBRw/nsuXLx/Iefb2\nZp5m7fp6B/vqYMhzGnM99F3b9vZ2dnd3bzp+9erVAz3PLEN6cIeU2diYPftwa2srly5dOrBzDbGq\ne67v+cyqh1Xd26u6f/rMej/03T8HbchzOujXaNbn4/b29kreEEZ2OOyuXLnynDflxsbGXN+rrff0\nsj+XF30PLnp9N5Y/yHbrvOdetUXPf+Nr19fGWcb5j5pZ7YLDqvX6HuT377x1d9i/+w/6c/PG9+6J\nEyeWXkEjew2+K8lOrfXVpZRXJXlfkjclSSllK8lPJ/nrSc4l+XQp5VeTvDpJaq3fVkp5IMl79su8\nPMn7a63vu5UTr33NmpG9EGuhDiZW+Z/QsVIHE+pBHVyjHuBgaXMsxmfScOpuMepvOHW3GPWX1yT5\n9SSptf5ukldMxV6S5JFa6+O11t0kn0pyf631Pyd56/7ffHOSJ/Z/fnmSv11K+WQp5UOllDN9JzYN\nCgBm0EBhbFrp2ABw2I2s/XVnkienfr9SSjlea73cEXs6yV1JUmu9XEr5hSR/J8nf24//XpJ/V2v9\nTCnlx5P8RJJ/MuvEo6oFAAB6XU/HTvIjmaRjAwDL8VSS6QyYjf2Omq7YmXw9iya11n+Y5MVJHiyl\n3JHk4VrrZ/bDDyd5Wd+JB2XWGNUB4HZgyggj9Jx07FLKK/r+eGNj46b7eHNzc3lXdxvY3t5e9yUc\nWidOnFj3JRxqOzs7676EQ0vdLWbV792Rtb8+neSNSX55f82aP5iKfT7Ji0op9yR5Jsn9Sd5bSnlL\nkm+stf7LJOeTXN3/9/FSyg/UWn8vyd9M8pn0GDoNauYiOwAALE1fOvZNblxMeHNzM1euXLnlk1lg\n+LnlV7nI67oX2D3oBYZPnDiRixcvruz8R83Ozk4uXLiw7ss4MKtcYHjeuhtZR8HcDvpz88b37m3Y\n6fpwkteVUn47ybEk31NKeXOS07XWnyul/FCSj2cya+mhWuufl1L+U5KfL6V8MslWkh+stT5bSvm+\nJB8opVxK8hf5+ro2nYZ21sw1qgMAh9Fhb7BxJPWlYwPAoTem9let9WqSt91w+AtT8Y8m+egNZc4l\n+fsdj/X7Sb7tVs89tLNmrlGd48ePN7fpvd2pgwkpirdlb3Un9eD9cI0pB/AcfenYAMARMbSzZq5R\nncuXZw/43Lhv+yLGnC7Z6qyaVQdDntOY66Hv2malKN6Ywr3oeWYZ0oM7pEzf6uZ9KcKr6mFe1T3X\n93xm1cOq7u1V3T99Zr0fVrk6/pDndNCv0awpB6vqwBnZbgSQdKRj9/1x1/t4+ti62wwHPU3poHU9\n/kGdc0wjx12WMQVt3ffbUbLs+2fsUxznfezpY+t+73kftGl/TQztrDGqAwCwYjPSsQGAI2ZoZ81c\nozoAcBite/QNAOB2o/01MaizxqgOAAAAwHIMzawBgCPPnGkAgNXS/poYbWeN1KeJVS16O5aFrrqu\nfcibdVXPZ5X36ZhfV+/XYVr1tu56XdVCy2P5/AEAgLHQZQUAAAAwIqPNrAGAdVt3dhMAwO1G+2tC\nZw0AAEux6DTHdZdfp3X/Z2XRuuu6/uljh/m1YfnWef+3zr2M98Y8lv3eab13WR2dNQAwg8YJAMBq\naX9NWLMGAAAAYERk1gDADLaOBABYLe2vCbUAAAAAMCIyawBgBnOmAQBWS/trQmYNAAAAwIjIrAGA\nGcyZBgBYLe2vCZ01ADBypZSNJB9M8tIkF5N8b631kan4G5O8O8nlJA/VWh+cir0gyWeSvK7W+oWV\nXjhrt7e313usK36YLHr9Q8pPl1lnqv7Ypwkc9nuL5Vrk/h1Sdp4y6/hcOUjLOP+6n9PtSpcVAMxw\n7Nixlf1r+K4kO7XWVyf5kSTvuxYopWwl+ekkr0/y2iRvLaXcNxX7t0meXUL1AAAcuFW2v8bc8b2S\nzJpWBXTFx1xpq+xZHHM9DNF6Pl0pb4uOei2zzCpfnzH3aB+1Hvxlj9jcapmu+JjvA5bqNUl+PUlq\nrb9bSnnFVOwlSR6ptT6eJKWUTyW5P8mvJHlvkp9N8qOrvVwAABYhswYAZtjY2FjZv4Y7kzw59fuV\nUsrxGbGnk9xVSvnuJF+ptX78oOoDAGDZVtn+GvP6OOO9MgDgmqeSnJn6faPWenlG7EySJ5L8oySv\nK6X8ZpK/luTfl1L+0gquFQCABVlgGABmGNFU1E8neWOSXy6lvCrJH0zFPp/kRaWUe5I8k8kUqPfW\nWj9y7Q/2O2zeVmv9i9VdMgDA/EbU/lornTUAMH4PZ5Il89tJjiX5nlLKm5OcrrX+XCnlh5J8PJOM\n2YdqrX++xmsFAGBBOmsAYIaxjOzUWq8medsNh78wFf9oko/2lH9gOVcGAHCwxtL+WjedNQAA3JZW\ntfvfEP6zcri1dm9c9PVd9+6Qy74/x3z/L7vuF338Mdcd87HAMAAAAMCIyKwBgBnGvJ0jAMBRpP01\noRYAAAAARkRmDQDMYN43AMBqaX9NyKwBAAAAGBGZNQAwgznTAACrpf01obOGJMNSzVa1ZeC6tybs\ns4xrm/WYY66HMVvVtqzLSNfsesyrV68e+HlmGXLPLePL1b0PAMDtRmcNAMxgzjSHXVdn5zwdoGPv\nLG29R9d9/X3Xd9g/X1p127r31v3aLGrR16/1/Nd9fyz7/Is8/rIH1Ra9N9d9b4/93roVh+EaV0F+\nEQAAAMCIyKwBgBmM7AAArJb214TMGgAAAIARkVkDADPYjQAAYLW0vybUAgAAAMCIyKwBgBnMmQYA\nWC3trwmZNQAAAAAjIrMGAGYwZ5qjbm9vb6mPv+7R0WWcf/pzYdn1t0yta19GfPrYuutu0Xtj2dff\nqr+Wsb/3Frm+Vtl131ssTvtrQi0AAAAAjIjOGgAAAIARMQ0KAGZYdxo5AMDtRvtr4kh11gx5Ucc+\np3HM1zfk2lplrl69eiDnGWIZz2dImStXrsz9mLOsar7nqj5Qj9p5VmnM76OjWN8AALCII9VZAwAH\nSUcSAMBqaX9NWLMGAAAAYERk1gDADLaOBABYLe2vCZ01AABHVNc6UqtcD2/RVPZFyy/juR6W9PzW\nc192vGsdwulj616XsfU6jvF1nq6z1vW16nfR8kf5P9OL3vuLPv6iln1vsDo6awBgBg0WAIDV0v6a\nOLpdogAAAACHkMwaAJjhKKd5AwCMkfbXhFoAAAAAGBGZNQAwgznTAACrpf01IbMGAAAAYERk1gDA\nDOZMAwCslvbXhM4aAAAGaaWqL5rKvuzye3t7C8WXadFrW3b86tWrc5efLtMqv6hF781l39tDTNdp\n6/VZ93vvMFv0uS/7tVm2rutf52fh7WwlnTWH+YvwoMqs8jFXVZ/LqLuuL+5lf5lfM+T5DLm2VplL\nly51Hh/ywT6kV3pzc3PuMkMMeT6rureHvK7L+OK9Xb4Yx/wdMfYGFQDAUTOm9lcpZSPJB5O8NMnF\nJN9ba31kKv7GJO9OcjnJQ7XWB0spm0keTFKS7CV5W631c6WUb03y4f1jn0vy9lrrzP94yC8CAAAA\nuNl3Jdmptb46yY8ked+1QCllK8lPJ3l9ktcmeWsp5b4kb0ySWuu3JXlnkvfsF3l/knfWWr89ybEk\nb+o7sc4aAAAAgJu9JsmvJ0mt9XeTvGIq9pIkj9RaH6+17ib5VJL7a63/Oclb9//mm5M8sf/zy5N8\nYv/njyX5jr4TW7MGAGYYUxouAMDtYGTtrzuTPDn1+5VSyvFa6+WO2NNJ7kqSWuvlUsovJPk7Sf7e\nfvxYrXXvxr+dRWYNAAAAwM2eSnJm6veN/Y6artiZfD2LJrXWf5jkxUkeLKXckeTqrL/tMjizppTy\n+/sXlyR/Umv9nqGPBQBjNLKRHQCAI29k7a9PZ7IGzS+XUl6V5A+mYp9P8qJSyj1Jnklyf5L3llLe\nkuQba63/Msn5TDpprib5bCnlgVrrbyZ5Q5L/2nfiQZ01pZSdTFJ4HhhSHgAAAGDkHk7yulLKb2ey\nKPD3lFLenOR0rfXnSik/lOTjmcxaeqjW+uellP+U5OdLKZ9MspXkB2utz5ZSfjiTLJvtTDp6PtJ3\n4qGZNS9NcqqU8hv7j/Fj+4vtdJ/k+PHe3rGtra2Bl3F0bG9vr/sSRuHkyZPrvoS1O3v27LovYRR2\ndnbWfQlr53NhYp33wshGduDIWfd7bMj55ymzt7c3KLaK+NWrM3eLPZB41/mvXLmyUPl5tF6njY3+\n1SBa5ReNDzFdZ4s+/qLlW6/PMt/bi94bi753WL51fzdM299a+203HP7CVPyjST56Q5lzSf5+x2P9\nYSa7Rt2SoZ0155O8N8m/S/KiJB8rpZSpuVvPcfly5+Ekk46aS5cu3XR8VS/QkDfjQb+Bt7e3s7u7\ne2CPt6oPmIOuu5MnT+bZZ5+96Xjry/ygDHk+Q66tr8zZs2fzxBPdUxeHvCdaDZEum5ubc5c56Mbu\nzs5OLly4MPdjHpRlN95vxazPhaHvh1V9Lhx0Pcy6F3TmAQBwlA3trPnDTLao2kvyh6WUx5L8b0n+\n14FdGQCs2ZhGdgAAbgfaXxNDd4P6R0nelySllG/IZMuqLx3URQEAAADcroZm1nwoyYdLKZ9Kspfk\nH82aAgUAh5WRHQCA1dL+mhjUWVNr3U3y5gO+FgAAAIDb3tDMGgA48ozsAACslvbXxEo6a8a8PdpR\n3JptDLsaDSlz8uTJzp3BVnVt09tJLrNMazeo8+fPd8aG7Oy0tbU1d5njx+f/WBhyba0P4a5d5Fa1\nI9YQQ65tlZ+NY/7SG/N3BAAArIPMGgCYYcydXHA7WPZ7cNmdwX2Pv2hHdWtgatnx1oBVV3x3d/eW\nyy/62rQGa1r31ubm5lIff0h8+jVZ9PFbVjXYtSrT99Oi99ayy2t7qINrdNYAAIxUKeWVSf5VrfWB\nUsq3JvlwJps7fC7J22ut86eSAgCjp7MGAGY4aiOLHC6llH+W5C1Jzu0fen+Sd9Zaf7OU8rNJ3pTk\n4XVdHwAsg/bXhFoAABinLyb5u1O/vzzJJ/Z//liS71j5FQEAKyGzBgBghGqt/7GU8sKpQ8dqrdcW\nO3g6yV2txzh+/PhNc/+HLD5/VA1ZF6G1lgmzveAFL1j3JRxqd93VfMszw/b29rov4VDb2dlZ9yXc\nlnTWAMAMFrhjZKbXpzmT5IlWgRt31tva2urceXGoZSyiOk98UfMu8ru5ufmchXEXWQT4dltg+AUv\neEEeffTRWy5vgeHnxu+66648+eSTB/b4retfxwLKy7K9vf2cxa0tMDyfnZ2dXLhw4Tm/L9vY6mBd\nTIMCADgcPltKeWD/5zck+a01XgsAsEQyawBgBiM7jMwPJ3mwlLKd5PNJPrLm6wGAA6f9NaGzBgBg\npGqt/zPJq/Z//sMkr13l+Q/TVId1WGQ6xKLTnFrTiBYtf+MUunnjXY9/8eLF6z+3rm/RaWKLTuNZ\ndJrQMqYZTU/laU3TWvbaSodtt57p+2XRaUzznGsZ5Y/65ypfp7MGAGbQIAIAWC3tr4nD1SUKAAAA\ncMTJrAGAGYzsAACslvbXxEo6axadc3pQ5zkoQ26eoXO6V/WcWvOGu7TmKndpbRc6PZf5VssMOc9B\nlWnN9+7SuheeeeaZzuND5gZvbW3NXebEiRNzlzl+fP6PklaZrvtrSB0s4/3aZch7dVWfjUOt6vxj\nrwcAAFg1mTUAMIORHQCA1dL+mrBmDQAAAMCIyKwBgBmM7AAArJb214TOGgCAI6przafpY8tuELfW\nnGqdf9lrVg1ZM2v6WKt835qArfUCW+viteKttQVb8dZafkMe/8KFC9d/bl3/kPUUpw1dL/Kazc3N\n3nhrHb1llN/d3b3lx2/FW2sHtsq37v3W81vks2fZawtaK4+x0FkDADMY2QEAWC3trwlr1gAAAACM\niMwaAJjByA4AwGppf03IrAEAAAAYEZ01AAAAACNiGhQAzCANFwBgtbS/JnTWAACwFuveInfI1tvT\nx1rbS/fFF92ae9Gttae3gR7y+K3yXfFz587d8uMvunV3y6Jbby+6NfaQ+PTW563yW1tbvfEh29ZP\nW+fW3q3H7io7XWbdnzstOiq4RmcNAMygwQQAsFraXxOj7awZ0uM5pMyQG2EZvbGzHnNV9dAavenS\nGhHp8uyzz84dnx5FuFXPPPPM3GWGnKc1atWlbyThxS9+cR599NHOWGuEpMvJkydXUmZ7e/tAy5w9\ne7bzXhhynlV92LdGiLosOqo1jyH1sKrPx9Z5fGEDAHC7GW1nDQCsm44iAIDV0v6asBsUAAAAwIjI\nrAGAGYzsAACslvbXhMwaAAAAgBGRWQMAMxjZAQBYLe2vCZ01AAAcSovuqteKX716tfdYV3xa326b\nrZ04W7tuLhrf3d3tjbd2yWzFL168eNOx6d06W9fX2nGzVfctm5ubS423drA8frz/v2Fd5c+dO3fL\nj9+6v1q7jLbqt3X9rXjrvdf3n/VFd+Bcxs6+05bd0bBI3Q19/GXXGd101gDADEZ2AABWS/trwpo1\nAAAAACMiswYAZjCyAwCwWtpfEzJrAAAAAEZEZg0AzGBkBwBgtbS/JmTWAAAAAIzIaDNrhvSmjXlL\nsaFbRw55TkO2Mmxtj9jl2WefnbvM9JaNtxp/+umn5z7P448/PneZJ598cu4yrW0vu/RtN/ia17wm\nf/qnf9oZ29nZmftcd95559xlTp8+PXeZU6dOzV3mzJkzvfGu+2vI+2FVnyWtLTy7tK5tzJ9psyxj\nJMToCgAAt5vRdtYAwLrpKIL1GjrYdavx1gBXV3z62JDy17QGylrxS5cu9cZbA0oXLlzojZ8/f36h\n8ufOnbvp2FNPPXX959ag38WLF3vjV65c6Y23XvvWIEsrvr293Rs/efJkb7w1CHfixImbjk0PanbF\np7XqpxVf5N5O2vV//Hj/f0P7Bjdbur67p6932d/tree+6GDdotd/GAYDtb8mTIMCAAAAGBGZNQAw\ng5EdAIDV0v6akFkDAAAAMCIyawBgBiM7AACrpf01IbMGAAAAYERk1gDADEZ2AABWS/trQmYNAAAA\nwIjIrAGAGYzsQL+9vb21Pn4rfvXq1QMvP33sypUrveUvX748KJYkly5d6o3v7u72xi9cuNAbP3/+\n/ELxp556au74Y489dv3nc+fO9ZZvXX+r7luvbevzfXNzszd+4sSJ3vipU6cWip8+ffqmY08++eT1\nn++4447e8q37a2dnpze+6HtnUVtbWwd67nnKrPu7f93nHwN1MCGzBgAAAGBEZNYAwAxjGdkppWwk\n+WCSlya5mOR7a62PTMXfmOTdSS4neajW+mApZTPJg0lKkr0kb6u1fm7lFw8AMIextL/WTWYNAIzf\ndyXZqbW+OsmPJHnftUApZSvJTyd5fZLXJnlrKeW+JG9MklrrtyV5Z5L3rPqiAQAYRmbNAMuYJzlr\nXmhrPm6X1hzVLhcvXpy7TGsuc5fpuba3Gn/00UfnPs+XvvSlucu0rq1Laz51l1ZP8Re/+MXO4ydP\nnpz7XF3znVvuvffeucs873nPm7tMay501+vRKnNQWvPQuwz5XGjNh+96vhsbq+tjX9WoRus86xxd\nGdHIzmuS/HqS1Fp/t5TyiqnYS5I8Umt9PElKKZ9Kcn+t9VdKKb+2/zffnOSJVV4wAMAQI2p/Dc1u\n3kryUJIXJjmR5Kdqrb9aSnlZkl9L8kf7xf9NrfU/zDq3zhoAGL87k0z3YF4ppRyvtV7uiD2d5K4k\nqbVeLqX8QpK/k+TvrepiAQCOiOvZzaWUV2WS3fym5DnZzX89ybkkny6l/GqS70zyWK31LaWUe5L8\n9yS/muTlSd5fa31fx3luYhoUAMxw7Nixlf1reCrJmanfN/Y7arpiZzKVRVNr/YdJXpzkwVJK//Yh\nAABrtsr21y20wZ6T3ZykM7u51rqb5FNJ7k/yK0nede3pZJJ1k0w6a/52KeWTpZQPlVKm22830VkD\nAOP36UxGabI/qvMHU7HPJ3lRKeWeUsp2Jo2E3ymlvKWU8qP7f3M+ydX9fwAA3JrO7OYZsaeT3FVr\nfabW+vR+Z8xHMlk7MEl+L8k/rbXen+SPk/xE34lvaRpUKeWVSf5VrfWBUsq3JvlwJjtLfC7J22ut\nGn8AsDwPJ3ldKeW3Mxmh+Z5SypuTnK61/lwp5YeSfDyTQZiHaq1/Xkr5T0l+vpTyySRbSX6w1vrs\nup4Ah9OQ9bhWef5Fr6+1FlpXfPpYa23BS5cuDYolye7ubm+8td7gs8/2v91b8dZafk880b8MVlf8\nscceu/7zuXPnesu3rq+1RmPrtW2Nph8/3v/fpNb6dq21Bu+4oz/RsWtdxK9+9avXf27dH611C1v1\ns+h7a9H16PrW6Buynsn082mVbz33Ma2nMsSQ6z/sz3lBg7KbSynflEn77YO11l/ajz9ca7324fhw\nkg/0nbjZWVNK+WdJ3pLJHKwkeX+Sd9Zaf7OU8rOZzNd6uPU4AHDYjKVxsj8o8rYbDn9hKv7RJB+9\nocy5JH9/+VcHAHBwxtL+2vfpTHbY/OW+7OYkz2SS3fze/V05fyPJ99da/8vU33+8lPIDtdbfS/I3\nk3ym78S3klnzxSR/N8kv7v/+8iSf2P/5Y5lsFaqzBgAAADhKhmQ3/0ySu5O8q5Rybe2aNyT5viQf\nKKVcSvIXSd7ad+JmZ02t9T+WUl44dehYrfVabtj1HSf6bG1t9aayDdkm96gZsi3zQTp79uzcZb7h\nG77hwK/j1a9+9YE/5mHz7ne/e92XMAove9nL1n0Jazdk6/WjaJ3fEavcKh0AgHG1vwZmN78jyTs6\nHu73k3zbrZ57yNbd0xMcn7PjxCx9c3JPnDjROed2SOpTa+5ll0XnPB5EmZMnT86cl9uaC92lNYe3\ny/nz5+cu89RTT81d5vHHH58Ze/WrX53f+Z3fuen4o48+Ovd5vvSlL81dpjU3u0vXfOKWvnvu3e9+\nd/7Fv/gXnbEhHXpD/rN/7733zl3mec973txl7r777pmxl73sZfnsZz970/EzZ3oXTO906tSpucsM\n6RxozW3vsrm5OTN2+vTpPPPMMzcdH/rlNeb5yX3nmfUdoZMfAICjbEir/7OllAf2f35Dkt86uMsB\ngPEY0bYvAgBIAAAgAElEQVSRAAC3hZFt3b02QzJrfjjJg/vbg34+k62oAAAAADgAt9RZU2v9n0le\ntf/zHyZ57RKvCQBGYcyjLQAAR5H218SQzBoAAI6AIevwrVLr+tYRnz7WWluwL963puOtxFvr5rXi\nrbUHn3iif1nKr33ta3M//vRjPv30073lW+sptuqn9dq0/jPYWotuZ2dnofis9Sqv6Vqv7atf/er1\nn3d3d3vLt+KttT4X/WxorbG3SHzR16713FqPv+zyLToybh8r6axZ9ItymQ56seBZWh+Is+JDFhge\nsujtkAWGhyzI+5WvfGXu+J/92Z/NfZ4vf/nLc5dpNUq6DKnr1v3zp3/6p53HW1/6Xe644465y7Qa\nT12G3D8tXYtRD/ly6lvE9yAtYwHfrs+FoV/QvtiHUW8AAKul/TUxnj2xAAAAADANCgBmMbIDALBa\n2l8TMmsAAAAARkRmDQDMYGQHAGC1tL8mZNYAAAAAjIjMGgCYwcgOAMBqaX9N6KwBAGAp9vb2jlz5\neR7z6tWrM2NXrlzpLXvx4sWF4ufOneuNP/30073xJ554ojf+5JNPzl3+8ccfv/7zU0891Vu+df27\nu7u98cuXL/fGW/8Z3N7e7o3v7Oz0xk+ePNkbb13/pUuXbjr2ta997frPrefXur9aNjb6J2Bsbm72\nxre2thYq3xdvlW29b3UEcFiYBgUAAAAwIjJrAGAGo28AAKul/TUhswYAAABgRGTWAMAMRnYAAFZL\n+2tCZg0AAADAiIw2s2bR1fuXeZ4hZfp2A+iLt1Z679LaHaDL+fPn5y7zzDPPzF3m0UcfnTveKtPl\nS1/60txlHnvssbnLtHZS6NK6f/7kT/6k8/ipU6fmPtedd945d5kh90/r/u7S2mXhy1/+8k3HWjsL\ndDlx4sTcZY4fn/+jcciuC0N2M1ilVY1qtM6zztEVIzsAAKul/TUhswYAAABgREabWQMA62ZkB/q1\nMgCXnSG46OMvWr6VVdmXedrKnt7d3e2NtzJhWxnQTz75ZG+8lT3cKt8Vnz42pPy0Z599tjfeqt+N\njf4x61Y27+nTp3vjd9xxR2/80qVLvfGue+uJJ564/vOi9+6iz7+Vhdwq34r3vbdaGd1ddTN9bNG6\na7UNWo+vbdGmjiZk1gAAAACMiMwaAJihNfIIy1JK2UryUJIXJjmR5KeS/I8kH06yl+RzSd5ea51/\n0TAAGDHtrwm1AAAwPv8gyWO11m9P8reS/Osk70/yzv1jx5K8aY3XBwAskcwaAJjBnGnW6FeSfGT/\n52NJLid5eZJP7B/7WJLXJ3m470GOHz9+0whlazc++p09e3bdl3BoffjDH173JRxqH/jAB9Z9CYeW\n9+1idnZ2Vno+7a8JnTUAACNTa30mSUopZzLptHlnkvfWWq+tXPl0krtaj3PjIqvb29vNhWvnMfYF\nhlvx1gLBN8bPnj37nEVeL1y40Fu+bxHg8+fP95ZtLfDbWkD4scce641/5StfWah8Kz5dT8mko+a7\nv/u7r/9+uy8w3Cp/5513Puf3D3zgA/mBH/iB67/fdVf/278Vv/vuu3vj995770Llz5w50xtvPf+T\nJ0/OjJ04caK37I0d0je+b1uvfaujYNGOhEUff9nXd6OdnZ3nfNatuuPmdmYaFADMcOzYsZX9gxuV\nUr4pyX9N8ou11l9KMr0+zZkkT3QWBIBDbJXtrzG3wXTWAACMTCnlviS/keSf11of2j/82VLKA/s/\nvyHJb63j2gCA5TMNCgBmGPNoC0fejyW5O8m7Sinv2j/2jiT/TyllO8nn8/U1bQZbdJrS2MsvGr96\n9ebNtqaPdcWn9U2zunTpUm/ZVrw1jercuXMLxRedhtVVfvrYvNOobtR6/q0pbi2ttZ1a529NA2pN\n0+q6N6enhrW+nzY3N3vjralErefXN00paU+VaT3/vvjx4/3/hW29b1c9jYj5eQ0mdNYAAIxMrfUd\nmXTO3Oi1q74WAGD1TIMCAAAAGJGVZNYMSUEdkhY7JF1qyHmGlGmlyc6KD0nhbKUVdmntZtCltUp/\nl1ZKa1f88ccfn/s8X/3qV+cu82d/9mdzl2k9ny6te+GLX/xi5/FTp07Nfa7nPe95c5dpXV+XVjpq\nl9YuAF3p0a2dB4acp0srNbhLa9eILqv6bEyGfT6uKgV1zOnI0nABAFZL+2tCZg0AAADAiFizBgBm\nMLIDALBa2l8TMmsAAAAARkRmDQDMYGQHAGC1tL8mdNYAABxRiy5UPnRR84MyZCH2g4x3Lbo/fay1\nKH/fRhGtDSEuXbrUG3/22Wd74+fOneuNnz9/fqmP37UJw/Sx1iYNrY0sWptjtOqv9Z/B7e3thR6/\ndW9sbPRPcOjauGD6NWttbNDaLOHkyZO98db90drAoVU/rfu/r/6GbNwyfaxV963PhTFvjHAQuq7/\nsD+nw0pnDQDMoHECALBa2l8T1qwBAAAAGBGZNQAwg5EdAIDV0v6akFkDAAAAMCIyawBgBiM7AACr\npf01IbMGAAAAYERk1gDADK3tPQEAOFjaXxMr6axp7VXfih/UeQ7KkPMMrYOrV6/Ofa7Lly/PXWZ3\nd3fuMufOnTvwMl3xZ555Zu7zPP7443OXeeyxx+Yu87WvfW3uMleuXOmNP/roo53HT506deDn6nLy\n5Mm5y9x5551zlxlyL5w/f37u8wy5t4fU26o+f8aeFjrk+lplxv6cYcy6PptW9Xl1K+dadrzVjhpS\nfvpY6/H7vk9a7bWLFy/2xi9cuLDU8s8+++xC5bu+s6ePDWkHzHP+S5cu9cZb/xlslW/dO63vrq2t\nrd74zs7OTcem66QrPq31+rfirfZTq35a741F4kM+F1b5udeyaLtn0XaRdtXhIbMGAGbQoAEAWC3t\nrwn5RQAAAAAjorMGAAAAYERMgwKAGaThAgCslvbXhMwaAAAAgBGRWQMAMxjZAQBYLe2vCZk1AAAA\nACMiswYAZjCyA+O2t7e3UPzq1atzx6ePtcpfuXJlUOxW4pcuXeqNX758eanxixcvzh2fPta6/t3d\n3YXO36q/RT/fW9d34sSJ3viy62/R17dVf8u+v/rOv+z37cbGevMZFr03l1F+1e0h7a8JmTUAAAAA\nIyKzBgBmMLIDALBa2l8TMmsAAAAARkRmDQDMYGQHAGC1tL8mdNaMXGsBrIMq01pErEtrYbUhZbri\nQ87z7LPPzl1myHkuXLgwd5lWXc96zCEfWkOub0jdtRaJG1KmKz7kPl3Ve6i1iOVBWdV5AADgdldK\n2UjywSQvTXIxyffWWh+Zir8xybuTXE7yUK31wVLKVpKHkrwwyYkkP1Vr/dVSyrcm+XCSvSSfS/L2\nWuvM/3iYBgUAMxw7dmxl/wAAWG376xbaYN+VZKfW+uokP5LkfdcC+50yP53k9Ulem+StpZT7kvyD\nJI/VWr89yd9K8q/3i7w/yTv3jx9L8qa+E+usAQAAALjZa5L8epLUWn83ySumYi9J8kit9fFa626S\nTyW5P8mvJHnX/t8cyyTrJklenuQT+z9/LMl39J3YNCgAmEHGC4dd19TJVU6nbJ1r3fHWtNeu+PSx\nIeVvtWxr6u8i515FvGs68/SxRZ/foq99S+v61nHvzvOcFn3+y67fRSy77hbVajss2rZYd/lVGNk1\n3pnkyanfr5RSjtdaL3fEnk5yV631mSQppZxJ8pEk79yPH6u17k3/bd+Jb6mzppTyyiT/qtb6QCnl\nZUl+Lckf7Yf/Ta31P9zK4wAAAAAcEk8lOTP1+8Z+R01X7EySJ5KklPJNSR5O8sFa6y/tx692/e0s\nzc6aUso/S/KWJOf2D708yftrre+bXQoADr+RjewAABx5I2t/fTrJG5P8cinlVUn+YCr2+SQvKqXc\nk+SZTKZAvXd/3ZrfSPL9tdb/MvX3ny2lPFBr/c0kb0jyX/tOfCuZNV9M8neT/OL+7y9PUkopb8ok\nu+YHa61P38LjAAAAABwWDyd5XSnltzNZf+Z7SilvTnK61vpzpZQfSvLxTNYDfqjW+uellJ9JcneS\nd5VSrq1d84YkP5zkwVLKdiYdPR/pO3Gzs6bW+h9LKS+cOvR7Sf5drfUzpZQfT/ITSf5J32Ps7Oxk\nY2P2WsanTp1qXcaRd/bs2bmO9/mGb/iGucu85CUvmbvMa1/72rnLtLznPe858Mc8bGzNPPHjP/7j\n676EtRvy/j+Kjh9f3/Jqfd9dAAAcvDG1v/a31n7bDYe/MBX/aJKP3lDmHUne0fFwf5jJrlG3ZEgL\n+OFa67W5VQ8n+UCrwIULF2bGTp06lfPnzw+4jPVpLXjWpWuRtWvOnj2bJ57onq42pG6efPLJ9h/d\n4NFHH527zCOPPNL+oxv88R//8czYe97zns7/oP/Jn/zJ3Of5oz/6o/Yf3eB//a//NXeZr371q3OX\n6Vuwbm9vb2ba38mTJ+c+17333jt3mb/8l//y3GVe/OIXz13mr/7Vvzoz9uM//uOdHXdDzvON3/iN\nc5e555575i5z5syZ9h/doK+jetbnwtCOiyHlhnxRDklb7Stz/Pjxzs/PdXbgAADAsg3psvp4KeVv\n7P/8N5N85gCvBwAAAOC2NmRo8vuSfKCUcinJXyR568FeEgCMw8gWuAMAOPK0vyZuqbOm1vo/k7xq\n/+ffT/JtS7wmAAAOAeusLU/rPyutqaqt8suOb25u9h5rlW9Nd110TYvW+be2tnrjXc9vnnjr+XXF\np4+1nn/r/Mt+/Q+zZX+uHeW642CZ9A8AM2hQAQCslvbXxHiWWQYAAABAZs2QNLchu0G1ysyK9+0c\nNMuqUpLH3OPZSv3sMmR3mRMnTsxdZnd3d9B1bG9vz32uVgpvlyHPaUjdDUn/HcPuROs29NrG/Jxa\nn1nrnGYx5noDADiKtL8mZNYAAAAAjMhtn1kDALMsuoAlAADz0f6aUAsAAAAAIyKzBgBmMGcaAGC1\ntL8mdNYAABxRXQt0r3LR7kUXEF/nAuMHoe8/HK3NEFoL8bc2EWjFhyz0P621IUFXfPrYzs5Ob/lL\nly71xlsuX768UPnWxg6t5996fkPi08da17fo/dG6P1vx1jSWVrzvvXPU/yO/7ue37u8Nvk5nDQDM\nsO4GEwDA7Ub7a8KaNQAAAAAjIrMGAGYwsgMAsFraXxMyawAAAABGRGYNAMxgZAcAYLW0vyZk1gAA\nAACMiM4aAAAAgBFZyTSoVhrTQaU5jXn/96tXrw6KD3lOQ+pzY2P+frudnZ25y5w8eXLu+B133DH3\nee655565y5w/f37uMkPq7dKlS73xe++9t/N4q+66PP/5z5+7zKzz9zlz5szcZVqva1d8SB1sbW3N\nXeb48fk/GofcC60yXfEh50mkkw41tL6BcVi0DdoVnz7WKr+5uTkz1vp8aX1/bW9v98ZPnDjRG2+1\n44a02aZ1tQ2mj+3u7vaWb7WBW/Vz5cqV3nhL32uXJKdOneqNnz59ujfeajt1tYOmj7Vev9brv2i8\nVf9D2jjzxG9nt0Obzus/oRYAAAAARsQCwwAww+0wegUAMCbaXxMyawAAAABGRGYNAMxgZAcAYLW0\nvyZk1gAAAACMiMwaAJjByA4AwGppf03IrAEAAAAYEZk1ADDDxoYxDQ63vb293mPrHr1c9vm7nv88\n5++KTx9rfUb0Pf7x4/3N8K2trd74zs7OQvGTJ08uFL/jjjt645cvX77p2NmzZ3vLTNve3u6NX7hw\noTe+u7vbG2/dG63Xp1W/red611139cZPnTrVe6wrPm3R17cVbz3/EydO9MZb9bu5uTkz1nrfdZWd\nPrbI+/Yg4rRpf02oBQAAAIARkVkDADMYHQMAWC3trwmZNQAAAAAjIrMGAGYYy8hOKWUjyQeTvDTJ\nxSTfW2t9ZCr+xiTvTnI5yUO11gdLKVtJHkrywiQnkvxUrfVXV33tAADzGEv7a9101qxIa5GkWfG+\nxbVmGVKmtQhYl9OnT89dprXYWlf8mWeemfs8XYvatQypt7vvvnvuMq1r+5Zv+ZbO462F5LoMub7n\nP//5c5e555575i7TWlivKz6kDobc261FHbu0FkLs0rrnuuJDv7x86R1635Vkp9b66lLKq5K8L8mb\nkmS/U+ank/z1JOeSfLqU8qtJvjPJY7XWt5RS7kny35PorAEAOAR01gDADCPq5HpNkl9Pklrr75ZS\nXjEVe0mSR2qtjydJKeVTSe5P8itJPrL/N8cyyboBABi1EbW/1kpnDQCM351Jnpz6/Uop5Xit9XJH\n7Okkd9Van0mSUsqZTDpt3rmqiwUAYDE6awBg/J5Kcmbq9439jpqu2JkkTyRJKeWbkjyc5IO11l9a\nxYUyLl2jk6scsWyda29vb6mPv+hzbdVfa5p733TX1lTY1pTcnZ2d3nhr6vAdd9zRG7948WJv/MqV\nK73xrtf2zjvvvP5z6/m3nt+FCxd6461p563rP368/79JrWnWZ86c6Y236r+r/PSx1nIErXjr/K34\nyZMne+Ot6eGteN/90XptuspOH1v258Zht+7vDb5OZw0AzDCixsmnk7wxyS/vr1nzB1Oxzyd50f66\nNM9kMgXqvaWU+5L8RpLvr7X+l1VfMIsppWwmeTBJSbKX5G1JLiT58P7vn0vy9lrr1XVdIwAsw4ja\nX2tl624AGL+Hk1wopfx2JosJ/1+llDeXUt5aa72U5IeSfDzJ72SyG9SfJ/mxJHcneVcp5Tf3//UP\nhTImb0ySWuu3ZTKF7T1J3p/knbXWb89kHaI3re/yAIBlklkDADO0pjisyn72xNtuOPyFqfhHk3z0\nhjLvSPKO5V8dy1Br/c+llF/b//WbM5na9h1JPrF/7GNJXp9JRx4AHBljaX+tm84aAIARqrVeLqX8\nQpK/k+TvJXldrfXaQiBPJ7mr9Rg7Ozs3NXpba5nQ77777lv3JRxaH/rQh9Z9CYfaz/zMz6z7Eg6t\ne+65Z92XcKi11hhiOXTWAMAM5kyzbrXWf1hK+edJ/luS6Wls1xeS7nPjIqynTp3K+fPnD/QaF9Fa\nYHjRBYivXu1f0qe1CO2lS5ee8/t9992XL3/5y9d/393d7S3ftwhu63VoxZ966qne+Ne+9rXe+GOP\nPbZQ+See6L/9nnzyyef8/qEPfSj/+B//4+u/nzt3rrd8K37UFxieXow5mXTUvOMd75gZv9HZs2d7\n43fffXdv/HnPe95C5Rd9/n3126r7GzsW7rnnnufcz63FrVtZHYvGF13geNVtk+3t7ed81q2i40b7\na0J+EQDAyJRS3lJK+dH9X88nuZrk/y2lPLB/7A1Jfmsd1wYALJ/MGgCYwcgOa/Sfkvx8KeWTSbaS\n/GAmO389WErZ3v/5I2u8PgBYCu2vCZ01AAAjU2s9l+Tvd4ReO8/jdE0jmj627gbxss/fevwh0xmm\nj7WmU/TFW1MJbpyCdaOTJ/s3d2tNM2lN4WpNI2pNUeuqu+mpM63nv7Oz0xtvXX9rmlNrilzr3tja\n2uqNt67/9OnTvfGuaUTTU59a04xa8dY0qtb1zTsV6UataWZ99TtkGtP0sbFPc1r24rqLTi9ldXTW\nAMAM6/6PLADA7Ub7a+JIddYMeVGHlBnSGzl05KXV69yl1dPfpdX736W1cFmXixcv9sa7dlhojex0\nGVJvrRGoLkMWaWyN9Hzrt35r5/Ehr2trVKVLa0G5Ln/pL/2lucs8//nPnzveGuXpMuTebo0WdRny\n+gz5XGiVOUirGnnxhQwAAM91pDprAOAgLTsVGQCA59L+mlALAAAAACMiswYAZjBFCwBgtbS/JmTW\nAAAAAIyIzBoAmMHIDgDAaml/TeisAQBgkGU3qFuP39q1rlW+FW/tytfa4bFvd8pW2dbOhK3yrR0U\nr1692htfdEfArrqb3km09fwuXLjQG7906VJvvLWbaOv5tRY4be082np+p06d6o137VJ65513Xv+5\ntePn9N92aZU/efJkb7y122br+W9vb/fG+3bZbO3A2Xrftl7bVrz1ubFo+UXji753V7UbKG2mQQEA\nAACMiMwaAJhBGi4AwGppf03IrAEAAAAYEZk1ADCDkR0AgNXS/pqQWQMAAAAwIjJrAGCG1o4OAAAc\nLO2viSPVWbOqdKmu7eCW9Zitbe+6DNlubUiZ1paPQ9x33303HRtS313bHbbcfffdc5c5f/783GVa\nW0m+8IUv7Dw+pB5a23Z2ueeee+Yu0/W6tbTquys+5Pm0tpbs0toSssuQ12fRrRkP8lzLPv8i5/GF\nDQDAOpRSNpJ8MMlLk1xM8r211kem4m9M8u4kl5M8VGt9cCr2yiT/qtb6wP7vL0vya0n+aP9P/k2t\n9T/MOveR6qwBgINkzjRHXWuwZtH3wLrfQ4t2BnfVz3SZ48f7m9J9A1utQa/Wa3PlypXB576Vx29p\n1d329vZNx6YHhJ599tne8ru7uwvFW/XT0rp3WoM0Xc9/WmswqWvg8ezZs9d/PnPmTG/5u+66a+7H\nn3by5MneeOv6W8+/NTDW995q1X1XfPpY695ddDBt0cdndHX0XUl2aq2vLqW8Ksn7krwpSUopW0l+\nOslfT3IuyadLKb9aa/1yKeWfJXnL/vFrXp7k/bXW993KiQ1XAgAAANzsNUl+PUlqrb+b5BVTsZck\neaTW+nitdTfJp5Lcvx/7YpK/e8NjvTzJ3y6lfLKU8qFSSm+vq84aAJjh2LFjK/sHAMBq21+30Aa7\nM8mTU79fKaUcnxF7OsldSVJr/Y9JLt3wWL+X5J/WWu9P8sdJfqLvxL25m/tpPQ8leWGSE0l+Ksn/\nSPLhJHtJPpfk7bXWg1+8BAAAAGB9nkoynQGzUWu9PCN2JskTPY/1cK31WvzhJB/oO3Ers+YfJHms\n1vrtSf5Wkn+d5P1J3rl/7Fj252sBwFEzolEdAIDbwsgyaz6d5DuTZH/Nmj+Yin0+yYtKKfeUUrYz\nmQL1Oz2P9fFSyt/Y//lvJvlM34lbCwz/SpKP7P98LJMVjl+e5BP7xz6W5PWZ9ArNdOLEid6FlloL\nWN0OhuxycxTdf//97T864n7yJ39y3ZcwCq985SvXfQlrN2Q3uKNoGTvwAQDALXg4yetKKb+dSZ/I\n95RS3pzkdK3150opP5Tk45kkwjxUa/3znsf6viQfKKVcSvIXSd7ad+Lezppa6zNJsr/wzUeSvDPJ\ne2ut15aPvz4nq8/Fixdnxk6ePNlcDf5WrWpkcsjq+X0r0p8+fTrPPPNMZ6y1zXOXCxcuzF1myBbU\nTz311IGWuf/++/PJT37ypuNf/epX5z7PY489NneZIc/noLfu/smf/Mn8xE90T108alt3953nla98\nZf7bf/tvNx1v7WzQZUgdDNn6fUjHSt9reuLEic7PzqHbWI85c6Pv2jY3Nzs/P1fVgWPbcACA1RpT\n+2t/yZe33XD4C1Pxjyb56Iyy/zPJq6Z+//0k33ar527WQinlm5L81yS/WGv9pSTT69O05mQBAAAA\nMIfWAsP3JfmNJN9fa/0v+4c/W0p5oNb6m0nekElHDgAcOWPOSIJb0ZUNPH3sqN/jrec3ZPR2nsy+\nIdnY1yz62iz63FvP8/jx/tUUujJOn/e8513/uZUJ3peZnyS7u7u98atXF9v/pFV/rfrZ2trqje/s\n7PTGu5aJuPvuu6//3MocbsVPnTo19/mntTKKt7e3e+Ot+umr39a92xWfPrboe6NVftF4yyKfK4fF\nUf9uulWtNWt+LMndSd5VSnnX/rF3JPl/9hfQ+Xy+vqYNAAAAAAtqrVnzjkw6Z2702uVcDgAAAMDt\nbTwr9wAAAADQnAa1NkPmqa2qzJB5gq3ztOb9zmNV2/0OqbvW/NTpuczXtOa8drnzzjvnLjNkN6gh\nO5m1dvj6lm/5ls7jQ+6RIbsanTlzZiVlWq9R185PQ55Pa054lyF1Pfa5tcv43FpVmbHXLQAAHLTR\ndtYAwLrpKAIAWC3trwnToAAAAABGRGYNAMxgZAcAYLW0vyZ01gAAHFFd61XNs4bV2BvMy76+jY2b\nk9C7ji3Dos+tVX5zc3OheGt9t5MnT9507N57773+88WLF3vLX7hwoTd+5cqV3nhrjcCWReuvtU5j\na03GrjUo77777us/t9bk66r/eeKt62vFW89/kfuv9R7seu2mjy26Vt6y40PWGeRo0lkDADOM/T+q\nAABHjfbXhDVrAAAAAEZEZg0AzGBkBwBgtbS/JmTWAAAAAIyIzBoAmMHIDgDAaml/TcisAQAAABgR\nmTUAMIORHQCA1dL+mtBZAwBwRHU1eKePtRrEi8Zb1t0gb51/b2/vpmMbGweTmN46d+s8i742rcff\n3NzsjR8/3v/fiO3t7ZuOnTlz5vrPOzs7veVPnz7dG798+XJv/OrVq73xRbXqr1U/Q+p3uk5OnDjR\nW75Vv12vzzzxra2t3viiz7+vfofc+6v83Fv35xpHx0o6axZ9Qx3UeQ7KkC/pri/7W3nM1gdhlyH1\nMOQ5tT5ku7Q+uO+8886bjg2pg9YXVJdWo6DL7u7u3GVa/spf+Sudx4e8Rq0v2oMqc+rUqbnLnDx5\nsjd+xx13zF2my5DnM6Sux/7FvKrPhaNm7K8rAMBRo/01oSUOAAAAMCI6awAAAABGRGcNAAAAwIjo\nrAEAAAAYEbtBAcAMFrgDAFgt7a8JmTUAAAAAIyKzBgBmMLLDYdd1D08fW/Y9vujjj/E9OH1Nm5ub\nt/y388QOIr6x0T8m27r2Vvz48f7/Rpw4ceKmY2fOnLn+86VLl3rLX7lyZaH41atXe+N7e3u98ZZ1\n1O/p06ev/7y9vd1bfmtra+7Hn7bo47fqpxXvu7+HlF3l517Lovfe7WDdr9FYyKwBAAAAGBGZNQAw\ng5EdAIDV0v6akFkDAAAAMCIyawBgBiM7AACrpf01IbMGAAAAYERk1qzIoiuiz2PICuNDyrRWie/S\n6iU9efLkTceG1E1rhftbPXdLa6eBLq06eMELXtB5vLVrwDrLDKnvrl0ipnW9HkPuuSHP5yDfj7eb\ng1OmtyAAABs5SURBVP782djY6Hyfreo1MrIDALBa2l8TOmsAAI6oro7NeTo7l7299Ngb5F3XN11/\nrQ7qvrpe99bdi2493Rqo6dpa+9SpU73xVvlprQGzZW/dvY76nx7IatV/K77o1t6Lbk2+SP0MeW/M\ns3X3Yf/c4ujQWQMAM2iQAQCslvbXhFx/AAAAgBGRWQMAMxjZAQBYLe2vCZk1AAAAACOiswYAAABg\nRHTWAAAAAIyIzhoAAACAEbHAMADMYIE7Druue3ie+7r1t8uOL2oZj7+x8fWxzr29vd6/7Yu3yq67\n7jc3N3vjx4/3/zfiypUrNx3b2dnpjU+7evXqQvFFXptbiS9av9P3UZeu+p+uv0Vfn1b5ReOt59eK\n99XfkLr3fX7rxlB/Xq8JmTUAAAAAIyKzBgBmMLIDALBa2l8TMmsAAAAARuRIZda05pYelCE9fYvO\n25xHaw7pQRlSD60y29vbNx0bUjdD6qA1d7rLkHuuVQdnzpzpPD7kOS3jXu0y5Npac6m77oUh51nV\n+2GIRefLz8MIxTDqDQBgtbS/JmTWAAAAAIzIkcqsAYCDZGQHAGC1tL8mZNYAAAAAjIjMGgCYwcgO\nh13XOmQHuU5e6z2yaHyM5rnmvr9trUt29erV3njrdVy07hc9f9fadCdOnLjlx2/FF62/Rcsv+97v\nqt+tra3rP7fW5GvFW6/fouWXeX+u+3NjVeukzrLs57/u55es/zUeC5k1AAAAACMiswYAZjCyw7qV\nUl6Q5DNJXpfkcpIPJ9lL8rkkb6+19g//A8Aho/01IbMGAGCESilbSf5tkmf3D70/yTtrrd+e5FiS\nN63r2gCA5ZJZAwAzGNlhzd6b5GeT/Oj+7y9P8on9nz+W5PVJHu57gOPHj9+0dsT29vbBXuVtprWW\nB7OdPXt23ZdwqD3/+c9f9yUcWt63i5leL2kVtL8mdNYAAIxMKeW7k3yl1vrxUsq1zppjtdZrKz8+\nneSu1uNcvnz5Ob9vb29nd3f3wK7zdltgeHNzM1euXLn++yILcS57gdxFH/+gFwA+e/ZsnnjiiaU9\n/qrLr3qB4ec///n5yle+cv13Cwzfuhvfty2H7XPpRot+NtwY39rayqVLl57zO6uhswYAYHz+UZK9\nUsp3JPlrSf59khdMxc8keaKrIABw+OmsAYAZDvvoGodXrfX+az+XUn4zyduS/N+llAdqrf9/e3cb\nI8l1Fnr8P+OdYXbt9Rp/wAKByIeYR/4AvsGG2Ipjr5TY9zpoZRMhPkQEri1fyxChRA4kJjiOhEAo\nUmxL+OJwvc5iXqW8sQhH8suHGyfYwUQ4QXKE/TgbQEKBSBB51y9re1+m+VC9pJnt6pqp7q46M/3/\nSSNN95lTderU6anTT51z6gngeuCL/ZROkqT5sf9VMVgjSZK0PXwQOBgRq8BzwOem3eC0HeJ5T3Mq\nscPeVZmappFMO9Vhq9NwNmqaJjTO6PSJvqd5Tbv9abWp/9H6m/b8zXua07TToKYxbtuj7837s9O3\neZSv9GPeqYoN1rT5B9lVI+qysbbZV9M/xz41Hc+uXWc3ya7qoE2no42msp177rkz21fJn4k2c61L\nbtvz6NSN22ZX7RTa1Xebeph23QJpp8vM/SMvr+mrHJIkqTvFBmskSZIkSZL6EhHLwP3ApcAbwC2Z\neWQk/QBwF3AKOJSZB0fS3gp8/MxNl4h4M/AQMAC+AbwvM2vvSpZ7m1qSpJ4tLS119iNJkqRu+1+b\n6IPdCKxl5pXAHcDdZxIiYgW4F7iOauTrrRFx0TDtQ8CDwNrItu4B7szMtwNLwA2TdmywRpIkSZIk\n6WxXAY8CZObTwOUjaZcARzLzxcw8ATwJnHlAwLeAd2/Y1mXAl4a/PwK8c9KOJ06DGkaKDgFvAr4P\n+G3gX4AvAN8c/tknM/PTk7YjSdJ25IgXSZKkbhXW/zofODby+nRE7MrMU2PSXgb2AWTm5yPiTRu2\ntZSZg41/W6dpzZpfAL6bme+NiAuBvwd+C7gnM++enFWSJEmSJGnbegnYO/J6eRioGZe2Fzg6YVuj\n69M0/W1jsOazfO+xkEtUi+ZcBkRE3EA1uuYDmflyw3YkSdp2CruzI0mStOMV1v96CjgAfCYirgCe\nHUl7Drh4OLDlFaopUJ+YsK2vR8T+zHwCuB744qQdTwzWZOYrABGxlypocyfVdKgHM/OZiPhN4GPA\nr03azurq6sRHwK6trdWmLYpZPpK46ZHIs8qzsrKy5TxN9uzZM/NtbjcXXHBB30UowiwfYb5d7d69\nu+8iFGF1dbXvIkjb1rgO7+h7fXeIp93/vMs/GAwmpjftvyn/NNueNn9T2aY99nFG+5tNfd/19doH\npADN5Zv39qfVVH/j0kePqen4pk1vU76tpM/TuHM37/M5S9PW3Tw+u31fK3p2GLg2Ir5CNYDlpoh4\nD3BeZj4QEbcDj1GtB3woM789YVsfBA5GxCpVoOdzE/62+dHdEfEjwwLen5l/HhEXZOaZ4TqHgfua\ntnHixInatLW1NV5//fWmTWxKV42oTWBj0j/E5eXl2gtGm2NquvjMKs/p06dnup89e/Zw/Pjxmexn\n1mWbpUlt4YILLuDo0Ymj4bakq8/EtB22jc4991xeffXVs95vE9Tsqg5m3QnYvXs3r7322lnvtw3s\ntsnXJk+bepj02VtdXR17DekqgLPgnRNJkqTOldT/Gj5a+7YNbz8/kv4w8HBN3n8Grhh5/QLVU6M2\nZWJPfPjYqceBD2fmoeHbj0XETw9/fwfwzGZ3JkmSJEmSpMmaRtZ8BPh+4KMR8dHhe7cD90bESeA7\nwK1zLJ8kSb0p6c6OJEnSIrD/VWlas+b9wPvHJL1tPsWRJEmSJElabI1r1kiStKi8syNJktQt+1+V\nHRWsaXNSu1p0s+2K6V0dU1ergjct/Dtu0dmu6qCUVeJ37Rr/seyqfKUsStz0BJNZ7WecUtrCOG3L\nNuuFf9WtiFgG7gcuBd4AbsnMIyPpB4C7gFNUTyI4OJL2VuDjmbm/00JLkiSptdk9L1qSJM3LjcBa\nZl4J3AHcfSYhIlaAe4HrqJ4wcOvwAQFExIeAB4G1zkssSZKk1nbUyBpJknaoq4BHATLz6Yi4fCTt\nEuBIZr4IEBFPAlcDnwW+Bbwb+JNui6tSTDtKcdGHojfVX9OoxXnW37Tbbir7tOnjjI58nrbumkZ/\nTjtKtu9Rtk1tr+2o/a7S+/xslG7exz6P7S/y+eqTI2skSSrf+cCxkdenI2JXTdrLwD6AzPw8cLKT\nEkqSJGlmHFkjSVKNgu4kvQTsHXm9nJmnatL2Ake7KpgkSdIsFdT/6pUjayRJKt9TwLsAIuIK4NmR\ntOeAiyPiwohYpZoC9TfdF1GSJEmz4sgaSZJqFHRn5zBwbUR8BVgCboqI9wDnZeYDEXE78BjVTZhD\nmfntHssqSZLUWkH9r14ZrJEkqXCZuQ7ctuHt50fSHwYersn7z8AVcyucJEmSZs5pUJIkSZIkSQUx\nWCNJkiRJklQQp0FJklTDOdPa7sa1Ydv15g0Gg//2emlp6az32ur7PDTtv+k429TD8vL37hOvr69v\nOX/dtsZp2v609T9tO2jaf9Nnt03+WWo6/mnb16yN7q/vz960Sqvbedju52hWHFkjSZIkSZJUEEfW\nSJJUwzs7kiRJ3bL/VXFkjSRJkiRJUkGKHVnTJprWNHd1VvvpMtLX1ZzDruqh6Ry1OYfjtClbm7qe\nx/mZVR3A9m4/s/qclXJex+lyvnnJ9dC0n2nXNZiGd3YkSZK6Zf+r4sgaSZIkSZKkghiskSRJkiRJ\nKojBGkmSJEmSpIIUu2aNJEmSyua6Au2VXndN5Wtab2xc+ug2pz3+pv13uS5cG23K13eZt2Lade9K\nPtZp216TeedvU/6Sz8dOZrBGkqQadk4kSZK6Zf+r4jQoSZIkSZKkgjiyRpKkGt7ZkSRJ6pb9r4oj\nayRJkiRJkgpisEaSJEmSJKkgBmskSZIkSZIK4po1kiTVcM60JElSt+x/VQzWSJIkaVsaDAZ9F0Et\nNX0Z89xONu2X2Z1cv9O2LQMFKkUnwZqmBj8uvc2HpE2e5eWtzwRrs58u/yGW/M+3TVtoc47a1MH6\n+vqW87TR5fnZaReb7dy2x2lq223a/ix19Zko2U77DEmSJJXO/lfFNWskSZIkSZIK4jQoSZJqeGdH\nkiSpW/a/Ko6skSRJkiRJKojBGkmSJEmSpIIYrJEkSZIkSSqIa9ZIklTDOdOSJEndsv9VMVgjSZK0\nQ43r8G6lEzwYDGZZHOm/lN625v1lsWn7TZ/daeuvaf9N229Kb3N8szKubKPvTXvsBhLUFYM1kiTV\nsEMmSZLULftfFdeskSRJkiRJKojBGkmSJEmSpII4DUqSpBoOw5UkSeqW/a+KI2skSZIkSZIKUuzI\nmpKjaV2uXt/VvtrsZx55+nwyQFdtru0xtsnX1TGV/Hlto83xzCPPLOu15PbTpJRySJIkSV1xZI0k\nSZIkSVJBih1ZI0lS3xzVI02nzxGzUPZnuO+6KV3p9VNy29qMpvrd7sc3ybTHPm3bnHb7056beR/f\nLOzk9rcVjqyRJEmSJEkqiCNrJEmq4Z0dSZKkbtn/qjiyRpIkSZIkqSAGayRJkiRJkgpisEaSJEmS\nJKkgrlkjSVIN50xLkiR1y/5XxZE1kiRJkiRJBXFkjSRJNbyzo0U37WegKf9gMJhq+6Wb5via8pb+\n/2lc+Uff6/vc911/bfa/lTx9H1+Trs9/3+1t1Pr6+sT0eZ+70tsGbI8ydsGRNZIkSZIkSQVxZI0k\nSTW8syNJktQt+1+VToI1TcPOmoZJlmYejWeW22xTd6XkmdV576r9lNxOoV35/OfYTun1Vnr5JtnO\nZZckSZLacGSNJEmSJEnSBhGxDNwPXAq8AdySmUdG0g8AdwGngEOZebAuT0S8BfgC8M1h9k9m5qfr\n9m2wRpIkSZIk6Ww3AmuZeWVEXAHcDdwAEBErwL3ATwGvAk9FxF8Bb6vJcxlwT2bevZkdu8CwJEmS\nJEnS2a4CHgXIzKeBy0fSLgGOZOaLmXkCeBK4ekKey4CfiYgvR8SnImLvpB03jqyJiHOAg0AAA+A2\n4HXgoeHrbwDvy8zJzyCTJGmbcb0cSZKkbhXW/zofODby+nRE7MrMU2PSXgb21eUBvgo8mJnPRMRv\nAh8Dfq1ux5uZBnUAIDPfFhH7gd8BloA7M/OJiPgDqiE9hzexLUmSJKkT83gQQCkPQ5g2/7RfhqZ9\ngMi0D53o+/jnbVz5SirztGVpyj8pfd4P+Cj9ASJNmuq26fhKb3s9eAkYHQGzPAzUjEvbCxytyxMR\nhzPz6PC9w8B9k3bcGKzJzL+MiC8MX/7ocOfvBL40fO8R4DomBGtWV1dZXq6fcbW2ttZUjIXV5oPR\nJs+k89Ol1dXVvovQuz179vRdhCLs3r277yL0bmVlpe8iFKHP/wsL3jlRzyLia1QdPoB/orph9hCO\nbJYk7WCF9b+eohrA8pnh+jPPjqQ9B1wcERcCr1BNgfoE1XV6XJ7HIuJXM/OrwDuAZybteFMLDA+j\nQH8E/Czwc8C1mXkmJHdmqE+tEydO1Katra3x+uuvn/V+m+DBOeecs+U8bfbTZePp6pHa6+tb7+vN\numyrq6sT28qs9jNLs663PXv2cPz48S3nm6WuAoST7N69m9dee22m29yKrgKek/azsrLCyZMnt7zN\nWSrhzlLd/wUDu9rpImINWMrM/SPv/RWObJYkqUuHgWsj4itUM4xuioj3AOdl5gMRcTvwGNV6wIcy\n89sRcVae4bZ+GbgvIk4C3wFunbTjTT8NKjN/KSI+DPwtMHrL+8xQH0mSJM3GpcCeiHicqr/2EaqF\nCTc9slmSJE1nOIL1tg1vPz+S/jDw8CbykJlfo3pS1KZsZoHh9wI/nJm/CxwH1oG/i4j9mfkEcD3w\nxc3uUJIkSY2OUw2lfhC4mCo4s7SVkc1QjeDbOIqvlKnP21WbkdyqOL13Oi4d0Z4jcrUdbWZkzV8A\nfxgRXwZWgA9Qzc06GBGrw98/N78iSpLUj8LmTGuxvED1ONAB8EJEfJdqZM0ZmxrZPBgM/tuUxuXl\n5VZTeOelhOmWW3HOOedw+vTpmWyr72PveoHhjdN7m9rhvBcYLt3G81O3dERfttMCw7NcamEzpjm2\nLtKblND3KaEMJdjMAsOvAj8/Juma2RdHkiRJwM3AjwO/EhE/RPUY0Mcd2SxJ0mLY9Jo1kiQtGu/s\nqEefAh6KiCepnipxM/AfOLJZkrTD2f+qLO30IYSSJLV17Nixzi6S+/bts2ciSZIWXpf9Lyi3D+YK\nc5IkSZIkSQUxWCNJkiRJklQQ16yRJKmGc6YlSZK6Zf+r4sgaSZIkSZKkghiskSRJkiRJKojBGkmS\nJEmSpIIYrJEkSZIkSSpIbwsMR8QycD9wKfAGcEtmHumrPH2JiK8BLw1f/lNm3tRneboUEW8FPp6Z\n+yPizcBDwAD4BvC+zFzvs3xd2VAPbwG+AHxzmPzJzPx0f6Wbr4hYAQ4BbwK+D/ht4B9YsLZQUw//\nwgK1BYCIOAc4CATV+b8NeJ0e24ML3GmnsN/Vjn2Vdry+T6fE6+F2ExE/ADwDXAucwrrbtI3fT4Hf\noeP6s/9V6XNkzY3AWmZeCdwB3N1jWXoREWvAUmbuH/4sUqDmQ8CDwNrwrXuAOzPz7cAScENfZevS\nmHq4DLhnpE3s6C/nwC8A3x2e9/8F/F8Wsy2Mq4dFawsABwAy823AnVSdg0VsD9I8LHy/a6vsq0zF\n6/t0vB5OYRgs/H/Aa8O3rLtNqvl+av31pM9gzVXAowCZ+TRweY9l6culwJ6IeDwi/n9EXNF3gTr0\nLeDdI68vA740/P0R4J2dl6gf4+rhZyLiyxHxqYjY21O5uvJZ4KPD35eo7nwsYluoq4dFagtk5l8C\ntw5f/ihwlMVsD9I82O/aOvsq7Xl9n4LXw6l9AvgD4F+Hr627zRv3/dT660mfwZrzgWMjr09HRG/T\nsnpynOqfyf+kGt74Z4tSB5n5eeDkyFtLmTkY/v4ysK/7UnVvTD18Ffj1zLwa+EfgY70UrCOZ+Upm\nvjwMRHyO6u7RwrWFmnpYqLZwRmaeiog/Au4D/owFbA/SnNjv2iL7Ku15fZ+e18N2IuJ/A/+emY+N\nvG3dbd5Z30+x/nrTZ7DmJWD0TvFyZp7qqzA9eQH408wcZOYLwHeBH+y5TH0Znfe4l+oOwiI6nJnP\nnPkdeEufhelCRPwI8EXgTzLzz1nQtjCmHhauLZyRmb8E/BjVfP3dI0mdt4elpaXOfqQ5s981vYW8\nPrXl9X16JV0Pt5GbgWsj4gngfwB/DPzASLp1N9m476cXjaR3Un9d9r9K7oP1Gax5CngXwHB41bM9\nlqUvNzOcMx4RP0R11+vfei1Rf74eEfuHv18P/HWPZenTYxHx08Pf30G1MNqOFREXAY8DH87MQ8O3\nF64t1NTDQrUFgIh4b0T8xvDlcaqO/d8tWnuQ5sR+1/QW7vrUltf36Xg9bC8zr87MazJzP/D3wC8C\nj1h3mzbu++nj1l8/+hz+epgq6vkVqrmsC7O47ohPAQ9FxJNUq2vfvMB3uT4IHIyIVeA5qiGzi+iX\ngfsi4iTwHb43X3mn+gjw/cBHI+LM3Pb3A7+3YG1hXD3cDty7QG0B4C+AP4yILwMrwAeo2oD/G6Tp\n2e+ann2VzfP6Ph2vh7PlZ3fzzvp+CvwH1l8vlgaDQfNfSZK0gF555ZXOLpLnnXdeueNwJUmSOtJl\n/wvK7YO5sJwkSTVKnscsSZK0E9n/qvS5Zo0kSZIkSZI2MFgjSZIkSZJUEIM1kiRJkiRJBXHNGkmS\najhnWpIkqVv2vyqOrJEkSZIkSSqIwRpJkiRJkqSCGKyRJEmSJEkqiMEaSZIkSZKkgrjAsCRJNVzg\nTpIkqVv2vyqOrJEkSZIkSSqIwRpJkiRJkqSCGKyRJEmSJEkqiGvWSJJUwznTkiRJ3bL/VXFkjSRJ\nkiRJUkEM1kiSJEmSJBXEYI0kSZIkSVJBDNZIkiRJkiQVxGCNJEmSJElSQXwalCRJNXwagSRJUrfs\nf1UcWSNJkiRJklQQgzWSJEmSJEkFMVgjSZIkSZJUENeskSSpcBGxDNwPXAq8AdySmUdG0g8AdwGn\ngEOZebApjyRJksrlyBpJkmosLS119tPgRmAtM68E7gDuPpMQESvAvcB1wDXArRFx0aQ8kiRJpeqy\n/1XyYsYGayRJKt9VwKMAmfk0cPlI2iXAkcx8MTNPAE8CVzfkkSRJUsEM1kiSVL7zgWMjr09HxK6a\ntJeBfQ15JEmSVDCDNZIkle8lYO/I6+XMPFWTthc42pBHkiRJBfMOmyRJNQqax/wUcAD4TERcATw7\nkvYccHFEXAi8QjUF6hPAYEIeSZKkIhXU/+rV0mAw6LsMkiQV6eTJk51dJFdWVmp7JiNPdvoJYAm4\nCfhJ4LzMfGDkaVDLVE+D+v1xeTLz+TkfhiRJ0lS67H/B5D5YnwzWSJJUo5RgjSRJ0qIwWFNxzRpJ\nkiRJkqSCuGaNJEk1nDMtSZLULftfFUfWSJIkSZIkFcRgjSRJkiRJUkEM1kiSJEmSJBXENWskSarh\nnGlJkqRu2f+qGKyRJEmSJEnaICKWgfuBS4E3gFsy88hI+gHgLuAUcCgzD9bliYg3Aw8BA+AbwPsy\nc71u306DkiRJkiRJOtuNwFpmXgncAdx9JiEiVoB7geuAa4BbI+KiCXnuAe7MzLcDS8ANk3ZssEaS\nJEmSJOlsVwGPAmTm08DlI2mXAEcy88XMPAE8CVw9Ic9lwJeGvz8CvHPSjg3WSJIkSZIkne184NjI\n69MRsasm7WVg34Q8S5k52PC3tVyzRpKkGi5wJ0mS1K3C+l8vAXtHXi9n5qmatL3A0bo8EbE+5m9r\nObJGkiRJkiTpbE8B7wKIiCuAZ0fSngMujogLI2KVagrU30zI8/WI2D/8/XrgryfteGkwGExKlyRp\nYa2vr3d2kVxeXi7qNpIkSVIfuux/weQ+2MiTnX6CalHgm4CfBM7LzAdGnga1TPU0qN8flyczn4+I\nHwMOAqtUgZ7/k5mn6/ZtsEaSpBoGayRJkrpVUrCmT65ZI0lSjcLmTEuSJO149r8qrlkjSZIkSZJU\nEIM1kiRJkiRJBTFYI0mSJEmSVBDXrJEkqYZzpiVJkrpl/6viyBpJkiRJkqSCGKyRJEmSJEkqiMEa\nSZIkSZKkgrhmjSRJ9Zw0LUmS1C37XziyRpIkSZIkqSgGayRJkiRJkgpisEaSJEmSJKkgBmskSZIk\nSZIKYrBGkiRJkiSpIAZrJEmSJEmSCvKfnaIMg2FpcqkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110206550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pylab as plt\n",
    "plt.figure(figsize=(20,10))\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(pacs100_psf[1].data[centre100-radius100:centre100+radius100+1,centre100-radius100:centre100+radius100+1])\n",
    "plt.colorbar()\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(pacs160_psf[1].data[centre160-radius160:centre160+radius160+1,centre160-radius160:centre160+radius160+1])\n",
    "plt.colorbar()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Set XID+ prior class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#---prior100--------\n",
    "prior100=xidplus.prior(im100,nim100,im100phdu,im100hdu, moc=Final)#Initialise with map, uncertianty map, wcs info and primary header\n",
    "prior100.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_ELAIS-N1-SERVS.fits',ID=XID_MIPS['help_id'][good])#Set input catalogue\n",
    "prior100.prior_bkg(0.0,5)#Set prior on background (assumes Gaussian pdf with mu and sigma)\n",
    "\n",
    "#---prior160--------\n",
    "prior160=xidplus.prior(im160,nim160,im160phdu,im160hdu, moc=Final)\n",
    "prior160.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_ELAIS-N1-SERVS.fits',ID=XID_MIPS['help_id'][good])\n",
    "prior160.prior_bkg(0.0,5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Divide by 1000 so that units are mJy\n",
    "prior100.set_prf(pacs100_psf[1].data[centre100-radius100:centre100+radius100+1,centre100-radius100:centre100+radius100+1]/1000.0,\n",
    "                pind100,pind100)\n",
    "prior160.set_prf(pacs160_psf[1].data[centre160-radius160:centre160+radius160+1,centre160-radius160:centre160+radius160+1]/1000.0,\n",
    "                pind160,pind160)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 3362 tiles required for input catalogue and 9 large tiles\n"
     ]
    },
    {
     "ename": "SystemExit",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mSystemExit\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "#from moc, get healpix pixels at a given order\n",
    "from xidplus import moc_routines\n",
    "order=11\n",
    "tiles=moc_routines.get_HEALPix_pixels(order,prior100.sra,prior100.sdec,unique=True)\n",
    "order_large=6\n",
    "tiles_large=moc_routines.get_HEALPix_pixels(order_large,prior100.sra,prior100.sdec,unique=True)\n",
    "print('----- There are '+str(len(tiles))+' tiles required for input catalogue and '+str(len(tiles_large))+' large tiles')\n",
    "output_folder='./'\n",
    "outfile=output_folder+'Master_prior.pkl'\n",
    "with open(outfile, 'wb') as f:\n",
    "    pickle.dump({'priors':[prior100,prior160],'tiles':tiles,'order':order,'version':xidplus.io.git_version()},f)\n",
    "outfile=output_folder+'Tiles.pkl'\n",
    "with open(outfile, 'wb') as f:\n",
    "    pickle.dump({'tiles':tiles,'order':order,'tiles_large':tiles_large,'order_large':order_large,'version':xidplus.io.git_version()},f)\n",
    "raise SystemExit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}