{
 "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": 8,
   "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": 2,
   "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",
    "SWIRE_MOC=pymoc.MOC()\n",
    "SWIRE_MOC.read('../data/ELAIS_N1/MOCs/DF-SWIRE_ELAIS-N1_MOC.fits')\n",
    "Final=Sel_func.intersection(SERVS_MOC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Final=Sel_func.intersection(SWIRE_MOC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Final=Final-SERVS_MOC\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in XID+MIPS catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "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_SWIRE_cat_20170725.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<Table length=10>\n",
       "<table id=\"table4687292624\" 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><th>flag_mips24</th><th>flag_mips_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><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><th>bool</th><th>bool</th></tr></thead>\n",
       "<tr><td>HELP_J160000.316+542000.690</td><td>240.001318675</td><td>54.3335249815</td><td>1928.47</td><td>1943.4</td><td>1913.42</td><td>-0.000955198</td><td>4.74435e-06</td><td>nan</td><td>1196.0</td><td>1.0</td><td>True</td><td>True</td></tr>\n",
       "<tr><td>HELP_J161145.788+525236.929</td><td>242.94078253</td><td>52.8769246635</td><td>244.531</td><td>531.299</td><td>80.7009</td><td>0.011866</td><td>5.16518e-06</td><td>1.00054</td><td>1760.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161138.168+525244.570</td><td>242.909032282</td><td>52.8790472905</td><td>239.355</td><td>531.635</td><td>65.6535</td><td>0.011866</td><td>5.16518e-06</td><td>1.00149</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161143.750+525251.314</td><td>242.93229323</td><td>52.8809206325</td><td>242.255</td><td>516.11</td><td>71.3392</td><td>0.011866</td><td>5.16518e-06</td><td>1.00016</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161131.547+525241.333</td><td>242.881447116</td><td>52.8781479465</td><td>262.442</td><td>560.661</td><td>81.3861</td><td>-0.00204203</td><td>4.91504e-06</td><td>1.00381</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161128.756+525312.827</td><td>242.869815896</td><td>52.8868964635</td><td>256.02</td><td>554.574</td><td>68.8019</td><td>-0.00204203</td><td>4.91504e-06</td><td>1.00221</td><td>1068.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161132.773+525325.277</td><td>242.886553205</td><td>52.8903547935</td><td>68.2279</td><td>171.721</td><td>16.9382</td><td>0.377187</td><td>4.94219e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161138.691+525423.576</td><td>242.911211156</td><td>52.9065487865</td><td>67.2236</td><td>176.488</td><td>18.5074</td><td>0.377187</td><td>4.94219e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161138.467+525405.234</td><td>242.910280018</td><td>52.9014539065</td><td>65.6982</td><td>167.994</td><td>17.4561</td><td>0.377187</td><td>4.94219e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161137.367+525312.879</td><td>242.905696275</td><td>52.8869108555</td><td>65.1079</td><td>157.618</td><td>17.9052</td><td>0.377187</td><td>4.94219e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=10>\n",
       "          help_id                 RA      ... flag_mips24 flag_mips_24\n",
       "                               degrees    ...                         \n",
       "           str100              float64    ...     bool        bool    \n",
       "--------------------------- ------------- ... ----------- ------------\n",
       "HELP_J160000.316+542000.690 240.001318675 ...        True         True\n",
       "HELP_J161145.788+525236.929  242.94078253 ...       False        False\n",
       "HELP_J161138.168+525244.570 242.909032282 ...       False        False\n",
       "HELP_J161143.750+525251.314  242.93229323 ...       False        False\n",
       "HELP_J161131.547+525241.333 242.881447116 ...       False        False\n",
       "HELP_J161128.756+525312.827 242.869815896 ...       False        False\n",
       "HELP_J161132.773+525325.277 242.886553205 ...       False        False\n",
       "HELP_J161138.691+525423.576 242.911211156 ...       False        False\n",
       "HELP_J161138.467+525405.234 242.910280018 ...       False        False\n",
       "HELP_J161137.367+525312.879 242.905696275 ...       False        False"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XID_MIPS[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGoCAYAAAAAZTE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XuYZFlZ5/vvWmvviMhbdVd3V0M3I9AOutFWwIMeGJGR\nQeUiHBWPwtMKfRAVvCDe4DyAcAAVbR104HgQFMXxAioKzAHUUQRpRkYPKuIN2SKIYoNQ0Le8RcTe\ne63zx1p7R2R1ZVZWdUZGZuTv8zzdVXHNtbOq4s13rXe9y4QQEBEROWx23gMQEZGTSQFIRETmQgFI\nRETmQgFIRETmQgFIRETmIpv3AFpnz67PvRzv9Ollbrtta97DOBS61sWkaz16zpxZM/Mew1GlDGhK\nlrl5D+HQ6FoXk65VjhMFIBERmYsjMwUnctDe9f5bAFhbHbC+Mezuf8SD7jWvIYnIFGVAIiIyF8qA\n5NhrMx0ROV4UgOTE2StgaXpO5PBoCk5EROZCGZAcG5pqE1ksyoBERGQuFIBERGQuFIBERGQuFIBE\nRGQuVIQgMmW3QgeVZ4scPGVAIiIyF8qA5Eg5bqXWyphELp0yIBERmQsFIBERmQtNwYnsw3GbGhQ5\nDpQBiYjIXCgAiYjIXGgKTuZCU1oiogxIRETmQgFIRETmQgFIRETmQgFIRETmQgFIRETmQlVwMlOq\ndhOR3SgDEhGRuVAAEhGRuVAAEhGRudAakBwIrfWIyMVSBiQiInOhACQiInOhKbgjaq8prd2Oez6o\n46E1nSYih0EBaM4u5cP+Yl9zvuevrQ5Y3xhe9NcWETkoCkAiM3BQ2ajIIlMAOiSa1hIR2UlFCCIi\nMhcKQCIiMheagjtgmmoTEdkfZUAiIjIXCkAiIjIXCkAiIjIXWgMSOUSX0uFCZFEpAF0CFRqIiNx9\nmoITEZG5UAASEZG50BScyBGh/nFy0igA7UFrPSIis6MAJHLEKTOSRaU1IBERmQtlQCLH1G4HDT74\nflfOYTQiF08BiMk/ZJ0SKovgUtYuNZ0n83CiApCKCkTO7zD+bSjIyblMCGHeYxARkRNIRQgiIjIX\nCkAiIjIXCkAiIjIXCkAiIjIXCkAiIjIXCkAiIjIXCkAiIjIXCkAiIjIXCkAiIjIXM23FUxTF84Cv\nAXrAz5Zl+Yuz/HoiInJ8zCwAFUXxCOBLgYcBy8Cz93p+UE8gEVlAxhiz3+eePbu+cJ+DZ86s7Xr9\ns8yAHg38DfBm4BTwnBl+rQNhjOGkxEFd62LStcpxMssAdBVwH+DxwHXAW4qiuH9Zlrv+jbmIHxRm\n5iiM4bDoWheTrvX4On16mSxz8x7GoZllAPoM8MGyLMdAWRTFEDgDfGq3F8z7p5mT9BOVrnUx6VqP\nnosJkrfdtjXDkczHmTNruz42yyq4PwYeUxSFKYriWmCFGJRERERmF4DKsnwb8JfAe4G3At9dlmUz\nq68nIiLHy5E5kO4oVMEdl5T+IOhaF5Ou9ehRFdzuVXDaiCoiInOhACQiInMx004IIiKX4v3vfz8/\n9mM/hnOOhz3sYTzzmc/c8fj6+jo/8AM/wNbWFr1ej5/8yZ/kzJkzMxvPcDjkOc95DrfeeisrKyvc\ndNNNXHHFFTP7eieFMiAROXJe/OIX87KXvYzXv/71/PVf/zUf+MAHdjz+5je/mc/93M/lda97HY99\n7GP5xV+cbZevX//1X+++3td+7dfyqle9aqZf76RQBiRyTL3pTW/iHe94B5ubm9x2221813d9F495\nzGN473vfy8tf/nKstdz73vfmJS95CaPRiBe84AXceeednD17lm/6pm/ihhtu4ClPeQpXXnkld9xx\nBy984Qv5oR/6IZxzhBB42ctexjXXXMNNN93E+973PgAe//jHc+ONN/Lc5z6XXq/HLbfcwtmzZ/nx\nH/9xrr/+eh75yEdy3XXXcb/73Y/nPe953Vif8YxnsLU12eNyv/vdjxe96EXnva6NjQ3G4zH3vve9\nAfiyL/sy/uRP/oTP//zP757zuZ/7uXzkIx/pnp/nOQA///M/z+d93ufx8Ic/fM/v06Mf/eju8X/+\n53/mBS94wY4xPP7xj+dJT3pSd/t973sf3/qt3wrAf/yP/1EB6IAoAIkcY9vb27z2ta/l1ltv5YlP\nfCJf8RVfwQtf+EJe//rXc+WVV/KKV7yCN7/5zVx//fV89Vd/NY961KP45Cc/yY033sgNN9wAwOMe\n9zi+6qu+ite97nU84AEP4NnPfjZ/8Rd/wfr6Oh/84Ae55ZZb+M3f/E3quuabv/mbechDHgLAtdde\nyw//8A/zhje8gTe84Q285CUv4ROf+ARvfOMbOX369I5x/tzP/dy+r2ljY4PV1dXu9srKCh/72Md2\nPOfyyy/nPe95D4973OO44447+LVf+zUAnv70p+/7+5Rl8ePvPve5D7/6q796wTGtra1141lfX9/3\n9cjuFIBEjrEv+ZIvwVrLVVddxalTp/jUpz7F2bNn+b7v+z4ARqMRX/qlX8qXf/mX8yu/8iu8/e1v\nZ3V1laqquve47rrrAPiGb/gGXvOa1/Dt3/7trK2t8f3f//185CMf4cEPfjDGGPI854EPfCAf/vCH\nAbqM5JprrukypNOnT98l+MCFM6Bf+7Vf4/d///cBuOmmm9jc3Owe29zc7D78W6985Sv5tm/7Np70\npCdRliXPetazeMtb3rLv79Ott97K1VdfDewvA1pdXe3GtLm5yalTp3b9WrJ/CkAix9jf/d3fAfDp\nT3+ajY0N7nnPe3LPe96Tn/3Zn2VtbY13vvOdLC8v89rXvpYHPehB3HDDDfzpn/4pN998c/ce1sal\n4He84x188Rd/Mc985jN529vexmte8xoe9ahH8aY3vYmnPvWpVFXFX/7lX/J1X/d1wPlbzLTvda4L\nZUBPfvKTefKTn9zdzvOcf/mXf+GzPuuz+OM//mO++7u/e8fzT5061QWlK664go2NjYv6Pl155ZXd\nY/vJgL7oi76Im2++mQc84AG8+93v5sEPfvCez5f9UQASOcbOnj3LU5/6VNbX13nRi16Ec47nP//5\nPOMZz8B7z+rqKj/xEz8BwEtf+lJ+53d+h1OnTuGcYzwe73ivL/iCL+C5z30ur3rVq2iahuc973lc\nf/31vPe97+VJT3oSVVXx2Mc+luuvv37m1/XiF7+Y5zznOTRNw8Me9jAe+MAHAvC0pz2NV7/61Tzr\nWc/qphqrquJHfuRHgPOvAe32fboYN9xwA8997nP5pm/6JvI852Uve9nBXOgJp04IU47LzuqDoGs9\n/t70pjfxT//0T/zgD/5gd9+iXuv5nO9a24zvoQ99aHff+b5Ph0mdEOZzHpCIyKG6//3vz7XXXjvv\nYcg+KQOactJ/elxUutbFdFyuVRmQesGJiMgRowAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJz\noQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAk\nIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJz\noQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzkc3yzYuieB9wZ7r5T2VZ\nfsssv56IiBwfMwtARVEMAFOW5SNm9TVEROT4mmUG9EBguSiKP0hf5/llWf7pDL+eiIgcIyaEMJM3\nLoriC4GHAr8AfA7we0BRlmV9vueHEIIxZiZjERGZo31/sNV1E7LMzXIs87Dr9c8yA/oH4B/LsgzA\nPxRF8RngGuBju71gVsFwv4wxcx/DYdG1LiZd69FzMT9Y33bb1gxHMh9nzqzt+tgsq+CeBvwUQFEU\n1wKngE/M8OuJiMgxMssM6BeB/1oUxR8DAXjabtNvIiJy8sxsDehihSMwkOOS0h8EXeti0rUePeYi\n5uDOnl0/+hd0kc6cWdv1+rURVURE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE\n5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIB\nSERE5kIBSERE5kIBSERE5iLbz5OKongY8IXALwEPKcvy3TMdlYiILLwLZkBFUXwv8KPADwCrwM8V\nRfHsWQ9MREQW236m4J4KPBrYLMvyM8CXAE+b5aBERGTx7ScANWVZjqduD4FmRuMREZETYj8B6Oai\nKF4GrBRF8XXAW4B3zHZYIiKy6PYTgJ4DfAj4K+BG4HcBrQGJiMjdsmsVXFEU9566+Xvpv9a1wL/M\nalAiIrL49irDvhkIgJm6r70dgM+e4bhERGTB7RqAyrK87jAHIiIiJ8teU3AvLsvyxUVR/BIx49mh\nLEuVYouIyCXbawruL9Kv7zqEcYiIyAmz1xTcW9Nvry3L8senHyuK4sdmOioREVl4e03B3QRcDXxN\nURSfc85rHgo8f8ZjExGRBbbXFNwbgc8HvoJYEdeqgR+Z5aBERGTxmRDuUl+wQ1EUl5VlecesBxIu\nNJBDYIzhCAzjUOhaF5Ou9egxxpgLPys6e3b96F/QRTpzZm3X69/PcQxfVxTFTwGn020DhLIs3UEM\nTkRETqb9BKAXAY8oy/JvZz0YERE5OfbTC+4WBR8RETlo+8mA/qIoit8G/oB4FAMAZVn+ysxGJSIi\nC28/AegyYB34D1P3BUABSERELtkFq+BaRVGcLsvytlkNRFVwh0vXuph0rUePquDuRhVcURQPBH4T\nWC6K4qHAu4EnlmX5voMbooiInDT7KUL4GeAJwGfKsvw48J3Aq2c6KhERWXj7CUDLZVn+fXujLMu3\nA/3ZDUlERE6C/QSgW9M0XAAoiuKbgVtnOioREVl4+6mC+07gl4Hri6K4HfgQ8OSZjkpERBbevqrg\niqK4B7ABOODqsiz/cT9vXhTF1cRzhb6qLMsP7vVcVcEdLl3rYtK1Hj2qgtu9Cu6CU3BFUTwL+L2y\nLDeJ/eDeWhTF0/fxuhz4OWD7IsYqIiInxH7WgJ4OPBygLMt/Bh4MfM8+XvcyYrXcxy95dCIisrD2\nswaUA6Op22NSQcJuiqJ4KnC2LMvfL4riefsdzEVkqjNzFMZwWHSti0nXenydPr1Mlp2cgwb2cx7Q\nTxDb8Lwh3fX1wHvKsnzhHq95NzFIBeBBwD8AX1OW5b/t9hqtAR0uXeti0rUePVoD2n0NaD8ByBE3\non45UAHvLsvyv+33ixdF8S7gO1SEcLToWheTrvXoUQC6ewfS/VlZlv8L8NsHNyQRETnp9hOAPlkU\nxcOB95ZlObrgs89RluUjLnpUIiKy8PYTgL4YuBmgKIqAjuQWEZEDcMEAVJblmcMYiIiInCz7OY6h\nBzwbKIj7f74PuKksy/GMxyYiIgtsPxtRXwmsEjeg1sD9gF+c5aBERGTx7ScAPbgsy+cDVVmWW8D/\nAXzRbIclIiKLbj8BKKRpuLY+/Sou0AlBRETkQvYTgF4O/CFwTVEULwf+HPgvMx2ViIgsvP0ex/D5\nwH8iBqyby7L864MeiDohHC5d62LStR496oRw945jyIFHAY8hBqGHFEWxWB0ARUTk0O1nI+ovAEvA\nzxMD1o3A9cRybBERkUuynwD0kLIs79/eKIrircDfzm5IIiJyEuynCOFjRVHcb+r2PYBbZjQeERE5\nIfZ7IN1fpTN+auDLgE8URfFOgLIsHznD8YmIyILaTwB60Tm3XzaLgYiIyMmyrzLsw6Ay7MOla11M\nutajR2XYd6MMW0REZBYUgEREZC72cxzDdcDjgc8BPPCPwFvLsvznGY9NREQW2K4BqCiKa4h94O4D\nvIcYeCrgOuANRVF8FPjBsiz/dfbDFBGRRbNXBnQT8JKyLD9wvgeLongg8OPAU2YxMBERWWyqgpty\nXKpqDoKudTHpWo8eVcHtXgW3nzWgzwaeQTwHqHujsiyfdiCjExGRE2k/G1HfSDwP6H+gg+hEROSA\n7CcAmbIsnzPzkYiIyImyn31A/7MoiicURaE9QyIicmB2LUIoisITp9zadZ/2iQYIZVm6gxyIihAO\nl651Melajx4VIVxCEUJZlrtmPEVR9O/uoERE5GTbz5Hcf3LObQv8+cxGJCIiJ8JenRDeCTwi/d5P\nPVQDb5ntsEREZNFdcCNqURSvKMvye2c9EK0BHS5d62LStR49WgO6xI2oRVE8mnQgXVEUNwL/K/Dn\nZVn+14McoIiInDx7FRq8HHg+MCiK4keAbwb+Dvj6oij+70Man4iILKi9ihAeBTyyLMt/Ix7H8DVl\nWb4KeALwVYcxOBERWVx7BaAt4Or0+08CK+n3K8RCBBERkUu21xrQS4A/K4riN4APAjcXRfGHwKOB\nnzyMwYmIyOLaNQMqy/KtwMOBjwM94E+AdeCpKkIQEZG7S+cBTTkuZZ0HQde6mHStR4/KsHcvw1aD\nURERmYu9OiF8HLjHeR6aSTNSERE5WfYqQvgS4J3AE8qy/MAhjUdERE6IvYoQbgG+n1gNJyIicqBU\nhDDluCxqHgRd62LStR49KkI4oCKEoiiefveHIyIicvFVcN8xk1GIiMiJs2c37PPYdypZFIUDXgMU\nxOO8v6Msy7+9yK8nIiIL6mIzoFdfxHP/N4CyLB8GvAB46UV+LRERWWAXzICKosiBrwSuArbTuUCU\nZfkre72uLMv/VhTF29LN+wC3382xiojIAtnPFNxvAdcAf0+cSiP9umcAAijLsi6K4peJRzh8w4We\nfxHFIjNzFMZwWHSti0nXenydPr1Mlp2cPf77OZL7g2VZ3v/ufJGiKO4J/H/A55dluXm+56gM+3Dp\nWheTrvXoURn23SvD/nBRFPe+2C9aFMVTiqJ4Xrq5Bfj0n4iIyO4ZUFEUf0Scarsa+Czgr5g6iK4s\ny0fu9cZFUawAvwTcE8iBm8qy/H93e74yoMOla11MutajRxnQ7hnQXmtAL747XzRNtT3x7ryHiIgs\nrl0DUFmWNwMURfEzZVl+z/RjqbDg5hmPTUREFthexzH8AvDZwBcXRXH9Oa+5fNYDExGRxbbXFNyP\nAvcFXsHOjtg1sSRbRETkku11HMNHgXcDzyTuA7oH0AB/XpblrYcyOhERWVh7TcF9KfDLwEeBf0t3\nXwN8TlEUTyvL8h2zH56IiCyqvabgfh746rIsPzR9Z1EU9wPeDHzhLAcmIiKLba+NqO7c4JN8hIvo\nii0iInI+e2VAbyuK4q3AbwCfSPfdE/hm4HdnPTAREVlse/aCK4rifwceB1xLzHpuAX63LMvfPuiB\nqBPC4dK1LiZd69GjTgi7d0K4YDNS6I5k+DxgXJblBw9wbB0FoMOla11MutajRwHoEpqRFkXxe+nX\nLyDu+/ll4DeLovibdJ+IiMgl22sN6B7p15cD31OWZRuQvhz4ReAhMx6biIgssP0cx3BZG3yg6xG3\nPLshiYjISbBXALpfURSvAkZFUXw7QFEUp4uieDaTqjgREZFLstcU3PXAlwC3EzsgQCzBfgjw1NkO\nS0REFt1eB9INyrIc7vXi/Txnvxa9Ci6EcKTOrz8uFUQHQde6mI7LtaoK7tIOpHtdURT/HfiNsizX\npx8oimINuBH4SuAJBzLKYy4E2Ouv2VEKPiIiR8FeAegbge8E/qwoituBfyUexXBf4EriMQ3fOOsB\nHh8BdSgSEdm//W5EfSDwOYAHPlyW5V8d9ECO+xRc+7Ljkugcl+mLg6BrXUzH5Vo1BXcJU3BFUbwH\n+BDwe8AfzCLoLJpAwKQs6Kit+YiIHDUX6gV3P+CxwKOAJeBdwO+VZfkXBz2Qo5ABnc+5gSSENNlm\n2qxnMmwzNQXXTsgd1Rh0XH56PAi61sV0XK5VGdDd7AUHseINeAQxID24LMsvO5DRJUclAF0oc2kD\nkDXgAxhCF4Js+gcxffvi3vvwsqbj8o/3IOhaF9NxuVYFoEurggOgKIrrgMczWQP6EPDTBza6I8dM\nBRaT7pkUGLQBpgnTlW0GCKkSLuZBIUwCSvvrhf4easpORE6SvdaAriH2gbsP8B7gH4EKuA54Q1EU\nHwV+sCzLf539MOclpGk00wWUaYbQPRZjVKANVzDJZnZO4WltSEQE9s6AbgJeUpblB873YKqM+3Hg\nKbMY2GGaDgomJjO0WQ0YGj+9zpNyIQONB/BYm4JOiIUIwbcrQ76bhjOmnawz3fu0b3uU14pERGZl\n32tAs3ZU1oDOp27a0BEDlU3BovYxcNgUq0LajWqIU3Tt89og1E7fufRAO9U3TWtAB0/XupiOy7Ve\nzBrQb739g+ERD7rXLIdz6O7uGtBnA88ArmJqp2VZlk87kNHNQVtIwNQqD8Rg4qfTEnZmJ22VW8x8\nIHhP3QSMNTgL3oN1BoLH+7hgZozB24CzMQPyAUITYlAybY6VAlT6n7IhETkJLhiAgDcCfwj8D+Do\n/7hxEdrP+XBOIArsbBM+mZabCEyVZBODUjs957uihAAmBpi4hjR5jbM73619v0l5t4oSRGSx7ScA\nmbIsnzPzkRyi6YASptZ62g/+qSUffACb0pK68efNULwPGAO1D/gQdry/bwLeB6w12DQ9hzF4TwpO\nO/cMhalsTERkke3nQLr/WRTFE4qi2M9zj7x2ztiYSYZhMPtrJJrWebwPOxKiNmsJAZompPebZFNt\nlVz8mmltyEzWhOLzwlT+FVD8EZFFt1cZtmfyw/l3AKEoCtLtUJalO5QRHrCupBroMp8dlWtxfWaS\nxYS456frgGDwIdA0oSsyaPwkePgA49pj2heQoo6PU3LdEpNpsx6DsWnvkA1dMOwCVxqPiMii2TUA\nlWW5a8ZTFEV/NsOZnbvuv2kn38I5z4PMpSo1AoQUHEyYFCOkbMc6Q9OE7nVtzAphsikVUogLcW2o\nDSrt4zbtWg2A8QZrQjdV1wY8ETk53vX+W1i0SrjdXHBarSiKPznntgX+fGYjmpHpD/KpGbVUgTYp\nhjYmBpW68YyrhqpqGFcNTe1pfMx8qrqhqhuG45q68fgm0HhP46FpPHXjCQHqxlPXnqbxqcLOM6oq\nxlVD3cT7myZ060TexwypTr8nBEJaV/IpqImILIq9puDeSez91k7HtWrgLbMd1sXbT5+1Vsw82hY5\nO6vO2nWZgKFuJqGpnqrOrn0sTnBTtdrex2zKWUPVBHJrYmAhkFmLweC9p2kAApkxWDspxSbtCWor\n5oydTMN11XrqoiAiC2SvKbhHAhRF8YqyLL/38IZ0aS62z1obdLoP+FQU4AMQYjbjLPgwtX5DTBmt\nNdgQ135CmD6GITBO2VNsxWNpGs+wCTjnU4YTp/6qqknrSAFjLZk1kFlMWqCqpg4YsobJviHtExKR\nBbFnGXZRFPcHPl0UxauJ+yo/zoyOYzgMOzOIc7KJqcKENPvV7etpiBcPKXsiPuCbEAsYArHogFgF\n1zSBPDOpYi5mNj5MrRVZqOpAjicYgyUQ7NRhDm1WZibdFGLgMhizc31JROS42qvQ4LuA3yB+9v4Z\n0AadXyiK4gcPYWwHrmscymS6La6thG4ajFRCHfOZyYd8CKTy6/j8xsf1nbr2VLVnXPtUkBCn4Xzj\nGY8bGh+oGt+t8cT3il93nNaHvA+M01rTuGq6daS45hS7Lfi2/LuNTLTFDoHj0I5ERORce2VA3wc8\nqCzLrek7i6L4aeB9wE/NcmCzM/kAn74H4t023Q5mat0oTZu1mUlb5dZmSr7xWGfBgrO2e13V+K6k\nu92MCnG9qE5l3JmzXWbWpCm6zDqaAD0bCyKsnZRkB2xXVdcNe6rX3PTtHde4x2MicrS86/23ACx8\nNdxeAagC8vPcv5QeO5bu8mGdygzaoAIx8zEWjE+3ps708e2mUmK5dt2AczZmOAbq2mMtaR0IhuMG\n5+Iaz6hqsNbgfYgByMbpucw5TMpurDVdJV5VN2TOYr3BOYs1cYovpDlA2/aSS8UUewUXBR4ROWr2\nCkAvBf6yKIp3AJ9I910DPBL4oVkPbFbOrYaL97W32ZFZWGu6cmqTCgF8CHHdJn34ZwbqOhYsNI1n\n3HiWnWNcp4q42pNnDh+gagI2BZ+68SwNMtqNsN7HogYX2mm3GOScjV8vpFYKcbowTRBOTSme7zoV\ndETkKNt1Dagsy9cDX0ZsQroFbKffP7wsy984nOEdvPZk0jabCanWuSs6SBtDvW/38MRAMt2NoF0D\nCiFQ17HizYd0bEPqhGCnvrObW2N8en7jPT5lOt4HqtozGsd9RhCbmg7H8TkBGI5rhuM67R2a7BcK\nqQrPt1OBTNoBtdcpInKU7VWEcK+yLD9eluWvAG8GNoFl4Fi24GndddG+bY8TP/CtoTtaoe3LZq3p\nzvDperl1RzPEzMQTYgaTnpdZi7WG4bgBa6h9XHvyPq4HZdbQzx2ByXlDbQeGce277bFNM1l3ahof\nR+R3rud0G2vbpnMiIsfAXp0Q3gpQFMU3AG8nHsX9ecDN6b5jKX5om64nW8CkbGUyHRcwZM6kKThD\n49ugM9U8NGUbmTU0tccZcG6SXYUA46phqedStjQZQ5sJxcynoWoaNrfHbA1TZwXv2R7FrAfoKuHq\nJlbcVc0kc6qbpsvGpt+/va0qORE5qvZzHMNzidNu/wxQFMWPAX8E/PYsBzZ7oWsGGsvfUnVc+rC2\n1mAwNN5jzaR7QpOyl9q3veFSR7kAvcymc4EC3sQCBWsNPeswBJwzBGupG4+zFh9iiXYbH8ZVw6Cf\nEUKgamJxAsQiB2gPy/P4YHAhBro8dVloq7PbibedrYc0HSciR89+jljYBD42dftTHPN5nthUoN1m\nurMwYapEAYAsi1VyTapSaz/KnU1ByZlYqWYMwbfvH0ure5ntSq/HdcxahuMaAmkfUZPKueP03bjy\nbA1rxpWn8Z7huGFUxUwovjb+2vhJt4W48dV3GVlXGu4n+5u6ta4pyopEjr62HHtR7ZUBXVYUxQeI\nQeqnge8riuLBwIuIxQh7KooiB14L3BfoAz9aluWR6SF3blLQFiC0H9htKDIYnCF1MojTdT6kjtXW\nEIIn69kXNAenAAAgAElEQVTYaHTsscSgNaxSFmUt9ShOpVV1DDTWQC/P2NiqyDOLr1PHBAPjypNn\ncX3I+4C3JgWc0J6bGsu+Q1xnCuli2i4Odrqxand1d70HVCknIvO1VxXcvydWwX0vccoN4J7A7wPf\nvY/3fjLwmbIsHw48Bvh/7t5QZ6edhmsziJ3aBf9JhVzbxDR2PUjZUfdmMfBkdpJp9XoOZ0n/xam8\nuvH08liebU1sVAqTo7rr1F17XDWMRnWsnpsq4W4a31XHjat6UhkXpjouTCbkODf4+MBd7hMROUx7\nrgGVZXkrMeC07lWW5Sv3+d6/xWSdyBC7aB9pgenMKE7S+biiEzeCuhgQfB0Y5I6qDhjjcTZQ13Fj\nqbMmrRG1a0IBE4Dcxa4GKTjkmU1dDgLDUYO1FgPdlJ01FmdNbNdDICcWOfgAwcf+SE2A0HgGfZc2\ntbZTa5Oqvelsru0hN73vacf1q1uCiByi/RQhTPsO4Of388SyLDcAiqJYIwaiF1zoNYfxwddVh6X/\n3XUqbtIZoS25DqmVjrGk9jkxg2mP3rbEgEITcw4fYKnnMKbdoBrfN24snXTRbrMaZyddrkMIjOuG\nuvZdVV1sxeMZV7Eyz5qazFpqG3DGMho3YGLga99nzKR4AsC0O2fjRULK9uzU9S9y4cKlXs/5Ni4f\ndcdlnAdh0a51ZbmHtTsnps6cWZvTaGbvYgPQRf1pF0XxWcQ9RD+bNrbu6VAXxu+y/jF19HVavHdt\nJ+ruKISQKtLi3iFDiM9JFW2YWBTgQ6CXWzIXM6Cq9pjG4Gzc7Tqu4zrQ9qjCWYsx0MstgcBo3LA9\n8iz10lHdIZBnsRqu8fGAurEx9JddXO9JEadpAq7tlJA6Z7s8vq49/nvSudvs2Du0V++4RdC2Urq7\njsP35KCu9Tg4Ltd6MUFyc2t8l/vOnl0/yOEcur0C6MUGoFfv94lFUdwD+APgmWVZvuMiv87MnfuX\n4tzb1sQMwaT9PU2TWt8QsxjLOZ0Hpv5zzmCJxQMA1oL1JvVwi0ErAP3cMa59N10WS68Ny31HnsVp\nNQI0tSeYWMTQZljGQGYhyx25jb3omsbiXJx+sxiaJlY2GADXZnftoNtNrGZHEOoeUoGCiMzYBQNQ\nqmb7SuAqYLsoihsBUoeEvTwfOA28sCiKF6b7HluW5fbdGO+hMMZ02U93n5tqd5OC0+QMobZzQiwo\nyPOYQm+PYhuD4FNQSpkKWTyKwTmLqT29nottdirPoOdwzjHoxfeo65hRjcY1Vd3uH/JkmaWqA6su\nI6S2Pt435MFic5dOW00XkabcMFNlCTu6J+zMhMz0E0644zgFJ3Jc7CcD+i1iE9K/Z/J5G4A9A1A6\nRfXIn6R6Pu3RB7udompMnBqb/pBup23b84WsMeSZoanBtMUIBqydFAkYLG6Q0QTo55aqjkUGzsAw\nlW6bqfWhWOkWg1TbMWFze8xSPyNPbX0AqrqJGZKLhQwBQ0PAtetbgA2TjbRtDV97NSG0h1YoC5q+\n/unNviphF7n79hOA7l+W5f1nPpI525EBdB0OdmrXUayxGBcntOrUM86aGFSa1BnUxZ2qGGJjURdS\nW5x2X4+HzIEP7fqOZ7mfM649mYOtjYrMQc+5ruForIDzDHo9xrVPe4ZcF/R6WVx32hzW9PKMqXwm\nZmtmKrhYdk67TX8v0vfDWn3AwjlrP+33T8FH5G7bTyeEDxdFce+Zj2TOLvSBMnm8nc9KfeSYdFXo\nFvdTwUK3HpTWYdo+cdYY+rlJQSuuN8XHQ2pIahj0bDrOoaFuPNYa+mlqb2O7ThtjA1vDinHdMBzX\nbGxXjMZ1V+DQdlJoGt91W5gc5hAzKt8t+pgUnCal2MdgfXemdvTUY+qHgAD1jhNuj3lrEJE52TUD\nKorij4j/rq4G/qYoir9iai9PWZaPnP3wDteFqmraM4GArhjAOtNN2cVAY2na220Nd/Bpn07ccFob\nQ55ZtoZ1XNNJGZFJgWlcewb9jKYJbI8qfAj0c8eglxFCxXBcc1mvz53jMWvLPZomZlBZ6pBgTeyo\nECvsTCxAaK8tQJbZNBUXOyfsqIjj/Ms/mnKK2sww+PhnD/H7aHZuIhORfdhrCu7FhzWI42LnekBa\nJZlqgTBZO0nPTeXQxgI+pjntvh4D5LmNHbUdOOsYYjAmQIhrRsEGlvpZ7HZQe/q5Jc8MYKnqhn6e\nxf1C1mCtZVw1bA1jcLPWdAfp1fWk9DpW8k1Of/UhxCyua+ezM8CadEUnMfi0WWAsdQdCyhjTkRoV\n8Qj2yXdJ5GCd2CO5y7K8GaAoip8py/J7ph8riuKXgZtnPLYjzbTlzVN8yiDagoSQNqlaLCEVEwQb\nMD72fnM2HtVtDJBDvxcYVQ2Zi1M/m8OapX467jvEhqO93OFc3HzqbAw205nb9rgGk7GcZWnTagx8\nPmVImds5tWansp/0m/MEm8n+oUV2viyv68DXTrd5aEz683M2fV/CVNgWkf3aawruF4DPBr64KIrr\nz3nN5bMe2HFkMISp6ufptQFj2vWCtgLOd1lHCJMPOGsgpI4Gg1SevbKUYUcNjY8fkrmz+MzjQzor\nqI7TQSa9fjiqIQR6vQyfBXJvCVlqU2pCF3RsqpA75yK6qjuI61zT03LTwWvR4tH5Amy7d8oYk7qO\nh27/V1U1WJO6U9i7lu6LyN72moL7UWIn61cAL5m6vyaWZMuU9qfn6c+g6SKyEFJ37PQJHvfzkPb4\nNCl4GPIswwWfgk0s2zbGYZ1je1jTpIDkN9MRDr6Jp6qmqTiXWbZHDZmz9ELcn+Rtu4geCKnFz/Qx\n4z5FwbBjzGZHdjSdAe183qKLAbvrIGEDvplkjO3BuSacmG+IyIHZawruo8BHgQe29xVFcUNZlr8+\n+2EdP/uvooN2702qnSPPHHXauxPXY2yc68lsqmCD3Blqx6S1j4nZj3MG3wSakI72DrG/XHv2UGwH\nZOMuoNzR7k+NP8XH6joTF4doq8+nRxo/Zyd7hbqNYCemKKEtvo/XXNdp/cyY1GjWp6k4j7Nu4bJC\nkVnaawruxvPc/cOpM8J+OiHIOXZWmcVP+rYO3vXcjrWcxhqqJlZctYfhDfoZVVOnAGQZVxWnVnvU\nTdwT1DSeOk3lOWepmwD4dLR4/Do+BGwwXTVXKsCeas668xN0x3Ri2kzbnSJ7IqRA6+NBgH4SgfHp\ndNy2kAOVIohclL2m4L4D+BzgrUz+Va0B/4l9dEKQu7prZ4W42GINNGGyVyg+OaRecAZMoCFOmWUO\nGmPIgNXlnKbxZM4wriY1bMaQihliD+4qZTykDg0+S6XEQEgnupr0wskaz9SH6XRqtGMtaPGzoMmf\nx2SvlvdQp/WgpolHZ3hvcNbiUs+9Bf+2iByIvQLQw4mnn34B8IyyLM8WRfGXZVl+y+EMbfFNf7hl\n7bQY7dqPxZhAMAEfDKaJP3H38oyq9iz1HFU6/qFqYkGDs3GPT/ve8QMxHbTnPaMqvm/tA/2ei925\np35ib4+g2C3D2U+3iEXVrgFlzjJsPMOxJ7NtVwQfT75tPNY67H62d4vInmtADfB/FUXxMOAtRVG8\nFG34vigXWifpChe65wTaM4p8W0FngJCCizOEVM9tjGF1KWdU1WAMvjFUoT0tNWY4Vd10xQbOxpLh\nPIs/sddNwFqPaUitgyZFCXFsaTzdGojpOjxMW/QMqNVeZjy5Ngb+UdXENktpw6+zAefac6ROxvdF\n5O644M9qZVm+B3g08ERiVwTZp4srTJjeLEoXcDJnu4Pl8ix2Q1ga5Az6GYN+xnI/xxpLr+dY6mfx\nGO+6nZbzbI8atscNoyr+165jBO/JrMOayRlC0wUG7fRaO7Fn7jLOk/UBa4xJBwRaGh8rGofjhqr2\nXS++xvtunU9ELmxfkwVlWd5JPNH0uUVR/PvZDulkaz/cp6e6smyyF6cNTiEEMhszl6V+Fo8Cd3Ev\nSm+q+3Zbbdd4z2hcMxzXjMYNdROPeBjV8fd1EybBKfU8837SNy6NjtB2fDgBG1OntaffGuI5TBCn\nO+smsLVds7k1Zmu7oqp9qkbUZIHIhewagIqi+IqiKD5WFEVZFMXXAX8MfA3wB0VRPPHQRnjCTH9w\n2dS41Jh4xlDsZGDp5Rl5lg6tIwacQc91P5kvDzJCiBlQP3f0cpsKEzyb21UqzzZ46E5w9T7gm7S3\nJQUgm7o0tIc4tL/GjZnnH/uifvC2wT/P4vcf4t4riI1JN4Y1G9s1de3VLVtkn/YqQvjPwGOBVeCP\ngAeUZfmhoijOAG8H3nAI4ztxzu03F3/qNqkqbnJ6aTwbCPLcMRo35JmjnzeEEMu5syywPMhomliy\nbQKp03bMbNa3Kvo9R5bFBql55rr3bYfQ+NQHoasAM2mZahJo2mk6c57xL6I6nUToTOxUsbldUdcN\n1hpGVZySaxqPT90RRGR3e03BZWVZ/i3wXuCOsiw/BFCW5Vku/ihvuQTT03AG4jRbZuMRDjauSVhr\nyHNH5ixL/R7Lg5xe7uhljrXlHoN+hvcBl8U1olhdFxfQrTHUdaCqQxd4fAipo3bMvrpahG5MUx3B\n21+5a3HCQgptM9n4fbHWMk5TbnUT0tRc2w9QwUfkQvYKJH9dFMXrgRXimUD/Gfgl4AnAPx7G4GRi\nOrMwqUMZqSyY0BBC3E9kDRhnaLJ4ZICzlswFmrrBYPHeM64Czlk2hxWZNdRpQ2qWOXITV30aD233\n7qm9l5ipRaE2G2u7JZyE1ff2SIteHpvEDnqO7VHN9sgzrgz9ND3qnCFz6owgspe9MqBvAf4QeBvw\nCGBMnHb7PODpMx+ZnFf703WeVsIzZ7BpbSjPHL08Y9CLGRDEjKaXOXq9LB6iFmBrFA+5AxjVgapp\nOyGE2HCzbYya1oC6H+ZDewjbpHw8Vs9NHXS3yJmQgSx9M/LMdSXuw3HDbXcOWd+quH1jFDsjpExT\nRHa3VwbkyrJ87dTtH0r/dYqiGJRlOZzJyOS8ptv1ZKlNT9xsGsgzEzetGkMvi4HCmrhR1TmDqVMF\nnYuv29yuYhcEHNvjGucseRaPbIhHiLdZlcHYdi0ofqqGHXNzZnqAO9aHFknsYR5SZhjPc+r3HKdW\nemxuV1RVQ9V41jfHDHLHoJ9pKk5kD3sFoNcVRfHfgd8oy3J9+oGiKNaAG4GvJE7JySHasXHVxDOH\nermLp506n47ejtNzTR4YVTXjynfrQXXjqRvY2K648rIBAHUD26MagiMES3BTR0dkZipVvutG1LbN\nT1f9tWCBp9WufwUbN502PjaAtcawMsjY3K7JrGVzWHOvnlPwEbmAvQLQNwLfCfxZURS3A/9KPIrh\nvsCVxGMavnHWA5S9OQuNbz/o4rRaPIiuvRUbk9rap5Y8qZs2saR4c6tieSkDH88LagONT+cOgcU0\nnoBNhQmTLGxHhnOe1jyL2DHbxIJEIO4HWlvpUdUNw1HG+taYYdXQ6znWN8esrfbI1JdHZFd7teLx\nwCuBVxZF8UBiY1IPfLgsy786pPHJHtqNkcbGkuwmmHQw3Tn7eVJW0mYot62PuccVS4Qq7WNpoJfH\n4LQ1qrHOQhMwfUO/Z+NO/yaAo+uc0J0NlKbbQirPbo/3hsWbgoO2HxxpqjOWY/dyR5471lb6NN7j\nvaffy3BGwUdkL/sqp04BR0HniNnZSy505/xgDP08tvAZV57aQz/PaJpAsIHV5Zyq8gQCTR17mzmb\nY0zs8myHVTrqOxYWDHoO7wPWOnzY0b70/NlQN7746wLGIeL3zNDvOZYbx2js2NqOa0QhENfXMksv\nUxAS2Y328xxjO0qzUzfsdnqo38vIc4CKngGf9qlsjxsuX+2zNawZDmsC4DLwwROa2EeuqsZcfqof\nuyT4wLj25M6mfmek6i9oQ9HOJqaTbanxd+fruH38I5I1BpvFq+vnGYNew/JSxuYwFnMM+lmawhSR\n3SgALaDu890H8txS1x7bbl61ce0nzwyDnqNqPFXlMbRTS7G786hq4ltQM/AW08txLrQNu2nYeUzQ\ndKjp9gu1xdk7jnE4/sEHJoHWpvY8S/2Mpb7reuRtblf0cks/d/MeqsiRpQC0QIyJzUi7c4WIe4Cc\niW1irLMs9XP6vbhzP88co6pme1THA9aI5dlLfRun63LY2q5wpsdyLJaj8bElkPdxCmrHKaoh4ImF\nETDJjGLz0sXIfFrGmHg6apoC7fczrrp8mVHqPN7LrTIgkQtQAFpAcS0o7gGKR28TS7K9p4mRiUHP\nUtUN3jtq61Nrnri7fzSu49k/PpBllu1hRZYZBr0sNUedHEhnbSxSaAsdTFoDmaxHxeY0U4Xj3ZlH\nx12bAWZt9/F0X5b2WY2rhkFf/8REdqN/HQvKmlgRF0KgbsA5WLI5deO5Y2PMylLO9sjifYUHrLMY\nawi1Z1w1VE2gzgPLJiPPMraGDePKszzI06bXyT6XLigxPf2WJuDatj2kTaxdp4TjH4FM2pgbTAzw\nee7wHjLi1Fuvp39ecuke8aB7zXsIM6c5ggXWrlM4254xFD/2Y4VcXPNp1y/axzI3OYuorj3DUc1o\nXFPVnnEVD1zz6eTVuk7/NelQtvONocsLpnrGLUDwmWZNe3igIc/i9zuEQF03i92aSORu0o9oC2x6\n4T9P80VtgLEGhuMaY+L5QcbAcNSwspQzHDWMa5+OXSD9vmFlkJM7G6fxfMDZgPftUQ1x7Wd6f9C5\nx1L7MDlyfPpz+bivDcUO2TDoZfQyS1WnVknOHvtrE5klBaATY1IoEI9dgOWeAwzjqqGfzhUK7cbV\n1BF7XDf4tLkyBN+tJzlnCD6eJxSCJWBjSbdlcopqmo8z3QgMIc3JBbhLgDrO0upXahYb2xg1Pp3F\ntEDXKXKQFIBOiEDoKtacs4wrz9Igp/YVzhq2hhX9PK7r9DLLJ2/bireNofaBZuzZHjXcsVFx+Vof\nZw1XnBqAsakA2xOCITjb7RM694O3zXx2Cz7HuVLOGHDGEIzDuTgd2W4KFpHzUwA6IdrFctLUWJa6\nWy/1M7ZHFYN+xqiKGQ4BLlvtMRrHYxsm+30MxpLuM2yPagKxS0KWOXLAmLRXyKbO2+lrn/s53J3s\nGkL3+HENPtOmpxfjt/v4BlWRWVMAOiHaIoO2e3bmXFeqbW2PUdVgjKGuPVvDilMrfe5kxHDU0ATi\ndBIw6LtUfAArS7HVT0hdTw0OawLjEBj0pvtnT7Kv84xs8rwF+bBue/QtWK2FyIFTADrhYofrQOYs\nmfN4D70sHqa21HP4JjCsmnhYXYjZz6hqcNZwx8YoBiTnABv7xKVMqao93k4yG0NsI71jP1DKjxrv\n4+F3TPYOnbtXqL3/OH2mL0pAlfl41/tv6X6/qCXZCkAnmDEG68BicXj6eZZ6ylmWgI3tmjx3uMyy\nPWoYVw0Gw52bY06t9NjYrtga1awOcurG0eu17X4MvgJrLZk1OAfG0mVKuTPUvt0DEE9ppe2sQKpd\nmOohF9oSb088GO8i7acU+iCCxbnvcTHvp2AlJ5EC0Ak1fbKqMxCswdsQ121sXJtZW87Y3Apsjz25\nM4xrqJsGQ6yOM21KEgLG9KjrBmehl7m4KO88jbW4xpDnDtouCcTuCZ7Jcd/phHBSdx8gBZ6pXqbp\nNPC7tDqd3v3adoGYaoc3NR1muscnh7m2nezOVxTRdrPbfQ3rnO/qVD+83QPK+bpBKPjISaQAdMLF\nEz4hS+sWmTVULrA1rFnq5zhryYYVm8O4j2d71IABZ2OT06r2NE2g8WPy3FL5wFIvo2oCl630sLaJ\np7UCzhkaD97HRqexbHsSOBoPDXE6sAmTzKX7aA7tCaTxZttqrQlggo9FDT5mVM7GXm0hlZ2b9IHf\nPg4xE2v3Ok3r9jJ1j50/gzo3s5p+6q7Bp6sCPP+fhzIhOUkUgCROxRFo0v4fZw25M/hgyILpzrWp\nG8u4auJR1LVnXHvG45pxZWiahqx2+NQhoe0Lt9SP5xA1jaeXxcBiaA/JS2tQFtrsxJh4dEQgNj6F\ngHNtfjI5/hsDNmVgPgWC7uC9EGiaNrjEgNVla0yyoqbZGQimN+6a9LyYTYVu6m86OExPFZquwcNU\nVtauiZk280pNYlP21y52ddnXVOXCdPbVPrZbWJpeH1PnBTlOFIAEiB+SmYVgbcoSMvLccfvGiOVB\nj36e0csrcudY3xpTNYFP3bYdg0rjcduGlaWccdUwrgNryznjKh5glznLnZuexoeuUeeplV5abzJY\nGwNL42Ow6GWWLLNsDmsya3E2HStuiMUNgbTPCAjxELh0FV2lXmwGGosbYs+6QPCxC3i/57rD+pyz\n6X3aDuJMHSI3yY7iBlNL6A46Z2pf1eSutvS8laU372b92t54OzKrFBE9XX+9GIDB2baj6+6H/8X1\nsfges8qelJnJLCgASSf+lB4/sDPbUHvDyiBje9RgiccLODvp9nxqOQYcnKWua7aG0POOUeXxjWeY\nO4ajmn7fkVlH3TT0MkcgUNeePLdxis+ZdBprWpOycZ9S1cS+c5mLveucjR/KsaVPoPZxT1P3k7+J\nQWJ77FIgbOLrUvBy1lDXgTs3A1nq4h0CZNmk/x3GkKXNUCF9qDsbT4e1No6jCQFnJutCtjugLwbK\n9kwgDIyJa2xtJnTu+UhtELJpyjDDpmAY72+aFHggdZEwmBDOyY7MJOjtOC49vv9BBA4Fn/marog7\nLvZTuacAJJ32PCEA23PkAXq5Jctq6jpgnWHceLbHDa4JXHvVCpvDmvWtMetbnvWtirBVsb5VcWql\nR+4MVeMZ9LJ4rHeA4bjhilN9zt62zTVXraQpOs9w3DAcN91YLlvJCRjquukCUTultbbSY3O7Yn1r\nzGUrfTaHFdvDml7PsdRzrC7lfPwzW4wrv2ND7cpSjrOG29dH3PPKZbZHcZ2rbhqaJtDrWXLncNbg\nnGFju2J1qcdSP0tZYZvZeU6f6jMaNzF7S2ckWWfwDawsZd3aVpOCpO365MWwsbKUpSnGmLw0Tczj\nVpfy7vykLIvnMkFIa3WG4LtNV2TtH1ZK07wP3XQqTHUgZ2o97dzsaR+ZTZjKvs73HiKXSgFIdpHm\npYgds33TkGXxmIGVpbz7ib+XG/o9S1U7hqMGCOQungjqbFyHGY5qlvs5Nt1f1w3b45p//dQGq4Ms\nfXAHtsY1dRPwTeCOdZum5kLqvBBofPx9LzOMK8/2qOKT1jIY5Gxsjgh4Brmjlxnu3BozruEztzny\nzDCuYyDMnGF9c8Qdd94JxuCsSyfDBvI8o9cdMxGPsfi03WZtuUeWWeq6YWOrDwTWt8ZkzmCNJcti\n6V6/Z6kbGNd1DGI2Zl4xCzNdJhPXl3w39WhNDEJ1E6sN2wyqDUohBDJr4nEZQPBT03Ht6pg1EGLv\nuTAVl4wJXWZrSL+nrcBLa2rep3Uq0z0euoA/nU2l10z9DWmDXtcBYo+1KpFzzTQAFUXxEOAnyrJ8\nxCy/jhy8trdZrCiz9HtgKlhb7nH5Wp/P3DEihMCoilnEIK/wHu7YHLM0yPnoJ+5keZAzHFWMKt/1\nh1vfqtjYrrh8tYe1luG4Js9ic9Ot7Zo7N0cA1HVNr9cDYGtri+FwSNM0NE1NM97CZgOq0QZ5byV2\ncBhvMR7eiW9SFmUsg9XTNNUQAGtzAoGmGhJ8g7GO3mAViB+2eW/A8uplZJnDGEvTNKyurKQP6zgt\nmTtDv5dxrzMr3Lk55h5XrHRjz1zc82Sd4bKVXlrbsqwuZeTOYqztjmmIFYL51Id+mvZ0hlEVr//U\nSo/hKHYrzzLbtU6arjHIs3ZqcpL1+MZ3pex2uqa9nYpL05c2bRKOhRakqbs4NdmuPxkCNv1daLtZ\nkFottT+etAUX04FqtwRpVhmU1qeOr5kFoKIo/k/gKcDmrL6GzJ41sczZWUOeWeq0JrHcd4yqhhAc\nFQ2DvqXfc/TGcdro1EqPEKCfO8a1x/tAHWJXbWdhfWMb5xyNb9jcrOPm1AaG21vU4xF106Q2QXF9\naXtri7oex4xovEnwAV9XbPFJXLZE8BWj7TsJvqEebQEwXD+LdTl1PSTLlwgBfD1KP917MJast4TL\ncozLqcdb3dlJPsD2ek7e62NtTtbrkznLqN8nsw1V3TAajVhZyun1cvpZRr9ncc5SVXFqLs8sW8OM\nteU8BbC45uVDYGObLqtpCzEaF19vgO1hHavlfIP3Fp+5VHZu0im3hiqtcQVncWmTrk8VGiFMTcnZ\nqf1PifcBk/r+tdrKu3bDcIDUq6/9umnzsNnxsvS8nbe7Fa6prOh8QeIggoeCz/E1ywzow8DXA786\nw68hM+bazgbpQy9zsUrO2n46viF05wpdcSqwPMi4Y2PMKB3x8Mlbtzm1YtkeVdx6Z4U1huVBzgc/\ndAunr7iSW8/eQlVV5P1ljHUMN27j9n/7EPlgDV+NqOshvq6px5tUww2wjtBU+HpMNdqiGm0RQkPw\nNU1dQfDc8amPcOrqz6a/fDmjjVuBuPeoGm3SWz7FYOVKhhufZrh5O6evKeKHs2/oLZ0CIOstkfVX\ncXkPazPywRr95VP0ly9n5fJ7cMu/xgBnsz6rq2tYa8nzjLXlAb1eRj93GGM4tZJz+dogno6aW5b7\n8bHV5ZyqjuteywOHjXXoOGsZ9BxL/djgdVz7FPgh9/EU2hBiFeDacs7tGyP6mWXQz/DG0NRxii0n\nFomEECv+DEyiQgoUbWBp72rPiKrjBqw47Zrm8vzUDuAmxPJ3e24Qajc1d90sJutQ7VTe+YLNzrJ2\nZTInzcwCUFmWbyyK4r4X85qj8JfvKIzhsOznWnduBp3s1cmcIWSxUq2fxw/N4dgCjnE/4/Rqj8YH\nrrqszx0bYzKbs7E1pkmlxmeuvJyqaVg9dZr122/DNzXOOLJ8if7KFQTfYPM+mYGaETYswWgTm/Vi\nFwUfp9ayXp+6ruLah/dgcvorV+Drino8xGY54+EGeX+FarQJIdDUI/LBGsb1qEYbLK1dha/GYB0u\ny52TfK8AAA22SURBVLGunR6zGJOl7wP4pqKpK6yN03TGxqm6Xp5D2itkzlnzqesmVQ3G7MKHkCqm\nA3lmqJtAz9J9X+O0WHyeTVNf3seD/1rWxH1YvSxmU03jsZnr/pym/8zqxmMzm+KH2RmHfCBMVfN1\nCZKJAcykEvJ2qs507xp3L6Wndtqqwum/W9NThnctgNi5Kfdi/u1N79laJCvLve4HkqPkMf/hvjN5\n3yNVhDDvTXTTlT6L7mKv1Zh4hLdPAaSXW6rcEXxcB8qcZamfM6o8aytjQliJveJStdqn7xgxqj2b\nWxWnT/U5tdLnM3dsMx43rKyssbGxwfb2kMHyGtbdl607P8Xq6WsZbt4KAapqSLV9NVlviTs+9RGq\n4Z24vM/qFffCuJxq607Gw3WyfInly65m87ZPsLR2BaPN2+gtrWFtzmVXX4frLRPqEVf8u+upRtuM\ntm6nv3SK5cuuxvWWU4CK60QuHzBYvZLeYBWX93FZTt5bYjAYMBgMqKoKay2nVgacWunRyx2Dfsba\ncly7Wl3uEbwnyyynVnqppN1zaqXH1rDm1ErO1rDGuZj5tHuWQiBNd8YqxDYLzZyhSVV2tYflzKZA\nARjI7aScvOdiEIxVdHFtKrTzZAbS2RxdeXgrbkCm28AL8RgPE9oMZZLl+Km9TCFMhaap95ve6Hvu\nGlBb7NA2sZ1+bLfXTN9/XP69XkyQ3Nwaz3Akl+7s2fVLfu2ZM2u7PnakApAcfVPr2HF9CMido8kC\nIcTpuUEvTiGF9MT2J92m8eTOxqO9s1idtr41ZnmcMegZ1tfBuRxf54yXDMZlLA/OMBpuAitsWsNg\naZme+XecWl1mOBqR95fpDdYYr6wx3rwNl/Wo105x6vKrsC4jXHUto607sC6P+3jyZawxrJy6Amsd\n9egKfGgYrF5Gf7CKyxy+HhOCJ8975L1l8sEKxvVwLqM/GLA0GNDv9zCmT+MbBr28q5Rb6mesLucY\nIM8MvTxmU7kzDHou7t4xJk69GcPaSk5Vhbg+5GKO0u436uXnlMbbSQVa24aonR4F03WCmP5zsjZW\nzaV2EJgw+XNsF3OMnWQjYepx302JtcFqUiV3bnei/XzInu85O5u37u81sjgUgOSi7PhAxGBt2syZ\nx+KDlcZzOsBwVONDYHO7Zlz7uHH0Pp47tyqsgVvObnJ6rcdw3DCuPHXTMKoa6jp0PxXXjefs7dtd\nIcPt6yPuccVy2l8U+NDHbueKtT6BuDF0aZCzvjViY7vmilMDxnVsC5Rnll7uaBrP7RuxyqyfO/Lc\nMsgzqqYhc5bTawOW+hnrW2O2Rw3XXLmcjtU29PJYBBA8LA1i1/ClfkbdeJb6GStLefppPMSzlDbH\n9FKQXV3u0c8dVeNTYUY8prst7HApmsQWQnEqL06vtR0YTOoyHjNQ1+3/ib+0GUTmLG3/uskajcFm\nk+mwWNjQVtBN3mY68Ngu6Ew2/Z57wu2kkexkzedig8UlNDaXBTPTAFSW5UeBh87ya8i8pXUPwJv/\nv717DbWsrOM4/l1r7X0ucxG1puyiGUHPGylf2D1LqIESDJNoojLHsPCFL/RNZKkgXaSsDLKLGUmG\nWXZFiDQQTSiiKIIk+kdFoNhluk2TNp7L3r141tp7nzPnjHOOc+Y5c+b7AWHW3uus/ewzzvqtZz3P\n81/D9ip82N5CGjA7nasI5BNnw+xCfvDd00+epanyyfjg3CJzC21IDPLTVodDWBgMmJlqODiXKygw\nzNOTZ6Zyzbn6jJOZm19k23SPg/OL7JidYttMw0lzi+yY7VNX1Whxa78trzMzVdPvN+3zj3Kvoarz\ng/h2bpti+2yPmak822/bdI+pdm1Pr87HGA5hZrpH09T0m/EZND+AL/8u+r2ak7b3qalo+rnUTw6N\nOi8g7dWjOnYw7o20Qy7trLU8tXqiU9P2QNrRl9E6oLaXM7FOp6qWzDcY/02Nekc55IdV1U5EGLdl\nhb/d0bqi5SanXNtT0XrYA9JTMnkV3e/VDJt8W2iuNxiVsllcHI5OllTbeWI+j6gvLAzyVf8gP3G1\nu59f1/V4Flc/F0BdWBxf+Td1xYHH59kx2+OJucUcBm3ALCwOmJsfMD2V68HNzw9yqR1y7beqqth1\nygz7H5tnfn6Rpsm9m36T5yQ3Td1OQ87FcJp21l+udtBut2f27rW6HvdoKnLvbaqXg25+IX/XYTux\noN+bKBzX6gKju8WVv0f7u6jaKgqjfSdmmTWHnvQng+DQTwKolrzR9dq6zz/Uyu/Ze9HRYADpqKi6\nwYFqPLNpMW+2g+HdSTavfxkC9KrcM6phfpDvLw2paRpoht203ZrZ6YqDc4NROZumrtg+26Oua2Zn\n8uB6ry0q2t2u6jd5QWXOlYp+L1/rd8fYNt0w3zB6dEO3/qduqrw+hmo0E63Xy2uUuoftjUqVtvs1\nbS+qrqs8xlIN20WjQ5qqnbpcVwwXx72dcSCPx2xG4zLkXsrSLszSPs1kED2V3seT/aw9G22karPM\nIhlugoYcL7NqjoZj/V278YduIH3QVr4eDLs1KeP1RsOJUe4uIKqqm5lVjapht099OOTu0XA4Lpi6\nsDAY3dLqAmW0VoWJ0/pwYoIF4+cJVctuZA0nPmP5uEi3kHOlHsNqQbGWADmSff1/ePOp1pDi+/Yd\n2PxfaI127dq56vfffBPOtXV1V/vt5mjge9lpvpsa3D2XZ/yE02rpuhIYBcfk7avRDC668ZHhkmPm\nfboQqiYONh7T6J4x1C3XHLdnvD1c1p7xgQ5/a2y97I1oq/EWnI6J5T2VZpVBhGbZSXZyv6Ya9zwA\nDh1OWfYhFeTFpO3Y0bK969FnLfvMw5znV+pxHcl7qx/PUNGJyx6QjitrPV97gpc2LwNIW54hJG1O\nBpAkqQgDSJJUhAEkSSrCAJIkFWEASZKKMIAkSUUYQJKkIgwgSVIRBpAkqQgDSJJUhAEkSSrCAJIk\nFWEASZKKMIAkSUUYQJKkIgwgSVIRBpAkqQgDSJJUhAEkSSrCAJIkFWEASZKKMIAkSUUYQJKkIgwg\nSVIRBpAkqQgDSJJUhAEkSSrCAJIkFWEASZKKMIAkSUUYQJKkIgwgSVIRBpAkqQgDSJJUhAEkSSrC\nAJIkFWEASZKKMIAkSUUYQJKkIqrhcFi6DZKkE5A9IElSEQaQJKkIA0iSVIQBJEkqwgCSJBVhAEmS\nijCAJElF9Eo3YDNIKdXA54AXA08Al0XE78u2amOllF4GfCwizivdlo2SUuoDXwbOBKaBD0fE3UUb\ntUFSSg1wK5CAIXB5RDxUtlUbK6X0DOAXwO6I+G3p9mjt7AFlFwIzEfEK4P3AJwu3Z0OllN4HfAmY\nKd2WDfZO4B8RcS7wBuDmwu3ZSBcARMSrgGuAj5RtzsZqLy5uAf5Xui1aPwMoezVwD0BE/BQ4p2xz\nNtwfgItKN+IY+CZwbfvnClgo2JYNFRHfA97bbj4P+HfB5hwLnwC+ADxauiFaPwMoOwnYP7G9mFLa\nsrcnI+LbwHzpdmy0iPhvRBxIKe0EvkXuGWxZEbGQUvoK8BngjtLt2Sgppb3Avoi4t3Rb9NQYQNl/\ngJ0T23VEbNmr5RNJSul04H7gqxHxtdLt2WgRcQnwQuDWlNL20u3ZIO8GdqeUHgDOBm5PKZ1Wtkla\njy17lb9GPybfQ78rpfRy4NeF26OjIKX0TOCHwBURcV/p9myklNLFwHMj4gbgcWDQ/rflRMRruj+3\nIXR5RPylXIu0XgZQ9l3yFdVPyGMFlxZuj46ODwCnANemlLqxoDdGxFYcuP4OcFtK6UGgD1y5Rb+n\nthAfxyBJKsIxIElSEQaQJKkIA0iSVIQBJEkqwgCSJBXhNGwVl1I6E/gd8Jtlb10QEQ+vsP9e4Dbg\n7RFx58TrVwI3Ac+PiD+llIYRUS07/hCYIpdwuTQiHkkpnQF8llzCpm73uyIi/naYNjfk2nLnkqfu\n3xoRn162z43ArojYe4S/CumEYgBps3g0Is5ew/6PAG8B7px47SJWr4G25PgppRvIJWveTC5qeXsX\nZimlq8l1xg5XL+9S4GnAi4BZ4OcppQcj4pftMV4H7AW+v4bvJJ1QvAWn49WPgHO6cjMppecBB1ha\n0+9wHiSXrAE4Ddg28d7NPHnl7IeA6yNiEBGPAX8ETm/bciq5GvVHj7At0gnJHpA2i2enlH41sX1H\nRNx4mP0XgHuB88lVr98K3AVc/2Qf1Jby30MuwQRwNXBHSul64D7gB+2xVtVWTe+O90rgpcDF7Uu3\nAB+kDSRJKzOAtFms9RYc5JB4DzmALiSH0WoBNBlw08DPyM9+IiLuSSk9BzgPeD3wceBt7TEPK6X0\nWuDrwDsi4l8ppcuAhyPivnasStIqDCAdz+4nV30+C/h7ROxPKa2274oB194uuzYiriI/E+qelNKH\ngD+nlHZFxL7VDphSugj4PLAnIh5oX94DPKsNu1OBHSmlm9rjS5rgGJCOWxGxSK52/UXgG+s8zH7g\nTSmld0289gLgr8A/V/uhlNJLyOGzeyJ8iIjdEXFWG3bXAXcbPtLK7AHpeHcXeezl7vX8cEQsppTO\nBz7V9nweJ0/RvqANuNVcQ/73c/tEr+u6iFhXO6QTkdWwJUlF2APSppRSugq4ZIW3Ho2I849RG/aQ\nZ8gdYh0TJiQtYw9IklSEkxAkSUUYQJKkIgwgSVIRBpAkqQgDSJJUxP8BtNwP9h8JCUAAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11763fe48>"
      ]
     },
     "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": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "good=XID_MIPS['F_MIPS_24']>20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "149329"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "good.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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": 13,
   "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": 14,
   "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": 15,
   "metadata": {
    "collapsed": false
   },
   "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": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x11b9c9630>"
      ]
     },
     "execution_count": 16,
     "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 0x11a861630>"
      ]
     },
     "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": 17,
   "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": 18,
   "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": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 8480 tiles required for input catalogue and 20 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_SWIRE.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_SWIRE.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 3",
   "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
}