{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mc741/anaconda3/lib/python3.6/site-packages/mpl_toolkits/axes_grid/__init__.py:12: MatplotlibDeprecationWarning: \n",
      "The mpl_toolkits.axes_grid module was deprecated in Matplotlib 2.1 and will be removed two minor releases later. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist, which provide the same functionality instead.\n",
      "  obj_type='module')\n"
     ]
    }
   ],
   "source": [
    "import pylab\n",
    "import pymoc\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import pylab as plt\n",
    "from astropy.io import fits\n",
    "from astropy.table import Table\n",
    "from astropy.table import Column\n",
    "from astropy import wcs\n",
    "\n",
    "import xidplus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook uses all the raw data from the masterlist, 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. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Sel_func=pymoc.MOC()\n",
    "Sel_func.read('../../dmu4/dmu4_sm_EGS/data/holes_EGS_irac_i1_O16_20190201_WARNING-MADE-WITH-Lockman-SWIRE-PARAMS.fits')\n",
    "EGS_MOC=pymoc.MOC()\n",
    "EGS_MOC.read('../../dmu17/dmu17_HELP-SEIP-maps/EGS/data/70101880.70101880-0.MIPS.1.moc.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Final=Sel_func.intersection(EGS_MOC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Final.write('./data/testMoc.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Masterlist\n",
    "Next step is to read in Masterlist and select only sources that are detected in mid-infrared and at least one other wavelength domain (i.e. optical or nir). This will remove most of the objects in the catalogue that are artefacts. We can do this by using the `flag_optnir_det` flag and selecting sources that have a binary value of $>= 5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "masterfile='master_catalogue_egs_20180501.fits'\n",
    "masterlist=fits.open('../../dmu1/dmu1_ml_EGS/data/'+masterfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "good=masterlist[1].data['flag_optnir_det']>=5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "134490"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "good.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create uninformative (i.e. conservative) upper and lower limits based on IRAC fluxes\n",
    "As the default flux prior for XID+ is a uniform distribution, it makes sense to set reasonable upper and lower 24 micron flux limits based on the longest wavelength IRAC flux available. For a lower limit I take IRAC/500.0 and for upper limit I take IRACx500."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_lower=np.full(good.sum(),0.0)\n",
    "MIPS_upper=np.full(good.sum(),1E5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#masterlist[1].header"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for i in range(0,good.sum()):\n",
    "#    if masterlist[1].data['f_irac_i4'][good][i]>0:\n",
    "#        MIPS_lower[i]=masterlist[1].data['f_irac_i4'][good][i]/500.0\n",
    "#        MIPS_upper[i]=masterlist[1].data['f_irac_i4'][good][i]*500.0\n",
    "#    elif masterlist[1].data['f_irac_i3'][good][i]>0:\n",
    "#        MIPS_lower[i]=masterlist[1].data['f_irac_i3'][good][i]/500.0\n",
    "#        MIPS_upper[i]=masterlist[1].data['f_irac_i3'][good][i]*500.0\n",
    "#    elif masterlist[1].data['f_irac_i2'][good][i]>0:\n",
    "#        MIPS_lower[i]=masterlist[1].data['f_irac_i2'][good][i]/500.0\n",
    "#        MIPS_upper[i]=masterlist[1].data['f_irac_i2'][good][i]*500.0\n",
    "#    elif masterlist[1].data['f_irac_i1'][good][i]>0:\n",
    "#        MIPS_lower[i]=masterlist[1].data['f_irac_i1'][good][i]/500.0\n",
    "#        MIPS_upper[i]=masterlist[1].data['f_irac_i1'][good][i]*500.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#np.savez('./tmp_mips_prior', MIPS_lower, MIPS_upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "npzfile = np.load('./tmp_mips_prior.npz')\n",
    "MIPS_lower=npzfile['arr_0']\n",
    "MIPS_upper=npzfile['arr_1']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([   51.6788325,   442.452371 ,   333.393812 , ...,  1654.2      ,\n",
       "        8876.5      , 21136.5      ])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MIPS_upper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Map\n",
    "We are now ready to read in the MIPS map\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename='70101880.70101880-0.MIPS.1.help.fits'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_Map=fits.open('../../dmu17/dmu17_HELP-SEIP-maps/EGS/data/'+filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "w_im:  WCS Keywords\n",
      "\n",
      "Number of WCS axes: 2\n",
      "CTYPE : 'RA---TAN'  'DEC--TAN'  \n",
      "CRVAL : 214.708480552  52.8132907943  \n",
      "CRPIX : 1451.259  1553.104  \n",
      "CD1_1 CD1_2  : -0.000680556  -0.0  \n",
      "CD2_1 CD2_2  : -0.0  0.000680556  \n",
      "NAXIS : 2902  3105 \n",
      " w_nim:  WCS Keywords\n",
      "\n",
      "Number of WCS axes: 2\n",
      "CTYPE : 'RA---TAN'  'DEC--TAN'  \n",
      "CRVAL : 214.708537  52.813316  \n",
      "CRPIX : 1451.259  1553.104  \n",
      "CD1_1 CD1_2  : -0.000680556  -0.0  \n",
      "CD2_1 CD2_2  : -0.0  0.000680556  \n",
      "NAXIS : 2902  3105\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: FITSFixedWarning: CRVAL1_0= 214.708537000 /CRVAL1 at Wed Jun 20 00:53:54 2018By: sjo \n",
      "keyword looks very much like CRVALia but isn't. [astropy.wcs.wcs]\n",
      "WARNING: FITSFixedWarning: CRVAL2_0= 52.8133160000 /CRVAL2 at Wed Jun 20 00:53:54 2018By: sjo \n",
      "keyword looks very much like CRVALia but isn't. [astropy.wcs.wcs]\n"
     ]
    }
   ],
   "source": [
    "w_im = wcs.WCS(MIPS_Map[1].header) \n",
    "w_nim = wcs.WCS(MIPS_Map[2].header) \n",
    "print('w_im: ', w_im, '\\n w_nim: ', w_nim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "MIPS_Map[1].header['CRVAL1']=MIPS_Map[2].header['CRVAL1']\n",
    "MIPS_Map[1].header['CRVAL2']=MIPS_Map[2].header['CRVAL2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "MIPS_Map.writeto('../../dmu17/dmu17_HELP-SEIP-maps/EGS/data/70101880.70101880-0.MIPS.1.help2.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAD+xJREFUeJzt3X+MXlldx/H3h64DsoQf2tVo27ElbSoNicGddBESsxGUNkspIWpa1CBptsFYxB+JFGNCjDGuCTGCVEyFWoybNs26gdYdLWaVFJOGtAsmttTNNhXo0EqLi1XRpBS+/jGz7KR2Zp4f99ln5uz79c/MPfPcc783nXx75nvPPSdVhSSpXS8YdwCSpNEy0UtS40z0ktQ4E70kNc5EL0mNM9FLUuNM9JLUOBO9JDXORC9Jjbtr3AEArF69utavXz/uMCRpRXniiSe+VlX3LPW5ZZHo169fz9mzZ8cdhiStKEm+1Mvnxlq6SbIjycEbN26MMwxJatpYE31VnaiqvS972cvGGYYkNc2HsZLUOBO9JDXOGr0kNc4avSQ1ztKNJDXORC9JjVsWL0zp+Wv9/se+8/0XH3pgjJFI7fJhrCQ1zoexktQ4a/SS1DgTvSQ1zkQvSY0z0UtS4zqfXpnkBcDvAi8FzlbVx7u+hiSpdz2N6JMcSnItybnb2rcleTLJxST755p3AmuAbwIz3YYrSepXr6Wbw8C2+Q1JVgEHgO3AFmB3ki3AZuB0Vf068EvdhSpJGkRPib6qTgFP39a8FbhYVZeq6iZwlNnR/Azw9bnPfKurQCVJgxnmYewa4PK845m5tkeBNyX5Y+DUQicn2ZvkbJKz169fHyIMSdJihnkYmzu0VVX9D7BnqZOr6mCSq8COiYmJe4eIQ5K0iGFG9DPAunnHa4Er/XTgEgiSNHrDJPozwKYkG5JMALuA4/104KJmkjR6vU6vPAKcBjYnmUmyp6puAfuAk8AF4FhVnR9dqJKkQfRUo6+q3Qu0TwPTg168qk4AJ6amph4ctA9J0uJcj16SGjfWHaYc0Ws+d5uSRsNFzSSpcZZuJKlxbiUoSY2zdCNJjbN0I0mNs3QjSY2zdCNJjTPRS1LjrNFLUuOs0UtS4yzdSFLjTPSS1DgTvSQ1zoexktQ4H8ZKUuMs3UhS40z0ktQ4E70kNc5EL0mN6zzRJ7k/yWeS/GmS+7vuX5LUn54SfZJDSa4lOXdb+7YkTya5mGT/XHMB/w28CJjpNlxJUr96HdEfBrbNb0iyCjgAbAe2ALuTbAE+U1XbgfcCv9NdqJKkQfSU6KvqFPD0bc1bgYtVdamqbgJHgZ1V9e25n38deGFnkUqSBnLXEOeuAS7PO54B7kvyNuBNwMuBDy90cpK9wF6AycnJIcKQJC1mmESfO7RVVT0KPLrUyVV1MMlVYMfExMS9Q8QhSVrEMLNuZoB1847XAlf66cAlECRp9IZJ9GeATUk2JJkAdgHH++nARc0kafR6nV55BDgNbE4yk2RPVd0C9gEngQvAsao638/FHdFL0uj1VKOvqt0LtE8D04NePMkOYMfGjRsH7UKStASXKZakxrnxiCQ1zhG9JDXOEb0kNc4RvSQ1zvXoJalxlm4kqXGWbiSpcZZuJKlxJnpJapw1eklqnDV6LUvr9z/G+v2PjTsMqQmWbiSpcSZ6SWqciV6SGmeil6TGOetGkhrnrBtJapylG0lqnIlekhpnopekxo0k0Se5O8kTSd48iv4lSb3rKdEnOZTkWpJzt7VvS/JkkotJ9s/70XuBY10GKkkaTK8j+sPAtvkNSVYBB4DtwBZgd5ItSd4IfAH4aodxSpIGdFcvH6qqU0nW39a8FbhYVZcAkhwFdgIvAe5mNvn/b5Lpqvp2ZxFLkvrSU6JfwBrg8rzjGeC+qtoHkOQXga8tlOST7AX2AkxOTg4RhiRpMcMk+tyhrb7zTdXhxU6uqoNJrgI7JiYm7h0iDknSIoZJ9DPAunnHa4Er/XRQVSeAE1NTUw8OEYdWGNeZl55bw0yvPANsSrIhyQSwCzjeTweudSNJo9fr9MojwGlgc5KZJHuq6hawDzgJXACOVdX5fi7uWjeSNHq9zrrZvUD7NDA96MWT7AB2bNy4cdAuJElLGKZGPzRr9FrK/Hr+Fx96YIyRSCuX69FLUuNcj16SGueIXpIa54hekhrnevSS1DhLN5LUOEs3ktQ4SzeS1DgTvSQ1zhq9JDXOGr0kNc7SjSQ1bqyLmkn9cIEzaTCO6CWpcSZ6SWqcs24kqXHOupGkxlm6kaTGmeglqXEmeklqXOfz6JO8CngPsBp4vKo+0vU1JOfUS73raUSf5FCSa0nO3da+LcmTSS4m2Q9QVReq6l3AzwJT3YcsSepHr6Wbw8C2+Q1JVgEHgO3AFmB3ki1zP3sL8I/A451FKkkaSE+JvqpOAU/f1rwVuFhVl6rqJnAU2Dn3+eNV9Trg57oMVpLUv2Fq9GuAy/OOZ4D7ktwPvA14ITC90MlJ9gJ7ASYnJ4cIQ5K0mGESfe7QVlX1aeDTS51cVQeBgwBTU1M1RBySpEUMM71yBlg373gtcKWfDlwCQZJGb5gR/RlgU5INwFeAXcDbO4lKzZk/HXKUfTvVUvr/ep1eeQQ4DWxOMpNkT1XdAvYBJ4ELwLGqOt/PxV3rRpJGr6cRfVXtXqB9mkUeuC4lyQ5gx8aNGwftQpK0BFevlKTGuR69JDXOEb0kNc7VKyWpcZZu1JT1+x8b6VROaSWydCNJjbN0I0mNs3QjSY2zdCNJjbN0I0mNM9FLUuM63xy8H651o1FxRUvpWdboJalxlm4kqXEmeklq3Fhr9NJzwXq9nu98YUqSGufDWElqnDV6SWqcNXo9r1iv1/ORI3pJapyJXpIaN5JEn+StSf4sySeT/NQoriFJ6k3PNfokh4A3A9eq6tXz2rcBHwRWAR+tqoeq6hPAJ5K8AvgA8Kluw9ZK4bZ+0vj1M6I/DGyb35BkFXAA2A5sAXYn2TLvI78993NJ0pj0nOir6hTw9G3NW4GLVXWpqm4CR4GdmfUHwN9U1efu1F+SvUnOJjl7/fr1QeOXBuZG4nq+GLZGvwa4PO94Zq7t3cAbgZ9O8q47nVhVB6tqqqqm7rnnniHDkCQtZNh59LlDW1XVh4APLXmy69FL0sgNO6KfAdbNO14LXBmyT0lSh4ZN9GeATUk2JJkAdgHHez3ZtW4kafT6mV55BLgfWJ1kBnh/VX0syT7gJLPTKw9V1fk++rR0o7Fb6IGsSySoFT0n+qravUD7NDA9yMWr6gRwYmpq6sFBzpckLc3NwaUFuACaWuF69JLUOBc1k6TGuZWgJDXO0o0kNc7SjSQ1ztKNJDXO0o3UJ1e91Epj6UaSGjfWF6aklcIRvFYy34yVBuSbs1oprNFLUuMs3UgdcHSv5cxEr85Zz5aWF2fdSFLjfBgrdcwyjpYbH8ZKUuOs0UvLkH8VqEvW6CWpcY7opTFy5K7ngiN6SWpc54k+ySuTfCzJI133LUnqX0+JPsmhJNeSnLutfVuSJ5NcTLIfoKouVdWeUQQrrTTPLGnsS2Qap15r9IeBDwN/8UxDklXAAeAngRngTJLjVfWFroOUWjJo0h+0nu9zAPU0oq+qU8DTtzVvBS7OjeBvAkeBnR3HJ0ka0jCzbtYAl+cdzwD3Jfle4PeA1yR5X1X9/p1OTrIX2AswOTk5RBjSytDFSF4axDCJPndoq6r6d+BdS51cVQeTXAV2TExM3DtEHFITTOgalWFm3cwA6+YdrwWu9NOBSyBI0ugNM6I/A2xKsgH4CrALeHs/HbiomdQeH/4uP71OrzwCnAY2J5lJsqeqbgH7gJPABeBYVZ3v5+KO6CVp9Hoa0VfV7gXap4HpQS/uiF6SRs9liiWpcW48Iq1AC83QmV8TdxaPnuGIXpIa54he0h1H/86YaYcjeklqnOvRS1LjLN1IDVnOD2AtD42PpRtJapylG0lqnIlekhpnjV56Hul6d6vl/ExAz7JGL0mNs3QjSY0z0UtS40z0ktQ4E70kNS5VNb6LPzvr5sGnnnpqbHGoW87EUJf6meXzXL9pO+i2iV1tt5jkiaqaWupzzrqRpMZZupGkxpnoJalxJnpJalznSyAkuRv4E+Am8Omqerjra0iSetfTiD7JoSTXkpy7rX1bkieTXEyyf675bcAjVfUg8JaO45Uk9anX0s1hYNv8hiSrgAPAdmALsDvJFmAtcHnuY9/qJkxJ0qB6SvRVdQp4+rbmrcDFqrpUVTeBo8BOYIbZZN9z/5Kk0RmmRr+GZ0fuMJvg7wM+BHw4yQPAiYVOTrIX2AswOTk5RBhaDnxJSqPSz+9Wl7+HC73INK54hjFMos8d2qqqvgG8c6mTq+pgkqvAjomJiXuHiEOStIhhSiszwLp5x2uBK/104JuxkjR6wyT6M8CmJBuSTAC7gOP9dJBkR5KDN27cGCIMSdJiep1eeQQ4DWxOMpNkT1XdAvYBJ4ELwLGqOt/PxR3RS9Lo9VSjr6rdC7RPA9ODXtw9YyVp9Fy9UpIaN9ZEb41ekkbPEb0kNc4RvSQ1bqxbCX4niOQ68KVxx9Gj1cDXxh3ECHl/K5v3t3INcm8/VFX3LPWhZZHoV5IkZ3vZo3Gl8v5WNu9v5RrlvbnomCQ1zkQvSY0z0ffv4LgDGDHvb2Xz/laukd2bNXpJapwjeklqnIm+DwvskduEJOuS/EOSC0nOJ3nPuGPqWpJVST6f5K/HHUvXkrw8ySNJ/mXu3/DHxh1Tl5L82tzv5bkkR5K8aNwxDeNO+3An+Z4kf5fkqbmvr+jqeib6Hi2yR24rbgG/UVWvAl4L/HJj9wfwHmZXWm3RB4G/raofBn6Ehu4zyRrgV4Cpqno1sIrZZdFXssPctg83sB94vKo2AY/PHXfCRN+7hfbIbUJVXa2qz819/1/MJoo1442qO0nWAg8AHx13LF1L8lLgx4GPAVTVzar6j/FG1bm7gO9OchfwYvrc5Gi5WWAf7p3Ax+e+/zjw1q6uZ6Lv3Z32yG0mEc6XZD3wGuCz442kU38E/Cbw7XEHMgKvBK4Dfz5XmvpokrvHHVRXquorwAeALwNXgRtV9anxRjUS319VV2F24AV8X1cdm+h7d8c9cp/zKEYsyUuAvwJ+tar+c9zxdCHJm4FrVfXEuGMZkbuAHwU+UlWvAb5Bh3/2j9tcrXonsAH4QeDuJD8/3qhWFhN974beI3e5S/JdzCb5h6vq0XHH06HXA29J8kVmS24/keQvxxtSp2aAmap65i+wR5hN/K14I/CvVXW9qr4JPAq8bswxjcJXk/wAwNzXa111bKLv3dB75C5nScJsjfdCVf3huOPpUlW9r6rWVtV6Zv/d/r6qmhkRVtW/AZeTbJ5regPwhTGG1LUvA69N8uK539M30NDD5nmOA++Y+/4dwCe76rinrQQFVXUryTN75K4CDvW7R+4y93rgF4B/TvJPc22/NbddpJa/dwMPzw1CLgHvHHM8namqzyZ5BPgcs7PDPs8Kf0N2bh/u+4HVSWaA9wMPAceS7GH2P7ef6ex6vhkrSW2zdCNJjTPRS1LjTPSS1DgTvSQ1zkQvSY0z0UtS40z0ktQ4E70kNe7/AKepyUv3ydfzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data=MIPS_Map[1].data\n",
    "plt.hist(data.flatten(),bins=np.arange(-1.0,10.0,0.1));\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3071634754641976"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(~np.isnan(data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7ff083b740f0>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAD8CAYAAAC8aaJZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnX+sZVd137/r3vfmh+3JPI/zQK5tFYdOIOaPgvsKlqiQktR+F/cPE6mtQCqMM280rmIkSFKpJvkDmqgSoSW0VEA6zIxsIoRBhQgrJW9wXSqUNoCfkTE2lvHE0HqwZU8YzzDGeOa9e1f/OHufu84+a+9z7u9f6yM93Xv32eecfe87a++111p7bWJmGIaxODQm3QDDMMaLCb1hLBgm9IaxYJjQG8aCYUJvGAuGCb1hLBhjF3oiahHRU0R0mojuGff9DWPRoXH66YmoCeCHAG4FcAbAwwDew8w/GFsjDGPBGfdI/1YAp5n5GWa+DOB+AHeMuQ2GsdAsjfl+1wF4Vnw+A+BtsgIRHQVwFACuvPLKf/TGN75xfK0zjBnmkUce+TtmXq2qN26hJ6WsML9g5mMAjgHA2toab21tjaNdhjHzENH/rVNv3Or9GQA3iM/XA3huzG0wjIVm3EL/MICDRHQjEe0C8G4AD4y5DYYxEK2VjUk3YSDGKvTMvAPg/QBOAXgSwJeY+YlxtsEYHa0DRybdhJGzvv8wAKB1zdGZFf6x++mZ+WvM/KvM/Hpm/vfjvr8xfNb33ZkJAzXQuubopJszcpgZ4A6A2ezoLCLPGIjWygao2QQRAe02wJ2ZFIQ6tFY2QEQgInC7AzQI6DDW993Z13f2WsO4MaE3+qa1spE9+AJud8CXtyfUotHROnAk/67MnHVynczxRM0m0OGe1P2ss2xMpIM0oTd6prV6V1fgO1wQfCIC7Vqeq9G+tbKRfc8OdwU+xP0Gdb53PgXqMHh7Z+wjvgm90ROFObsb6UDBY+QEZJ4EHwDQoKLANwh5GLvrELjdQWtlI2rbaK1sAO12fh4tL4GIxir4JvRGbdb33ZkZsLZ3igfa7e77BnX/MJuGLolX671AA+hqNp3iqO/n+2gQ+PLlkuBL9d/XQ4eBZnOsqr4JvVGL1oEj+dy1QDCn96O8r8ftzsQMVoPSWtnIBZ2aja6Ay98g/P4OajbzawCi8xMdotcO0G53taMxuAFN6I1K1vcfBoejO6DO6UPy0W/GWN93Z7Gg2URhRequ5exVdAChqp+79bxNwCPe5yO+MBKOWvBN6I0kuZvKjVyqgIejf1jeoJlT86nZzObwzUb2PbZ3up1Xg4DAQ8HMoGajO1ePdXTKb8XtTtcTQAQ0myMVfBN6Q6V14Eim3ob5FmIC7mlQ8RXdh3pWItjydvppijfYaR3e8lLXwNeJWPZ9hxkhH+09zkYyqt/LhN4osb7/cFe4O53CsUInEKr2Depa8kMV1jHtEXtaR8ftTsEvX+j4/Hw8pQFJQ+dyZGFrpDMdheCb0BsF1vfdmUWc+Qe/kT0iUbU1nKuGlvyQdntqR3zfrvA7Fn6PGM4KX0Ab4TXbSM12DQsTeiPHh9QC3Qefmo3CK4BqAZhxtO9XUsGBTO1vdzUh3tkpns+dpJGzl3YMU/BN6A0AzkIvH3bhVpKv0Yi0kAYVDHnhqDdto33Jhw6ov0c4vZGdYRi4I12X8hrJjkA55u85LGOoCf2C01rZyOPACw97+LC6h1FVdcMHNfTnh2q/O2daLPqtA0fK83itc3OGupJdAyiFIhfKIx1oCdFRyvvL+IBh2ERM6BeYwmgbMbyFx0qhqOFxQBXwEp3hjl794u+vzeMBlA1v3o3n3udogiyFvY6aH/wPfOfCzN2OZgg2ERP6BcU/ONH5eUwVrfMAJ46XRslJx+hTo7tMNqRBXcObNkqnohNDga9ydWpNEyO97JQGDeAxoV9ApJCV5qH+NRTu1MMOpH3RcgRrKo9cHwIxDForG+CdnaI6LlFU90J5eEzWl/P52HViZQFhx1zLm5DAhH6B8PP3AtpDWzX3DM8DMpXeq8Kp8yLHSmGvI6bwO0QET4bV5p937y5WSs3PU5+rzhfnhFOPQQ17JvQLQmv1ru6HcERKURViK0n5oCtGNFpeGptFv5T8I6KlhCMqEQGXLqUvHv6m3pBXY5GOxBtTNQMjOlmw0Oa545XXUZvY11nGTNE6cKQokH5ErggP7RvNHlCzcxmbK0/GugeGR2k4i7on6/x23lUpp0o15/cly70sbzSwef5E9f1jzer7TGMmaB04UowfB8rGqRDhV681d5RuKXndlEU7fC/KRin46rWVdoar30r4wBt/XHYCYRx9rwE64X3F/47bbZy6cLK364WXH+hsY6oJXXJ+BVeOH+HCkdn71RtU/fD7+imrfg2DoBoYNGRkQowY3j+exyqkpjfyuP8t/cge6/jqGC3D+wqbwqmL91afX4EJ/RzSOnBET1zRoGyEClXT8AH1x2V5lcEpNW+V5X4dekAhQcUIVuTlLkq/eEauApTzdh92HKr1dd2UfmSPuTwbNTpRBSIaeITPmzCUqxhTQ67OA7qrqeaiGAD1563+urK+tHhLdnZqC1A/1mltE4rWgSN527TsN4VOIOZ7r7LShx1nmDcQAJZ3lY/J30IJBPIMMocPGfcGlsYIaV1zNBc+8oJbx4AWMy75a+zeXW21bjbzTDE5DQIFRVjelb5W0JbWykatBz4X9AYB3hPQbIIvX86ST16xF3zx5cIimTxBBoQGEBJMQSqjEfOydlnb2b7sLiRy7cnfXybMFLEDwxR4wEb6uaF1zdGucamuwHuq3EkpIfUPLXfK9wut1vJaMT+2ol1U+fDzxUJeeLyhkjt5Xj+++HLeJr/OQK6Qq2Olj9bRjJOy84rN74NpVKEzQtbJDFvggQFHeiL6MYCLANoAdph5jYgOAPgigNcB+DGAf8nML1H2i/1nALcDeAXAncz83UHubxRHOG53sgem2Syq6NJglxrZ67iTtGvkqn27XL+uQctfU7oWXRkhPs3IpjOdctJO7bsEkXK1VguG6wg0mk2gs1N8X/W9lc5NTj2YeWhz+JBhjPS/zsxvZuY19/keAA8x80EAD7nPAPBOAAfd31EAnxnCvRcavzrMp2fOVUL5oMr32oPYiy9dnhNa66vsBHXQbBHefx+b33cYtGtX/FpBmZr2KmW70I6F58uOKgxQihny2u1kMNOoBB4YjXp/B4D73Pv7ALxLlH+OM74FYIWIrh3B/ReC1oEj+cie769W1+0VsaAXiAkCNeKjmJIhViVQ5aN5+Fw9uYmEJE9tFdoSEqije+i61I5p7esF6eGomE6NQqUvNGXA8xnA14noESLyC31fy8zPA4B7fY0rvw7As+LcM66sABEdJaItIto6e/bsgM2bT7yFXhp7SoKjWY+BTNBkJtfYAxw+7F6gU+pubA+7cM7r7+mEtSpzrFxl5kd8H0obXSzj76UEuvQUWizP27XcuxYTiUuIMWqBBwYX+rcz883IVPe7iegdibrar1X6FZj5GDOvMfPa6urqgM2bPwo51MUDSM1GcQVbbPSLrXWPGdb8q3S1aX5o/7p3T/k6sQe/Ssg0OmIprjxf00y0+ILIPWWnqU4BOpx1av2M8kDWCWudRqcb8jsOgQcGFHpmfs69vgjgLwC8FcALXm13ry+66mcA3CBOvx7Ac4Pcf2EJfcN1DXAaVQYnLXBHRuCF5/3iVf06vZIavVPuxaprRX6HcHuqyg6lDv685aWsE97ZKbXduwHHJfDAAEJPRFcS0T7/HsBtAB4H8ACAQ67aIQBfde8fAPA+yrgFwAU/DTDqsb7vTj2ENOYrrmOk60VdDY2E8h6aEU6LPIvE29dti3drdbP19hDd1m/nWMeCryFdqLKjFH+0tDRWgQcGG+lfC+Cvieh7AL4D4L8z8yaAjwK4lYieBnCr+wwAXwPwDIDTAD4L4HcGuPdCcurivaClpeKDnxLsToVg9CIEVZqC9rmueh3TNpS2aHnj+opHGKROXXzbtA7RsfnTY8O7X01omtMZr62t8dbW1qSbMXXIlXN9PaRuy6bcZdTvdUJ68fNXHeulTf22P+XL7/X+QPy88Hu5usMe4YnoEeE6j2IReTOIT56Q53ZLqcxaWbut++/dFk3RZaJV1B1JNat6eH4vanikLl15hX6edJ+Fx8LFRim0TEEpTca/p8HWww+KCf2Mkj80ofoYE9LlJf3hDIVHzttllFsoqL3Mo0NCtVerm7q+FNiYka1B4J+/Er+3Rofr70DTbOqBOBEbRt4BNWgiKr3EhH6GUXdd8YQP9vaOLrCyXujmC4NWpMBoQhcT8LqEnVLqu0mtJXbfmsk6+0KLQNTm7z6OwnVA/aa4GiYm9DPM5vkTZQt66mGuMphpri0ZRaaNxFLo+hCkqE0pjKWX7N2jew5CYlb3sDPoNeAmPC8h7PJ3nKRKLzGhn3E2z5/orhZLqLo9o3UQMSNbalROtENdquq/QyrNlI8FqNPJaB2VXNraQ3tLaJ2k4irdPH9iKkZ4jwn9HJAvzkipuoB+LDbqperKB7sXw1dQJxkHPyy0TksLvIlF9Gn0YM+YltFdYkI/J2yeP1Eegf2D7K3SdSLZQsGQ2Vw0i39qC6uYcNQVmpjQyTaFnU+v+MCZOqsEU9MJRfOZptFdYkI/R2yeO158YGUyjVALqLLC+/LtGumtYqRsB3XvrxGuuQfi30/hlbe/QW9fbD4OZVPPRJuZeWoFHjChnztKgu/R5uHyNaTK9xwTsCpXWur+8jqxPHIa2pLemJ2hQbjifz+VbpPyWd2Oy+O/q2vzKNfCDwMT+nmkw2Cfb02q9uHIJYn4l1ViHYimZVTRbIKu2Fu8jsz2I8tjxJb0VqnhGlXfWytzWgbv7EzlHD7EhH4O2Tx/AtRs6pl0XNRdbkTT5qn9hrWmzokJU7sNfuUXxbIekmIU8HP9MKssULRvpNpTR/OR53s/fLsz9SO8x4R+TvEjTmlE395Jq+5anHhILx1CP3P2mMoP6ALt8XN9LaounPMP2qkFv9GsCDxgQj/XlFRN6mZZjar3nWBKUIfU1CBmO6h7bS2yMNaGOn71fiPxGlTqcMad/GJYmNDPOYURyKn5Mv2USjglqEKbGvRr8e+XmEAPqx3UKHU4g+wcO0lM6BcAbSQqJdEcNBZ9kOuEghmm3NLq+uXBVcQ6oNQ0QUN0grM6wntM6BeEPHjHUSvnuyemOkujWJhht9+FNkA55ZZHCOqr//j1/V+/Qfo0oWaADzUnuzR2UEzoF5TCSF8loFUW7Q4Xk0amOogqZD0phH4pq7vHnm/9cDDtRGtPxZQmH+FnUKWXmNAvEPnD2qBisEm/anIv1v0qj4B2rhTC2CKZmlF4hTo9dhayg5zlEd5jQr9g5Mtx6wqhpMICX5l6LUzk0c+9U/EEHY7P1cN4hB6iDP1UaJbccilM6BeQfLSKrZWPUeFzD+0EpU4gtrjHU9c45++pCW7MpafVDebwckNLWX9UG0lOChP6BaX0EGsCiCDPXEw7iKjLtbZ1lvgVbzHB79XiXkW7Hd+VVrRhXkZ4jwn9ArN57nh5NJZJLIB4nrkaVMb3x4h0QOrqul6umzLeyX31hBYxTyO8x4R+wfGjmBdQ3qmZGLIGPY/0QG9Wde9K7DW6z0ccyutwB1heGvne8NOACb3RTVvVoKiKOxCjis7rNaw2tgrQX0e4BOdNpZeY0BvlB9znv4+p0DEhVoJbaiWfkNRd+hq2oU5gTVUOAces++GrMKE3AHQTbDJ3R7zCwhyf1CKlTitqeE+Rf/4cv75eIya4qcCaca8DmHIqhZ6IThLRi0T0uCg7QEQPEtHT7vVqV05E9EkiOk1EjxHRzeKcQ67+00R0aDRfxxiEcMQnEkE8MvVWDfo24gHgS5f1AyMWXt/mfCvsOaXOSH8vgFZQdg+Ah5j5IICH3Gcg26f+oPs7CuAzQNZJAPgwgLch2876w76jMKaL1Fw2KsiKMNYy4iUSa5QIA3t6IZE6S5JvHtJhtFY2+rvXDFAp9Mz8TQDnguI7ANzn3t8H4F2i/HOc8S0AK26P+nUADzLzOWZ+CcCDKHckxpSwef5EyXUFoBy6K7LG9EUvQizddaEAh4t9atyrYGuQC4pEvXkV/H7n9K/1e8u719e48usAPCvqnXFlsfISRHSUiLaIaOvs2bN9Ns8YlFMXTubzeikQ+WgvBKafeXuSXuLzgXiOvMT1C22mRjQd1jyq+sM25Gn/LU6UlwuZjzHzGjOvra6uDrVxRm9snj8BIiqM5H353nsllo+uV6oy9nrknn11YvRnnH6F/gWntsO9vujKzwC4QdS7HsBziXJjBiDqrpcP5/WVi2yGYXzr140XO69qdWBwvfX9h2s0cnboV+gfAOAt8IcAfFWUv89Z8W8BcMGp/6cA3EZEVzsD3m2uzJhy8qg0p0ITUUHlr1TtY4KnxdH300FIu0Oqc+iwOm8vaQM+UEfkCKRmY64Ev47L7gsA/gbAG4joDBFtAPgogFuJ6GkAt7rPAPA1AM8AOA3gswB+BwCY+RyAPwbwsPv7I1dmzAAFw54jZryrHPk92r57o0y7XSdcN7aWv5NFLM7L/J5q/5MmwNraGm9tbU26GYbDj3Z+dFd3nZ1mBskF2KA8QGnzp8eG2KjhQUSPMPNaVT2LyDNqc+rCyYLbrnL9/CjQ3HNVG1h4+k3g4c9tt/vfiGOKMKE3ekJdjuvoa9TvNd2V5p7T5vSxdQDavD7WHo05CNwxoTd6Jo9cQ3F0j470u3fHL5babadH8vvvWo6vA6jKElTTRjDLgm9Cb/RE/rB7y7boAEr743kuXer9Rv3uNw+kg3Vi8/qqHXLkeS6t16wa9kzojdq0Vja6o2kq2GWQubOn7g47Vbn863oIvKtu9269zfI87mTeixkN3DGhN3oiFyw/0u+7Kn3CsHPTA7W0gFK4cHi9cBdbj9RKwgw7nk7XazGLar4JvVGL2MPNF1/u74K9JL0Iqdi7PupKlBlufVqwVI5+b7H3+OQiQOF11gTfhN6opO9otNTGFL1skNkjdbwItHtXuoKc4/tXuc130FnMkuCb0BtJWisb/efNq9gcI3nNQfLwQyT6jHkUvLFP5toP7RQ+W5C3+GuhwzO4FNeE3uiNflbA9bpUNlYWXq8TF2yiTEild0FtS2qqIFffhaq+wjRHt0pM6I0ordW7qitNwoId3LOkzkt1PKwf62BSiTZDlT7i3iOajfm9Cb2h0lrZALZ3hjN6DbNj6NUA2O8GG6lRPbVF1gz4703ojRL5aFWVFUdzgaXQ1Oy9e+J1NKTKXQcpoDdcWx2EE1KnfoNKiUOnWfBN6I0C6/vu1H3cGrGlqHXq+/e/eDVep+61YobAsPzZ59MbZMj6odEu9OtXrd3H9Br2TOiNIo0RPRIp912qrM5xrQPwwl21dTVQrKPZAzrcXV2nGTJnLDLPhN7IKbnn6lJHXa4jGFVpsiPXKNgd6s7NpVHOJ/SQWoAWvlsVnqscn8bR3oTeAJAF4Hjh6Vnwtfl8OFcP69R1/dXIj5e75aqEMnZ9LRuuJvyagIf39Z+l/37K5vcm9AbW9x/OdrPpNwvOtrLT7S9e1SPw+hXuANXu0Kua7dvSbuvTAJEnryTk0o0npxGRzmyacuyZ0C84XuBDCvvYpUhZ7WNz7UGIJeSMXTfcklpLcw0A2zvlPfTC4BwAWIrstBOG6AYaADUbaF1zNPXNxoYJ/YITG92jI3+oyqZGcC3AJSzvFU3gUpZ2Oaf3arymujcI/GrgSZCdgxdiba2+FqcftrVGRN+4MKFfYLzKWRmnLomp05owxh5yLc9cr/5zyfZOUei0zTa1OXjY/rAzkBl/Uq6+utpLszkVhj0T+gVFWurD14HQRvEqI54v60X11wJ9QvdZKLQxt5s/N+x0whgCzWgnr+9fY+5I5yWYtOCb0C8gMuIufEAHDrvVRndZJle1efrx2cdGaXlumKbL14tZ+cMY+zC3XzhF0Yx84X0indkkLfom9AtGYZTpcGnTCjna1+4A6qrly0tZBxAIQWnjjH4DXiJr3UvXTn2WaLn9/LU1W0bK7be8q3SvSQm+Cf0CobmNSktP0RX2wtbUKWIBLSGaaw89WOLrEo684Txes0X4Or0aGVPag2zD9uVyGU1G/EzoFwi/D5003BX2aRf1APQufP0Ka4NUraJQpk0LUu3QVG7lvqXOoMrCXqcN0pvgr9vhsruPOxMZ7evsZXeSiF4kosdF2UeI6CdE9Kj7u10c+xARnSaip4hoXZS3XNlpIrpn+F/FSCHVejmyjnxbqgoh2Tx/ApvnjqvtKEw1dnaqO5VQ1ZZ/Wr3gWJ7kU3O/yevK6/jOyI/aoZsQ6GoPOzvFe1MDoMbYBb/OSH8vgJZS/glmfrP7+xoAENFNAN4N4E3unE8TUZOImgA+BeCdAG4C8B5X1xgDmsCPbB86LWbdI1RnZu7uiAsU3mtptntqay9GQAFffLlooEtpB97u0G5nwivdhGHHwJ2iUY8amTbAHYA72Dx3vP53GwKVQs/M3wRQd4fZOwDcz8yXmPlHyHavfav7O83MzzDzZQD3u7rGiIm5h6Tw568JYalt1EuNxmIEPHXhZOmwF/xcwJcrklfG7j2ITSAcpeUUwHcIoWC329muOruWVTsCXXVV8XrtNtBug7d3xi7wwGBz+vcT0WNO/b/alV0H4FlR54wri5WXIKKjRLRFRFtnz54doHmGVBtL208FG0H6LZ+q9qkbxKXnz5WjekjhWGpnHC9cdZbOptT1MDQ3PCYFPLV55eXt7E9xJfLLL3fbKuwNpy7eG7/eCOlX6D8D4PUA3gzgeQAfd+XaUMGJ8nIh8zFmXmPmtdXV1T6bZ7RWNnJXWEGV96vS3LFQba5Soyt3qk1oC5qnQGPz/Im05iFdetIjoFntgaLP3M/DZacXCwwKrfKpwBv/GnY2y0t6WzFjLjtmfoGZ28zcAfBZZOo7kI3gN4iq1wN4LlFujAD/MKmRdsLPXGkxTwlozYUvBeFtNmurs6cunOwt7FUKtVTBNet9u50Z1bRrxeLoww4kbI/siLztosOZoEvBD9ozldZ7DSK6Vnz8LQDesv8AgHcT0W4iuhHAQQDfAfAwgINEdCMR7UJm7Hug/2YbMVorG0k1XaKp7IXROCZcsWOCwtp8V3fzp8cq2xSlrnEuIlylz6kOResINGHXNAsf7uvtEdKarwj/VM7piegLAP4GwBuI6AwRbQD4GBF9n4geA/DrAH4XAJj5CQBfAvADAJsA7nYawQ6A9wM4BeBJAF9ydY0h4hNhaJb5ylFdlqUEuiMMf9o1nM89vH9qHh/Dn1PZpkQ7fZtKwi9X4jUoS/oRBuiIzS7oir3Fa0oBDo17HdbtEYGRcBICDwA0zQn619bWeGtra9LNmAn8uvh8hN21nLuUSiM5RuiyU+hH4CV5ko8rrwD//JXBGlOlxYSqeuo6sTraMXltagym9UQgokeYea2qnkXkzQGtlQ3QUmbFztfBb++Ad7K89dRsFNbHS4EPO/2qZba1BglhER9U4IFsfs/M4F/8onv9OqQW1nirfyp2vupavrzqmOgEuN0ZicD3QsTfYcwCt+39V6Bd2dyRhOrohToMp9Vy4HntoO4y27C+ihvRhqm+nrpwspvlp66qn6oXrgOoGt219QXhe6neK9cbRgc4DGykn1Far/nXoKWl6KKYkl8eyOpH5vreFlCKy0fZol/oNK68Qr3/KOarpy6cjBvTYp9TxEb51AIaxQ+fXHQzRI1nWJjQzyDr++4ELm+DXIRYuBw2VN3zhTZO3Zd/XiNo/NK+/H1e7jqUwkKd0CUXzLHD8Nphw+1OOuV1yuPg2bsnM9aFq9xCK3yK2KIexUswKYNdDBP6GWN9/+FM2KGP0gWhDebw+Z+b4xeu8fOflzoJ+RDXSZ7JzGp47TA5deFk7SCfnNCSf3k7D4VV69UReB+SG1PxvUtuikZ4jwn9DOEDOUrC6Qg/y7ryvR8tw6AdOaLHlteG9yuo+rt7jJXvk4IgxYR/17K+Nt4vkpGfU4a4GMGIHv7um+eOT6XAAyb0M8P6/sN5WG1saWwo4CWDnkK4xj5Wv9J41yBsvvhntb/PoGyeO54OmvERd9IvH7O0h5F0sRz2ku2dbFnw+RPdDtP9f6ZV2D0m9DNAptI3SkIXE2zNAh+ORKWY+2ajW0cZIWNW/XwRzQTmrV7wS+m2RKix+h7Qhd9vb5VKpKGo7ZvnT2QGzqWlqRd4wIR+6tGWxobprMKR2r/GfO6aoa8gEGI1WThvD6E9eyb6oJcScKRi5Osa6lJprFyIc/h/2Tx3fOL+97qY0E8xcn+5cDRLhabG5vwhUvhj8/iUuo/du7H5wqfrf6ERUeh0fOBNajmtZnn3733svIYL0aVmA2hMPpV1v5jQTyl56GkQTQd0BVELqpHvtfBbT5W6L+vF1P5pEHhPLvhecJd36XvKhxZ6aW3XrPcyN5/vEDoui3AvHoQpwoR+CimlqRZqeiqaLkaqXjgFCF99x+P9/J5Ru+YGhXbvKi+EkUKqbcChjf7SFiA6BqLhRhyOExP6KWN9350AilZ1KbSp+bk6V/fHrthbmvuHFnvtHn4aIe0H02qsytvVbmfZajgYjcM19aE1P+WnDxJlTOtvUAcT+inCB95ovva6o3p+bhhU47ZoCjuTUPgL58j3zi017SN8YfT1QTSA7q6rE8Pvz3GZa0GNmR3hPSb0U4If4YF4CqvCSO46h1AbyOkUY+nD64Uaghae62HmmRrdvO8cQHSDjaSVX6vLHdCe3TNjoU9hQj8FrO8/DDTKobE5DSoJsJ9fS/+9jJfXrqGVF0Z1LTy3hnYxjXjfeY5cUhuGy3rB14TfbW+9ee44/uon/2U8jR8xJvQTRu4eG6K50EKLvQ+pDRfEFGLtnR8+XJVXsu477YCaje5iG2fEm0WVttBmL9DbYtMMn2gEKBn88t+m3Z4ZDacuJvQTREtvFarbdUdZqbJH3XPttu7vDwlcW9M+j0+RC6yW5uvydnasQaWkGv43nDeBB0zoJ0a464xU36sCajTVO7bQxo+mKWDqAAAUU0lEQVTcBXweO7G+PmYbmJeH3tslfGBNgQ535/6d+fvuISb0E8CP8AVXWTg3D4TRl8skFnKJaWiMK9RZ6iZI8lMComKMOREVUm4B8/PQy9Ge251uDoLQziFSV8/Ld9cwoR8zcvFMYXQVy121SLq8vlhpl/vQfR0fvRfG5F/e1g1zYncYZs6TaE6zL75ffJ697u/Yzpe/5poQz8YquUExoR8jXuDD7C+pubv0u8f86jlBcElhnh9kwQGQq7RapN88ItNtyS2lNs8dz42V8y7wgCXGHAvr++4ENZuFNdeeUIhVv/rSUiGVtXZObBltfrwTtxWEQj7PD37MCzHLxspesZF+xLRW7wItL5Wi4ICiwIbH6Yq9XWH0IaMoqvuFTmLXcrcj0EZ1FEd+oJwME1ish39RsZF+hHgLfWwkjlnrmRl4Jcvxnhv0pIbQbOYjfy7Yzv3kR3VJ6t5yumECvxjYSD8iWgeOVLrhwlE/VNVLlntByQ7gpg8xQ6BsQ8E74Cz5JvCLQ5297G4gom8Q0ZNE9AQRfcCVHyCiB4noafd6tSsnIvokEZ2mbP/6m8W1Drn6TxPRodF9rcniN5EsRcUhENZG0ede8Kk3iudIeGenHMDj3W8NXX0vTC0agTEw2KvemG/qjPQ7AH6fmX8NwC0A7iaimwDcA+AhZj4I4CH3GQDeiWy32oMAjiLbyx5EdADAhwG8DdnW1h/2HcU80VrZKAkbxJzeI416hdHYq/Gd8ggfRtxpI3vugw8o1AvCeseZ0NKYPJVCz8zPM/N33fuLyHadvQ7AHQDuc9XuA/Au9/4OAJ/jjG8BWKFsa+t1AA8y8zlmfgnAgwBaQ/02EyZMn5QLpZtvhyO3NpKHEXcy2UPo5tM+ax6AmIoPmOFuEelpTk9ErwPwFgDfBvBaZn4eyDoGAK9x1a4D8Kw47Ywri5WH9zhKRFtEtHX27NlemjdRNIHXBNMfA6DnZQ+jxJzanqvxbkmtNl0IbQK5Or8cbG7p7mMCv5jUFnoiugrAlwF8kJl/lqqqlHGivFjAfIyZ15h5bXV1tW7zJopcOOOJGeAKo+5Oea23VM9Lqjt3U1bl19+zpzTCy4AeAGAXkRdOA4zFpJbQE9EyMoH/PDN/xRW/4NR2uNcXXfkZADeI068H8FyifKZZ3384f6/NpTVLulTbY+dIX3vhNYjRx6VLpfvk1wl2srFR3gDqWe8JwAkATzLzn4pDDwDwFvhDAL4qyt/nrPi3ALjg1P9TAG4joqudAe82VzazeJU+NjcPV7cVXHjKKC/PDTPf5K8ihFdzwYXX0bSFWVwbbwyPOsE5bwfwXgDfJ6JHXdkfAPgogC8R0QaA/wfgX7hjXwNwO4DTAF4B8NsAwMzniOiPATzs6v0RM58byreYAKU5/O5d4EuXi5U6rHYIgD7nL7nhlPql62k7vHiazULnYgE4BlBD6Jn5r6HPxwHgN5X6DODuyLVOApj5p87npJcjLjnfebYtUjk3XWhBjwm7rJOfl9iRRRP4/HqBNsGvXurxmxrziIXh9ogXeCCInffC19Yj7zzqaK3UL6jwbX1BTqgNhBpCeN+vX/p87e9pzC8WhtsDPp8d/dI+APrGECHhvL5KcP31wnk6LS1VxvD799IQ6D+bWm94TOhrImPpOxcyj2VoFVct9UL9DgVSnhNa50vqfuiqC2L1Q61B3tcE3pCY0Ndgff/hUoYaKXwpNR0oRstp4bcFK71ws6Us8v48/xoG5OTt6pg/3ihiQl+BXB4rc6jlNMor5Xx92v9LpetJQY11FmFd7/7TOokwMs8f8xqFzBBjGIAJfRK5VTSAUkhsNpLq8e8A0Dl/oXA9L4jawhmJGonntIOSGh/Uk9c0f7yhYUIfwUfaxSzv4TGpVsvlrSULu7ayLiHI+bXFlCKFb+M8p7wyBsNcdgGt1bvyVXEhany7cgztboplzY+uBd+U7uUy2uT3Utx2MU+ACbyRwkZ6wfr+w4UND73gSvdXSlhL6nrFrqjJ68WCcZQAHvmZRS57w9AwoXcUDHYoutzkyreYTz6Mky8Y/hRSlvewXKLN/+WfGe6MKkzo0fXBp9xbmrW8itiimpS2UOUKVEd3V2ZqvVGHhRf61spGPpJTkKWmSuhivnR5XKPKcBcrS93LBN6oy0ILvdxTrrA9cYRQDdcWzkhCP3vVdasEvY52YRhVLKzQh0tjNVea9hlQVsmJ0FntXOln9+WlyDmUDYGpaD95LQuzNXphIYVeZrsBiplqNNU7Gg7rz/EdRjiaLxd3gZXX9K8xa3ylS891GKbWG72ycH76ML2VjHmPEVvNRp1gLXywrJbFrjPaNanZKEX0he81a753zZml3uiHhRrppUpfFfteZaTTrhHmtZPXKnwONIQwlj5/H2gfeb1OB7D5vdEnCyP06/sPFxbHaIRGszoW9MI5wUYTUd9+YiOLwj07EXW/0cCpl++DYfTDQqj3cl/4FCkfuWatT10jfB/W0cJp5WvsWjaPNwZl7kf6XKWvCImlpW7/l/Kbh6NzzGJPRIVr+mtELf1iQU1Uy2g2TOCNgZlroZcr5SpH6J0g5r6CcFQO3W/Yuye/ZiG8VqSwLqBseNlrLL9h1GFuhd7ns6sdcCOMcLFz/GicC7wyOufCLfaX10ZuLSmGGr8v2mmjvDEM5lLo5TZTYWqrkFywOvqS2UJdYW0nokIKrdL1Ip8rR/wA32mYwBvDYu6EXrrliCi6dXNIryGuvcbGy/NU1R0AuX3iw2tbxJ0xTOZK6EvprVCee/uyfmmIvHcx/7qGF2h5XukcsRbeq/82whvDZm6EPtyEAtDzwYflnlgwTljm01+HdaoCfcJsPHmCjsDPn9evqf4bRq/U2cDyBiL6BhE9SURPENEHXPlHiOgnRPSo+7tdnPMhIjpNRE8R0boob7my00R0z7C+hNxmKhRAbSlsnUU0vn4dQpW9ity4FxFsfz1T641RUCc4ZwfA7zPzd4loH4BHiOhBd+wTzPwfZWUiugnAuwG8CcDfA/A/iOhX3eFPAbgV2bbVDxPRA8z8g0G+QLivnGsDACVePbF4JhkfX5HzrlBfuVbYNi2uPsTUemNU1NnA8nkAz7v3F4noSQDXJU65A8D9zHwJwI+I6DSAt7pjp5n5GQAgovtd3b6FPkxx5d8XFtJEBM5/DuuE9bwhUFtQo9ZPCHhoA4hdw0Z4Y5T0NKcnotcBeAuAb7ui9xPRY0R0krI954GsQ3hWnHbGlcXKw3scJaItIto6e/ZstC2h0S5Ur8OgGaCcj06mqg7Pk58r59bKdfz9wvdhRyNfBzEwGkZdags9EV0F4MsAPsjMPwPwGQCvB/BmZJrAx31V5XROlBcLmI8x8xozr62urqptaa1s5NbwXFCazaiKHotxj6Wn1jQHSdjByPX0eaeyVFaiNNU+vLeN8saoqbXghoiWkQn855n5KwDAzC+I458F8Jfu4xkAN4jTrwfwnHsfK69F68CRrt+93S4KsQijDUf0lNU+/w5iQU44f69rRS+cL9qXs7wEXN4uaR6m1hvjpI71ngCcAPAkM/+pKL9WVPstAI+79w8AeDcR7SaiGwEcBPAdAA8DOEhENxLRLmTGvgd6aSxv75TWogP6Ypg60Xjy3NLusoEWEHYiWlitrKfVwfaOfn6nYwJvjI06I/3bAbwXwPeJ6FFX9gcA3kNEb0amov8YwF0AwMxPENGXkBnodgDczcxtACCi9wM4BaAJ4CQzP1G3obnRTvqyFSt5QegiHUPVyJ2ywIejs9cQ6mgFJSNig8DbO0BjbsIljBmApjkAZG1tjbe2tvLP6/vuzNNUS+Gr4wbL6wbr6uussw+vUVlPuabmOWDOsuBY2itjGBDRI8y8VlVvpoaYUxfvBXYt59bycO4ec4tpLjhPlcCH1v86xAyEsk3+vQm8MW5mSugBYPPFPyukm+J2W58/Q88uG1rVCy68AFpail47FuRTSbOZd0IWgGNMgpkTeiCLVssFTJkPh4JciLvf2VHDZtVR3C2AibnX6pSVzvWLamyjSWNCzKTQA5ngp4RMC9QB9OQVsWvEouxiVHUO8rgltjQmxcwKPVD2a2sRcP59LshukUs0Dl+4+YC4kMulsmEnoUUD+uOm1huTZqaFHnBz9GC3WU2gowtqwnqKEU4dwYMgG+1emgfB/PHGpJl5od/86bHSbrPa+4KKHyS0SKrsV+wt+f+rrPmlICEi58azebwxeWZe6IHyMlTpv/cCV2C7uG98yhXX+fkrJR973cUxubC7QCFzzxnTwFwIPeDm951iVJyfQ4fJKmoF2FRY9qsi78JFPabWG9PC3Ai9R1sRF/ucOpZcnFNjlM/Pt7h6Y8qYK6E/dfHe7gYTLggnFZ6bHwvqVlG1iEceM5XemDbmSugB5//udErBL0l1fGenVBb18yvaQOiay4+Z4c6YQuZO6D1+cU2vpDSC1L0KdLL4fgvAMaaRuRT6UxfvzYRX2YHGE5vP11HdY9fIy5tNC8Axppa5FHogc+PFQnWB9Fxf+6yVa8E3/t6GMa3MrdB7CpZzEbmXWncfooXVRhfsdOqtzTeMSTH3Qg8IwRcr8upqAGFd1fovOhCz1hvTzkIIvbbllWZxl8dTATlhph7/3vzxxiywEEKPTqeUVgvoPdIOQMkjYAJvzBq1UmDPMutXHSpkwKmV3cYhY/g9anotm8cbM8Rcj/StA0fyaDsiQuPKK/JjWgcQBuGE2oGWcQewqDtjtphboV+/4r3Z9tCNbrLMzs9fyY+H83IiyjajEMgMup4w8aap9casMbdCf+qVPwftWs6z32phsiVVX1lymyXfLK7SIyLQ0pL5442ZZG6FHgA2zx0XCTTjefGYObmZZb4uX3YSkU0rDWPamWuhB1x0XEMPrin42TtFQ18hg67rFLojfxubZ//rBL6NYQzO3As9kI34m+eOAyKtlkdqAlomXani+89muDNmmTobWO4hou8Q0feI6Aki+neu/EYi+jYRPU1EX6RsU0pQtnHlF4notDv+OnGtD7nyp4hofVRfKsbmT49leeo6HbDbVTZX+4N5u2hzwWpvhjtj1qkz0l8C8BvM/A+R7UXfIqJbAPwJgE8w80EALwHYcPU3ALzEzP8AwCdcPRDRTch2qn0TgBaATxNReegdMacu3gtqNrvJNP1uOSIfvvwsVX4z3BnzQKXQc8bL7uOy+2MAvwHgv7ny+wC8y72/w32GO/6blEnNHQDuZ+ZLzPwjAKcBvHUo36JH/Dy/YMDrcMHY5zuD/LNhzAm15vRE1KRsm+oXATwI4G8BnGdm7+M6A+A69/46AM8CgDt+AcA1slw5R97rKBFtEdHW2bNne/9GNdk8d9wl0xTW/EDQpbDbKG/MC7WEnpnbzPxmANcjG51/TavmXjVfFifKw3sdY+Y1Zl5bXV2t07yB2Dx/opSTPnTV2TzemCd6st4z83kA/wvALQBWiMiHsF0P4Dn3/gyAGwDAHd8P4JwsV86ZKJvnjgPMpc0obAsqYx6pY71fJaIV934vgH8K4EkA3wDwz121QwC+6t4/4D7DHf+fnOnJDwB4t7Pu3wjgIIDvDOuLDMqpl+/Lsukmkl0axjxQZ6S/FsA3iOgxAA8DeJCZ/xLAvwXwe0R0Gtmc3Q+JJwBc48p/D8A9AMDMTwD4EoAfANgEcDczT1262FMXTuaCb6O8MY/QNFum19bWeGtra9LNMIyZgIgeYea1qnoLEZFnGEYXE3rDWDCmWr0noosAnpp0Oxy/DODvJt0IwTS1x9qiM+62/H1mrvRzT3u6rKfqzFHGARFtTUtbgOlqj7VFZ5raIjH13jAWDBN6w1gwpl3oj026AYJpagswXe2xtuhMU1typtqQZxjG8Jn2kd4wjCFjQm8YC8bUCj0RtVxardNEdM+Y7vljIvo+ET1KRFuu7AARPejSgj1IRFe7ciKiT7r2PUZENw9475NE9CIRPS7Ker43ER1y9Z8mokPavfpsy0eI6Cfut3mUiG4Xx9Q0aMP4HxLRDUT0DSJ60qVr+8CkfptEWyby2/SNTxYxTX8AmsgSdfwKgF0AvgfgpjHc98cAfjko+xiAe9z7ewD8iXt/O4C/QpYn4BYA3x7w3u8AcDOAx/u9N4ADAJ5xr1e791cPqS0fAfBvlLo3uf/PbgA3uv9bc1j/Q2QLvm527/cB+KG759h/m0RbJvLb9Ps3rSP9WwGcZuZnmPkygPuRpduaBDL9V5gW7HOc8S1k+QWu7fcmzPxNZHkHBrn3OrJVkOeY+SVkWY5aQ2pLjFgatKH8D5n5eWb+rnt/Edmy7uswgd8m0ZYYI/1t+mVahb5Waq0RwAC+TkSPENFRV/ZaZn4eyP7pAF4zxjb2eu9Rt+n9TmU+6dXpcbaFsszKbwHwbUz4twnaAkz4t+mFaRX6Wqm1RsDbmflmAO8EcDcRvSNRd1JtTN17lG36DIDXI8uI/DyAj4+zLUR0FYAvA/ggM/8sVXXU7VHaMtHfplemVegnklqLmZ9zry8C+AtkatgLXm13ry+OsY293ntkbWLmFzjLldgB8Fl0MxmPvC1EtIxMyD7PzF9xxRP5bbS2TPK36YtxGQ96NJgsITO03IiuoeNNI77nlQD2iff/B9mc7z+gaDD6mHv/z1A0GH1nCG14HYrGs57ujcxI9SNkhqqr3fsDQ2rLteL97yKbqwLZPgbSWPUMMkPVUP6H7jt+DsB/CsrH/tsk2jKR36bv52xcN+rjn307Muvo3wL4wzHc71fcj/89AE/4eyJLBfYQgKfd6wHxAHzKte/7ANYGvP8XkKmG28hGgo1+7g3gMDKD0WkAvz3Etvy5u9djyPIdygf9D11bngLwzmH+DwH8E2Sq72MAHnV/t0/it0m0ZSK/Tb9/FoZrGAvGtM7pDcMYESb0hrFgmNAbxoJhQm8YC4YJvWEsGCb0hrFgmNAbxoLx/wEJTHKNQ1jOagAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(np.log10(MIPS_Map[1].data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in PSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_psf=fits.open('../../dmu17/dmu17_HELP-SEIP-maps/EGS/data/output_data/dmu17_MIPS_EGS_20190204.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADXlJREFUeJzt3W+snvVdx/H3R/bPTD1zo9OFUtulhKyayZITtmQ+mJPMAjswzYxUH6A2NEvEzMREIRgTHxgxJmqILOboCHtAQETZ6OzCELfwBDdgY5Omq+uQhVqyOueO/xKQ7euD+2o4HO/T3of7vs913b++X0nTc//uf9+Ui0+vfq/f9fulqpAktet7+i5AkjRfBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpca/quwCACy+8sHbv3t13GZK0UJ544olvVtWOc71uEEG/e/duHn/88b7LkKSFkuTrk7zO1o0kNc6gl6TG9Rr0SVaSrK6trfVZhiQ1rdegr6rDVXVoaWmpzzIkqWm2biSpcQa9JDXOHr0kNc4evSQ1bhA3TKk/u2/625c9fubWq3uqRNK8GPR6GYNfao9B37CNoS3p/GTQ66zG/WXhWb60WJx1I0mNc9aNJDXO1o22zAu20mLxzlhJapxn9I1who2kzXhGL0mN84xeU3MKpjRsTq+UpMY5vVKSGmePXpIaZ49+QTnLRtKkPKOXpMYZ9JLUOINekhpn0EtS4wx6SWpcr7NukqwAK3v37u2zDM2BK1xKw+ENU5LUOFs3ktQ4b5haEN4gJemVMui1LVzhUuqPrRtJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOKdXqjcukyBtD4N+gLw5StIs9dq6SbKSZHVtba3PMiSpaS5qJkmNs3WjwbBnL82Hs24kqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjfOGKQ2W+8xKs+EZvSQ1zqCXpMYZ9JLUOHv0A+D685Nz4TNp6zyjl6TGGfSS1Li5BH2SDyT58ySfSPK+eXyHJGkyEwd9kjuSnE7y1Ibx/UmOJzmR5CaAqvp4Vd0A/BLw8zOtWJK0JVs5o78T2L9+IMkFwO3AlcA+4ECSfete8tvd85Kknkw866aqHkmye8Pw5cCJqnoaIMk9wLVJjgG3Ap+qqi+M+7wkh4BDALt27dp65RLePStNYtoe/UXAs+sen+zGfg24Avhgkg+Ne2NVrVbVclUt79ixY8oyJEmbmXYefcaMVVXdBtw25WdLkmZg2jP6k8DF6x7vBE5N+ZmSpBmaNugfAy5JsifJa4DrgAcmfXOSlSSra2trU5YhSdrMVqZX3g08Clya5GSSg1X1InAj8CBwDLi3qo5O+plVdbiqDi0tLW21bknShLYy6+bAJuNHgCMzq0iSNFO9LoFg60aS5q/X1Sur6jBweHl5+YY+61BbXOFSejkXNZOkxhn0ktQ4e/SS1Dh79NvM3aQkbTdbN5LUOINekhpn0EtS47wYK0mN6zXoXetGkubP1o0kNc6gl6TGGfSS1DiDXpIa56wbSWqcSyDovDNuGQqXMlbLeg16aTu4vpDOd/boJalxBr0kNc6gl6TGGfSS1DinV0pS41zUTJIa5/RKif8/BdN59WqJPXpJapxBL0mNM+glqXEGvSQ1zoux0hgufKaWGPRz5oJakvrmDVOS1DhvmJKkxnkxVpIaZ9BLUuO8GCtNyGUStKg8o5ekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1rtcbppKsACt79+7tswxpJlzaWEPVa9BX1WHg8PLy8g191jFLLkt8/vC/tRaFrRtJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaN/OgT/LWJB9Nct+sP1uStHUTBX2SO5KcTvLUhvH9SY4nOZHkJoCqerqqDs6jWEnS1k16Rn8nsH/9QJILgNuBK4F9wIEk+2ZanSRpahMFfVU9Anxrw/DlwInuDP4F4B7g2hnXJ0ma0jR7xl4EPLvu8UngnUneBPwe8I4kN1fV7497c5JDwCGAXbt2TVGGtNgm2XvWTcY1jWmCPmPGqqr+DfjQud5cVavAKsDy8nJNUYck6SymmXVzErh43eOdwKnpypEkzdo0Z/SPAZck2QP8C3Ad8Atb+YAkK8DK3r17pyhDGq6NbRlbMOrDpNMr7wYeBS5NcjLJwap6EbgReBA4BtxbVUe38uVVdbiqDi0tLW21bknShCY6o6+qA5uMHwGOzLQiSdJMuQSCJDWu16BPspJkdW1trc8yJKlpvQa9PXpJmj9bN5LUOINekho3zTz6qTmPXuebSZY7eCXv2875+d4bsHjs0UtS42zdSFLjDHpJapxBL0mN82Ks1IBxF2u36yJpn9+tyXgxVpIaZ+tGkhpn0EtS4wx6SWqcQS9JjXPWzRRe6e3s0qLwGG+Ds24kqXG2biSpcQa9JDXOoJekxhn0ktQ4g16SGuf0yi1wqpkWmYuPnb+cXilJjbN1I0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS47xhStLMbbw5a+ONWZPcvOUNXrPjDVOS1DhbN5LUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIa51o30nls0Ta8P9caOhrPtW4kqXG2biSpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjZv5DlNJXg98BHgB+GxV3TXr75AkTW6iM/okdyQ5neSpDeP7kxxPciLJTd3wzwL3VdUNwDUzrleStEWTtm7uBPavH0hyAXA7cCWwDziQZB+wE3i2e9l3ZlOmJOmVmijoq+oR4Fsbhi8HTlTV01X1AnAPcC1wklHYT/z5kqT5maZHfxEvnbnDKODfCdwG/GmSq4HDm705ySHgEMCuXbumKGM+Nu42Ly2aIR3Dk9SyXfWO+55nbr16W767L9MEfcaMVVX9N/DL53pzVa0CqwDLy8s1RR2SpLOYprVyErh43eOdwKnpypEkzdo0Qf8YcEmSPUleA1wHPLCVD0iykmR1bW1tijIkSWcz6fTKu4FHgUuTnExysKpeBG4EHgSOAfdW1dGtfHlVHa6qQ0tLS1utW5I0oYl69FV1YJPxI8CRmVYkSZoppz9KUuN6DXp79JI0f70GvT16SZo/WzeS1LhU9X+vUpJ/Bb7+Ct9+IfDNGZYzb9Y7f4tWs/XOV8v1/khV7TjXiwYR9NNI8nhVLfddx6Ssd/4WrWbrnS/rtXUjSc0z6CWpcS0E/WrfBWyR9c7fotVsvfN13te78D16SdLZtXBGL0k6i4UN+iR/mOQrSb6c5P4kb1j33M3dPrbHk/x0n3WekeTnkhxN8t0kyxueG1y9sOmewIMxbi/jJG9M8lCSr3a//2CfNa6X5OIkn0lyrDsWPtyND7LmJK9L8vkkX+rq/d1ufE+Sz3X1/mW3eu1gJLkgyReTfLJ7PPR6n0nyj0meTPJ4NzbTY2Jhgx54CPixqno78E/AzQDdvrXXAT/KaJ/bj3T72/btKUYbpz+yfnCo9Z5lT+AhuZMNexkDNwEPV9UlwMPd46F4EfiNqnob8C7gV7s/06HW/Dzw3qr6ceAyYH+SdwF/APxxV++/Awd7rHGcDzNaUfeModcL8JNVddm6aZUzPSYWNuir6tPdUskA/8BL+9ReC9xTVc9X1T8DJxjtb9urqjpWVcfHPDXIetl8T+DB2GQv42uBj3U/fwz4wLYWdRZV9VxVfaH7+T8ZhdFFDLTmGvmv7uGru18FvBe4rxsfTL0ASXYCVwN/0T0OA673LGZ6TCxs0G/wK8Cnup/H7WV70bZXNLmh1jvUus7lh6rqORgFK/DmnusZK8lu4B3A5xhwzV0b5EngNKN/RX8N+Pa6k6yhHRd/Avwm8N3u8ZsYdr0w+svz00me6PbShhkfE9PsGTt3Sf4O+OExT91SVZ/oXnMLo38S33XmbWNevy1Tiyapd9zbxowNYSrUUOtaeEm+D/hr4Ner6j9GJ53DVFXfAS7rroHdD7xt3Mu2t6rxkrwfOF1VTyR5z5nhMS8dRL3rvLuqTiV5M/BQkq/M+gsGHfRVdcXZnk9yPfB+4KfqpXmive1le656NzHUvXeHWte5fCPJW6rquSRvYXQmOhhJXs0o5O+qqr/phgddM0BVfTvJZxldW3hDkld1Z8lDOi7eDVyT5CrgdcAPMDrDH2q9AFTVqe7300nuZ9Q2nekxsbCtmyT7gd8Crqmq/1n31APAdUlem2QPcAnw+T5qnNBQ6516T+CePABc3/18PbDZv6S2Xdcv/ihwrKr+aN1Tg6w5yY4zs9mSfC9wBaPrCp8BPti9bDD1VtXNVbWzqnYzOl7/vqp+kYHWC5Dk9Um+/8zPwPsYTdyY7TFRVQv5i9FFy2eBJ7tff7buuVsY9RKPA1f2XWtX088wOkt+HvgG8OCQ6+3quorRjKavMWo/9V7ThvruBp4D/rf7sz3IqCf7MPDV7vc39l3nunp/glHb4Mvrjturhloz8Hbgi129TwG/042/ldHJyAngr4DX9l3rmNrfA3xy6PV2tX2p+3X0zP9nsz4mvDNWkhq3sK0bSdJkDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhr3f5BU0nJ9qHz0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dat=MIPS_psf[1].data\n",
    "plt.hist(dat.flatten(),bins=np.arange(-20.0,50.0,1.0));\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "centre=np.long((MIPS_psf[1].header['NAXIS1']-1)/2)\n",
    "radius=5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7ff09141be80>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAD8CAYAAAD9uIjPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFoBJREFUeJzt3X2MXdV97vHv4xmbF5vXONhgQyHFyY0vzUs1JVDUl9SIEppC21uuIEqumwZZkUqbppESqKWmansletMmrdSojZXQooYmRTQIFBzAJEVRdAXBDjRgTBLHJTCYN4cQSAi2Z+bpH+cMnZkzM+fMnDXn7LP9fKQtn33Ont/+MS8/1t5r7bVkm4iIOlnW7wQiIkpLYYuI2klhi4jaSWGLiNpJYYuI2klhi4jaSWGLiCUl6WJJ35K0V9I1s3x+lKR/bX5+n6Qzuz1nCltELBlJQ8AngXcAG4ErJW2ccdj7gB/YPhv4BPCX3Z43hS0iltK5wF7b+2wfAj4PXDbjmMuAG5qvbwY2SVI3Jx3u5osXasXwsT5m+QnF4o2tLJv+xPJysVz1/2V09Wuz9DReNt7wy+WesHnDhrXFYgF8e8/+YrFeOfgCh8Ze7uqn+6tvX+nvP9/ZD2DXNw/uBl6Z8tY229um7K8DnpiyPwq8bUaYV4+xPSbph8BrgAMLTP1VPS1sxyw/gfPPfl+xeAdGTi4WC+DHp5X7ax8/puyjaqULpYfKxqPwk3nLXypbedfsOlQs1j13fKRYLICLfu5Pi8W695Ft7Q9q4/vPj/P1O8/o6NihU7/ziu2ReQ6Z7Qc587elk2MWpKeFLSKqz8AEE6XCjQKnT9lfD8xsok4eMyppGDgBeL6bk1b9gikiesyYwx7vaOvA/cAGSWdJWgFcAdw245jbgM3N178NfMVdzs6RFltEtCjVYmveM7sauBMYAq63vVvSnwE7bd8GfAb4Z0l7abTUruj2vClsETGNMeMFpzOzvR3YPuO9P5ny+hXg8mInJIUtImYxUbo3qMe6usfWbkRxRAweA+O4o62qFl3YOhxRHBEDaAJ3tFVVN5eir44oBpA0OaL4kRKJRUR/GDg84EsGdFPYOhlRjKQtwBaAo5cf38XpIqIXXPHLzE50U9g6Gi3cfLxiG8AJx5w62N+tiCOBYXzA/1K7KWydjCiOiAHTePJgsHVT2F4dUQw8SWNQ3buKZBURfSTGqz5LQhuLLmxzjSgulllE9EWj8+AILWww+4jiiBhsjXFsR3Bhi4h6mjiSW2wRUT9psUVE7RgxPuAzmqWwRUSLXIouwIb/uY47dv5FsXi/dMn/KxYL4IU3lvt2HF5TbipqgKEVgz6yaGFePlh27vInTiy3oMXPbf54sVgAJ+4qOJjAr7Q/pl0IxKHic8f3VlpsETFNY4BuLkUjombSeRARtWKL8cqvHzm/FLaIaDGRFltE1Emj82CwS8NgZx8RxaXzICJqaTzj2CKiTvLkQUTU0kR6RSOiThoPwS99YZN0MvCvwJnAY8D/tv2DWY4bBx5q7j5u+9J2sQe7LEdEcUYc9lBHW5euAb5sewPw5eb+bH5i+y3NrW1RgxS2iJjBhnEv62jr0mXADc3XNwC/0W3ASSlsETGDmOhwA1ZL2jll27KAE62x/RRA899T5jju6GbseyV1VPxyjy0ipjEspDV2wPbIXB9KuhtYO8tHWxeQ0hm290t6HfAVSQ/Z/u58X5DCFhEtSnUe2L5wrs8kPSPpVNtPSToVeHaOGPub/+6TdA/wVmDewpZL0YiYxogJd7Z16TZgc/P1ZuDWmQdIOknSUc3Xq4ELgEfaBU6LLSKmaSy/15PScB1wk6T3AY8DlwNIGgHeb/sq4I3ApyRN0GiIXWc7hS0iFqo3Cybb/j6waZb3dwJXNV//f+BnFho7hS0ipjF58mBBvvPwKBe//sPF4j39rlOLxQI4fEq5dQqWH3O4WCwADfYzyQs2NDxeNN7h4XJrRjzvo4rFAlj2rvOKxRr/0r1l4mQ+toioE1tpsUVEvTQ6D7JKVUTUStY8iIiaaXQe5B5bRNRMJpqMiFqZfPJgkC26LEs6XdK/S9ojabekD5RMLCL6Z4JlHW1V1U2LbQz4kO1vSDoO2CVpRyePO0REddlweKK6RasTiy5szfmTJudSeknSHmAdHTygGhHV1bgUPUIL21SSzqQxlch9s3y2BdgCcPTw8SVOFxFL7Ih/8kDSKuDfgD+0/eLMz21vA7YBnHD0Wnd7vohYWkf8cA9Jy2kUtRttf6FMShHRX0fwpagkAZ8B9tj+eLmUIqLfJo7gS9ELgPcAD0l6sPneH9ve3n1aEdEvjV7RI/RZUdtfgwEv6xHRog4DdPPkQUS0OJIvRSOiho74XtGIqKcjtld0MQ6+Zjn73lNuOu+DZxwsFgtg+Khy01EvW1Z2yN7GtU8XjXfGyh8Ujbf/JycUjbfnuTVF47ngj+Pwa8reWH/h7OXFYo0XmLXcFmMpbBFRN7kUjYhayT22iKilQS9sg30hHRHFTY5j62TrhqTLm3M5TjRXf5/ruIslfUvSXknXdBI7hS0iWkygjrYuPQz8FvDVuQ6QNAR8EngHsBG4UtLGdoFzKRoR09gw1oOJJm3vAdD8q4GfC+y1va957OeBy2gz72MKW0S0WMBl5mpJO6fsb2tOVVbKOuCJKfujwNvafVEKW0RMs8BnRQ/Ynu/+2N3A2lk+2mr71g7iz5ZI21GJKWwR0cKFekVtX9hliFHg9Cn764H97b4onQcR0aJHnQeduB/YIOksSSuAK4Db2n1RCltETGPTq+EevylpFDgfuF3Snc33T5O0vZGLx4CrgTuBPcBNtne3i51L0YiYQYz3plf0FuCWWd7fD1wyZX87sKAJbFPYIqJFqXts/ZLCFhHT5FnRiKgfl53mqR9S2CKiRaYGj4hacY86D5ZSCltEtMilaETUTnpFF8DD5uApY8XirTzhlWKxACYmyv0wz1n7VLFYAFetnXNml0X56eVl1zz47uGTisa7Xr9QNN6DT64rFmvo2HK/wwCHTiq3hoIL/EXbKWwRUUMZ7hERtZN7bBFRK0ZMpFc0IupmwBtsKWwRMUM6DyKilga8ydb1hbSkIUkPSPpiiYQiov9sdbRVVYkW2wdoTAB3fIFYEdFnpuyYzn7oqsUmaT3wa8Cny6QTEX1nwOpsq6huW2x/A3wYOG6uAyRtAbYADJ18Ypeni4heGPRxbItusUl6J/Cs7V3zHWd7m+0R2yNDq1Yu9nQR0UvucKuoblpsFwCXSroEOBo4XtJnbb+7TGoR0R/V7hjoxKJbbLavtb3e9pk0lsT6SopaRE0cwS22iKgjgwe8V7RIYbN9D3BPiVgRUQUpbBFRNxW+zOzEYD/CHxFLowf32CRdLmm3pAlJI/Mc95ikhyQ9KGlnJ7HTYouI6SYH6C69h4HfAj7VwbFvt32g08ApbBHRohcDdG3vAZDKF9GBLmxS2e/+8PBEsVivX/VssVgAG1eUXaNg/fCqovGgbH5rj36xaDyp4JoHw+PFYgEcXlHw97hUjei8V3T1jMvDbba3FcpikoG71PiD/1Qn8Qe6sEXE0lhAm+GA7fnuj90NrJ3lo622b+3wHBfY3i/pFGCHpEdtz7u6UQpbRExXcPCt7QsLxNjf/PdZSbcA5wLzFrb0ikbEDB3O7NGDDgZJKyUdN/kauIhGp8O8UtgiolVvhnv8pqRR4Hzgdkl3Nt8/TdL25mFrgK9J+g/g68Dttu9oFzuXohHRqlw/2pxs3wLcMsv7+4FLmq/3AW9eaOwUtoiYrnfj2JZMCltEtCg8kqrnUtgiotWAF7Z0HkRE7aTFFhEtcikaEfViFvJIVSWlsEVEq7TYIqJucikaEfWTwhYRtZPCFhF1IudSNCLqKL2iEVE3abFFRP2ksC2EUMEm7tjYULFYAMuWlZur5ZmDxxeLBfC9sWOLxpvgR0Xjfefwa4rGe+xHZeOVXJyk9CrpOlgwXolf4dxji4haSmGLiLpRDyaaXEqZ3SMiaicttoholUvRiKiVdB5ERC0NeGHr6h6bpBMl3SzpUUl7JJ1fKrGI6KMeLL+3lLptsf0tcIft35a0Aig72Coiek4Mfq/oogubpOOBXwR+B8D2IeBQmbQiom9qcI+tm0vR1wHPAf8o6QFJn24uQT+NpC2SdkraOf6jsqPdI2KJ9GYl+I81b2N9U9Itkk6c47iLJX1L0l5J13QSu5vCNgz8LPD3tt8K/BhoOantbbZHbI8MrVrVxekiomd6c49tB3CO7TcB3waunXmApCHgk8A7gI3AlZI2tgvcTWEbBUZt39fcv5lGoYuIATc5J1u7rRu277I91ty9F1g/y2HnAntt72ve7vo8cFm72IsubLafBp6Q9IbmW5uARxYbLyIqpPMW2+rJW03Nbcsiz/i7wJdmeX8d8MSU/dHme/Pqtlf094Ebmz2i+4D3dhkvIvrNC+oVPWB7ZK4PJd0NrJ3lo622b20esxUYA26cLcTsGc6vq8Jm+0Fgzv+oiBhQhXpFbV843+eSNgPvBDbZs04uNQqcPmV/PbC/3Xnz5EFEtOjFcA9JFwMfAX7J9stzHHY/sEHSWcCTwBXAu9rFzuweEdGqN72ifwccB+yQ9KCkfwCQdJqk7QDNzoWrgTuBPcBNtne3C5wWW0RM16PHpWyfPcf7+4FLpuxvB7YvJHYKW0RMIwb/yYPeFrZlxseOtT+u03AF1ygobdczsw3J6cYvFI225qgXi8Z79KU1ZeM9c0rReBMF1yk4/NKKYrEAVj1d7o7QskJ/XilsEVE/KWwRUTspbBFRKzWY3SOFLSJapbBFRN0csRNNRkR95VI0Iuql4usZdCKFLSJapbBFRJ3kyYOIqCVNDHZlS2GLiOlyjy0i6iiXohFRPylsEVE3abFFRP2ksEVErSxslapKSmGLiGkyji0i6mnWlfAGRwpbRLRIi20hBBoq9x07dHB5sVgAKriGwvh42ZUN793/U0XjlTY2NtTvFOY1drDcr/qqvWV/73Zf98FisfSXf7Sr6yAZoBsRdZTOg4ionRS2iKgX05POA0kfA34dOAR8F3iv7RdmOe4x4CVgHBizPdIudtkbQRFRC3JnW5d2AOfYfhPwbeDaeY59u+23dFLUIIUtImbjDrduTmHfZXtyied7gWKrjHdV2CR9UNJuSQ9L+pyko0slFhH9MTlAt8MW22pJO6dsWxZ52t8FvjTHZwbukrSr0/iLvscmaR3wB8BG2z+RdBNwBfBPi40ZERVgL2SiyQPzXR5KuhtYO8tHW23f2jxmKzAG3DhHmAts75d0CrBD0qO2vzpfUt12HgwDx0g6DBwL7O8yXkRUQaG+A9sXzve5pM3AO4FN9uw9Frb3N/99VtItwLnAvIVt0Zeitp8E/gp4HHgK+KHtu2ZJfMtkM3X8xR8v9nQR0UO96DyQdDHwEeBS2y/PccxKScdNvgYuAh5uF3vRhU3SScBlwFnAacBKSe+eeZztbbZHbI8MHb9ysaeLiF4xMOHOtu78HXAcjcvLByX9A4Ck0yRtbx6zBviapP8Avg7cbvuOdoG7uRS9EPhP2881k/kC8PPAZ7uIGRFV0INHqmyfPcf7+4FLmq/3AW9eaOxuCtvjwHmSjgV+AmwCdnYRLyIq4oh9CN72fZJuBr5Bo0fjAWBbqcQion+O6OX3bH8U+GihXCKiCjK7R0TUTWOA7mBXthS2iGiV2T0iom7SYouIesk9tgU6LJY9c1S5eIWby+NHlftpjp10qFgsgIkJFY2nsuGYKDwV+tjLZX81j3lsRbFYJafyrqYFPStaSWmxRUSrXIpGRK1kweSIqKW02CKidga7rqWwRUQrTQz2tWgKW0RMZzJANyLqRTgDdCOihlLYIqJ2UtgiolZyjy0i6ii9ohFRM86laETUjElhi4gaGuwr0RS2iGg16OPYyk6iFRH1YHe2dUHSn0v6ZnOx5LsknTbHcZslfae5be4kdgpbRExnw/hEZ1t3Pmb7TbbfAnwR+JOZB0g6mcZKeG8DzgU+KumkdoFT2CKiVQ9abLZfnLK7ktnnFPlVYIft523/ANgBXNwudu6xRUSrHt1jk/R/gf8D/BB4+yyHrAOemLI/2nxvXj0tbD+zZg07P/RHvTzlgpx35V8Xi/X9cwqu7QCMrSw3Zz+UnyF1+UtlG/+vv/7xovG+9L1PFI1XawY6X/NgtaSdU/a32d42uSPpbmDtLF+31fattrcCWyVdC1xN6wLss63O0Ta5tNgiYgaDO/4/3wHbI3NGsi/sMM6/ALfTWthGgV+esr8euKddsNxji4jpTE86DyRtmLJ7KfDoLIfdCVwk6aRmp8FFzffmlRZbRLTqzT226yS9gcZw4O8B7weQNAK83/ZVtp+X9OfA/c2v+TPbz7cLnMIWEa16UNhs/6853t8JXDVl/3rg+oXETmGLiBnyEHxE1I2BAZ+2qG3ngaTrJT0r6eEp750saUfzEYcdnYwEjogB0oMBukupk17Rf6J1pO81wJdtbwC+3NyPiFro2SNVS6ZtYbP9VWBmL8RlwA3N1zcAv1E4r4joF4M90dFWVYu9x7bG9lMAtp+SdMpcB0raAmwBOOOMMxZ5uojoqc6fPKikJR+ga3ub7RHbI6997WuX+nQRUcKA32NbbIvtGUmnNltrpwLPlkwqIvrIrn+v6BxuAyYnfNsM3FomnYiohLq32CR9jsZDqKsljdJ4SPU64CZJ7wMeBy5fyiQjopeMx8f7nURX2hY221fO8dGmwrlERBUsbNqiSsqTBxHRqsJDOTqRwhYR0xhwWmwRUSte0ESTlZTCFhEtBr3zQO5hl62k52hMKNfOauDAEqezWFXODaqdX5Vzg2rn12luP2W7q5Hwku5onq8TB2y3XTWq13pa2Dolaed886j3U5Vzg2rnV+XcoNr5VTm3KsqaBxFROylsEVE7VS1s29of0jdVzg2qnV+Vc4Nq51fl3CqnkvfYIiK6UdUWW0TEoqWwRUTtVKqwSbpY0rck7ZVUqXUUJJ0u6d8l7ZG0W9IH+p3TTJKGJD0g6Yv9zmUmSSdKulnSo83v4fn9zmmSpA82f6YPS/qcpKP7nE8WUOpSZQqbpCHgk8A7gI3AlZI29jeracaAD9l+I3Ae8HsVyw/gA8Ceficxh78F7rD9P4A3U5E8Ja0D/gAYsX0OMARc0d+ssoBStypT2IBzgb2299k+BHyexqIxlWD7KdvfaL5+icYf5rr+ZvXfJK0Hfg34dL9zmUnS8cAvAp8BsH3I9gv9zWqaYeAYScPAscD+fiaTBZS6V6XCtg54Ysr+KBUqHFNJOhN4K3BffzOZ5m+ADwNVfHr5dcBzwD82L5U/LWllv5MCsP0k8Fc0Jkx9Cvih7bv6m9Wspi2gBMy5gFJUq7BplvcqNxZF0irg34A/tP1iv/MBkPRO4Fnbu/qdyxyGgZ8F/t72W4EfU5FLqea9qsuAs4DTgJWS3t3frKJbVSpso8DpU/bX0+dLgpkkLadR1G60/YV+5zPFBcClkh6jcQn/K5I+29+UphkFRm1PtnBvplHoquBC4D9tP2f7MPAF4Of7nNNsnmkunEQWUGqvSoXtfmCDpLMkraBxA/e2Puf0KkmicY9oj+2P9zufqWxfa3u97TNpfN++YrsyrQ7bTwNPSHpD861NwCN9TGmqx4HzJB3b/BlvoiIdGzNkAaUFqMx8bLbHJF0N3EmjZ+p627v7nNZUFwDvAR6S9GDzvT+2vb2POQ2S3wdubP5Pax/w3j7nA4Dt+yTdDHyDRs/3A/T58aUsoNS9PFIVEbVTpUvRiIgiUtgionZS2CKidlLYIqJ2UtgionZS2CKidlLYIqJ2/gsbkGIpfwom/gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.imshow(np.log10(MIPS_psf[1].data[centre-radius:centre+radius+1,centre-radius:centre+radius+1]/np.max(MIPS_psf[1].data[centre-radius:centre+radius+1,centre-radius:centre+radius+1])))\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set XID+ prior class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: AstropyDeprecationWarning: \n",
      "Private attributes \"_naxis1\" and \"_naxis2\" have been deprecated since v3.1.\n",
      "Instead use the \"pixel_shape\" property which returns a list of NAXISj keyword values.\n",
      " [astropy.wcs.wcs]\n",
      "WARNING: AstropyDeprecationWarning: \n",
      "Private attributes \"_naxis1\" and \"_naxis2\" have been deprecated since v3.1.\n",
      "Instead use the \"pixel_shape\" property which returns a list of NAXISj keyword values.\n",
      " [astropy.wcs.wcs]\n"
     ]
    }
   ],
   "source": [
    "prior_MIPS=xidplus.prior(MIPS_Map[1].data,MIPS_Map[2].data,MIPS_Map[0].header,MIPS_Map[1].header,moc=Final)\n",
    "prior_MIPS.prior_cat(masterlist[1].data['ra'][good],masterlist[1].data['dec'][good],masterfile,flux_lower=MIPS_lower,\n",
    "                     flux_upper=MIPS_upper,ID=masterlist[1].data['help_id'][good])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "prior_MIPS.set_prf(MIPS_psf[1].data[centre-radius:centre+radius+1,centre-radius:centre+radius+1]/1.0E6,np.arange(0,11/2.0,0.5),np.arange(0,11/2.0,0.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([215.08957607, 214.97755187, 214.86830064, ..., 215.03025421,\n",
       "       214.31734531, 213.82151274])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prior_MIPS.sra\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Cannot determine equinox. Assuming J2000. [aplpy.wcs_util]\n",
      "WARNING: Cannot determine equinox. Assuming J2000. [aplpy.wcs_util]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([<aplpy.core.FITSFigure at 0x7ff0913cb240>],\n",
       " <Figure size 720x720 with 2 Axes>)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAJdCAYAAADJBYiHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XlcVOXiBvBnAFkaxB1wl7yYpmRZmvuCgZoioKlpWu6ZGXldSjMJya30mlpqZZpWam6JFRgiJpoLmqbyQ8sNd2UUBXVYBmbO7w9iZJlBkJl5Z+Y838/nfpIz2zPHg+e573nPOQpJkiQQERERkZ6D6ABERERE1oYFiYiIiKgYFiQiIiKiYliQiIiIiIphQSIiIiIqhgWJiIiIqBgWJAP27t2LHj16ICAgAF9//bXoOHbjxo0bGDZsGHr16oXevXtj7dq1oiPZJa1Wi5CQELz55puio9iVe/fuISwsDD179kSvXr3w119/iY5kV9asWYPevXujT58+mDRpEnJyckRHsmnTp09Hu3bt0KdPH/2y9PR0jBgxAoGBgRgxYgQyMjIEJrR+LEjFaLVaREZG4ptvvkF0dDR+/fVXnDt3TnQsu+Do6Ihp06Zhx44d2LhxI9avX891awbfffcdGjduLDqG3ZkzZw46deqE3377Ddu3b+c6NqHU1FR899132Lp1K3799VdotVpER0eLjmXT+vXrh2+++abIsq+//hrt2rXDzp070a5dOw4APAILUjEnT55Ew4YNUb9+fTg7O6N3796Ij48XHcsueHp6onnz5gAAd3d3PPnkk0hNTRWcyr7cvHkTe/bswSuvvCI6il158OABjhw5ol+vzs7O8PDwEJzKvmi1WmRnZyMvLw/Z2dnw9PQUHcmmtW7dGlWqVCmyLD4+HiEhIQCAkJAQ7Nq1S0Q0m8GCVExqaiq8vb31P3t5eXEnbgZXr17F6dOn0bJlS9FR7MrcuXMxdepUODjwV9uUrly5gurVq2P69OkICQnBjBkzkJmZKTqW3fDy8sLIkSPRrVs3dOzYEe7u7ujYsaPoWHYnLS1NXzw9PT1x584dwYmsG/8VLcbQnVcUCoWAJPZLrVYjLCwMH3zwAdzd3UXHsRu///47qlevjhYtWoiOYnfy8vJw6tQpDB48GFFRUXBzc+PhCRPKyMhAfHw84uPjsW/fPmRlZWH79u2iY5HMsSAV4+3tjZs3b+p/Tk1N5VCvCeXm5iIsLAxBQUEIDAwUHceuHDt2DLt374a/vz8mTZqEQ4cOYcqUKaJj2QVvb294e3vrRzx79uyJU6dOCU5lPw4cOIB69eqhevXqqFSpEgIDAzkJ3gxq1KgBlUoFAFCpVKhevbrgRNaNBakYPz8/XLx4EVeuXIFGo0F0dDT8/f1Fx7ILkiRhxowZePLJJzFixAjRcezO5MmTsXfvXuzevRuLFi1C27ZtsXDhQtGx7EKtWrXg7e2NCxcuAAAOHjzISdomVKdOHZw4cQJZWVmQJInr10z8/f0RFRUFAIiKikL37t0FJ7JuTqIDWBsnJyeEh4dj9OjR0Gq16N+/P3x9fUXHsgtHjx7F9u3b0aRJEwQHBwMAJk2ahC5dughORvRoM2fOxJQpU5Cbm4v69etj3rx5oiPZjZYtW6JHjx4IDQ2Fk5MTmjVrhkGDBomOZdMmTZqEw4cP4+7du+jcuTPeeecdjB07FhMnTsSWLVtQu3ZtLFmyRHRMq6aQDE26ISIiIpIxHmIjIiIiKoYFiYiIiKgYFiQiIiKiYliQiIiIiIphQTJi48aNoiPYNa5f8+G6NS+uX/PhujUvrt/yYUEyghuSeXH9mg/XrXlx/ZoP1615cf2WDwsSERERUTG8DpIRrw7ph3v3HgAouXqKLCll9UkG/mR8ScHbleX9yvieUuE/lvV9Sz7yqC2ktCRGXyqV9rmGvkvZ16HRx6TSXlfK2i32AtN+rpHXl/Lk0j7jkc8z8mDp2Ur+ZZXnvY0+vyzfsTz/Opngc4095/Hfp9y/PA8XP9Z3L++GUNrzCv0emzVLaSuhlA95jF+oR3+8iTbKsqywQs8p9d+k0v4Rfaz3lgz+0dDCDp3bYtWqVcY/x87xStpG3Lv3AIu+mAGdpNMvKygvuiL/cOQ/rpMKfgEfPv/hsqJboU6Sirxv4efkP1Z8Z130M4y9t65YPsnA83UGshQvZbpCWUq8j/6xou9R/DuUfKwoQ9+zYI1IklTk/Q19RuEMhd+n4I86lPKYgc8tsQ6kh68p8fchFX1Ofq5iOYusHxhYZrh8Fn5uaY89zFLy70OSDP/XWHZjn/fw9RKKba6lvrfBx3SPfn6RLDrD6wcSypRFv2/Rf66Bnb2h99G/zsDfefGcusIPFl9W8nNhYB2U2GD1i6VCzy++IRR/k8LfqeRGqf8uhV+iX78GN+ISQQ2+R8FzSqyrYhuZgV8GydCGYGCjLPFeZflcg798upKfa/AXqtjzDPxiScU3QF3JjfJhFl2JbQlS0UzFn/9wma5IhhKfK0kGlwFF90NSsV8+ydBnFPm7fvj43bt3IWc8xEZERERUDAsSERERUTEsSERERETFsCARERERFcOCRERERFQMCxIRERFRMSxIRERERMWwIBEREREVw4JEREREVAwLEhEREVExLEhERERExbAgERERERXDgkRERERUjJPoANZKIVXCB5OWiI5ht+7evYtq1aqJjmGX7H3dKh7xs7nZ+/q1uIK/QEeuW8MKraAKKu/6lfvfBQuSES4uLvjpp59Ex7Bb/fr14/o1E65b8+L6NR+uW/Pi+i0fHmIjIiIiKoYjSEYcPXoUCoWlB+/lo0GDBly/ZsJ1a15cv+bDdWtej7N+JUkyUxrrx4JkREZGhvANIyJuaf5/A8JM8jxrsnHjRgwaNEh0DLvEdWteXL/mw3VrXuVdv3IvqwpJdAuwUgqFQmhBKig9ADC+5au4fPkyGjRoAE9PT6PPs6WCRERE1k30flA0zkGyQoVLz8m4w2j2XHO8+uqr6Nq1K7Zu3WrweYZ+JiIiosfDESQjRDXngpKjycrB6T+OI/rzrYACcHJ2RD3PupAkCVFRUVh1citclW5wdnMp+R4cSSIiogqS+wgS5yBZkYJylHrhGuK++RkXj5+FNk8LAMjNAs7dOwcHBwe0af8iajepj6pe1dG6byd4+dQRGZuIiMjucATJCEs358IjR+tmrMCV05cgaXWlvsbTpza8nqyLPmEDS4wkcRSJiIgqQu4jSJyDZAUKzx06/cdxXP6/lEeWIwBQpdzA/+05hmO/HTRnPCIiItlhQbIiKcfPIGbZ1kc/sRBJq0P8ql+QcuJskeWcsE1ERPT4WJAEKygyZxKTsXHWKmgyc8r9HnmaPOz+NhqarBxodfIdDiUiIjIVFiSBIuKWQpOVg/2bdmHjrG+Qo85+7Pe6+ncKTsQfgTZHU+T9iYiIqPx4FpsgEXFLkXrhGuJXR+PskWSgogM/OmDHsq04e/gUuo/owzPbiIiIKoAjSBamVqsxaUMkHty9j0PbEnDr0s2Kl6N/SVodzh/9G3u+i4EmK/9QHUeRiIiIyo+n+Rth6tMb1Wo1du/ejemLPoKjkxMkSUK66i7SrqQiM/2ByT4HAFyUrhjz+WTUqPfwtiQ87Z+IiMpD7qf58xCbCTzqZrGJiYn48MMPkbB/HxwcFFBWrQznJ1yQev4aJDNMqs5RZ+Ov2IN4pnsbuLq7wVXpZvLPICIismcsSBUUvnMpHEq54fGhQ4cwePBgXLp8Kb8MKYCcrBzAzGeb7d+4G/s37kbNBl6o38wHd2+mYcXY+Wb9TCIiInvBOUgVEBFXtBwVn++jVquxePFiXL1RaKRIgtnLUWG3L6fin8PJ2PtDLNRqtcU+l4iIyJaxIJmJWq3G6dOn8ef5E8jLyRWaJfPufVxOvoDr168LzUFERGQrOEnbiEdNTivt7LCeldtg2bJlyM7ORnRsDLLuZ5ojYvkogFdmvIHNH68RnYSIiGyA3CdpsyAZUdqGUVo5Opt4CgdX70ReXh5u3kpFbrbG6HMtSeHkiGYdnkHwpMGY23eq6DhERGTl5F6QOEm7nEorRwc278butb9Cq9FaMFEZKSRkP8hCtjpLdBIiIiKrxzlI5VBaOYpb9TPiVm63znIEQMrV4UryBaRdu1Xmi0fyIpNERCRXLEhlVFpZOLXvOA5sjLdgmseTm63R39T2UViOiIhIzliQysBYWXhw9z6O/nYIUQvXWTjR47v6dwpSL1znKBIREVEpOAfpEYwVhJO7DmPPd78h/fZdSHk6C6eqAB2wed63GPjhSCBAdBgiIiLrxLPYjFAoFPho5xKDj6X8dQbrP/oaeblaQGtD5agQV6Ur+k4ego2zVhVZbqwQ8l5uRETyIvez2HiIrZw0WTmIWrQBedm5NluOACBbnY1Nkasxbtw4/RW2C5cjC17sm4iIyOqwIJVTysmzuJd6R3QMk/nqq6/QqlUrDI4cC01WDrT/NqPi95fjXCQiIpITzkEqh7OJp/Dz/zaIjmFyZ86cwdlZZ9HArzG6DO0Jn5a++sd0UsmyREREZO84glRG+9bHYv3Mr5CZ/kB0FLOQdBIunTiHDTO/xuGf9+mXl3YzXiIiInvFglSKgonJd2+kYffaGMFpSvLwrAKFo2n/CnOzNfht+VbsXR9r8HpJLElERCQHPMRWBldPXwCsbNJy4FuheK5nO9w4cxmxX21D6oVrgInmjEs6Cb+vicGl/0tB4Ji+8PKpY5o3JiIishE8zd+Iwqc3vjJjOLbOXSs40UMenlXx9soP4OzmAiD/zLq067dx/uhpnIw7jFuXUk3yOQ5OjmjcqilemfGG/rMK8LR/IiL7xtP86ZGq1/OEYyXrGWxb+9W3cHLNLywRAWFwdnNB7cZ10XHgSxi9dDLaD+wOJ5dKcHR0rNDn6PK0uPDXP0g5fkZ/dhsREZEccATJiMLNWa1W4+n2LXEl+QIkrdjV9frrr2Pt2pKjWQVzg7Q6CY4OCty9kYYrp1OQ8MNvuHP1VoU+s0YDL4z9fDJHkYiIZIQjSPRISqUSnYf0gJu7UnQUvPfeewaXRwSEQScBjv+edlatdg38NO877ImJR7t27Sr0mWmXU/HHZuu/GS8REZGpsCCV0f9GfGQVh9lmbVpk9LEip+T/O7rj5+eHuLg4LFiwoEKfu29dLFJOnC2yjGe0ERGRvWJBKqOMjAxUdn5CdAxc/edSuV+jVCrx/PPPV+yDJWD7wvUGT/0nIiKyNyxI5VDRSc+mcGRbAt5Y+E6J5YZGc1QqFfbt24fvv/8e48ePr/BnZ6Tewdkjpx/5uURERLZO/DEjG1GnTh20bNkS169fh1arFZYjT5OHg1t/xwdNFpSYNA3kn/Kfrc5Cv/eH4vjmA0hNTUVmZqbJPn/LvLXQ5ubime6tTfaeRERE1oYFqYyUSiXmzp2LO3fu4MCBA0Kz3Lp8A6kp1+Do5AS3yko4VnKEq9INfx84iUNRCdBJEm5fvAmFQoG8nFzTfrhWh20Lf0BujgZPd30ebk+4mvb9iYiIrABP8zfC2OmNKpUKI0eORHRMtNCrazs4OcLJ2QmAAjXre+LujdvIume6kaKy8GpcF33eHYhv3jE+cZyIiGyT3E/zZ0EyorQNIykpCX1eC8HlpAsWTmV9XJRuGBQxCj4tfXldJCIiOyL3gsRJ2o/Bz88Ppw6eRIchgaKjCJejzsKmWatxfNdhqNVq0XGIiIhMggXpMSmVSrw0vDeC/vsqnJ5wFh1HqOwHmdj+6Tq06PAskpKSRMchIiKqMB5iM6KsQ4sRcUtx90YaVoz/BLlqmV8jyAFo0aUVDv2yF0ql+KuOExHR4+MhNqqwnMws058tZot0wLmjf2P37t2ikxAREVUIC1IFRQSEISczGwrFo58rB9n3MjFo6GAcPnxYdBQiIqLHxoJkAjXqeeEJj8qiY1iNrHtq/Pe//4VKpRIdhYiI6LGwIJnAwoEz0Gt8Pzg/UfLK1nJ18HAixowZw0nbRERkk1iQTOTpTs/inW9nIvCtEMCRx9ukPC12796NSZMmlXr6P+/lRkRE1ogFyUQiAsLgXq0ynu/ZHnWeaiA6jlV48OAB9u3bh3Pnzhl8vKAcsSQREZG1EVKQPv/8c3Tq1AnBwcEIDg5GQkICAODkyZP6ZX379kVcXJzB10+bNg3+/v76554+nX+HeUmSMHv2bAQEBCAoKAjJycn616xZswahoaGIiYkx2/cquJK0u4e72T7D1uTk5GDp0pIFiKWIiIismdlvVpuYmIht27Zh/vz5RZYPHz4co0aNKrLM19cXW7duhZOTE1QqFYKDg9GtWzc4OZWM+d5776Fnz55Flu3duxcXL17Ezp07ceLECURERGDz5s1Qq9VISkrC5s2b8c477+Dll182/Rf9V7Y6CwBQydUZudkas32OLVm9ejW8vb0xZ84cAIbLUUTcUt6qhIiIrIZVHWJzc3PTl6GcnBwoynnufHx8PEJCQqBQKPDss8/i3r17UKlU+gtdlff9HsfHfSfDVfkET/sv5tNPP8WOHTtKHTniqBIREVkLYQVp3bp1CAoKwvTp05GRkaFffuLECfTu3Rt9+/bFrFmzDI4eAcBnn32GoKAgzJ07FxpN/khNamoqvL299c/x9vZGamoq3N3d0aRJE/Tv39+so0dA/i1I2g/0R/V6nmb9HFuTl5eHzz77DJosmV9tnIiIbILZbjUyYMAAaDQaZGZmIiMjA7Vr1wYATJkyBc2aNUO1atWgUCiwZMkSqFQqzJs3r8jrz58/j/fffx/r1q2Di0vR0+dVKhVq1aqF3NxczJw5E/Xr18eECRMwduxYjB07Fi+88AIA4I033sDUqVPRokWLcuev6CXWk5OT8eyzzyIvL++x38PeuLi7YcCHI9C41VOlPo+H2oiIxOOtRsxk8+bN2L59O2bPng1/f39s374d27dvR6dOnVCzZk04OjrCwcEBAwYMMHitnMaNG8PNzQ1nzpwp8ZinpycUCgWcnZ3Rr18//eu9vb1x8+ZN/fNu3rwJT08xIzlVq1ZF5ZpVhHy2tcp5kIVfFv+Ia/9cMvocnXx/F4mIyIoIOcRW+ArLu3btgq+vLwDgypUr+hGXa9euISUlBXXr1jX6ekmSirze398fUVFRkCQJx48fR+XKlYUWpJoNvYR8tjXLuHkH34QtQuxX2ww+7qDgXCQiIhLP7GexGbJgwQL8/fffAIC6desiMjISAHD06FGsXLkSTk5OcHBwQEREBKpXrw4AGDNmDGbPng0vLy9MmTIFd+/ehSRJaNq0KWbNmgUA6NKlCxISEhAQEAA3NzfMnTtXxNcDkD8XqfOQHjj752lAy2GRIiTg0NY9eK5HW3g2qi06DRERUQlmm4Nk60x17LXj4Jew/8d4EySyP0+2fhqDPhwOZzfDt2jhXCQiInE4B4nM6vleHaBw5Go25MKRU9gQ/jWunr4oOgoREVERHEEywpTNueuwnkj4IdYk72V3FEDlGh4YGD4K9Zo2KvEwR5GIiMTgCBKZ3Z7vf0PX4b1Fx7BOEnD/9j1ELfyB10giIiKrwYJkIe1Cu+Dpzs+KjmG10i7fwsFtewDkn+qv/fd8f57RRkREIrAgWcjcvlPR+bUeqN2kvugoVuvI9n14cPc+HBSAowPv1UJEROKwIFmQl08dDPn4TVT1riY6ilXKvKfG6f0nSyyfGbtEQBoiIpIzFiQLiggIg3u1yuj/wXA4u7uKjmN1JK0OO7+Owsn4I0WWW+Imw0RERIWxIAlQr2kj9H13EBwrCblOp1XL0+Qi7pufcfPiDf0yXl2biIgsjQXJwgpOW/dt0xyNX2gKZXUPwYmsjE7Cg7R7+P69L3Bq3/EiD7EkERGRpbAgCeLs5gL/4b3h2ZC32jAkM/0Boj/fjAd374uOQkREMsSCJEDBKFKtRnXQ861QVHKtJDiRdcpMf4DjcYlFlnEUiYiILIEFSZCIgDA4KADPRrXRaXAgwHnIBu35LqbEpG0iIiJzY0GyAp0GB+LVj99EjfqeAM/YKkKr0SLum5+LHGrjKBIREZkbC5JAhe8z5tv6aYz9Ygpe+eANeHhVFZjK+qjT7+Ps4VOiYxARkYywIFkJB0X+xO3mXZ5D9xFBouNYFUkr4diOg0i7qtLfr42jSEREZE4sSIIZult9/WY+cHKW18RtHx8fvPrqq0Yfv3oqBetnfoXozzcjNeW6BZMREZEcsSBZoWq1a6BJu6dFx7CoQ4cOYeDAgaU+586120j6/U/EfLEFmqwchO/kKBIREZkHC5IVMDSK1OW1XgKSiDNkyBCMHz/+kc+TtBIuJ53Hbyt+Qmb6fR5qIyIis2BBslKejWqjRfdWomNYTHx8PG7evFnm5//12yEsHzOPlwAgIiKzYEGyEoZGkYLCXoVX4zoC0tiGrHtq/LJ4I4YteEd0FCIisjMsSFbM2c0Foe8NQ40GXqKjWK28nFzsWRMNtVotOgoREdkRFiQrYmgUycunDoYveAeevhxJMubq6RSM/t9E0TGIiMiOsCBZGUMlyb1aZQS9Mwgu7q4CEtkACYhZ9hNSUlJEJyEiIjvBgmQj6jVthH7T3oBjJUfRUazSPdUd/KeJL9atWyc6ChER2QGFJEmS6BDWSKFQQOSqCd+5FA4Gbst2bMdBxK2IQnZ2tuVD2QAHBwdER0ejZ8+eoqMQEdk00ftB0TiCZKUMlSMA+HnRBuzYsQNt27a1bCAbodPpMHToUBw+fFh0FCIismEsSDakYH5S165d8eWXX8LDw0NwIuuUlpaGCRMm8Mw2IiJ6bCxIVqr4ZO3iPz/xxBMsSKU4cuQIVqxYIToGERHZKBYkK1ZQigyd2VanTh00b94cCoWRY3GEefPmQaVSiY5BREQ2iAXJRimVSkRGRqJBgwZwdeXp/4bcuXMH06ZNEx2DiIhsEM9iM8JWZu8nJiZi7ty52Lt3L9LT00XHsUqRkZGYOXOm6BhERDbFVvaD5sIRJBv34osvYuXKlXjppZfg5FpJdByrFB4ejlmzZhl9PCJuqQXTEBGRLeAIkhG21pyTkpIwadIkJPyRgNzsXNFxrNKFCxfg4+NTZNnM2CVw/PeaCobmehERyZWt7QdNjSNIdsLPzw9RUVEYOv9teDauLTqOVVq1alWRnyPilurLUa5Wvv8IEBFRSSxIdkSpVGL1u5/B/40+cKzkJDqO1VmzZo3BayNpdRIqOSp4qI2IiPRYkOyQT0tf+LZuBhclz24r7Nq1a4iOjoZarS5ShhyNXbaciIhki3OQjLD1Y69vfT0NsV/+hJS/zoqOYlUqV64Mj3rV0W1EEBo/62vwOZyLRERk+/vBiuIIkp1aMXY+ugztBYUTR0cKu3//Pq6dvoT1HyzHyfgjRR7TyfffASIiKoYFyY7VqOeJap41RMewSro8HaK/2IIHd+/rlxUcaeNcJCIiYkGyYwsHzkC7Af5wcuH1kQzRqLNx+OcE0TGIiMgKsSDZuWf8X8DTHZ9Fbd96oqNYpX3r4rBvw84iyzgHiYiIWJDs3Ny+U9F+oD/qNGkAz8Z1RMexSru/j8HdG2n6n3mIjYiIeBabEfY0ez8ibik0WTnIVmchZvlW/PPHSdGRrE67V7oicGxokWVanYSPe7wrKBERkVj2tB98HCxIRtjbhlEwKqLJysHXby9A2tVbghNZn+deboeeb4bC2c2lyHIeciMiObK3/WB58RCbzDi7ucCv2/OiY1ilv2IO4ut3/ofUlOuioxARkWAsSDJReBTkyVZNBSaxbmmXU/Hr0k3IzsrRXxeJc5KIiOSHBUmGvJ6sA08f3tDWmKvJKbhz/TZ4BxIiIvliQZKRglEkZzcX9Jv2OqrXqyU4kfU6eygJQP7VtXUSR5GIiOSGBUmmvHzq4M1lU9F3ymtQOHEzKO7w9r24evoiHBTgSBIRkQxxzygzheciObu5oFmnlniyZROAJaCIzHQ1Yr+KgiYrBwBHkYiI5IYFSeZc3VzQ9Y1eeKKqu+goVufqqRSc2p9/zSgHBW9mS0QkJyxIMlT8uj71mjZCyOTX4OjsJCiR9dq+aD1STpwFkF+SOIpERCQPLEgyVbwk+bZ5Gn3/+yqcXJ0FJbJSeTqsm/GlviQBLElERHLAgkR6z3Rvjddmvyk6htXRavKwYebXRUoSERHZNxYkGTN0C41Gz/wHfi+9ICCNdcvN1mBD+NdIOZ5fkjiKRERk31iQZE5rYOZxv/eGoWnnlgLSWLfcLA02zvoGV/++KDoKERGZGQuSzBm6W31EQBhOJxxH2wFdLR/IyuWosxHzxRY8uHufo0hERHZMIcn5Vr2lkNNdjEvb0a8YPx+qczcsmMY2ePvWRciUoVgxdr7oKEREZiGn/aAhHEEig3ORCrQN7mLBJLbj5tlr+OmT7zFl0xzRUYiIyAxYkMiggqlJjZ5pAgcnR7FhrJTqwnX88MEKJCUliY5CREQmxkNsRshxaLH4obaCkaXAN0OQsGYHNBqNiFhWz/M/dZG0/xg8PT2NPqdg3ZY2WkdEZE3kuB8sjCNIZFRE3FJExC3Fzq+i8Pfff2PixImiI1kl1blr6NGjBxITEw0+zsncRES2hyNIRsi1ORfemeskIDLw4YjHtWvX0K1bN5w9ywsmGtK4cWOsX78ebdq00S8zVI44ikREtkCu+8ECHEGiIgrvvAuXIwCoWrUqGjRoYOlINuPixYuIjIyEWq0GwJEjIiJbxoJEJegkw6McSqUS48aNE5DINmi1Wuzfvx+HDx/WlyNDF+JkcSIisn4sSFRC8ZGjwnr16oW2bdtaMI1tSU9PR+DLPXBgy24AgKODosLvWTAXjIiILIcFicpFqVRiyZIlqFatmugoVisvJxcJP8Ti7o00o88pa+Ep/LzTKXsqGo2IiMqIBYnKrU2bNli1ahUcHLj5GCQBmqxsnD5wstSnPaokFTyu+HcQatP50t+PiIhMh3s4eiyhoaGYM2cOnJ2dRUexThKw57sYXD190XRvKQHhO3mojYjIEliQ6LFNmzYNf//9N2bMmAFUfKqN3cnN0uDXzzci7aoKmqycEo/rJOOjSIWXa3UPlzsogNnxn5s8KxERFcXrIBkh9+s/lFeBp3V0AAAgAElEQVTHwYHY/2Oc6BhWyee5JniiihJtQ7uiXrNGJR4vfsZg4XL04O593E29AxdXZ9xLywAkoGqdmuhd3xf/qd8WderUgVKpNPdXICIZkvt+kAXJCLlvGOWVkpKCxr6NIWm5zopzcqkER0dHOLlWQuveHfDk80/By6cunN1c9M+JCAiDWq3GzJ//B1elGwDgyC/7kBi1D5qsbORkZgMFq1YBVHJ1RrXKVVGvXj1MnDgR/fr1Y1EiIpOS+36QBckIuW8Yj6P7yD7YvSb64Y6c8ilQdJ04KODZqDaC3h2kH1HqWbkNXp84Cln3s+Dk5IQHGfdw+1JqmT/i+eefx7fffgs/Pz+TRici+ZL7fpAFyQi5bxiPq/fEV/Hbsi3Q5WlFR7F6lWtWQej7w1CpkhN++ODL/FGiCggICMC2bds4kkREJiH3/SALkhFy3zAe1wc/L8CGiK9x6fh5rr8ycHSpBDelGx7cuVfh93J2dsaWLVsQFBRkgmREJHdy3w/yLDYyqbl9p6L7iCC4V68sOopN0ObkmqQcAYBGo8Fnn32mvxccERE9PhYkMrl6TRshdNrrLEkCHDhwADExMaJjEBHZPBYkMrmIgDA0fMYX/aa9UeRMLTK/nJwcjBo1Cjt27BAdhYjIpnEOkhFyP/ZqChFxS3H174v4c9VunDhxQnQcWXF2dsaSJUswbtw40VGIyEbJfT/IESQyq2/eWYSdO3eiTZs2oqPIikajwVtvvYXFixeLjkJEZJM4gmSE3JuzqUTELUVEQBiSkpLQNegl3LmkEh1Jdj777DNMnDhRdAwisjFy3w+yIBkh9w3DVArfNkN18QZWjJ0vMI18zZgxA7NnzxYdg4hsiNz3gzzERmZV+D5jno1qo/NrgQLTyNecOXMwaNAg0TGIiGwGCxJZVIeBL6HxC01Fx5ClTZs2ITY2VnQMIiKbwIJEZld4FMnZzQVdh/WCa2U3gYnkKzg4GFu3bhUdg4jI6rEgkcXVa9YI3Uf0ER1DlnJycjB58mSoVJwsT0RUGhYksojCo0gA0LRDS7i6cxRJhEuXLmH37t2iYxARWTUWJLKYgpKk1Ulwr1YZXUf0FpxIvsLDwzmKRERUChYksiidBDg6KAAALwZ1QuuQzoITydPZs2fRr18/JCUliY5CRGSVWJDIoiID80eRCkaTXh7fH091bikykmzt378fQ4YMgVqtFh2FiMjqsCCRxemKXXes3+TX4F6zipgwMvd///d/GDp0qOgYRERWh1fSNkLuVxA1t4IrbOskQJIk3L50Az/NXwtVyk3ByeRp/vz5eP/990XHICIrIvf9IEeQSBitToKDIn9OkpdPHQybPwG1feuJjiVL06ZNQ0JCgugYRERWgwWJhIgICNNP1i7gXq0yOgzsDkWx5WQZAwYM4HwkIqJ/sSCRVfFt0xw1GniJjiFLt27dwttvvy06BhGRVWBBImGKXzwSyL8VScDoYAFpCADWrl2LxYsXi45BRCQcCxJZnSZtnoY/LyIpzOTJkzkfiYhkj2exGSH32fuWVHBGWwGtToKjgwIX/voH6z/8CtpcraBk8tW1a1f8+uuvUCqVoqMQkSBy3w9yBImEK36orWDydr2mjeDz3FNwcHQUEUvW9uzZg9WrV4uOQUQkDAsSWS1nNxe8NCoI9Zv7wMGZJcnSwsLCsGPHDtExiIiE4CE2I+Q+tChC4YtHAkDB2f6arBwk7TmG6KWbIGl1gtLJk4uLC7Zt24ZevXqJjkJEFib3/SBHkMjqOCjy/6eT8v/n7OaCuk3q8/pIAuTk5CA4OBjr1q0THYWIyKJYkMhqFJ+LVNCH0h5oUbVOLbRq+ZyAVJSbm4uJEydCpVKJjkJEZDEsSGTVHBRANaUj5vedii+++AKVK1cWHUmWbt++jePHj4uOQURkMZyDZITcj72KVPy0/8IjSzt27EBISAg0Go2lY8nek08+ifPnz4uOQUQWIvf9IEeQyKoVP+zWq1cvrF69Gu7u7oISydeFCxfw3nvviY5BRGQRHEEyQu7NWbTwnUsRGVjyViQFYmJi0Ls3r7YtwvTp0zFjxgxeRJLIzsl9P8iCZITcNwxrp1ar8cwzz+DChQuio8hSq1atsGbNGvj5+YmOQkRmIvf9IA+xkU1SKpXYsGGD6BiydezYMcydOxdqtVp0FCIis2BBIpsVk3EI3UcHiY4hW3/++SfS09NFxyAiMgsWJLJJBWe6/ef5plA4cjMW4dy5c/jhhx9ExyAiMgvuWcjmFL4MQPU6tVCjXi2BaeRt2rRpvKktEdklFiSyKRFxS/X3aitQuUZVMWEIADBhwgSkpKSIjkFEZFIsSGQzil9AEgCy1Vm4demGgDRUICsrC23btkVSUpLoKEREJsOCRDahcDkqfM/apL3H8SDtnoBEVJhKpcKYMWN4vzYishssSGQTil9RGwA2zV6NXV9uE5CGDElMTMTrr7/OkSQisgssSGSTrv1zCaf3nhAdg4r5888/sXDhQl4fiYhsHgsS2YzCo0jn/jwtMAkZk5aWhq1bt+LgwYOioxARVQgLEtmkKrV45pq1UqvVmDBhAkeRiMimsSCRTSkYRfpP6+aCk1Bp/vnnHyxZskR0DCKix8aCRDYpQ3UHUDz6eSTOggULkJycLDoGEdFjYUEim+N+QoM1Uz4H5HuTaZuQnp6Ojh07Yv369aKjEBGVm0KSJO5mDFAoFOCqsT4pKSlo1aoVMh7cg5SnEx2HykChUCAmJgY9e/YUHYWIykHu+0GOIJFNOX36NDJzsiDpWI5shSRJCA4ORmJiougoRERlxoJENiX2dmL+kTX2I5ui0WgwdepUntlGRDaDBYlsRkTcUlSrXQMt/J8XHYUew759+3D48GHRMYiIyoQFiWxC4Xux1XuqgcAkVBEvvfQSFi1aJDoGEdEjsSCR1StcjgCgkquLoCRUUTqdDu+//z6OHDkiOgoRUalYkMiqzYwtebFBZzcWJFuWl5eHrl27YuvWraKjEBEZxYJEVs3RQQGtLv8003//A11ursBEZAqZmZmYOHEiUlJSREchIjKIBYmsWkRAGBwdHl4yWycB3o3rQ+HETdfWXb16FSNGjEBSUpLoKEREJfBCkUbI/QJZ1mpm7BI4OiiQsD4We9bEiI5DFeTk5IQ+ffrghx9+gFKpFB2HiAqR+36QBckIuW8Y1qb4RG2dBBzfeQi//G+DoERkKlWrVsWOHTvQtm1b0VGIqBC57wd5nIJsQkRAWJGfHRRAqx5t0Tqkk6BEZCrp6enw9/fnpG0isiosSGQzFA+nIuknbL88/hU0addcTCAymaysLEyYMAEqlUp0FCIiACxIZEM+eunhKFLBvG2dBPSf9gaeqO4uKBWZys2bNxEVFSU6BhERABYksnEOivzrIvV6e4DoKGQC48aNw44dO0THICISU5A+//xzdOrUCcHBwQgODkZCQgIA4OTJk/plffv2RVxcnMHXX7lyBQMGDEBgYCAmTpwIjUYDIP+GmBMnTkRAQAAGDBiAq1ev6l/zySefoF+/frwXlI0rPhepQLMOLdGk4zMWTkOmJkkSXn75ZaxevVp0FCKSObMXpMTEREybNq3E8uHDh2P79u3Yvn07unTpAgDw9fXF1q1bsX37dnzzzTcIDw9HXl5eidcuXLgQw4cPx86dO+Hh4YEtW7YAADZv3gwPDw/ExcVh+PDhWLhwIQDg/PnzAIB169Zh3bp15vqqJJCjgwL9pw5FrUbeoqOQCbz99tv6/+NERCSCVR1ic3Nzg5OTEwAgJycHisKzcv8lSRIOHTqEHj16AABCQ0MRHx8PANi9ezdCQ0MBAD169MDBgwchSRJ0Oh0cHBxkf8qivYgICIOjQvHwytr//tfZzQX9p7+BBi2eFBeOTCI7OxszZ86EWq0WHYWIZEpYQVq3bh2CgoIwffp0ZGRk6JefOHECvXv3Rt++fTFr1ix9YSpw9+5deHh46Jd7e3sjNTUVAJCamoratWsDyL8AXeXKlXH37l34+voiOzsbQ4YMweDBgy30DcmccnWSfqK2g+JhSfLyqYNeYQPhWMnJ+IvJJhw5cgQxMbwYKBGJYba9yIABA6DRaJCZmYmMjAwEBwcDAKZMmYLBgwdj/PjxUCgUWLJkCebPn4958+YBAFq2bIno6GicP38e77//Pjp37gwXl9JvTlow0mRodKjgsZkzZ5ry65FgkYFhRS4eWehuJPBuVButQzrh0ObfBSQjU8nOzsY777yDhg0bok2bNqLjEJHMmG0EafPmzdi+fTtmz54Nf39//XyjTp06oWbNmnB0dISDgwMGDBhg8F5MjRs3hpubG86cOVNkebVq1XDv3j393KSbN2/C09MTQP5o0o0bNwDk3zH8/v37qFq1qrm+Igk26D/GJ2W/0Ku9BZOQuaSmpiIsLIyH2ojI4oQcYit8Mbhdu3bB19cXQP7ZaQXF59q1a0hJSUHdunWLvFahUODFF19EbGwsAGDbtm3w9/cHAPj7+2Pbtm0AgNjYWLRt29bgPCayD818uhY5q01XaABRB8DDk+XYHiQmJmLNmjWiYxCRzAgpSAsWLEBQUBCCgoJw6NAhTJ8+HQBw9OhR/Wn+EyZMQEREBKpXrw4AGDNmjH6u0dSpU/Htt98iICAA6enpGDAg/xo4r7zyCtLT0xEQEIBvv/0WU6ZMEfH1yMIiAsKgk4oeZqtcowrq+NYH2I/twoQJE/Dhhx+KjkFEMsKb1RrBM95sy+mUPdh47qS+KOkkIDXlOmJX/IRLJ86KjkcmMmvWLISHh4uOQSQLct8PsiAZIfcNw1ZFxC3VlySFAsjJzMHG2d/iwpHToqORibz77rtYvHix6BhEdk/u+0Grug4SUUVFBITpD7VJUv61kUKnvCY2FJnUkiVLMGPGDNExiMjOsSCR3YkICIO20IztB3cyOBfJznz66adITk4WHYOI7BgLEtmlj3u8CycHBRQKICszW3QcMrG8vDz07t2b91YkIrNhQSK79WH3d5CnlVCrnhdc3Z8QHYdM7NKlSwgNDUViYqLoKERkh1iQyK593ONdVK3hge5j+oqOQmZw/fp1TJ06lReSJCKTY0Eiu/dh93fwXGBbBLwZIjoKmcG+ffswdepU0TGIyM7wjp4kCx/3eBcLXZ2gzc3D7tW/io5DJrZixQpkZmZi2bJlUCqVouMQkR3gCBLJxpQub6PDoAD4vdRadBQyg7Vr16JVq1YG7+1IRFReLEgkK5GBYQgcEwz3Gh6io5AZnDlzBn369Clyv0ciosfBgkSys3DgDAS+1Q+OrjzCbI8uX76M4OBgjiQRUYXwViNGyP0S63IwctG7+HbyUtExyEy6deuGX375hXOSiB6T3PeDHEEi2Vo9aQleeeUV0THITA4cOMALSRLRY2NBIllbs2YNGjZsKDoGmUFOTg4GDx6M1atXi45CRDaIBYlkTalUYuPGjahbt67oKGQGqampGDduHFatWiU6ChHZGM5BMkLux17l5vDhwxgyZAjOnz8vOgqZQY0aNXDq1Cl4enqKjkJkM+S+H+QIEhGANm3a4MSJE6hf31t0FDKDtLQ0rF27VnQMIrIhLEhE/1Iqlbh8+QZaPNNEdBQyg/feew/z588XHYOIbAQLElExSSf+QWBgoOgYZAbTp0/HmDFjRMcgIhvAOUhGyP3YKwFNmjTB2bNnRccgM1ixYgXGjRsnOgaRVZP7fpAjSERGxMbGio5AZvLWW2/hvffeEx2DiKwYCxKRET4+Ppg1a5boGGQmCxYswIcffig6BhFZKR5iM0LuQ4v0UEREBIuSnVIoFEhKSkLz5s1FRyGyOnLfD3IEiegRpk6dilatWomOQWYgSRJeeeUVJCYmio5CRFaGI0hGyL05U1GHDx9G165dkZWVJToKmUHNmjURHR2NNm3aiI5CZDXkvh/kCBJRGbRp0wYxMTHw8PAQHYXM4Pbt2xg9ejTUarXoKERkJViQiMqoa9euWLNmjegYZCZJSUk4efKk6BhEZCVYkIjKITQ0FHPnzoWDA3917NHy5cuhUqlExyAiK8A5SEbI/dgrlS45ORkdO3ZEenq66ChkYp6enpg3bx5GjhwpOgqRUHLfDz6yIKWlpeHYsWNQqVRwcXFBkyZN0KJFC7v/f9By3zDo0fbs2YNu3bqJjkFmEhkZiZkzZ4qOQSSM3PeDRgvSoUOHsHLlSqSnp+Ppp59G9erVodFokJKSgitXrqBHjx4YOXIk3N3dLZ3ZIuS+YVDZzJs3Dx988IHoGGQmo0ePxsqVK0XHIBJC7vtBowXpk08+wbBhw1CnTp0Sj+Xl5WHPnj3QarXo0aOH2UOKIPcNg8pGrVajU6dO+Ouvv0RHITNZv349Bg8eLDoGkcXJfT/IOUhGyH3DoLJLSkpCcHAwUlJSREchM6hVqxY2b96MLl26iI5CZFFy3w+WWpD27duHXbt2ITU1FQqFAp6enujevTs6d+5syYxCyH3DoPJRqVTo0KEDzp07JzoKmYGPjw82bNiAF198UXQUIouR+37QaEGaM2cOLl68iJCQEHh5eQEAUlNTERUVhYYNG9r9TR7lvmFQ+e3ZsweBgYHIzc0VHYXMoGnTpkhISICnp6foKEQWIff9oNGC1KNHD8TGxpZYLkkSevTogZ07d5o9nEhy3zDo8Sxfvhxvv/226BhkJp6enti1axf8/PxERyEyO7nvB42eq+/s7GzwqrJJSUlwcXExaygiWzV+/HgsXrxYdAwyE5VKhcmTJ/OWJEQyYHQEKTk5GREREVCr1fD29gYA3LhxA+7u7vjoo4/QokULiwa1NLk3Z6qY5cuXIywsDFqtVnQUMoPIyEhMmjQJSqVSdBQis5H7fvCRZ7HdunULqampkCQJ3t7eqFWrlqWyCSX3DYMqbseOHQgNDUVOTo7oKGQG3bp1w5IlS3i4jeyW3PeDpRak+/fvY9++fUXOYuvYsaMs7mgu9w2DTGPbtm3o378/tyU79cILL2DPnj0cSSK7JPf9oNE5SFFRUQgNDUViYiKysrKQmZmJQ4cOoV+/foiKirJkRiKbFRoaijlz5sDV1VV0FDKDY8eOISEhQXQMIjKDUs9i27x5c4nRooyMDAwcONDgGW72RO7NmUwrJSUFc+bMwapVq0RHIRPz8PDArFmzMHHiRNFRiExK7vvBUu84q1AoSr7AwUHWK4zocfj4+GDJkiVo37696ChkYvfu3cN///tfjB8/XnQUIjIhoyNI27Ztw7Jly9ChQwfUrl0bAHD9+nUcOHAA48ePR79+/Swa1NLk3pzJPPbs2YPevXsjMzNTdBQyg4SEBFncaYDkQe77wVInaWdkZOCPP/4ochZbx44dUaVKFUtmFELuGwaZz7p16zBmzBhkZWWJjkIm5u7ujjVr1qB///6ioxBVmNz3g488zf/27dtFzmKrWbOmpbIJJfcNg8wrOTkZzzzzDHQ6negoZGLVqlXD0aNH4ePjIzoKUYXIfT9otCCdPn0aH330Ee7fvw9vb29IkoSbN2/Cw8MDH330EZo3b27prBYl9w2DzG/dunUYOnSo6BhkBk899RTWrl3Lm9uSTZP7ftBoQQoODkZkZCRatmxZZPnx48cRHh6On3/+2SIBRZH7hkGWsWPHDvTt2xd5eXmio5CJubu7Iz4+Hm3atBEdheixyH0/aPQstqysrBLlCACeffZZzp0gMpFevXph165dsrj4qtw8ePAAkyZN4n3biGyU0RGk2bNn4/LlywgJCdHfi+3mzZuIiopCvXr1EB4ebtGglib35kyWdfjwYQwfPhynT58WHYVMbNiwYfjuu+9ExyAqN7nvB0udpJ2QkID4+HioVCpIkgQvLy90794dXbp0sWRGIeS+YZDlqdVqLF26FJGRkcjOzhYdh0xo6tSp+PTTT0XHICoXue8HH3kWm1zJfcMgcZKTk9G+fXvcu3dPdBQyobfeegvLly8XHYOozOS+HzQ6B0mr1eLHH3/E4sWLcezYsSKP8ZecyHyaN2+ORYsWiY5BJrZixQokJyeLjkFEZWS0IIWHh+PIkSOoWrUqZs+ejXnz5ukfi4uLs0g4IrkaNWoUVq1ahRo1aoiOQib05ptvctI2kY0wWpBOnjyJ//3vfxg+fDg2bdqEzMxMTJgwARqNRtZDbkSWMnLkSBw5cgRt27YVHYVMZP/+/WjXrh2SkpJERyGiRzBakHJzc/V/dnJywscff4ymTZvi9ddf532kiCzEx8cH8+bNg6urq+goZCJJSUl46623oFKpREcholIYLUgtWrTA3r17iyybMGEC+vfvj2vXrpk9GBHl69q1K5YtWwYHB6O/rmRj9u/fj759+3IkiciK8Sw2I+Q+e5+sz4YNGzBkyBDRMciEOnfujJiYGCiVStFRiEqQ+37QaEHauXOn0Rc5Ozujfv36aNy4sdmCiSb3DYOs0/z58zFjxgze5NaODB06FN9//73oGEQlyH0/aLQgTZ8+3eiL8vLycP78ebRq1Qoffvih2cKJJPcNg6zX3r170b17d96/zY6sXr0aI0aMEB2DqAi57wcf+xCbTqdDUFAQoqOjTZ3JKsh9wyDrNm/ePMyYMYPbqB2JjIzEzJkzRccg0pP7ftBoQdq+fTuCgoKMTgy9fPkyVCoVXnjhBbMGFEXuGwZZv23btmHQoEFFzjgl29anTx/8+OOPnJNEVkHu+0GjBWnt2rXYunUrmjdvjhYtWqBatWrQaDS4dOkSjhw5gmrVqmHy5Mlo1KiRhSNbhtw3DLINv/32G3r37s05SXakZ8+e+PTTT+Hn5yc6Csmc3PeDpR5i02q1OHToEI4dO4Zbt27BxcUFjRs3RufOnVGnTh1L5rQ4uW8YZDuWL1+Ot99+W3QMMqFu3brhl19+KdNIUkTcUkQEhFkgFcmN3PeDPM3fCLlvGGRbli9fjnfffZcTt+3Id999h2HDhpX6nJmxS+DooAAAliQyObnvB3nlOSI7MH78eGzatAmVKlUSHYVMZPz48dizZ4/Rx8N3LtWXIyIyPRYkIjsRGhqKpUuXQqHgTtMePHjwAIMGDUJiYmKJx2bGLkHhbqST8gsTEZkOCxKRHRk3bhxWrlzJe7fZCZVKhfDwcCQnJ+PPP//U37+t+MgRB5KITO+RBWnRokW4d++e/ueMjAx89tlnZg1FRI9v1KhROHXqFJo2bSo6CpnAzp070aJFC7Rv3x5NmjTBS2ODkZWZDa3u4dwQnZRfkjiKRGQ6jyxIe/fuhYeHh/7nKlWqlLiJLRFZFx8fH3z55ZeoUqUKnJ2dRcchE8jNzUVGRgbiV/6MH95fhtuXbkAnPSxHQP5/I+JYkohM4ZEFSavVQqPR6H/Ozs4u8jMRWacuXbpg2bJlqF69uugoZGLX/7mMX5ZsRF52DoD8kkREpvXIgtS3b1+88cYb2Lx5M7Zs2YIRI0YgJCTEEtmIqIJee+017Nq1C1WrVhUdhUzs2qmLOHskGUDJOUgcRSKquDJdB2nv3r04ePAgJElChw4d0KlTJ0tkE0ru138g+7J161aMHz9eP8mX7EONep4Yu2wKnFxd4KAoeriN10WiipL7frBMBenatWu4dOkS2rdvj6ysLGi1Wri7u1sinzBy3zDI/kRHRyMoKIjbtZ3pM3EQnn+5vcHHWJKoIuS+H3zkIbZNmzYhLCwM4eHhAIDU1FTe1oDIBnXt2hWtW7cWHYNMLPqLLTi24yAA6CdtE1HFPbIgrVu3Dhs2bNCPGDVq1Ah37twxezAiMi2lUolvvvkG3bp1g5OTk+g4ZCJSnhY7lm3FmSOn4KAoOh+Jc5GIHt8jC5Kzs3OR04R5ryci2+Xn54cff/wR7dq1Ex2FTChPk4vNs79FyvGzoqMQ2Y1HFqTWrVvjyy+/RHZ2Nvbv3493330X/v7+lshGRGbg6emJ6dOni45BJpaXpcF373+Bw7/+UWQ5R5GIHs8jJ2nrdDps2bIFf/yR/0vXsWNHDBgwwO7v9yT3yWlk/6ZOnYqFCxeKjkFm0LBlYwxf8HCCtlYn4eMe7wpMRLZI7vvBMp3FVjDnSE4XnJP7hkHyMHr0aKxatUp0DDKD515uh8A3Q+H87yUAeEYblZfc94NGC5IkSfjiiy/www8/6H92cHDA0KFDMWHCBIuGFEHuGwbJg1qtRqdOnfDXX3+JjkJmUNW7BnpOeAVPtXkaAEsSlY/c94NG5yCtXbsWx44dw5YtW5CYmIjDhw9j8+bN+Ouvv7BmzRoLRiQic1Eqlfjyyy9RrVo10VHIDNJvpuHHD7/Cnu9iREchsjlGR5BCQkKwevXqEofV7ty5g5EjRyIqKsoiAUWRe3MmeUlISED37t2h1WpFRyEzCZk2DH7dXkBkIEeRqGzkvh80OoKUl5dncM5R9erVeao/kZ3p0qULR4YtrHp9L4t+XtT873Hst4MW/UwiW2a0IFWqVMnoi0p7jIhsU2hoKDp27Cg6hixU8a6OFl2eQ8NnGlv0c6MX/4g3Fr5j0c8kslVGD7E1a9YMbm5uJZZLkgSNRoPk5GSzhxNJ7kOLJE9JSUkYO3YsDh06JDqKXXOtokSrnu1wJfk8rvxfikU/u5KbC+7eSoNSqbTo55Ltkft+sEyn+cuR3DcMki+1Wo1p06bhiy++EB3FrlWuURW5eRpkZ2Ra/LOf7vockn8/ZvHPJdsi9/0gC5IRct8wSN7UajWaNm2Kq1evio5CZhITE4NevXqJjkFWTO77wUfeaoSI5EepVGLlypWiY5AZDRk1FCkplj28R2RLWJCIyKCePXti4MCBomOQmaTfuIPnnnsO69atEx2FyCrxEJsRch9aJALyD7XVr18fd+/eFR2FzMTV1RW//fYbunTpIjoKWRm57wc5gkRERimVSiQkJKBFixaio5CZZGdnY/To0VCr1aKjEFkVFiQiKpWfnx9+/vlnVKlSRXQUMpNz585h+vTpomMQWRUWJOXbLo4AACAASURBVCJ6JB8fH8yePVt0DDKjL7/8kpO2iQphQSKiMhkxYgQ6d+4sOoZNs+ZRuNzcXGzcuFF0DCKrwUnaRsh9chqRIYmJiejduzfS0tJERyEzqFGjBuLj4/Gf//yHV9om2e8HOYJERGX24osvIiYmBrVr1xYdhcwgLS0NAQEBGD58OJKSkko8Hr5zKSLilgpIRmR5LEhEVC5t2rTB/v378cQTT4iOQmZw69Yt7Nu3D4sXLy5xZpuDIv+/LEkkByxIRFRuPj4+WLBggegYZCapqan4+++/kZ6erl/GUkRyw4JERI9l/PjxmDVrFipVqiQ6CpnBgQMHcPjwYQBFy5Hu3ykpLExk71iQiOixhYeHY//+/XB2dhYdhcygf//+6Pp6L2h1DyfqFhxmI7J3LEhEVCGtW7fGjBkzRMcgM5AkCQnf/4ZNs/JvXKwrdkJT+E6OIpH94mn+Rsj99Eai8vL19cW5c+dExyAz6TdtGPz8XyixXCcBkYFhAhKRucl9P8gRJCIyiaioKNERyIziVv4CTVZOiVEkHnIje8WCREQm0bx5c4wePVp0DDKT+2npSLt+m4WIZIMFiYhMZuXKlXj11VdFxyAzyUi9U2IECeAZbWSfOAfJCLkfeyV6XCqVCn5+flCpVKKjkIlVesIFoz6biFqN6uiXFR5RigjgXCR7Ivf9IEeQiMikPD09sXz5cri6uoqOQiaWm5mD9TO/xu1LN+CgyC9HhkaUiOwBR5CMkHtzJqqo5ORktG7dGllZWaKjkIk5ODmg24ggtOvfDQDg+O8wEs9osy9y3w9yBImIzKJ58+b44osvRMcgM9Dl6RC/ajvupd6BQqGAVidBJ0HWO1OyPyxIRGQ2I0eOxKJFi0THIHPQAevDv4SDIn8EqeC/nLBN9oIFiYjM6r///S+GDRsmOgaZwe1LKhyLPVRkGeckkb3gHCQj5H7slciUVCoV6tSpA61WKzoKmVi1OjUxbsV7cHZzKbKcZ7TZPrnvBzmCRERm5+npibVr14qOQWaQmf4AGWkZ+hvacgSJ7AULEhFZxGuvvYYVK1aIjkEmlpOZjex7D6BQKKCTHl4XiTeyJVvHQ2xGyH1okchcXn75ZezYsUN0DDIhhaMCXYf1QpvQrnB1cylSlHiozXbJfT/IESQisqiYmBh06NBBdAwyIUkr4fc1MVgVthCpKddFxyEyCRYkIrK4n376CW5ubqJjkIndvqTCL4s3QpOdo5+LxNP+yVaxIBGRxXl6euKTTz4RHYPM4Nrpi7hz7RYATtgm28aCRERCjBw5Ek899ZToGGQGB7fE6/+skziKRLaJBYmIhFAqlfjss89ExyAz+L/dx3D4l328oS3ZNBYkIhKmV69emDVrlugYZAYHfozDg7v3OReJbBYLEhEJFR4ezusj2aHM+2rcTb0DSZL0p/wT2RIWJCISbty4cQgJCREdg0xIm5OHxG2/Q6F42I44ikS2hAWJiKzCDz/8gAYNGoiOQSaU/PtfuHbmkv4wG+cikS1hQSIiq6BUKrFq1SpUrlxZdBQyofUffoWrf18EkD9hm7cgIVvBgkREVqNdu3YIDQ2Fh4eH6ChkItkZaqydvERfkgAeaiPbwIJERFZDqVRiypQp8PPzEx2FTEiXp0PcN79Ak5XDCdtkM1iQiMiq+Pn54fvvv4dSqRQdhUzo/p0MZD7I0v/MUSSydixIRGR1fHx8sHQpd6D2xL1qZbgoXUXHICozFiQiskojR45EZGSk6BhkIu0G+MPtifyCxLPZyBawIBGR1Zo5cybeeOMN0TGoApTVPfBUez80fraJfhnnIZEtcBIdgIioNGvWrIG7uzuWLVsmOgqVU62GXqjb1AetQzrD2c2lyGMRAWGCUhGVjUKSJA52GqBQKMBVQ2Q9unbtioSEBNExqIyqeFXFa3PfhnsND7i4uRYZNWI5sg1y3w/yEBsR2YQ9e/Zg2LBhomNQGVRydcYrM0agRj1PuD3BckS2iSNIRsi9ORNZI5VKhRYtWuDWrVuio5AR7jWrYNDMkajXrBF0EliObJjc94NCRpA+//xzdOrUCcHBwQgODtYPm+/fvx/9+vVDUFAQ+vXrh4MHD5br9QDw1VdfISAgAD169MC+ffv0y6OjoxEaGoo1a9aY9bsRkfl4enpi9uzZcHXl6eLW6tWPRunLUWEsR2RrzD5JOzExEdu2bcP8+fOLLB8+fDhGjRpVZFm1atWwYsUKeHl54cyZMxg1alSRkvOo1587dw7R0dGIjo5GamoqRowYgdjYWDg6OiI6OhpbtmzBlClToFareRE6Ihv12muv4dNPP8X58+dFR6Fi6j7dEHWfalhiOcsR2SKrmoP09NNPw8vLCwDg6+sLjUYDjUZT5tfHx8ejd+/ecHZ2Rv369dGwYUOcPHkSAPTDhHIfMiSydWFhYSxHVsqn0Kn8QP7hNZYjslXCCtK6desQFBSE6dOnIyMjo8TjsbGxaNasGZydncv8+tTUVHh7e+uf4+XlhdTUVABAYGAg+vfvjxYtWsDd3d0M34iIzC05OZmHya1Y03Z++kNrvNYR2TqzHWIbMGAANBoNMjMzkZGRgeDgYADAlClTMHjwYIwfPx4KhQJLlizB/PnzMW/ePP1rz549i4ULF2L16tUG39vY6w2NDCkU+b+loaGhCA0NNcM3JSJLOXr0KHQ6negYZECzTs+UOLzG0SOyZWYrSJs3bwZgfA5SgQEDBmDcuHH6n2/evIkJEybgk08+QYMGDQy+pmbNmgZf7+3tjZs3b+ofS01NhaenZ4W/CxFZh+effx4KBwUk3qvCqrTq3R69wwYBeHgbkchAliOybUIOsalUKv2fd+3aBV9fXwDAvXv3MHbsWPx/e/ceF1W1/g/8s2eQiwMooQiI91BRKc+xvPxKTWwwRES8pBUdb6VddFJK00wOYaWmnRLL0vJ2UvuWmjfwhrcyj5p6zFtq1JHUVCgFQW4DzPr9QTPOAENiwN4z+/N+vXwJe/bMPHs5sh6etfZacXFx6NKlS7WfHxYWhpSUFBiNRly6dAnp6em47777aukqiKiudezYEW17dJI7DCon53q2zZAakyNyBrJsNTJv3jycO3cOANC0aVPLhpSrVq3CxYsXsWjRIixatAgAsGzZMvj6+mLGjBkYMWIEQkND7T4/ODgYERER6N+/P7RaLeLj46HVamW4QiKqDfE7k6AfOxDnD5ySOxSy8tOhH5CZfhWNWgQwOSKnwYUi7eDdbkTKkpCaBAAoLCjCv6csxNUfL8kcEVlz9XTH0f98h9DQULlDoRqi9n5QUbf5ExFVZuaOBTCJsvkt7h5ueOzFodC4sDqsJMZbhegz4FHk5eXJHQpRjWCCRESKdvbCPmg1EjTS7VvHm4e0xMh3DWjcuLG8wZGN6xczsWHDBrnDIKoRTJCISNG++OmkpXoE3P57meFfOH36NAIDA+ULjip4dsJzrCKRU2CCRESKlqA3WCpH5s1Pzevr+Pn5Yfv27QgKCpIxQrJWeDMPg6eOlDsMor+MCRIROQTzEFv5xQdDQ0Nx7NgxDrcpyK5PN2L79u1yh0H0lzBBIiLFS9AbKuwOb83Pzw/vvvtu3QVEVTIZSzFgYJTd3RCIHAETJCJyGFVtXTF48OAqF5ilulVaXIKJEyfi66+/ljsUorvCdZDsUPv6D0SO6NSpUxg/fjwOHjwodyj0h969eyMlJQU6na7Sx83rW3HfNuVRez/IChIROY3Q0FCsXr0anvd4yx0K/eH48eO4cuVKpY/F70yq42iI7hwTJCJyKq1atcLHSYsgafnjTQlycnIsm5dbi9+ZZFOdMFeSiJSCQ2x2qL20SOToRs6fiF2fbsaVHy8C/K8sq3r13XD5wkX4+fkBKEuGzJPuNdLt5RsADrUpidr7Qf6KRUROaeUrCxE790X0GRUpdyiqV5xfhIULF1q+t06OiJSKFSQ71J45EzmLhNQkfPTiPGSmXZY7FNUzGAzQ9Q2Cq4ebzXHrCpJJAInhrCIpgdr7QVaQiMjpZfx4CR07dpQ7DNVLSkrC3KHTceFEmuWYdXIEsKpEysEKkh1qz5yJnM2JEyfQuUtnoFTuSAgAIiYMRdeBPe0+zrlI8lN7P+gidwBERLUtITUJxoIiNGvfCpfOXJA7HAKw7YN1yMvKRY/H+8L9jyE3zk0iJWEFyQ61Z85EziZ+ZxKunE/H+rdXIvvaDbnDoT80aRuEAS8OQ1BIS5SaBCRJ4h1tCqH2fpBzkIhINQLbtcSQ6SMR2L653KHQHzJ+vIxVr3+MH749Aa1GYvWIFIMVJDvUnjkTOaP4nUnQSICxoAjfJX+L/at2wFhQJHdY9Ieesf0Q9o/+lqE2IQRm9XtJ3qBUTO39ICtIRKQqJgG4erjh4WF98dzHr8KzEbclUYr9q3Yg7bsfkBhugEYq66CJ5MIEiYhUw9zxmisUPgG+iH3rebh71Zc3MLJIWfAl8vLyYBJlFaTXti2QOyRSKSZIRKQqpnIjBk1aBeLFT19Do5Z+8gRENm7+loVHnoywfM8t9Ugu/OgRkapYV5HMfzx9vDB42mh4Nmogd3gE4Ojm/fjpyA+QJAlajcSNbEkWTJCISHUqW28noHUgRs43oJ6HqzxBkY1vv0hFUUGh3GGQijFBIiLVMe/1ZU6UzH83CmyEfs8Plikqsnb5h3Sc3nPM8m8Tv5NVJKpbTJCISLXMFSTrSlKXx3qgx+Nh8gREFqLUhF1LN+NWdi6Asn8jDrVRXWKCRESqZK4iWc9HAsr+7j4kDJ738PZ/uRnzCrF5/ioAFSfXE9U2JkhEpFrlO13z994+XtA/O7DuA6IKfj5yDme/PWGp8nGojeoKEyQiUi3rKpL5b/PXLe8PhpuXh0yRkbWtH65DQX4hq0hUp5ggEZGqmTdEFULYdMDuOg808L9HpqjI2q3rOchMv2pJYDkXieoCEyQiUj2TqLithauHG3oN18sUEZX3v2PnKtx1SFSbmCARkeqZh9rKC36wA5qHtqnjaKgyP3x93HJHG8C5SFT7mCAREcF2/pGZi7sb+k8Yinu7dkD9hp6Qyp9Adeb3ixlYM+MjXD6bXum/FVFNk4QQLFZWQpIksGmI1MV6botJ3O6ECwuKkPv7TeTn3sKORV/h2s+/QpSaZIpS3TQuGgx/41m0fbADgNtzyKjmqb0fZAWJiMgOY0lZ5+Dq7gbfID+06NAa0S8/iVZ/a1uj76Pz9YKLW70afU1nZSox4fPXF+M/6/cCAGbuWCBzROSsmCAREf0hQW+wWTRSq7ldSTJXk5q0CsTBlK9xT1DjGnnPe4Ia4+m3X0DEC0M4hHenBJC6eCNSPliLnFzu10a1w0XuAIiIlEgjAbC6sy1Bb0D8ziTLhO7HZ47BmvglyMnIuqvX923eBD1HPIqWXUPh5eWBJq0CIQDsXrYFBTfzauAKnN/Rzd8i/eSPuJF+GZ+9slDucMjJcA6SHWofeyVSM/NcJHMlSSNVPtfltc3zsPffW3Fo/b5qvX7DAF+MfX8yPH28KjyWmX4Vy+IWoOhWQbXjVrPZs2dj2rRpcofhVNTeD3KIjYioHOtkqKoO4u2BU9BvfAyiJo9Avfpud/TaklbCwMlPVJocAYBfywAMnDwCHg09qxe0yk2fPh2jRo2SOwxyIhxiIyKywzzM9md3Sv09ogfadu+EvZ9tw3+TD9g9T+fjhfBx0Wjarjlyfs+Gu84Drh4VE6sOPTujeac2uPbzrziwdhfSj6f91UtRhZUrVyIkJASvvvqq3KGQE+AQmx1qLy0SUdlQm0nYX0jS+jyzrKvXsfezrTi15yggAI1WC11DT/gG+WHg5CdQmF+IIxu/QVFhEdzc3dBtcG80aRVoeX6pSUBrNVk76+p1JI1MrPmLc2J79+7FI488IncYDk/t/SATJDvU/sEgorLE507W2alsztKFE2k4tH4fJA3g4aVD98GPwMffF5vmr8bvlzMhaTQQJhMaBTVB9CtPVlpJKv3jBQ9/tRepSzbV3IU5uTZt2uDEiRPQ6XRyh+LQ1N4PMkGyQ+0fDCKqnsoWmTQWFKEwr8AylPZUywj0CHsYHl71oXXRorSkFAW5+Xj8jWfg16zJ7dfSGypsyPrr+V+w8tUPUZxfVGfX5Mg+++wzxMbGyh2GQ1N7P8hJ2kRENcS6ggSUbVXi6dsQrh5uVpUoAZg7HSEgICBBsqy/ZF2xsl6TSZIkuLnf2URwAuLj45GXx+US6O4xQSIiqgEJekOF/cHM35uTnsDAQPTp3hvGomIU5hfCWFSMJq0C0aBRAwAV5zpZL1Dp3dgH9dxda/UanEl6ejpOnjwpdxjkwJggERHVAnP1xzrp0el0iI+PR2DbZvAN8kNg22boHfsYXD3cKiRH5ec+efp4IWxMFNw8PeokfkcnhMCFCxfkDoMcGOcg2aH2sVciujvWc4fsTfDOy8vDzM3vWuYm2Tuv/LwmAPj9l6tYOvl9GPO4xcafGTZsGL788ku5w3BYau8HWUEiIqoFVd39ptPp4N2o/Nykyl/DXIkyD7f5tQxA5ISh0NbT1kbYTmXt2rWIi4tDZmam3KGQA2IFyQ61Z85EpBzm9ZiAsiTJJIDj2/6D5Pe/kDcwB9G0aVMsWLAAQ4YMqfRx6z326Da194OsIBERKZz1BHBzohTapwsCO7SQLygHcvXqVbz88suVVpISUpPK9tort6wCERMkIiIHYX2XnIu7Gwa+NAJNmST9KZPJhMzMTJw/f77iY8J2OQUiMyZIREQOwiTK7s4yJ0qNWwbi/HdnsGDBAmi1nJNUlYKCAqSnp9scM1ePgLLkM34nq0h0GxMkIiIHYB5m02okS6KUGG6ATqfD2LFj0b17d7lDVDyDwYDPP/8cgO2QmnWSRGTGBImIyMGYEyUznU6Hjz76CAEBATJGpXzZ2dl48sknEdy9E25l5QK4PbxWahJlK5lzLhL9gQkSEZGDsL7tv/zyAKGhoUhLS8O7774Ld3d3mSJ0DD8dPoP3nk7A4c37LVUjSWL5iGwxQSIicjD2bknX6XSIi4vD77//jr7joqF1c6njyByHyViC7R+sQ/IHa1FSWGSzrQvnIhHAdZDsUvv6D0Tk+Cb++59YGbcAuddvyh2KojUPbQ39s9EIbNcSQFmiVNUCnmqh9n6QCZIdav9gEJFzeDx+DNbN/TeEsVTuUBTNu7EPhr0+CkEhLS3H1J4kqb0f5BAbEZETa/9wZ/zj7RfkDkPxcn7Lwp6VKTAWFAEAJ2wTEyQiImeWGG5A89B70frBELlDUbz079Pw8/dpAHjLPzFBIiJShZhXnoJUjz/yqyJMAoe/2mupIgGcsK1m/N9CROTkEsMN8PTxwtBpI1HP3VXucBTtl1M/4dSeo5bvWUlSLyZIREQq0aFnZzw9bwK0LtyWxC4T8M3qnZaFJAHORVIrJkhERCqQoDeguFSgadsW6Dasj9zhKFpudi4unvkfjAVF3MRWxZggERGphFZTdtt22MgBeGDgw3KHo1iipBR7lm3B9kVfIePCFQCsIqkREyQiIpVIDDdAkiRoNRL0YweizQPt5Q5Jsa5f/g1nD3yPfcuTUWg1aZvUgwkSEZEKuXq44cEoVpGqUnirED9+dwbHkr8FwCqS2jBBIiJSEfM+biYBtOjcFj4BjWSOSOEEsGdZMn787gcATJLUhAkSEZHKlJoENBLg6u6G8BcGyx2O4plKTdjy3v/h8tl0TtpWESZIREQqM6vfSwDK1vhp27UjWndpJ3NEyldqNOLQhn0wFhZx8UiVYIJERKRCCXoDTALIz85Fzm/ZgJYrIlalILcA/zv+I84fPqPqDVzVhAkSEZFKaaSyTVpNJaXo0bW73OEoXsHNPGx65zOc2nOUc5FUgAkSEZFKJegN8GzsA42LFum/X4ZvkB8kLbuFqogSE1I/2Ygcq5W2yTnxfwIRkYrpGniiz8j+KC4yQphM8Gnii4ef1MsdlqLlZ93CseRvWUVycpLgYGqlJEniODMRqUJCahJuZeUiOzML3n4+8GzohdN7jmDD3FVyh6ZoD0b3xHcbv5E7jFqj9n6QFSQiIoKnjxeC2jWHt48XNBJwX98H8cDAh+QOS9GObNqPR8dFIzMzU+5QqBYwQSIiUrkEvaHCMZMAHhzQU4ZoHMvuTzajffv2WL9+fZXncWkAx8MEiYiILEyibCFJAPBrGYD7+3WVOSLly8rKgsFgqLKSpJG4CrejYYJERESWdZE0UtncE7NBLz+FzhE9ZIzMMVy5cgUbNmyo9DEmRo6JCRIREQEoS46svzZvq/Ht2lRERkbKE5QDSUxMxOHDh22OlU+OmCw5DiZIREQEoOJcJHPCpNPpMHPmTPj6+soQleMoLCzEBx98gLy8PAAVkyHu4+ZYmCAREZFFgt5gSYxM4nbS1K1bN2zduhUuLi4yRqdsN2/loKSkBNnZ2VVOyuaEbcfABImIiCqVGG5bUeratSt+/fVXeHl5yRSRspUaS3Dq1Cm8s3uJzXClmfmYEIJDbQ6AC0XaofYFsoiI7Nm3bx/69u0Lk8kkdyiKVE/nhr4jI/G3ft3h6uFmOW4eYjMnSpUtr6Akau8HmSDZofYPBhFRVZYuXYpnnnlG7jAUrWXnYDz2/GA0aRWIUpOAJEk2lSWTqFilUxK194McYiMiomobO3Ystm7dKncYipb+fRq+eGMpbmXl2iydAMCypAIpFytIdqg9cyYiuhPbtm1D//795Q5D0dwb1Mfjr49Bq/uDbSpJ5iRJqUNtau8HWUEiIqK7FhERgffee0/uMBSt8GY+Vk1fhP9uOwitpuIwGykTEyQiIvpLJk2ahF6x/eQOQ9FMJSZsWfAFTuw5Wvb9H/OPuAWJcnGIzQ61lxaJiKpr1LsGrHxlodxhKJskoddT4fj6s+1yR/Kn1N4PsoJEREQ1onnovXhoRF+5w1A2IfDtmlQsX77c5jAXj1QeJkhERFQjEsMNeHTMQHQb2kfuUBTNZDJhzJgxmDp1qs3xGdsXyBQRVYZrxhMRUY0xCUD/TDRKS0pwdON+ucNRtHnz5qFevXqo90gTaCRA8LZ/RWEFiYiIakxiuAGSJKF7TB+4urrKHY7ivT3nbWRdvQ6gbM4PJ2wrBxMkIiKqUYnhBvgG+GL27NnQarVyh6NsJuDCiR8BcOFIpWGCREREtSIuLg4HDx6UOwzFS/10My6fS7d8zyqSMjBBIiKiGpegNyB+ZxJSsg+iz6gIucNRtMKcfHyRsBQXTqTJHQpZYYJERES1wryNxpYP1sLdu77c4SjarRs5WPfWCksliVUk+TFBIiKiWpGgNyAhNQnz/rMUES8OkTscxcvPvoVdSzfDWFAkdygEJkhERFSLzHuN3dfnAbTt0VHeYBzApTPpSPvuBwCsIsmNCRIREdWaxHADSk0CJgHETBuJex8MgaYe72yzx1RSiuQPvsTlc+ncyFZmTJCIiKhWSVLZDvbuHm4IGzsQLUPbQHJh92NP4c187FmWgpLCIm5BIiN+QomIqFYlhhss1ZAmrQIx/J/PYOQ8AzwaesobmIL9+uNFZP+WzbWRZMQEiYiIap25o9dIgIu7G5p1aIUhr42SNSYlM+YX4tS+YzAJzkWSCxMkIiKqdQl6A4CySdsaqexPm87BiJgwVObIlOvbNTvw320HeVebTLhZLRER1QlzcmSt68CeyMnOxYFVO+QJSslMQMr7/4f/bj2ArGvX8dG4OXJHpCqsIBERUZ1IDL9dRbLWa1hfNA9tAzc3NxmiUr6rP17C//3zU7y2eZ7coagKEyQiIqozCXqDpYpkTpRcPdyQvHoD+vXrJ19gCpd97TqWvbIAZ86ckTsU1WCCREREdcokKg63hYaGYs2aNRg+fLh8gSlcRtqv6NSpE8aMGYO8vDy5w3F6khCCS1FVQpIksGmIiGqHeX0f835tZnl5eYiIiMDBgwdRUlIiV3iK17NnT3z44YcIDQ2ttfdQez/IChIREcnGOjkCAJ1Oh/nz56NFixYyReQY9u/fj+HDhyMzM1PuUJwWK0h2qD1zJiKS0+HDhzF58mQcPHhQ7lAUrWPHjvj8889rpZKk9n6QCZIdav9gEBHJLS8vDxs2bMDTI58GTHJHo1xPPPEEPvnkE+h0uhp9XbX3gxxiIyIiRdLpdIiNjcXHiz6WOxRFO3HiBLKzs+UOw+kwQSIiIkWLiYmBt29DucNQrPT0dBiNRrnDcDpMkIiISNH8/PywbPGncG9QX+5QFEkIgevXr8sdhtNhgkRERIo3ZMgQjJ5ngNaVO2SVV1BQgKNHj8odhtNhgkRERA5h0bOz8eJzL8gdhiJNmjQJq1evljsMp8IEiYiIHMakSZPg7u4udxiKU1RUhJdffrnKdZHMi3PSnWGCREREDqNVq1aIj4+HJEl/frLKZGRk4OOPP650GxJzcsQk6c4xQSIiIocyffp0rF+/Xu4wFGnWrFmIjY3FqVOnLMfidyah1CSgkWz3v6OqMUEiIiKHc1x3CQMmcWPb8kpKSnDs2DG8//77yMvLQ/zOJAghoNVIMP2x5mNCKqtId4IJEhEROZzEcAM69ukCv3sD5Q5FcTIzM5GRkYGpX8yGRipbEZuVo+pjgkRERA7J1d0Ng6c8jZb33wtNPd7+b1ZUVIQjR45gz/JkXDybDgCW6pFJlP3hXKQ/xwSJiIgcUmK4AY1aBCDmtVFo160j6vt4yh1SBa4NGsjyvpmZmTh/6BTWv7UCl8+lyxKDo2OCREREDk3XwBO9n34M/q2DAIXd3Wa8eVO29xYlJuRkZmHfyq0wFhZVeJxVpKoxQSIiIoc1q99LkCQJTVoFInpqLFr8rS28/bhvm7XfL2Ug9/fbiZoQQsZoHAcTJCIicmgaqWxejWdDL/QbNwjNO7aG1lUrd1iKUZCbD1gV1qzXkGIVyT4mSERE5NAS9AbL1wGtA9ErNgJN2gTJGJGylBQa5cAdwAAAFfVJREFUUVpSCiGEpXpk/loIgZk7FsgcoTIxQSIiIoeXGH47SfJq1AD3BDSCi0c9GSNSloyff7X5niuR/zkmSERE5BTMQ22u7m7oNqgXvH04F8ns+I7D+PX8LzbHJEmy/GEVqSImSERE5BTMQ20aCQhq3xIxr8YisF1zmaNShl9OpuHz1xcj7cgPNseth93IFhMkIiJyGuYqEgAEhbTE8MRxaHHfvfIGpQQmgcLcfHwRvwSHt3xT4WFJkjhhuxwmSERE5DQS9AYIISxJkrePF3oMfQQaLbs7AIAAUj/8Ct+s2Q5jYZFliI0q4ieGiIicirbcxmOt7m+LpiEtZIpGeYQQ+Prf2/H564srrLLNKtJtTJCIiMipJOgNNpuzunq44dFnB/GuNmtC4OKpn/HlG0tx4USa3NEoEhMkIiJyOtabswJA85CWGDpjDLTc1NZG3o0cbJz7b1w+l26ZsM0qUhlZEqSFCxeiZ8+eiI6ORnR0NL7++msAwIEDBzB48GBERUVh8ODBOHjwYKXPz87OxujRoxEeHo7Ro0fj5h973Qgh8Oabb0Kv1yMqKgpnzpyxPGfFihWIiYnB1q1ba/8CiYhIVtbrIpmTpOAHO2Bg3Ag0DPCFu3d9mSJTnoLcAhz+ah+Ki4xyh6IotZ4gHT58GNOmTatwfNSoUdi0aRM2bdqE3r17AwB8fHzw0UcfYcuWLZgzZw6mTp1a6WsuWbIEPXr0wM6dO9GjRw8sWbIEAPDNN98gPT0dO3fuxKxZs5CQkAAAyMvLw6lTp7B27Vps2bKldi6UiIgURyPBMtymkYD7+j6I0e9PxlNvPY8xC+PgE9hY3gAVoNRYjKwrv6PwVoFlwjbXRVLYEFuHDh3QpEkTAEBwcDCMRiOMxooZ7e7duzFo0CAAwKBBg7Br1y6b45IkoXPnzsjJyUFmZqZljQfO1CciUo8EvcHyByirJJlE2Z1tgW2bo1m7FoiaPFzmKJXh98uZKCkusfSXs/q9JHNE8pMtQVq9ejWioqIwffp0yxCZtR07diAkJASurq4VHrt+/Tr8/PwAAH5+frhx4wYAICMjA/7+/pbz/P39kZGRAU9PT7Rt2xZDhgxB//79a+mKiIhIqczDbOZ1ksxVpVb3B+Ott96SLzCFcKvvhsLcfAC2w5NqVmuz1YYNGwaj0Yj8/HzcvHkT0dHRAIBXXnkFTzzxBF544QVIkoQFCxZgzpw5mD17tuW5aWlpmD9/PpYtW1at96xsNVBz1Wj8+PEYP378X7giIiJyRAmpSZbEyDo5AsqqTHn/Lw9btmzBoUOH5AtSZqUlJng1bsjKkZVaqyCtXbsWmzZtwptvvomwsDDLfKOePXuiUaNG0Gq10Gg0GDZsGE6dOmV53rVr1zBhwgTMnTsXzZtXvkS8r68vMjMzAQCZmZm45557AJRVjK5du2bzWuZKExERqZN5iA2wXWnbfFyn02HJkiXwC24qR3iKYCou4TSUcmQZYjMnNwCwa9cuBAcHAwBycnIwbtw4xMXFoUuXLnafHxYWho0bNwIANm7ciL59+9ocF0Lg+++/h5eXFxMkIiKyoZFskyYACA0Nxf+On8e93TvKFJW8SktKEBPQS+4wFEWWBGnevHmIiopCVFQUDh06hOnTpwMAVq1ahYsXL2LRokWWJQCuX78OAJgxY4al0jRu3DgcOHAA4eHhOHDgAMaNGwcA6N27N5o1awa9Xo+ZM2fin//8pxyXR0RECmM9Ubt8cmSm0+nw/a7DaNzav9LHnZm7qzu8vb3lDkNRJMFtfCslSRJ3OCYiciLxO5PuaAJy9JQnsXn+53UQkTJIkoSHHnoI27dvh06nszmu5n6QCZIdav9gEBGpVWZmJtp0aItb1yveYe2Mmjdvji+//BLdunWzOa72flBR6yARERHJzc/PD33HRkHr6vzbkoSEhGDt2rUVkiNigkRERFTBxrmfYXj8WEha5+wmQ0JCkJycjCNHjqBr165yh6NIzvkvT0RE9BcFd+2AMe+9hAb+DeUOpUYFBARg+fLliIyMtJlzRLY4B8kOtY+9EhFR2SKTxoIibHxvDc7u+17ucP4yn6aNsP2rlDuqGqm9H2QFiYiIyA6TAFw93PD4a6MxZOZo+LZognoebpA0jrWootbNBUNnjsal8+kcUrtDrCDZofbMmYiIyiSkJgEoS5ZKCotw8/ebKDYacXzbf3B08wGZo7sDEhA5aQSS/1W9pQvU3g+ygkRERFQF64UlXdzd4Bvkh8A2QWjRqY2MUd25vuMGVTs5IiZIREREf8q8ya2xRFg2u215f1to3ZS7FIC7pwci40Zg18cb5A7FITFBIiIi+hPmFbjd65VlRyYBePp4YdCUpxXXk/o2bYyIicPw/NIZeOCxHnKH47CUm/oSEREpkLmaBACdenVGo+bTsPSlf6GkwChvYBIQM/1ptO0WCld3t0o35aU7p7C8l4iISJnMyYZGKkuSzPxbBiDmlaeg8/GSKbIyLULvRXDXTkyOaggrSERERNVgXUEy69CzM4I6tcGJ1MPYuzwFotRUpzE169gKES8Ogau7GwAmRzWBt/nbofbbG4mIqHLxO5MqJEjWDm74GqlLNtZZktS4VQBi57wIz4ZeEEJgVr+XauR11d4PcoiNiIioGszJkUnc/mP+HgC6PNYdwV071HgPK2kk1C83jOfhrcOASSPg2bDseE0lR8QKkl1qz5yJiMg+68UjzayrShHe3fD444/j4sWLNfae9X08MWDSCJzZewy3sm7BTeeGniP0CGrfEiZx+067mqL2fpBzkIiIiKrJPA+psqE28/yftWvXImbk47hy7pcaec+Ae5uh5f3BaNO5LfJvFaC+pwfeHjgFQNmwH9UsVpDsUHvmTEREVatsLlL5ydF5eXkY8upI7Phw/V96r0Yt/DF4+kgEtA60HLOeLF4bk7LV3g9yDhIREdFdSAw32AyxVZak6HQ6dB3YC2GjI4G73N82oG0zjHxnAgJaB1Y6pMc71moHh9iIiIj+oqqSlMRwAxIkoNMjXXBg7W4c33kYJmNJla+ncdHinqaNERDcDD2GhsHTx8umYlRqEpAkqcbnHdFtHGKzQ+2lRSIiujMJqUl3VMVJSE2CSQD52bk4s/97/HjwFCBJKCoohIuLC1zdXVFaUop2D4UiuEsHaOtpUa++O9w83G3unKurypHa+0EmSHao/YNBREQ1y3oitUYCjAVFKMwrgKvOAwBQeKsAb0W/jHn/WYpSk4BWI1kqRdZznWrjjrXKqL0f5BwkIiKiOlA+qXH1cIN3o4Zw93CD+x9f63Q6m8qQViPZbG1SV8kRMUEiIiKqMxoJEEJUWFwSsE18JEmyWYjSXEFiclR3mCARERHVkQS9wVIVAm5vfFt+PpE5EbJea8mk3tEuWfAuNiIiIhmYE547rQqxelS3WEEiIiKqQ+ZqkUaqOukxr7NUWYWJah8TJCIiIhncSdLzZ0kU1R7e5m+HJN3lkqdEREROQs0pAucg2aHmDwUREZHacYiNiIiIqBwmSERERETlMEEiIiIiKocJEhEREVE5TJCIiIiIymGCRERERFSO0yZIV69exdNPP42IiAhERkZi5cqVAIBt27YhMjIS7du3x6lTp2yec+7cOQwfPhyRkZGIiopCUVFRhdfNzs7G6NGjER4ejtGjR+PmzZsAypYFePPNN6HX6xEVFYUzZ85YnrNixQrExMRg69attXjF8qhuO1++fBn33XcfoqOjER0djfj4+Epfd+HChejZs6flvK+//try2OLFi6HX69GvXz/s37/fcjwlJQUxMTFYsWJF7VysTGqrjflZtnU3PzMA4MqVK/jb3/6GpUuXVvq606ZNQ1hYmOXf4+zZswDYzjXdzpcuXcKwYcMQHh6OSZMmwWg0AgCMRiMmTZoEvV6PYcOG4fLly5bnzJ07F4MHD8Z3331XC1dKDk84qYyMDHH69GkhhBC5ubkiPDxcpKWliZ9++kn8/PPPIjY2Vpw8edJyfnFxsRgwYIA4e/asEEKIGzduiJKSkgqvO3fuXLF48WIhhBCLFy8W77zzjhBCiH379omxY8cKk8kkjh8/LoYOHSqEEOLWrVsiLi5OFBcXi+eee65Wr1kO1W3nS5cuicjIyD993aSkJPHpp59WOJ6WliaioqJEUVGRuHjxoujbt6/l3+n5558XJSUlYtKkSeLWrVs1dIXyq6025mfZVnXb2WzChAli4sSJlX5ehRDi1VdfFdu2batwnO1cs+1sMBhEcnKyEEKImTNnitWrVwshhFi1apWYOXOmEEKI5ORk8dJLLwkhhPjpp5/EnDlzRH5+vjAYDDV+neT4nLaC5Ofnh44dOwIAPD090bp1a2RkZKBNmzZo3bp1hfMPHDiAdu3aoX379gAAHx8faLXaCuft3r0bgwYNAgAMGjQIu3btsjkuSRI6d+6MnJwcZGZmWhacdNaVuavbzn/V7t27ERkZCVdXVzRr1gwtWrTAyZMnAcCmrYUTLfRZW23Mz7Ktu2nnXbt2ISgoCMHBwdV+P7ZzzbWzEAKHDh1Cv379AAAxMTHYvXs3AGDPnj2IiYkBAPTr1w8HDx6EEAImkwkajcbpfl5QzXHaBMna5cuXcfbsWdx///12z7lw4QIkScLYsWMRExODTz75pNLzrl+/Dj8/PwBl/9Fv3LgBAMjIyIC/v7/lPH9/f2RkZMDT0xNt27bFkCFD0L9//xq8KuW5k3Y2nzdo0CDExsbi6NGjds9bvXo1oqKiMH36dMvwT/l2btKkCTIyMgAA4eHhGDJkCDp16gRPT88auCLlqck25mfZvjtp5/z8fHzyySeYMGHCn77ee++9h6ioKLz99tuWoR+2c821c1ZWFry9veHiUrY5hLktgbJ2DggIAAC4uLjAy8sLWVlZCA4ORmFhIZ588kk88cQTNXhV5CycfquRvLw8GAwGvPbaa1V2mqWlpTh27BjWrVsHDw8PjBo1Cp06dUKPHj3u6H0q+w3E/Bvg+PHjMX78+Lu7AAdxp+3s5+eHvXv3wsfHB6dPn8aLL76IlJSUCs954okn8MILL0CSJCxYsABz5szB7Nmzq2znmJgYy2+Kzqim29gefpbvrJ0XLlyIkSNHQqfTVfl6cXFxaNy4MYqLizFz5kwsWbIEEyZMYDvXcDuXZ27Lqtp55syZ1XpNUhenTpCKi4thMBgQFRWF8PDwKs/19/dH165dcc899wAAevXqhTNnzlRIkHx9fZGZmQk/Pz9kZmZazvf398e1a9cs5127ds3y27mzq047u7q6wtXVFQDQqVMnNG/eHBcuXEBoaKjNeY0aNbJ8PWzYMDz33HMAKrZzRkaGKtq5NtqYn+WKqtPOJ06cwI4dOzB//nzk5ORAo9HAzc0NsbGxNueZ287V1RWDBw/GsmXLALCda7KdfXx8kJOTg5KSEri4uNi0pb+/P65evQp/f3+UlJQgNzcXDRs2rNXrI+fgtENsQgjMmDEDrVu3xujRo//0/Icffhjnz59HQUEBSkpKcOTIEdx7770VzgsLC8PGjRsBABs3bkTfvn1tjgsh8P3338PLy0sVP+yq2843btxAaWkpgLK7TtLT09GsWbMK52VmZlq+3rVrl2XuQVhYGFJSUmA0Gi3Pv++++2roapSpttqYn2Vb1W3nNWvWYM+ePdizZw9GjhyJ8ePHV0iOgNufZSFEhc8y27lm2lmSJHTr1g07duwAAGzYsAFhYWEAytp5w4YNAIAdO3age/fuTju/i2qWJJx0dtrRo0fx1FNPoW3bttBoyvLAuLg4GI1GzJo1Czdu3IC3tzdCQkIst41u2rQJS5YsgSRJ6NWrF6ZOnQoAmDFjBkaMGIHQ0FBkZWVh0qRJuHr1KgICArBgwQI0bNgQQggkJiZi//798PDwwNtvv13hN3ZnVN123rFjB5KSkqDVaqHVajFx4kTLDzLrdp4yZQrOnTsHAGjatCkSExMtncdHH32E9evXQ6vV4rXXXkPv3r3lufg6UlttzM+yrbv5mWG2cOFC1K9fH2PHjgUAPPvss3jzzTfRpEkT/OMf/0BWVhaEEGjfvj3eeOMN6HQ6tnMNt/OlS5cwefJk3Lx5EyEhIZg/fz5cXV1RVFSEKVOm4OzZs2jQoAHee++9Sn9hICrPaRMkIiIiorvltENsRERERHeLCRIRERFROUyQiIiIiMphgkRERERUDhMkIiIionKYIBFRtYWEhCA6OhoDBgzAc889h5ycHJvHV6xYgdDQUOTm5tp9jczMTMtq0YcPH7Z8vXnzZkRFRSEqKgojRoywLPcAAN988w369esHvV6PJUuWWI5Xdyf38+fPY9q0aTXTGETklJggEVG1ubu7Y9OmTUhOTkaDBg2wevVqm8eTk5MRGhqK1NRUu6+xfPlyDBs2rMLxoKAgrFq1Clu2bMHzzz9v2Q6itLQUiYmJ+PTTT5GSkoLk5GT89NNPAID58+dj1KhR2LlzJ7y9vbFu3ToAwNq1a+Ht7Y3U1FSMGjUK8+fPBwC0a9cO165dw5UrV2qkPYjI+TBBIqK/pHPnzpaNQQHg4sWLyM/Px6RJk5CSkmL3eTt37kSvXr0qHP/73/+OBg0aWF7bvB3HyZMn0aJFCzRr1gyurq6IjIzE7t2772ondwDo06dPlfERkboxQSKiu1ZaWoqDBw9aVuoGyqpHkZGReOCBB3DhwgVcv369wvMuXbqEBg0aWPaMs2fdunWWJCojIwP+/v6Wx5o0aYKMjIy72skdKNun7tixY3/h6onImTFBIqJqKywsRHR0NLp164abN2/ioYcesjy2detWREZGQqPRQK/XY/v27RWe/9tvv8HHx6fK9zh06BDWrVuHV155BUDVu7LbO17Vc8yb9RIRVYYJEhFVm3kO0t69e1FcXGyZg3Tu3Dmkp6djzJgxlo2Fk5OTK32+eSJ1Zc6dO4fXX38dixYtsiRS/v7+luE2oKw65OfnZ7OTO4BKd3IHUGEn96KiIri5udVAaxCRM2KCRER3zcvLC6+//jqWLVuG4uJipKSkYOLEiZbd17/99ltkZGTg119/tXley5YtKxwzu3LlCiZOnIh33nkHrVq1shwPDQ1Feno6Ll26BKPRiJSUFISFhd31Tu7p6ekIDg6u8TYhIufABImI/pIOHTqgffv2SElJQUpKCh599FGbx/V6fYXJ0PXr10ezZs3wyy+/ACiby2Sej/Thhx8iOzsbb7zxBqKjozF48GAAZXOI4uPj8cwzz6B///6IiIiwJDhTpkzB8uXLodfrkZ2dbbk7bujQocjOzoZer8fy5cstw3VA2dICjzzySK20CRE5PklUNkhPRFTLUlNTcfr0aUyePBkrV65ERkYGpk6dWifvbTQaERsbizVr1lgmdxMRWWOCRESyWbt2LY4fP460tDS8//77aNq0aZ28b3p6OjIyMtCtW7c6eT8icjxMkIiIiIjK4RwkIiIionKYIBERERGVwwSJiIiIqBwmSERERETlMEEiIiIiKuf/A5l35UEe2bXJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xidplus.plot_map([prior_MIPS])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[215.08957607 214.97755187 214.86830064 ... 215.03025421 214.31734531\n",
      " 213.82151274] [53.00814137 52.94286013 52.7669635  ... 53.12972201 52.42888554\n",
      " 52.20276907]\n"
     ]
    }
   ],
   "source": [
    "print(prior_MIPS.sra,prior_MIPS.sdec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate tiles\n",
    "As fitting the whole map would be too computationally expensive, I split based on HEALPix pixels. For MIPS, the optimum order is 11. So that I don't have to read the master prior based on the whole map into memory each time (which requires a lot more memory) I also create another layer of HEALPix pixels based at the lower order of 6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 705 tiles required for input catalogue and 4 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,prior_MIPS.sra,prior_MIPS.sdec,unique=True)\n",
    "order_large=6\n",
    "tiles_large=moc_routines.get_HEALPix_pixels(order_large,prior_MIPS.sra,prior_MIPS.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='./data/'\n",
    "outfile=output_folder+'Master_prior.pkl'\n",
    "with open(outfile, 'wb') as f:\n",
    "    pickle.dump({'priors':[prior_MIPS],'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": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prior_MIPS.nsrc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (herschelhelp_internal)",
   "language": "python",
   "name": "helpint"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}