{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Selection Functions and Number Counts\n",
    "\n",
    "Using the depths maps for ELAIS-N1 you can calculate the probability that a source of true flux f will be detected in each healpix\n",
    "\n",
    "In the field with an associated error calculated in the depth map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy.table import Table\n",
    "from astropy import units as u\n",
    "from astropy.modeling import models, fitting\n",
    "from astropy.modeling.models import custom_model\n",
    "from astropy.coordinates import SkyCoord, search_around_sky\n",
    "from IPython.display import clear_output\n",
    "import scipy\n",
    "from scipy.optimize import curve_fit\n",
    "import scipy.stats\n",
    "import pickle\n",
    "import os\n",
    "from pymoc.util import catalog\n",
    "from pymoc import MOC\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as patches\n",
    "from matplotlib import gridspec\n",
    "#import utils\n",
    "\n",
    "from herschelhelp_internal.utils import flux_to_mag, mag_to_flux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "DMU_LOC = '../../../'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_center(bins):\n",
    "    \"\"\"\n",
    "    Get the central positions for an array defining bins\n",
    "    \"\"\"\n",
    "    return (bins[:-1] + bins[1:]) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaus_cdf(val, mean,sig):\n",
    "    \n",
    "    return(0.5*(1 + scipy.special.erf((val - mean)/(np.sqrt(2)*sig))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaus_prob(errors, confidence, hist_errors):\n",
    "    '''\n",
    "    Returns the probability of a source of a given true flux being detected to a confidence level\n",
    "    The fluxes used for this are in the range 0,98\n",
    "    This is done assuming gaussian errors\n",
    "    \n",
    "    Parameters\n",
    "    -----------\n",
    "    Errors: a list of errors found in the field\n",
    "    confidence: a integer, the confidence level you are working at eg 2/3/4 sigma\n",
    "    hist_errors: a list of the number of regions of your field that have an error given in errors\n",
    "    \n",
    "    Returns\n",
    "    ---------\n",
    "    Prob: the probability that a source of given flux will be detected in the field averaged across\n",
    "    all the regions in the field\n",
    "    '''\n",
    "    \n",
    "    prob = np.zeros(len(hist_errors))\n",
    "    cutoff = np.zeros(len(hist_errors))\n",
    "    true_flux = np.arange(0,len(hist_errors),0.1)\n",
    "    \n",
    "    cutoffs = confidence * errors\n",
    "    #print(errors)\n",
    "    #print(hist_errors)\n",
    "    #print(cutoffs)\n",
    "\n",
    "    for n in range(len(errors)-1):\n",
    "        prob = prob + (\n",
    "            1 - scipy.stats.norm(np.array(true_flux[:(len(cutoffs)-1)]),errors[n]).cdf(cutoffs[n])\n",
    "        )*hist_errors[n]\n",
    "    return(prob/sum(hist_errors))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load in the depth map and data for Lockman-SWIRE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "depth_elais = Table.read(DMU_LOC + 'dmu1/dmu1_ml_ELAIS-N1/data/depths_elais-n1_20180216.fits')\n",
    "SERVS_elais = Table.read(DMU_LOC + 'dmu1/dmu1_ml_ELAIS-N1/data_tmp/SERVS.fits')\n",
    "SWIRE_elais = Table.read(DMU_LOC + 'dmu1/dmu1_ml_ELAIS-N1/data_tmp/SWIRE.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "servs_coords = SkyCoord(SERVS_elais['servs_ra'],SERVS_elais['servs_dec'])\n",
    "servs_elais_moc = catalog.catalog_to_moc(servs_coords,10,11)\n",
    "elais_servs_area = servs_elais_moc.area\n",
    "\n",
    "swire_coords = SkyCoord(SWIRE_elais['swire_ra'],SWIRE_elais['swire_dec'])\n",
    "swire_elais_moc = catalog.catalog_to_moc(swire_coords,10,11)\n",
    "elais_swire_area = swire_elais_moc.area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/rs548/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/matplotlib/font_manager.py:1297: UserWarning: findfont: Font family ['serif'] not found. Falling back to DejaVu Sans\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR4AAADdCAYAAABtwh14AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXdcVfX/x58fLhsERYaKgoLiwg0oqDjR3LMcYVaO0mz9srKhOdL6VpY5sqxsqJnl1pypuRdq5hYXbpaK7HH5/P44eAVFGQKX8Xk+Hudx7/2czznnfS6X13l/xvv9EVJKFAqFoigxMbYBCoWi7KGER6FQFDlKeBQKRZGjhEehUBQ5SngUCkWRo4RHoVAUOUp4FApFkVOihUcI4S+E+CdjOyuE+MrYNikUipwRpWUCoRDiZ+AnKeV2Y9uiUCgeT6kQHiGEGXAU8JZSphvbHoVC8XhMjW0AgBBiDPA80ABYLKV8PtM+B+BHoBMQBbwnpfztgVMEAVtyIzqOjo6yevXqBWO4QqEwcOjQoSgppVNu6hYL4QGuAx8DnQGrB/bNAVIAF6Ax8JcQ4qiU8kSmOk8DP+XmQtWrVyckJOTJLVYoFFkQQoTltm6xEB4p5XIAIYQPUPVeuRDCBuiH1oSKA3YJIVYDQ4BxGXXMAF9g2KPOL4QYCYwEcHNze3KDU1Ph8mVITNQ2V1eoXBmEePJzKxRlgGIhPI/BC9BLKc9mKjsKtMn0uSOw9XHNLCnlPGAegI+PT947te7ehStXoE4d0OngmWdg5cqsdRwcoEED8PbWXhs2BH//PF9KoSgLFHfhsQViHiiLAcrd+yClXA+sL1QrfvsNRo2C69c1z8bHB6pVAysrMDODyEi4dEnzgg4c0LwgJyeIiNCOnzYNbG3htde0z+npYFKiZzIoFE9EcReeOMDugTI7ILbQr3zlChw9ClJqTah33oFdu8DSUvNmHoWUmuDcuQNr1mhlK1aAjQ3UqKF9HjNGE6aWLSEgQNuqVSv0W1IoigvFXXjOAqZCiFpSytCMskbAiccc8+RcvQqBgRAXB7NnQ5Uq2pYbhAAXF227x4QJmiAB6PXg6wtnz8J338HMmVp5tWrQvj0MGgQdOoBpcf/TKBT5p1j8uoUQpmi26ACdEMISSJNSxgshlgOThRDD0Ua1egEBhWrQqFFw7Rp8+ilYWxfMOe91POt0MHSo9j4tTWuinT4Np07B0qXwyy+aN/Tdd9CnD8THax5U9eqq81pRaigWwgN8CHyU6XMwMAmYCIwG5gMRQDQw6oGh9IIlNRV27gQ/P6hVq9AuA2heTc2a2ta9u3btQ4dg924IDdWaaiEhMHky7NgBrVvDv/9qZU2bQv36YGFRuDYqFIVAqZi5nBd8fHzkY+fxrFsH3brBBx9A8+ZFZ9ijiIrShCYwUPO+Fi/WNtA6tmvX1vqOatSAdu2gSxclRgqjIIQ4JKX0yVVdJTwPMHgwrF0LP/+s/WMXN9LT4eZNuHABzp/XOsEjIiA8XBtNK18e+vaFH37Qmma3boGdneozUhQ6eREe9WvMTGysNj+nbdviKTqgDcPf6+xu1ep+uV6vjcJt366Nvq1dCz16aEJ665Y2zA/w11/alAA3N6hYUfUbKYyCEp7MrFiheQ1t2uRct7ih02n9Pk2b3i9bswaaNNH6jtas0byl557TJkSCNjXAyUnzklxctPvu3BmaNVPzjBSFimpqZaZTJ61T9+uvS6cnICXcuKGNpEVFadvdu/dHzi5e1Oq9/DLMnat5gP/3f9oQf/v2kJSkeVUNG2qTJxWKTKimVn4ZNEgTnNIoOqDd1+PmJN25o42aOTpqHlJsLCxfrolMfLw25P/uu5o3VLeu5k01aHB/c3Utvd+dokDJlccjhNiRy/MlSSk7PZlJhcsjPZ57s4wVjyYuDo4d0zq2723R0ff3V6igCdCMGZoohYZq0wD69dOac3FxmjDZ2BjvHhSFRmF4PL7AyzldF/g6l+crfpw8Cc7O2tNekT22tlrga+bg17g4rekWFnY/Xm3nTm329+bNMGsWdOyoCc/s2fDee9r8o7Ztwd1dmzx55ozWfHv55ceHoyhKDbn1eLZIKTvkot6mEunxSKl1rnp6wvvvG8ew0khqKty+rY2e6XRamMiRI5rInzql9RmVL6+Nsp0/DykpWnP3t4w8b4MGaZMrp0zRPqvg2mJNgXs8uRGdjHrFWnQey/jx2g9bUXCYmWle5D28vLQNtHCRpCTNiwKtP2nrVi1U5V6zNypKm3+0Zo0mYqNGaU25e8G1qanw33+al+Xjo3lWKrtkieCxHo8QQocWG5UOrJFS6jPKn5ZS/lk0JhYsqo+nhBIXBwsXap5SWFjWh4S1NSQkaO+//x6GD9eEbPduLczExkabJmFmpondkSNaE696dWjUKKs4KvJNQXo8vwJhaKlH3xZCDJVSngNGASVSeLJl9mzt6dqsmbEtUTwKW1utDwg0kQkNBXNzrZ/Iykqbwf3ff1qzec0abSLl9OnwzTdQtSqsWgW//qp5Wg8+bN3cIChIC5Xp00cr27JFa9a1a1e091lGyMnj+UdK2TbjvTuaEE0Exksp2xeFgQXNQx5PaqrWoezvD6+8YjzDFAVLeLg2Y9vDQ4tdO3VKm71tYQGenhxITWX57t04JiXheucODW7coLapKWa3b8OaNZz96CMuX7uG5bJlVKhQAddJk7BPTETUravFx3l5aa9OTpoQHj6s9Wf17GnsOzcaBenxmAshLKSUyVLKMCFEd+B3tNUgSgcHDmiT6Bo3NrYlijxy49Ytftm6lQ2HD2NraYmzvT2D27ShY+PGXNXpmHf4MP8tX86liAhMdTpMTUyYO2oUTTw9ObZpE1/u3UtqWprhfDoTEyIXL6aCrS0/u7ryyZEjWlMtA2shiFi/Hhu9niXAEaCyhQUXkpPZDkRaWHD1jz8QPXqwpEULbnl7Y+bvT73atWkmBBaNGmkjf+fPE29nh02bNlqnexkkJ+F5HSgPhANIKWOFED2BQYVtWJGxebM2t0QN4xY7xs6fz5WoKBKSk0lITiYxJQX/2rWZPmwYd+LiqDFiBMmpqTT19ORuQgL/XrhEs5r16NgY4pOSmPrnn3i4VMXKvDouDgmkpKVyOTKSJp6eDOvUib7+XfhuQ0XCY27QxGMXZ65dJiU1FYARPfvQtVUrEpOTuRUXx9WoKCJiYrAZMoToi3fZ+Ov3LPxvD6nJyVjpdARUc6Nvs44IIWDBAj45epSjBw7A/PkAWABvAVOBNMANqKHTEVijBkmpqdy8e5d2/fvz6rx5ACTHxrL74EHOnTtHJ39/qje4/6xPT0rCxNKySP8WBc1jhUdKeTCbMj2wsNAsKkpSUrQ+gLp1oVy5nOsrCozQ6zaEXrfBxGQLM1YtJS4pCSd7e+pVq8bkwUOYsdqDlfuuY2Z6HWsLC6zMzbG1tMQh4+9U3taWGcNeQWcaQOi1Jmw75siNO/ZM+C0Vj0pH6Nwkne3TNtB9cltiEsywMLvDtql7sLXScyncis9XePLL1mrEJ2n/Ai3rBrJ0XAjxyTr6f1qXdSEujB9wlrF9zmNmKom4Y86v26rS+oNK7DntgCCYT547wdD2RwE7XvrGl4+WVSY88SJfDTvB1vnzSUlL42qUZO6aSP46EMrMFHuadPGmi/813vhrHQsPX2PGuUvYCDPczNN50U7L8ntm2jQafziepEzrFwRaW7PN35+bx2/RNhwsLcN4tlVTGrVsybn//sOyRg2GT59O2KkE3uvej9qtKlOrc0dcHR0pf/sOF2Lq0HlAXazti0fwc55itYQQzlLKiEK0p9DJ0sfz7bfaEO1HH6mO5SJASsmR8+H8vKUZ36xvjD79c2ActpYueFWpSlJqFMmp6TiX383e01oO6q4+4Sz8vyPcjjPjs+We/P2vE0Jov9mL4dbo000wNUnHyzWOBu6xHAgtz8VwG17ocJk/dlfBxkJPx8aR/L7DlcYeMTTxiOHnLdUQAtp4R9PdN5wbtyyZ/Vd1bCz1xMSbIgR4VYnnWJgdjWvEUN8tliW7qpCmN6GGSzx+te5wOcqKvacd6OF3kxOXyxEWYUUzzxgOhFagVb1oXn4qjN93VmHDYWfS9CZ4VIpHAOdv2vBix8uE37HgrxAX7KxSuJtoTk+/m7zZ6wIh5+z5bfXfHIk+BwRhbloDnVyJh+kMZlUMYFDEz4Sn7caEiaRzxPDdeli7Ub1mCFv/cwICgZ0PfPtP4c8E1lV7iaGxVzgTH8/ttHQcLF2xd3meZNuXeef9ynSqdYiNP81nT3pt/liWSPdOXnwzux2W5ctrf8OM0USRMZ8qKgoWLdL67MeOLaR8PEKIrSW1U/keBuFJTtYyDFpbw2efqRijJyAhWcf1aAuquyRiqsv6ezpx2ZbxCyuw78xGbsUtJDn1OLCUjo0C8XL9j/UhBwiLfJ10aY3OJB0zU4mZLp3hnS6TmKLjx81uVLBJJSrWHBMh8al5BwvTdNKlwKV8Mg2q36VetTgszLR/iORUwcYjzvywyZ3GNWJ4o+cFHMqlEnLOnukrPUnTC4IaR9LP/wYV7VINdl6KsGLGag+qOycwpO1VKtqlsudUBb7d4E5yqo72DaPo6hNONcckQBsYW7m/Er9urUplh2Re7XaRutXi2H7cgW/W1SAxRYdDuRTaekfRrkE07s6JpOkFi3e4snRPZSzN0glue5WuPuGsOVCJRdtdSU7V+nsqV0iiZd1bdG4aQarehEmLvYhNNEWfLrC3TuOVbhdZvMOV01dv413pAGciW5Cqd6VS+WTaN4rCwTaVHzY5YmF2HjPdVW7FJhLofIvdkcOoa3URe9GGm/rqXEiqjZ79wCnsLF/jbtLXdKi6mi1Xe2X5G5pjxp/1vfBJcuD5C+5slUuxNnHG3MSJW2kWSGKBL4GgQhOebVLKEj2+aBCeuXNh9GiYNEmLK1Lkiz2nKjDoi6ZcjrTG3FSPl2s8Hi4JuDklcjUqhZX7nwMOAhJby8ZUcRjAM6388POyNZzjboIpZ67ZcvaaDXfizegXcINKFZIBOH3Vlnkb3fF2v0uv5jepWC41e0Me4Mw1G9ycErEyv99ciYwxx1QnqWCbu3MApOkF6RLMTbP/P7kWbYm9TSq2lnpD2c3bFkTdNadutVh02Uy0Pn/TGlvLNFzKpxjKrt+yICzCmtqucTg8cI+340z5dGktzM3SGdv7PPY2aaSmCX78243/LtrRrGYMLetG4+Uaj0nG8/NqtCWfL/ckJdWEt/qcp2blBA6ft+d/y2qSmKJDZ5JOYP1bPNXkJhVi14KLF3tDffhtuzPJaZfo0iSFAYG3mLkklMM3j9DOqhV7k/uRmn4ee9MvSNHfIE1G4qBLoqZFPPXbrGPu+nrK43kUPj4+MmT3bm0qvpubFm2tvJ1cExJqz7kbNpiZSo5etGPanzVxtEuhT4ubnL95mVNXthKXlEaqfiLp6QKHcr1p4eVKmwaBVHeuaWzzyxTpGf/aJpl+3qev2rD/bAU6N4mgUoWUh465eduCuCQdNStrEzL16TBzjQf/HKtI2wbRDAq8ZngoPEjPj5urtBiPZd06LTufm5sSnUeQpheEXrehmmMitlZ6EpJ1jPulDrPWemSp51NzL1Udv2fNwa1ci9aWzm7i0ZyJg3ogEZiI8cYwX0FWwblHnarx1Kka/8hjHhQVnQm80fMCw4IuY2ed9oij8k5ehad0/JeammpJ0RVZiLhjzsJ/qrLhsDN7z1QgLtEUU106PjXvEHXXnHM3bOnhdxO/WkcpZ+WCuakJG4/MZvX+JTRwb0p336fxqxWIk722pljp+LEohKBARQfyLjwDC/TqimLB0Yt2TFlSi1X7K5GmN8HdOYHA+tHUrBTP9VuWHL9cjqTUW/Tw/ZzzN1ax5sBRpjw7h0Y1fOjT4ln6+gdTwbaisW9DUYLIk/BIKcMftU8IMU5K+emTm6QoKiJjzBm/qDbfb3THxlJPN59wghpH4eaUaKiTkBxHXNIcNv+7iqtReqo51iC47SiqVnQDwKGcyl+kyDsF2ccTCCjhKSHsOO5A3098iUkwpbtvOANbX8PWSv9QPZ2JKccuhdC5aR86N+lNdeea2uxcheIJKDDhkVJ2LahzKQqX+Zur8fLchjjbJzNp8JksHk5yahJb/vuLf45t4OPg2ViYWfL1yIWYm6pFAhUFR9kc1SpDJKeacDysHEcu2HP4vD2Hzttz4GwFghpH8kLHy4b5J6HXT/H30TXsOLGJ+KRYvFzrcys2mkoVqijRURQ4eRYeIcTkR+xKBq4CGx7XF6QoGlLTBC983dgw1R/A2iIND5cEgtteoV/ADcPktjPXjvP2T8MwN7XAv047OjfpTX23xqpJpSg08uPxeAF9gAPAFaAa4AesAXoA3wgh+kkpNxSYlYo8oddD8JdN+GOXK119wvF2u4tnpQRcKiQb5nak6dM4ffUUdao2oLarN+/0nUpjj+bYWqpgWUXhkx/hMQEGSilX3CsQQvQCBkspWwghhqJ1MivhMQLp6TBidiP+2OXKCx0u08f/5kN1Tl89xjfrPuVa9GW+e2UpjnYutKrX0QjWKsoq+RGezjycj2ctsCDj/UJg9pMYpcgfxy6V49V53mw/7siEgWfwqRmTZX/4ness3vED2/5bh0M5R8b2mULFcirfsKLoyY/wnEfLuZxZXF7OKAdwBB49J1tR4KSnw7hf6vLlKg+sLdJ5petFmnlmFZ3YxBjGfDeI9PR0ercYzIDWw7C2UAvrKYxDfoRnOLBcCPEucA1wBfRA34z9tYEiCdARQlRHC30+kVH0tJQysiiuXZwY8mUTfttRlaDGEQxtfzXL9PaL4Wep4eJFOSt7Rnd5F2/3ZoaQBoXCWORZeKSUh4UQtYAWQBXgBrBXSpmasX8HkNsljwuC7VLK/kV4vWLFqn0u/LajKu0aRDGm2yVDzGtM/G3mbZzOzpOb+eS5b6nv1oR2DdVUK0XxIL/zeNqi9fM4Sym7CyF8hBB2UsqtBWdarmkphNiJlnLtA5mXPB8lmNQ0wfbjFQn+qim1qsTxSreLCKFl+dtxYhPfb/qShKQ4BrcZiZert7HNVSiykOf1YIUQrwJzgbPAvRT8icDH+TVCCDFGCBEihEgWQvz8wD4HIcQKIUS8ECJMCDE40+4bQE20cA1n7jf3Si1Rd83pM82HCoOfImiCP+a6dN7vH2pIVDXtz3eYvnIClSq4MmPEAga2HoaZrnjk2VUo7pEfj+cNoIOU8lJGPw/AabS+nfxyHU24OgNWD+ybg7agoAvQGPhLCHFUSnlCSpmMNnERIcRytObfsiewo9jz/IzGrD/kzFNNI2hY/S4Na9wlJTUcKSsihMC3ViuaefoT1KQXOpOyuXSKoviTZ48HKIc2cRDgXrPGDE0c8oWUcrmUciUQnblcCGED9ENbQDBOSrkLWA0MydifebZba+Bcfm0oCZy7bs2Gw050ahLBy13CCKh7m38vbGL0t8+w4bA2rapTk1481ayvEh1FsSY/wrMDGPdA2WvAtic35yG8AL2U8mymsqNA/Yz3rYQQhzL6eFyB37I7iRBiZEZTLiQysuQOen24qA6mOsnA1tdJTk3im3X/47Pl71PVsTpNPFoY2zyFItfkp6n1KrBGCDECKCeEOAPcRQuXKGhsgZgHymLQvC6klOuB9TmdREo5D5gHWs7lAraxSAgJtWfJTleeaXWNhORQJi7+gEsR5+jjH8yQtqMw1al4X0XJIT/D6TeEEL6AL+CO1uw6IGWm1ccKjjjA7oEyOyC2EK5VbDl91ZYRsxvhaJdMX/8bnLp6g9tx0UwY+CU+NVsa2zyFIs/k6zGZMWR9IGMrTM4CpkKIWlLK0IyyRtyfMFiqSUoxYcoSLz5f4YmZLolezRdjbVGPZp7+zBuzHCtza2ObqFDki1wJz2NSYWRBSjkhP0YIIUwzbNEBOiGEJZAmpYzPGK2aLIQYjjaq1QsIyM91ShrPTm/C8r1VaFn3P8LvDGHJzhO08V5CFYdqSnQUJZrcejzVMr23RBtpOgiEoa0/78eTDWN/CHyU6XMwMAmYCIwG5gMRaKNeo6SUpd7jiYk3Zf0hFxq4r+T45RdJTklibJ/JVHGolvPBCkUxJ1fCI6V84d57IcTvwCAp5bJMZX2Bp/NrhJRyIprIZLfvFtA7v+cuqXy3wZ3ElDkcv/wmrg5uTBsyl2qONYxtlkJRIOSnj6cL8OwDZauAn57cHAVo6UpnrPGgcoUY3J0DeaPnBBVJrihV5Ed4zgGvADMzlY3mfloMxRMQcecOM9cIbtyyZOKg/jT2CMJE5Ge6lUJRfMlvWowVQoh3uJ8WI40yECdV2Kw9eJBhM2cSHWtJA/f2NPGIRSjRUZRC8jOP58jj0mIo8o5er2fsTz8xY/VqnO3roE9fSqcmCWpZd0WpJb/zeFLR0lAonpDk1FQGfPYZq/bv59VuPVix71eqOZrTqt5xY5umUBQayo83Muampjja2TFr5Eia157I1ejyDAq8Zlh6RqEojagAHyNx9to1hBDUqlKF78eMQZ9uQr1XvKjunEBA3dvGNk+hKFTUc9UI7D55Ev933mHojBlIKRFCsPAfV0Kv2zK4zVXD2lcKRWlFCU8Rs+nIETpOmEDFcuVY8OabCCG4FG7Fa/MaUKtKHM297hjbRIWi0MmX8AghnitoQ8oC60JC6Pnxx9R2dWX3//6HZ+XKJKWY0P9/PujT4a3e59VIlqJM8Ng+HiFEveyKgZeAXwvFolKKlJJPly7F292dTZMm4VBOS574+vfeHDpXnvefPksVh2QjW6lQFA05dS7vA5aiiU1m3AvHnNLJvX6c1R9+CEB5W1sA1h9yZt5Gd8b1C6VFbdXEUpQdchKeU8DbUsoHcyH/VXgmlS42Hj7M7L/+Ysk77xgE5x7jfqmDY7lk/FS/jqKMkVMfTxDw0H+FlLJb4ZhTuthy9Ci9p03jalQUSSlZc+EfvWjHf5fs6eYbjqmuRGZjVSjyzWM9Hinl3cyfhRDOUsqIwjWpdPDPsWP0mDKFWpUrs3nKFEOfzj2+Xl0DCzM9nZsUXPJ5E3Oo6GuGmZ14uHGsUBQEElLvSqIPppKe73Vl8j6B8Hegff4vVzbYdfIk3SZPpoaLC39PmYKjXda00eG3zVm03ZWOjaKwtdIX2HUr+ppR2cMJO+vyCDU8pigEpJTcjb8DRBK5O//hmXkdTle/5lxga2lJcy8vtnz8Mc7lyz+0/9sN1UlJ09HdL7xAr2tmJ5ToKAoVIQR2NuU1r/oJyKvHozojHsPdhATsrK1p7OHB1qlTs62Tphd8u8Gdrj7hVK2YVLAGCJToKAodIZ68Ka9mLhcQ0Xfv0nzsWD76Lds1BQ1sOuLEzduWNK7x4HJhCkXZIa8ej3qcZkNCcjI9Pv6Yi+HhdGjY8LF1F2yrSjmrVJrVLH3Cc+t2NL0GdAcgIjIcnYmOihUdAdiydjvm5uYFcp24+Dhee/sVTp85hURS3r48yxatwtzMHBcPB+rVqW+o+3SfAbw26g2e6hNEdHQUFhYWWJhb8PXnc/Cu583IV4fROiCQIYOGGo5ZtXYFvy9bzOKf/uCzrz5hxZplmJjo0Ol0zPjfLJo2blYg95Fbduz+BytLa3yb+RXpdQuTvArPwEKxogSTptcz6PPP2XfmDH+++y6B3t6PrBsTb8rK/ZVo3zAKs1I4hO5QoSI7N+0F4NPpU7GxseXVl1/PUkdKiZQSE5P8O9tzv59N1SrVmP/NLwCcPXcGM1MzAGxtyxlseJD5c3+hQf2G/LLoJyZOG8/ShSvo1/tpvv1hThbhWb56Kf169WfP/t1s3bGV7Rv2YG5uTlR0JGlpBTcYkFt27N5ORQfHUiU8efrrSynDhRBfCiEaF5ZBJY1Xv/uO1QcOMHPECPoFPH65r2V7KpOUoqNdg6gisq54cOHiefw7+PLmuNdo81RLrl6/ins9V8P+Zav+5LWxrwCapzRk+CDadW1Nh25tOHjo4TUjb0bcpEqlyobPXjVrY2Zmlmt7fJv5cePmdQDaB3bgxKnjREZps0Ti4uPYtXcnXTp1IzziJo4OFQ2emmNFJyq5VHrofOcuhNLzma60CmpBm6dacvlKGOnp6XwwaRz+HXwJ6ODHqrUrAPhn5zaeHXb/+f3muNdYsux3AOr7ePHp9KkEdg6gZcfmnLsQysVLF1iw+BdmzZ1B607+7A/Zx7JVf+LfwZdWQS3o8XSXXN93cSI/+XjMgI1CiEhgAbBISnm1YM0qObRt0IDKDg6M6d49x7q/bquKl2scXlXii8AysO3/VI51Ujt2ITnDK7Ht/xQpzwST8kww4lYUNiODs9SNW7oh37acOXuaOdO/5atPZ5KWlvbIeuMmvM1ro97Et5kfl6+EMeD5/uzdcjBLnSEDh9I/uDfL1yynTas2DOr/LB41PDUb42Jp3cnfUHfsa+/Qq3ufLMdv+edvunXW/l5mZmZ07dydlWtXMOL5l/hrw1raBrbHxtqGDm2D+OLrz/ANbEKb1m3p27M/Ac0fXjJ6+Csv8O7/vU+XoK4kJSWRLtNZuXY5Z86eZtemfURFR9G+WyABLXJebtrJyZkdG/fw3Y/fMGfeLL76dCZDBg2looMjo4Zr4vza2NGs+XM9zk4uxMSUzFnv+cm5/KoQ4g3uL3PzoRBiP1rQ6HIpZVwB21gsSU9Px8TEhAGtW+eqfliEFduPOzLl2dNlMgK9hrtHrvpG/tn5D6HnQw2fY+7cITExESsrK0NZ44ZNOLLnONu2b+GfXdto1y2QLWv/obpbjcc2tV4cNZSEhHiklPyzYZehvH+vp5n6xRRGPP8Sy1cvZeiz2jJyduXs2L5hN3v372bnnh08/1Iwkz+cxsD+gwzH3rlzm+hb0XQJ6gqApaUlAPsO7KV/n2fQ6XS4OLvQws+fI0eP5NjP1aNLTwAaNWzC5q2bsq3T3MefUW+MpFe3Pob6JY385lzWA2uBtUKI+sBvwM/ANxkL/n0kpbxWYFYWM5JTU2n7/vuM6NSJF4OCcqwvJUz7sxYAwW2vcizMLocjCoa8eiiZ60sHxyfycB7E2vr+kssmJiZIeb+PKzn5flS+lDJO3/uvAAAgAElEQVRXHdHlbMvRs1tvenbrjZSSv7duYvjzLz32mPlzf6F2rTp8NPVD3h0/lp++XQBAQItWXL5ymeMnj3Pk6CEWfH9/ZNLU1JTWLdvQumUbanvVYcXqZVmEB7KfwpD5/jJjqtORnp6e7b0DWJhbAKAz0ZGmz94z/Prz2YQcPsjGLetp1cmf3Zv3Ub58hcfee3Ejv/l47IQQw4QQ24AdwH6gNVAXiAPWF5yJxY/xCxey78yZbCcHZseUJbWYt9Gd3i1uFJnoFGdMTEwob1+e8xfOkZ6eztoNawz72rZuyw+/zDN8Pnbiv4eO33tgj6GJkZyczNnQM1Sr6para5ubmzN+3ET27N/NuQuhBnt6d+/DqDdG8FTHrgbROxN6mgsX7y8Xd/zk8YeuU758BSo6VGT95nUAJCUlkZCYQECLlixbtRS9Xk9EZDj7D+6jSaMmVKvqxumzp0hJSeHOndvs2L09R5ttbcoRFxdr+Hwp7CK+zfz44O0JlLcvz/WM/qqSRJ49HiHEUqAzmuB8C6yUUiZn2v9/QOkbK87g1JUrfLV6NcOCguju65tj/ZlravDRb3Xo0DCS5ztcKQILSwYT359M/yF9qFqlKrW96pCSrAX+fD71S9567w0WLVmAXp9Gq4BAvpj6VZZjL1w8z1vvvQGAlOl07tiVrp27o9frH+rj6dShM+PfnZjleGsra0YNH8Ps72Yy43+zAOjX62m++X42U8ZPM9SLj49n3IS3uRsbg4mJjlqetQz1MzNv1o+8+e5rfPy/SZibm/PLvEX06taHkMMHadWpBQLB1Amf4OToDEC3zt1p2bE5nh41adQg53Garp278fxLQ1izfjWfT/2SGXO+JOzyJSSS9oEdskwfKCmIR7mEjzxAiLHAQinlzcfUsZZSJjypcYWBj4+PDPnoo3wdK6XkqYkT2X/2LGfnzs3R49n2X0XafxhAnxY3eK79lUJfOaJKV3M83WoV7kUUCuD85VCur8saJdrz4+aHpJQ+uTk+1x6PEOLFjLe3gK7ZtGslEA0ckVKWykf7vxcusOnIEb4eMSJH0YlL1DHoi6ZUrpDE4DZquRqFIjN5aWoNyUUdO6COEOIdKeWcfNpUbGni6cmhr76igXvOCRjfX1CHiDsWTHvuFBZm6TnWVyjKErkWHillu9zUyxjl2gCUKuG5FB5OdRcXmnp6Ep+kI+quKZUfkSN5x3EHZq31oIfvTeq7lYnZBQpFnshzA0AI8dhxTinlCWBRvi0qhuw8cYKaL73Esj17AHj7p3rUfaUdV6Mss9STEn76uxo9PvbDo1I8Q9qV2XmVCsVjyU/PQ5wQ4j8hxAIhxFtCiI5CiJpCiJ/uVZBSjitAGx+JEMJFCLFHCLFdCLFVCFE556PyRlxiIs9//TXuzs50btIEgJX7KhETb8bIOQ251zcfGWNOr6m+vDizMdUcE3mvfyiW5qqJpVBkR34mEDoDjTO2RsAYwA24UIB25ZYooJWUMl0I8TwwDPi4IC/w3q+/cjE8nO3TpmFrZcX1aAtu3LakunMC6w+58OvWqjT1jKHjeH9ux5kxLCiMHn7hajVQheIx5Cdk4g7wT8YGgBDiY+BGgVmVe1syhwqXA04UxHmTU03YecKByLt3+Hb9Rp5r15XW9bW5EjtOVATg1e4XWbGvEq//4I0+XWCuS+fToaeoVURxWMURY6fFmDh1PDU9azHyhZcB6DWgGx41PPnq05mAFgdWw70GQe07M/SlYHZu2ss/O7fx3Mhnca/mTnJyEl0792Di+5MB+PW3n5n86UdUrlTFcO35c3+hlqdXgdxHbpkzbxbDnhthCMcoDRTUIO/HQL6bV0KIMUKIECFEshDi5wf2OQghVggh4oUQYUKIwQ/sb5wRKzYGOJxfGzIzdn49gib4M/gLHWnpNuw9M9Gwb/vxiliZ6/GoFE9w26skpZhQqXwy0188UaZFB+6nxdi5aS8vBA9j1Igxhs/3REdKmSVkID/cS4uxZ8sB9m45yNefzcbM1Aw/n+YcCNkPgF6vJ+ZuDCdO3X8WHTi0n+a+/g+dr3VAoCZCG3azZt0qQg7fD0p9us8Awz3s3LS3yEUH4Jt5s0hOLuBslUYmP53Lc4QQI4QQvkKIexJcGXiSCYPX0cRrfjb75gApgAtaUOrcjJEzAKSU/0opmwPjgfeewAYDiSna1zJ5sBvdff/l7DUfzt/QYo12nKhIvWqx6EygikMy8145yv+eP0lFu/wnvi7tFFVajOa+/hw4pAnPiVPHaVCvIVaWltyNvUtiYiLnL57Hu16DR9ppbWWNdz1vQ8qM3LJoyQJadmxOq6AWvPJ/mrcVdvkSPZ7uQsuOzekzqAfXrmuhiyNfHcZfmUJEqnq5AFq6jF4DujFk+CB8A5vw8usjAPjm+9lERkfSpW8neg/sTlpaGi+9NpyADn74d/Dlux+/yZOtxYX89PFcQ1tp4nXAQwhxDa3fZ6MQog/aIoChDzSDHouUcjmAEMIHqHqvXAhhA/QDvDOi3ncJIVajzSkaJ4Qwl1Lemz4ZQy7F78JNaxbvuP/D9699i/aNoh+odRvPynFUc0xkXYg78/92442eFzh5pRzPtbs/P7KC7aNTPBiT9z4y59iJgp212KB+Op9Myt+aJkWRFqOaazX0ej03bt5gf8g+fJv54ejoRMjhg1iYm9PIuxGmpo/+yd++fYtLly/Rwu++V/TniiXs2rvT8PnBJuOxk8f4+psv2bhyCxUqOHD79i0A3nr/TYYMep5n+g7g54XzeX/iO/wy7/GDvUePH2Xf1hCcHJ0I6tmOg4cOMHrEGOZ8N5P1yzdhb1+ekMMHib4VzZ4tmiCXpbQYhmAWIYQZUA9okLGNyHh1AgqiQeoF6KWUZzOVHQXaZLxvLIT4AtADScCLZIMQYiQwEsDNzY2F21z5aHEdw343pwTCftzywFGD+HRZDFODv6apZww/balGo4w8yd7usSjyRlGkxajpUYvmPs05ELKPAyH7+b9Xx+Lo6MSBkH2YW1jg59M822vu3LODlh2bE3r+LG+//q4hpgq0ptYnkz57pL07d2+nT49+VKjgAGB4PXQkhCW/LAVgYP/BTPtiSo737tPE15BorEG9hly+GvZQ1kGPGh6cuxDKuAlvE9S+M+3bdMjxvMWRfKXFuIeUMhVNCI5mLhdC5C5sO2dseTjgNAatIxkp5QEgMKeTSCnnAfNAi9Vq5HGX5e9pT9EvVnpyOcIqS/27CZeAjdSuoqVZCGocySdLazFxcW2szPV4Vi7+fTn59UwKi6JIi1HToxZ+zZpz4NB+zoSepnatOjg5OvHjL/MwMzNn2HMjsj1X64BAFv34O2fPnaFr30507dyd+nUfncI2M1LKPK3sYWpqaujj0uv1WVJfWFhYGN6b6EyyTbPqUKEiuzbv4+9tm/hu/lzWrFvJjM9m5/r6xYVc+eJCiJzlWqs3CQwjXwVBHFoYRmbsgCdyOUwEmOokpjqJiXg4SPb8Ta0N3q5hbwB8a93B2T6ZU1fKEVD3VqnMl1yUFGZajOa+/qzb9BdOjk6YmJjgWNGJyKhIDh05iF8OOYu9atbmtVFvMnPujFzfS5tW7Vi+eqmhiXXv1aepLyvWLAPgj+W/GzIXulV1499jRwBYu341en3OPRK2trbExmsz4KOiI5FS0rt7X9576wOOHj+aw9HFk9x6PG8IIeaT8yoTrwH5C/3OnrOAqRCilpTynv/diAIaNn8U4XcOAN5UsHUC9JjqJM+1v8IXK2oSWP/BviBFfiiMtBgADeo3JDIyggF97+c1ruNVl9TUFOztc3bEhw0dgU/rRly5pvXjPdjH89WnM/Fpej8dinc9b14b9SZd+3fGVGdK4wZNmDX9Gz7/eDqvjh3NV7On4+zkzOzp3wLwfPAwnn1xAP/s2Ea7Nu0Nib8ex9BnX6TPwB64VnFl0gcf8+rY0QZPa+L7ufIJih25SoshhEhHiz7PSXiSpJTWOdTJ7vymaCL4EVrn8gggTUqZlpHRUALD0SYtrgMCMkIz8oyPj4/8qMv9kYDPlntyKdyaKz/9DWirRlj3DyZV/xy/jR2GraX2RLpxy4KPfqvN+8+EUt05MT+XLnRUWgxFUfGkaTFy1dSSUppIKXUZr4/b8iw6GXwIJKLNBQrOeP9hxr7RgBUQASwGRuVXdHKDlJLW9b/gwX7qyg7JzBvzX7EVHYWiJPFEncsFhZRyIjDxEftuAb2LyhYzU1NquHQB3IFDRXVZhaJMUSyEpzix9uBBomN1aMKjUCgKAyU8mZBS8tKcOViZ/wd0NbY5CkWpRSXkROu5Bi3Z1/Vbt3Apr/WPqQBzhaJwyEvO5W3c/x/NDimlLJnTKDPYdeoUAC7lc149QqFQ5J+8eDwL0TILPrj9AzQEHg77LQFk9mp2nTyJvY0NFWyLPgK5tFDRzY7Wnfzx7+DL8y8Fk5CY/9jhXXt2MGBofwDWbfqLr2ZPf2TdmJg7WSYe5pZPp09l1rdfP1Qeev4s3fs/RetO/jRv25Q33hkDQGDnAMOkxrS0NFxrORvWPgdo26UVR4/9y29/LOTtD/7PcI16zWrRupM/Ldo1Y+nKPwz1R7/5Eo3869O6kz+tO/nTqVfRP7vz+909CbkWHinlj5k3YCXaAn5vAcvR4qpKNHtPn6Zl3boIoVqg+cXK0oqdm/ayd8tBzMzM+WnBj1n25zctRtdO3XhzzFuP3B9zN4Yff/0+z+d9FOMmvM3ojLQe+/85zIiMHD9+Ps3ZH7IPgOMnj1HTsxYHDmmf4xPiuXT5UrYR8PdShCz6cQn/N+51UlPvZzOY/MHHhrQbm1Y9GDNY+BT0d5cb8pMWwy4jhOIcWqqKplLKkVLKEp9geO/nnzN31Chjm1Fq8G8ewIVL57l8JYzmbZvy1vtvGNJibN2+hU4929PmqZY8/1IwcRkhAX9v24xfmyY81SeINetXG86V2YOIiAwneNhAWgW1oFVQC/aH7GPiJxO4dOkirTv5M37KBwDMnDuD9t0CadmxOZ98cT8x5RczP8M3sAm9B3Yn9EIo2XEz/CZVKt/PYHAvdqu5j78h58+BkH28EDyMYyeOAXD4SAiNvBuh0+ke+Z14etTEysqKOzG3c/096vV6xk95n4AOfrTs2Jx58+cCsH3XNgI7BxDQwY8xb40yxLw1bFGP6FtRABw5epju/Z8CNM9rzFuj6N7/KRoHeBtSajz43d0Mv0nXfp0Mnuue/btzbWtuyUsfjxXwBpqH8w9aytFCDV0oamwsLbEpRVne7v3gMtO7R1+GDx1JQmICzwzp+9D+wc8EM/iZYKJvRTF0ZHCWfWvzsJZ6Wloaf2/bTIe2HQEIPR/K7OnfMn3aDKJvRfHFzP+x4vc12FjbMGPOl3wzbxavjXqT198Zw+olf+FRw5MXRz2X7bnHTXiblv6tWfjj79rqofFxTHxvMqfOnGTnpr0AbN2+hQsXz7Fl7XaklAx64Rl279uFjbUNy1ctZfvG3aSlpdH2qVY0btDkoWuMHjGGngO64desOe3bdODZZ4Kxty9Pc5/mTP1cy1C4/9B+3n3zPZat+pPYuFj2H9pPc58Wj/1ejh77F48anlki4CdM/ZAvZmoR8HW86vL97KxpqX5eNJ+wy5fYsXEPpqam3L59i6SkJEa/+TKrlqylpkctXn59BPMX/MCo4a889vpnz51lzR/riIuPxTewKS8+N+Kh7272dzNp36YjY197B71e/0TN5UeRl+H0i4AO+AwIAVyEEC6ZK0gptxagbUXKir17OXD2LB8HB+dcWfFIEpMSDUsI+/sFMGTgUG6G36BaVTdDioeDhw9y5uxpnuqtiVJqagq+TbW0FO5u7nh61ATg6b4D+WXRTw9dY8fu7cydoTUNdDod9nb2D+Wl2bZjC1t3bCWwcwCgLUd84eJ54uLj6P5UD6yttEn2XTplP23i2QFDaN+mI1v+2cy6TX/x88If2blpH27V3ElNSSE8IpzQc2ep5elFk0bNOHQkhAMh+w1pVx9k7vez+fW3n7l0+SJLF67Ism/yBx/Tq3ufR36n23du44Uhww25hCpUcODYyWO4u7lT00MLkRn09GB++HlejsLTqUNnLCwssLCwwMnRkYjIiIfqNGnUlFffGk1aairdnupBg/oNH3vO/JAX4UlCG9V6VFtEAh5PbJGRWLZnD/8cP84nQ4ca25QC43EeirWV9WP3V3RwzJOHc497fTwPXS9TWgwpJW0D2/PjnJ+z1Dl24j9EAU1ikFLy5pi3eCF4WJbyuT/MyXUai8qVKhM88DmCBz6HfwdfTp05SeOGTfBt5seqv1bg4lwJIQS+TX3Zf3Avh/8Nwbdp9iOio0aM4dWXX2fNulWMen0kh3cfy3UO5WxTbzwmxtJUdz/1RtIDKVMzB6Wa6HTo9Q8nZWvZohV/LdvIpi0beOm14bw26g0G9h/8UL0nIS+dy9WllDUes5VY0QE4cuECTTxK9C2UGLR/1H1cuHgegITEBM5dCKWWpxdhV8K4eElbsGTZqj+zPT6wVVvmL/gB0Po/7sbexdbWlri4+4sntm/TkUW/LzD0HV2/cZ3IqAgCmrdk7YY1JCYmEhsXy4bN67O9xt/bNhs6gMMjwrl9+5Yh6XtzX3/m/jDHkGbDt1lzfl+2GGcnlxwj4Ht07UXjRk1ZvDT3S8+1a9OBnxb8YMjcePv2LWp5enH5ymXDd7hk2e+0bNEKALdqbvz7378ArFm3KsfzP/jdXb56GSdHJ4Y++wJDBg3l6LF/c21rblHDN0C6TOD0tWsG4ZFSTR0sTBwrOvHNV98yfMwLtOzYnKAe7Th77iyWlpbM+N8sBgztx1N9gqjm6pbt8Z9O+oyde3YQ0MGPtl1acfrsKRwqVKS5Twv8O/gyfsoHtG/Tgf69n6ZTz/YEdPDTOrDj4mjUoDF9evYjsHMAQ0c+i3/zgGyvsW3HFgI6+NEqqAX9nu3FpA+m4uKs9Sw092nBpbCL+DbTMhpWcqmEXq9/ZIbDB3nnjXF8M2+2wSuZMPVDw3B6607+pKRkjfp+btDzVHWtRqsgLa/znyv/wNLSkjlfzuX5l4cQ0MEPExMTXhgyXDv/m+/x3kfv0KVv0GM7uu/x4He3e+9OWnf2J7BzAKvXreTlYaNzdV95IVdpMUoTD6bF+Hy5J6evHifybluWv/ceffz9GTazEfP/duP3t0Owtig5i/KptBiKoqJI0mKUdvTp4djb2NDU09PYpigUZQIVJApYmffgyk+/GdsMhaLMoDyeDIQQeUraXSyRUNaazoqiR0r5+KjNXFDmhSddphIR04olO3fmXLmYk3pXcjf+jhIfRaEhpfYbS737ZL+xMt/USkwOJVV/GH16O2Ob8sREH0wFIomyi1I5PRSFg9QecNpvLf+UeeGJS9IijUvDHJ70FIjcrZZSVhR/ynxTKy7pGAJrvKpUMbYpCkWZocwLT0zCHsxMm+ZqopVCoSgYSrTwCCHGCCFChBDJQoif83q8Pl2PjUV9rMyfLgTrFArFoyjRM5eFEH2BdKAzYCWlfD6nYxwdHWX16tUL2TKFouxx6NAhKaXMlTNTojuXpZTLAYQQPmgrkOZI9erVCQkJKVS7FIrcIKUkKSmJ2NhYYmNjSUhIICUlhdTUVFJSUgxbdp/T09MN2Rwzv+bmfeay3DgeuXVODh06dDi3916ihSe3CCFGAiMB3NyyDzxUKJ6UpKQkLl++TFhYGNevXycyMpKIiAjDa3R0tEFk7m16vd7YZhuFMiE8Usp5wDzQgkSNbI6iBKPX67l48SLHjx/nxIkTnDhxggsXLhAWFsbNmzcfqm9hYYGzszNOTk5UrFiRqlWrUq5cuWw3a2trzM3NDZuZmVm2n83MzDAxMcHExAQhhOE1N+8zl937XFDkZeZ/mRAehSI/SCm5dOkSe/bsYe/evezbt48TJ06QlHQ/uZa7uzu1atWiW7duuLu7GzZXV1dcXFywtbUt+aE4hYASHoUiE2FhYaxfv57NmzezZ88egxdja2uLn58fo0ePpn79+tSvX5969epRrlw5I1tcMinRwiOEMEW7Bx2gE0JYAmlSyofzOSoU2ZCamsru3btZt24d69at48QJbf0Cd3d3OnbsSEBAAAEBAXh7e6u5XgVIiRYe4EPgo0yfg4FJwESjWKMoEUgp2b9/P4sWLWLJkiVERkZiZmZGYGAgw4YNo2vXrnh5eakmUiFSooVHSjkRJTKKXHLu3DkWLFjAokWLOH/+PJaWlvTs2ZMBAwYQFBSkmk1FSIkWnrwihHBo1qyZsc1QFCHp6emsX7+e2bNns2HDBoQQdOjQgQ8//JC+fftiZ2dnbBPLJGVKeIA5xjZAUTTcvn2bn376iTlz5nDhwgUqV67MpEmTGDZsGK6urjmfQFGolBnhEULYAP2MbYeicLl+/TrTp0/nu+++Iz4+ntatW/PJJ5/Qp08fzMzMjG2eIoMyIzyAF6AH1K+vFHLp0iU+++wzfvzxR/R6PYMGDeKtt96icePGxjZNkQ1lSXhsgRig9CyOruDixYtMnjyZhQsXIoTghRde4N1338WjFCR2K82UJeGJA1RPYikhOjqaqVOnMmfOHExMTHjllVcYO3YsVavmKlZYYWTKkvCcpWzdb6kkMTGRWbNmMW3aNGJjY3nhhReYNGmS6jAuYZToRGB5QUoZDyw3th2K/KHX6/n111+pXbs27777Lq1ateLo0aP88MMPSnRKIGVGeDIo+EWgFYXOpk2baNasGUOHDsXZ2ZmtW7eydu1avL29jW2aIp+UKeGRUt4ytg2K3HPkyBE6depE586duXv3LosXL+bAgQO0a1fylyIq6xSp8AghmgohHul1CCFGCyHU+GcZJywsjOeee45mzZpx6NAhvvrqK06dOsXAgQMLNH+MwngUdWfrZLS5NN88Yn8QWv7kXkVmkaLYcPv2baZNm8asWbMAeOeddxg3bhzly5c3smWKgqaoHx8+wOPWCt4B+BaRLYpiQlJSEtOnT8fT05Pp06czcOBAzp49y6effqpEp5RS1MJTHoh/zP4kwKGIbFEYmfT0dBYtWkSdOnUYO3Ysfn5+HDlyhJ9//lnlxi7lFLXwhAEBj9nfCrhaRLYojISUks2bN+Pj40NwcDAVKlRg06ZNbNiwgUaNGhnbPEURUNTC8ycwWAgx4sEdGStBDASWFrFNiiJk//79dOjQgU6dOhEdHc2CBQs4dOgQQUFBxjZNUYQUdefyNKAT8K0Q4j3gBCABb8AdOAJMKWKbFEXA8ePH+fDDD1m1ahVOTk7MmDGDl19+GQsLC2ObpjACRerxSCkTgEC0dKWxQAegY8b7CUCrjBnGilLCyZMnCQ4OpmHDhmzbto0pU6Zw4cIFXn/9dSU6ZZgij12SUiYBH2dsilLK4cOHmTp1KsuXL8fGxoaxY8fy7rvvUrFiRWObpigGqKBJRYEhpWTXrl1MmzaNDRs2YG9vz/jx43n99deV4CiyoIRH8cQkJiayePFiZs+ezZEjR3BycmLatGmMHj0ae3t7Y5unKIYo4VHkm0uXLjF37lx++OEHbt26hbe3N99++y3BwcHY2NgY2zxFMUYJjyJPxMTEsGzZMhYtWsS2bdswMTGhd+/evPrqqwQGBqq1qBS5QgmPIkeSkpLYuHEjixYtYvXq1SQnJ1OzZk0mTJjAsGHDqFatmrFNVJQwlPAosiUsLMywrO/WrVtJSEjAycmJkSNH8uyzz+Ln56e8G0W+UcKjAODKlSvs2bOHvXv3snnzZk6ePAlAjRo1ePHFF+nevTvt27dXS8QoCoQSLTxCCAfgR7TZ0FHAe1LK34xrVfFGSsm1a9c4ceIEx48fZ//+/ezdu5erV7UQOSsrKwICAhg+fLhaQ1xRaJRo4UFbGTQFcAEaA38JIY5KKU8Y1yzjotfruX79OmFhYYbt4sWLnDx5khMnThATE2Oo6+bmRqtWrQgICCAgIICGDRsqr0ZR6JRY4cm0Mqi3lDIO2CWEWA0MAcYZ1bh8IKUkNTWVlJQUUlJSDO+TkpKIjY3NdouLi+P27dtERkYSERGR5TUtLS3L+Z2cnKhbty6DBw/G29ub+vXrU79+fRwdHY10x4qyTIkVHjJWBpVSns1UdhRok9OB+/bt49VXXwW0f/hHbYW5X6/XG0QmJSXlIaHILeXKlcPZ2RknJyfc3d3x8fHBxcUFNzc33N3dqV69Om5ublhbW+fr/ApFYVCShefeyqCZiQHKPVgxI+XGSNCaFubm5jg5ORn6LoQQj9wKa7+pqSnm5uaYm5tjZmaW7ft7W7ly5bLdbG1t0el0hfPtKhSFSEkWnuxWBrVDi3TPgpRyHjAPwMfHRzZt2pR169YVvoUKhSJbSrLwnAVMhRC1pJShGWWN0HL8PJJDhw5FCSHi0UbBSguOlK77gdJ3T2Xhftxze7C41xdREhFC/I6WSGw42qjWOiAgp1EtIUSIlNKnCEwsEkrb/UDpuyd1P1kp6YsUjQasgAhgMTCqrA+lKxQlgZLc1Lq3MmhvY9uhUCjyRkn3ePLLPGMbUMCUtvuB0ndP6n4yUaL7eBQKRcmkrHo8CoXCiCjhUSgURU6ZEh4hhIMQYoUQIl4IESaEGGxsm54EIYSFEOLHjHuJFUIcEUJ0MbZdBYEQopYQIkkIsdDYtjwpQoiBQohTGb+780KI1sa2Kb8IIaoLIdYJIW4LIW4KIWYLIfI8SFWmhIes0ezPAnOFEPWNa9ITYQpcQYtPswfGA38IIaob0aaCYg5w0NhGPClCiCDgf8ALaOE8gcAFoxr1ZHyDNn2lMtrcuTZo01ryRJkRnkzR7OOllHFSyl3AvWj2EomUMl5KOVFKeUlKmS6lXAtcBJoZ27YnQQgxELgDbDG2LQXAJGCylHJfxt/ompTymrGNegJqAH9IKZOklDeBDUCeH95lRnh4dDR7SfZ4siCEcM1sfD8AAAZBSURBVEG7zxI7iVIIYQdMBt4yti1PihBCB/gATkKIc0KIqxlNEytj2/YEfA0MFEJYCyFcgS5o4pMnypLw5DqavSQihDADFgG/SClPG9ueJ2AK8KOU8oqxDSkAXAAzoD/QGq1p0gT40JhGPSHb0R7Wd4GrQAiwMq8nKUvCk+to9pKGEMIEWIDWfzXGyObkGyFEY6Aj8JWxbSkgEjNeZ0kpb0gpo4Avga5GtCnfZPzONgLLARu0QNEKaH1YeaIsCY8hmj1TWY7R7MUdoSX9+RHt6dpPSplqZJOehLZAdeCyEOImMBboJ4Q4bEyj8ouU8jaaV1BaZuk6ANWA2VLKZCllNPAT+RDSMiM8Usp4NKWeLISwEUL8f3tnF2JVFYbh5+1v1BQ0SGQEpQwkdXQSjKAuRI3mIlGEfkglL4TqIjOxhCyQMvq5mupCSSHFEiUNwcKLICSCSPCHIsrUnDHUMcVMZBw1+bpY69jyeI4zZ85pHz37e2DD3muvvb7vDMzLWt85+10PA7MIM4WbmVXA/cBMMzvfW+cbnI+BMYQlSSuwGvgKeKyeSVXJJ8CLkoZLGgYsBr6sc079Is7YDgMvSLpN0lDgWUKttCJyIzyRhnqbXdJo4DnCP2mXpHPxmFvn1PqFmXWbWVfhICyPe8zsZL1zq4K3CD8L+A34BdgLvF3XjKpjDtAGnAQOAv8AL1c6iL+r5ThO5uRtxuM4zg2AC4/jOJnjwuM4Tua48DiOkzkuPI7jZI4Lj+M4mePC4zhO5rjwOI6TOS48OUJSh6QZGcRZJ2nl/xzjZ0lTazSWRXfAmv+iWNI30Unxu1qPfTPjwtOARIE5n7xCcU5Sc73zqiVmNt7MdkLNBHWSmS3vS0dJOyUt7EtfM5sGPF9VZg2IC0/jMtPMBifHsXon5DgFXHhyTFxi3Jdcr5O0UtIYSaclTY7tzZJOlVvaSHpA0p5oOL8ZGFB0v1nSVkknJR2WtCi51yFpqaQfJf0tabOkAcn9ZZKOxrH3S5qePDdD0gZgFLA9zuxelfSKpK1FOXwkqb2Cv81ySauS62GSLqW5xfaqY+URFx7nGszsELAM+EzSIIK1w7rC0iZF0h0EB7oNBL+Wzwne1oX7twDbCdYJI4HpwGJJqdXFk4Q3nu8BJgIL4rNjCcZmU8xsCMEeo6Mo1/nAEf6b4b0PfAq0RdsG4i4IT1GZBUoLsC+5bgX2m1lPUb9axModLjyNyzZJZ+JRsTWlma0BDgA/EHYUKFf/eIhg79luZpfMbAtX7w4xBbjbzN40s4tm9juwBng66fOhmR0zs9MEkWqN7ZeBJmCcpNujqf2hPuR+HPgWeCI2tQGnzGx375/8CqWE5xrfmRrFyh0uPI3LbDMbGo/Z/RxjDTCBYN15oUyfZuCoXe2v0pmcjwaaExE8A7xGcEws0JWcdxP8sTGzgwTjrBXAn5I2VVAkXw/Mi+fzqGAGEmdxY4CfkuZJXC1ENYmVV1x48k03MCi5HlE4kTQYaCfYqq6QdFeZMY4DI6MFa4FRyfkfwOFEBIea2RAz65NdppltNLNHCAJmlPb3LWUqtQ2YKGkC8DjBCL+vjCOIaTdcsZedSnmnvWpi5RIXnnyzD3hG0q2S2gibsxX4ANhtZgsJ9qOry4zxPcGFblG0w5wDPJjc3wWcjUXigTHWBElTektO0lhJ0yQ1AT0E8/TLJbqeAO5NG2ItZguwEdhlZkd6i5fQAgyPRfaBBBfB0RTVl2oUK5e48OSbl4CZhM3z5hK3KZE0i1CrKPz+ZAkwuZSlqpldJNhhLgD+IhRWv0juX44xWgl+vaeAtYSdT3ujCXg3PtMFDCcs04p5B3g9LuWWJu3rCSJS6dKnhbCbwg6CvecJwu6faZ2reJbV31i5xK1PnYZF0ijgV2CEmZ29Tr8e4AKhyP2GpB3AWjPbWqb/HsLuoNuStpKxJH1NKMDvMrPptfhcjUDFm607zs1A/Bp/CbDpeqIDYGYDippaCMbspcYdT9jVY29fYpnZo5Vn3/i48DgNh6Q7CcujTsKSsZJnhxGWdAdK3HuP8K3VMjPrrDZWnvGlluM4mePFZcdxMseFx3GczHHhcRwnc1x4HMfJHBcex3Eyx4XHcZzMceFxHCdzXHgcx8mcfwG1mu0ZJ/glewAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a4c9256d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = np.arange(0,20,0.1)\n",
    "bin_cent = get_center(bins)\n",
    "\n",
    "\n",
    "#fig = plt.figure()\n",
    "#fig.figsize=[5,3]\n",
    "#ax1 = fig.add_axes([0.1, 0.6, 1.5, 1.2])\n",
    "#ax2 = fig.add_axes([0.1, 0.1, 1.5, 0.5])\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1]) \n",
    "ax1 = plt.subplot(gs[0])\n",
    "ax2 = plt.subplot(gs[1])\n",
    "fig.subplots_adjust(hspace=0, wspace=0)\n",
    "\n",
    "\n",
    "mask = ~(np.isnan(SERVS_elais['f_servs_irac_i1']))\n",
    "plotting_servs = np.histogram(SERVS_elais['f_servs_irac_i1'][mask],bins=bins)\n",
    "ax1.plot(bin_cent,plotting_servs[0]/elais_servs_area,c='red', linestyle='-.', label='True SERVS counts')\n",
    "weights = np.zeros(len(SERVS_elais[mask]))\n",
    "weights = weights + 1/elais_servs_area\n",
    "ax1.hist(SERVS_elais['f_servs_irac_i1'][mask],bins=bins,color='red',alpha=0.3,weights=weights,log=True, label=None)\n",
    "\n",
    "mask = ~(np.isnan(SWIRE_elais['f_swire_irac_i1']))\n",
    "plotting_swire = np.histogram(SWIRE_elais['f_swire_irac_i1'][mask],bins=bins)\n",
    "ax1.plot(bin_cent,plotting_swire[0]/elais_swire_area,c='blue', label='True SWIRE counts')\n",
    "weights = np.zeros(len(SWIRE_elais[mask]))\n",
    "weights = weights + 1/elais_swire_area\n",
    "ax1.hist(SWIRE_elais['f_swire_irac_i1'][mask],bins=bins,color='blue',alpha=0.3,weights=weights,log=True, label=None)\n",
    "\n",
    "errors = np.arange(0.5,20.5,0.1)\n",
    "mask = ~np.isnan(depth_elais['ferr_irac_i1_mean'])\n",
    "irac_error_elais = depth_elais['ferr_irac_i1_mean'][mask]\n",
    "irac_error_elais = np.histogram(irac_error_elais,errors)\n",
    "probg_swire = gaus_prob(errors,5,irac_error_elais[0])\n",
    "cdf = gaus_cdf(bin_cent,4,1)\n",
    "\n",
    "probg_swire_interp = np.interp(bin_cent,np.arange(0,19.9,0.1),probg_swire)\n",
    "#ax1.plot(np.arange(0,19.9,0.1),probg_swire_interp*plotting_servs[0]/elais_servs_area,c='red')\n",
    "ax1.plot(bin_cent,cdf*plotting_servs[0]/elais_servs_area,c='black', linestyle='--', label='Predicted SWIRE counts')\n",
    "\n",
    "ax2.plot(bin_cent,cdf,c='black')\n",
    "#ax2.plot(np.arange(0,19.9,0.1),probg_swire,c='red')\n",
    "ax2.set_xlim([-0.05,8.2])\n",
    "ax1.set_xlim([-0.05,8.2])\n",
    "ax1.set_ylim([10**1.7,10**7.1])\n",
    "ax2.set_ylabel('c',fontsize='xx-large')\n",
    "ax2.set_xlabel('Flux density [$\\mu$Jy]')\n",
    "ax1.set_ylabel('N [$\\mu$Jy$^{-1}$ deg.$^{-2}$]')#,fontsize='xx-large')\n",
    "\n",
    "\n",
    "ax1.legend()\n",
    "\n",
    "plt.rc('font', family='serif', serif='Times')\n",
    "plt.rc('text') #, usetex=True)\n",
    "plt.rc('xtick', labelsize=12)\n",
    "plt.rc('ytick', labelsize=12)\n",
    "plt.rc('axes', labelsize=12)\n",
    "\n",
    "column_width_cm = 6.9\n",
    "width_cm = 1.4 * column_width_cm\n",
    "hieght_cm = width_cm / 1.3 # 1.67\n",
    "width_inches = width_cm/2.5\n",
    "hieght_inches = hieght_cm/2.5\n",
    "fig.set_size_inches(width_inches, hieght_inches)\n",
    "\n",
    "\n",
    "plt.savefig('./figs/completeness.pdf', bbox_inches='tight')\n",
    "plt.savefig('./figs/completeness.png', bbox_inches='tight')\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/rs548/GitHub/herschelhelp_internal/herschelhelp_internal/utils.py:77: RuntimeWarning: invalid value encountered in log10\n",
      "  magnitudes = 2.5 * (23 - np.log10(fluxes)) - 48.6\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([22.41929121, 22.04124374, 22.15134019, ..., 23.62154175,\n",
       "       18.94781666, 25.58594778])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flux_to_mag(SERVS_elais['f_servs_irac_i1'][mask]*1.e-6)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/rs548/GitHub/herschelhelp_internal/herschelhelp_internal/utils.py:77: RuntimeWarning: invalid value encountered in log10\n",
      "  magnitudes = 2.5 * (23 - np.log10(fluxes)) - 48.6\n",
      "/Users/rs548/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/matplotlib/font_manager.py:1297: UserWarning: findfont: Font family ['serif'] not found. Falling back to DejaVu Sans\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR4AAADYCAYAAAA9D4zLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd0lEUXh5+bThKSEELoIXTpAQIoTXov0qsFUVCx8KkoKkoRuzQVQRQsIE1EpUvvIiRIkd4h9JqEAKnz/fFuYgiksi2bec7Zk+y8ZX672b2ZuTP3XlFKodFoNNbEydYCNBpN3kMbHo1GY3W04dFoNFZHGx6NRmN1tOHRaDRWRxsejUZjdbTh0Wg0VkcbHo1GY3W04dFoNFbHxdYCrE1AQIAKDg62tQyNxiEJDw+/opQqlNl5ecbwiEhHoGO5cuUICwuztRyNxiERkVNZOS/PTLWUUouVUoN8fX1tLUWjyfPkGcOj0WjsB214NBqN1dGGR6PRWB1teDQajdXJM4ZHRDqKyLTIyEhbS9Fo8jx5xvDoVS2Nxn7IM4ZHo9HYD9rwaDQaq6MNj0ajsTq52vCIyCMist70OCwiE2ytSaPRZE6WYrVEZGMW73dHKdXqAfRkC6XUX0ATABH5AfjdWn1rNJqck9Ug0TrAc5mcI8CkB5OTM0TEFagLPG2L/jUaTfbIquHZqpT6MbOTRKRvTkSIyIvAU0A1YI5S6qlUx/yB6UAr4ArwllJqdppbtATWKKWSctK/RqOxLlkyPEqp5lk8L6fTrHPAWKA1kC/NsclAHFAYCAGWishupdS+VOf0AL7PYd8ajcbK2IVzWSm1UCn1O3A1dbuIeAHdgHeVUjeVUpuBRcDjqc5xxZgKbraiZI1G8wDkyPCIyBPmFpIOFYBEpdThVG27gSqpnrcA1mY0zRKRQSISJiJhly9ftpBUjUaTVTKcaolI5fs1A4OBnyyi6G68gbTBVZFA/uQnSqnlwPKMbqKUmgZMAwgNDVVm1qjRaLJJZj6ebcACDGOTmlKWkXMPNwGfNG0+QHR2b5Q69alGkyFnzsCBA5CQ8N9DKShTBipVgthYOHQIqlcHDw/juasrONmF5yJXkJnhOQAMU0ql9b0stZykuzgMuIhIeaXUEVNbDWBfBtfcF6XUYmBxaGjos+YUqLFT7tyBGzcgMBBE4OJFOHgQdu+Gffvg6lUYPRqqVoVFi+Dll2HdOihdGt5+G2bNuv99nZwgyTSrnz0b+vSBUaPgq68gMtI4fvw4FC0K+dKuk2iSyczwtARi0jYqpdqbU4SIuJi0OAPOIuIBJCilYkRkITBGRJ7BWNXqDNTPQR96xJNX+Osv6NjRMC6uruDmBjGpPsb580OBArB8OZw4AYcPG6OZDRvg33+haVOoUuW/UYyzszHiOXcOTp40DErx4uDpCYsXG88fewyWmv4fv/02HD0K7dsb96pTB2rWNO6nAUCUyrrLQ0QClVKXzC5CZBQwMk3zaKXUKNM+nhkYRvAqMPw++3iyTGhoqNJVJhyAhATDEAQFGc+nTQM/P+jZE6KjoV07qFzZMD5xccYIpFgxCA42jI6k9R6Ykd27YcsW2L4drl0z2goVgr59YdAgQ5eDIiLhSqnQTM/LpuFZq5Rq9kDKbESqEc+zR44cyfR8jR1y9Kgxqli3Dtavh4IFYeJE49iQIYZRGTbMlgrvRim4csXwB23ebBiipCR4/nl4/33DADoYljI865RSTR9ImY3RIx7rERtrzFxu3YLEROP57dvG85s3jdnPrVvGTKVDB6hQ4b9rt3z9NfOmTiX2zh18PT05c/o0p65fZytAkSJ85O3NBV9f+vbuTb2KFY2bu7vb6qVmjagowy+0YoXxs1cvYzTm4uIwjmk94kmDHvFYlps3ISzMmGUcO2YYnL/+Mny86RMFrAT2APnwz3eDIqXewsnFD6fLL3Dw4rf4CESrBAJcXGhYsiQ/DBuGR4kSDJk6lemrVhEbH89HTzzB8O7drfI6zUJEhOEjEoGffjKc3Tt3GitkuRxteNJBj3genKQkuHwZTp2C1avhjz8Mo5O82JPse630UBJB+bZwO2I1V879w42ok3TzKEUDgeO3TtH3xl4S0ty7eqGf8fBry8GIRKJu5wfcCfSNpWXIZZ5qfoam1a7g7GycG3XrFs9PmcLsDRv4avBghrQ365qHdVi3zlgF+90xEivoqVY6aMOTM6Ki4LffjEWclSsN/20y5ctDrVpQsUISwcWjKFjUj8h9Gxn9bguOJcYDxhb5EuLM6IJBtC5YiqOJCcy6fYMqtbsRVLUVKiEOFR+Lq5fh90hScOpSPg5GeHMwIj/bj/gRc8eFMkVieKfHER5vGoGriyI+IYEen3zCmj17OPT11xQrWNAG746ZOH4cSpSAbt1srSTHWMrwFFZKXXwgZTZCT7WyT3Q0bNwICxcaLok7dwx/bu3ahh83IAAKF7rAvgO/sGvXak6GLaVViUr0fGI6SQl3mDa5OzVLhhBcoREFyz2Ci0favaBZJy5B+PtQAX7bVoSj570pUySGUX0O07dxBPGJcRy/cIHKyStcuZWxY40p12+/GU6vXIhFDI/pxo2VUvckBhORPkqpOdm6mQ3QI56MuXQJfv0VfvkFNm0yVq3d3aFxY2jZEsp6/EupU3vxPbGL0et/Zu61cySgKOZfkjpuntQvXYdKzV+ymD6lYMcRP2ZvLM7xC15ULhnN96/som6FGwB8v3o1zWvUIKhQIYtpsBi3bsGIEXD6NKxdCw0a2FpRtrGk4bmMsa9mhFIqXkT8gG+AmkqpChlfbXu04bkXpWDNGpg61fDXJCQYI/569SAkBEq6hrF95mssO7iVO4kJnAOSnF0Z5u5NlG8RQtq/RakiFa2qOUnBXwcLMGNVENdvujJu4H56NtxJxeefp3LJkmz6+GNckp1BuYnISBg+3FjyW7fO+APkIixpeIph5L4pDHwJjAKWAa8qpe7Z5Wwv6KnWvcTGwpIl8NFHEB4OPj7QrBk0b27sywuf+CQT1v1EDKCApvn8qBsUQsv6j3OryEMo56zmkbMc0bedmbSoDNuPFKB/kwha1fyCJyZ8xui+fXmvd29by8sZFy/Cm28aIR99+8I778BDD9laVZawmOEx3Twf8DdGeorpSqlB2ZdoG/SIx1jq/uQTI0QpKgqKFIHu3aFx40Ri9q4gPvoaRbwrcib8V8L++QNPN0/qNXoGv9KZfp5sglIwf3Mxft5QgkaVrxLo153ft21k66efUreC3Q/C78/164avZ+VK4z/ECy8YsWV2vunQkiOeEOBn4AjwHTAR2AE8r5S6kQOtViUvG57r143P7pdfGkvejzwCNUqEcevgxyw7f5iIiAMkJiZQ1tmVz4etxdnFzdaSs8XGff58uaQ0wYHniL5dm8J++QmfMAGxZHiEpYmMNDz769cbm6TKl7e1ogzJquHJyVh5DfCmUuo7U0frgC+AvUDJHNxPY2Gio2HSJPj4Y8N/2bo1PNfkAKFLR/PaD/OYB1QpVoXH6vWljFJUK9cg1xkdgMZVruHvHc/78yrg6zWVwa33oJTK3YbH19cIsejb14iuL1/eCMMICLC1sgciJyOeMkqp4/dp76SUWmQ2ZRYiL414bt2Cr782DM7VqxBa4zxVEl5j37HFTLlzk5B83myt1IpTVVvhF1TT1nLNxqGzXnwwvzy+Xgls+mgLJQtluH06dxERAe+9B9u2QdmytlZzD1kd8WQ7QEQpdVxEWorIDBFZbOosFCNpl90iIh1FZFpkZNqEho5HXBxMnmx8LocNg5Ilkxj88Ouc3F2cH/fNIS4hnr01O7P6hYVEtnvToYwOQMXiMYzsfZhLNyKpOXQd+07b7ZpH9nF2NpbZd++2tZIHIicjnpeAVzD8O28ppXxFpArwrVIq23lyrI2jj3hWrDD8kCdOGNkXHn8cVn1UiLVRV6jj5sngDm8TWLmlrWVahfV7LzP+jw4U83+HU9MfxsXZwbLePvqokVvIjqaSFhvxAEOBFkqpj4HkBOsHAetu5NDcxZEj0K8ftG0Ld6IP0KpoC8Z23EyVmO2EFqnMe9Xa887rq/KM0QFoUq0QpQrV59y1H/luZVFbyzEvN24YycU++MDWSnJETgxPfuCM6ffkfyGuGLWvNFbm6FHo1AkqVoR582KpWnUMV65XZ8v5NVw9tA6Ahj0/I7Tzezg52X7fjbXp++hjQARv/XSMm7dz4YbC9PD1NTZbvfuusc08l5ETw7MRGJ6m7WVg3YPL0WSHsDAIDYW1q+KpW/g5inqX5t9/R/JwxaZMe3YWhau2trVEm1O3QgN8PYtwI2Y64363P2dsjhGBl14yks8/8QTs2GFrRdkiJ4bnJaCLiJwE8ovIIYxKnq+aU5gmfZSCuXOhSRNF/tgr7LxThZgL31Dw9g1G9vqcYV3Hkr+wfe/3sBbOTi50qteF4v4ufLowiPPX7DxZWHZwdTXyO3t7G/PsHGwGthU53bksGNU7S2FMu7bbe91yRwmZ2LIFXn0Vtm+/gJ/zQFYl7sDvkdYcCOkMBXN5dLYFOX/NnSHfVKPfo2f5YeguW8sxL3/+aSxjrlhhbNKyIZZ0LqMMtiulflFKbbN3owNGeRul1CBfX19bS8kRMTHwyivQsOEmdu9ogzPFuJW4jM31WnCg+Uva6GRCUf9YmlXfyo9rA/j7kJ+t5ZiXZs2MDYXvv59rRj1Z8jaKyJisnKeUeu/B5Gjux9mzRkqKAwf24UIL/EikU7HKNGj8LP7lHrG1vFzBsfMHWfnPk/h5/cDL33bir083O0qaY2PK1a0bfPONUaKnSRNbK8qUrL71JVM9ymM4l5sD5YBmpufaqWABTp0ycuGcPAmPlR5FQYlnwsAf6Pj0DG10skHpIhUI9C1K8YJT2H64APM3F7O1JPPSooUR7TtyZK4Y9WRpxKOUGpD8u4jMBfoopX5N1dYVw8GsMSNr1sDTT8ONG4mMGeNMxYq/0GvxbLwK66KE2cVJnGhRowNzNn5H0QIHGf5jJXo2POc4ox53dxgwwMhtYkcbCtMjJ297WyBtZuo/gHYPLif7iEgTEVkjIutEpIstNJibK1eMmMAWLRQ3ry/BL8YP/+2fwfbt2ug8AM2qG8ngSwVO4dRlT5aGFbaxIjPzyCNGBVSw+1FPTgzPUWBImrYXgGMPLid7mEodvwa0VUo1VUr9Zm0N5mbVKihf7gbz5r6Gp2dprkV3xNUpgYRYB4o3shGBfkWpUboOpy//TiGfO3wwv7y9fz+zj1IwdCi8/rqtlWRITrayPgP8JiJvAGeB4kAC0NWcwrJIfeA2sFhEbmHkBLpgAx1m4Ysv4JVXFN7SG1ErqVysPo881J/GVVrh7pr7ay7ZA4PbDCN/Pl+2HDjPlOWlWbsngOY1rthalvkQMQL14uMNI2Sn066c7uNxBR4GigHngb+UUvE5FiHyIvAUUA2Yo5R6KtUxf2A60Aq4ghGYOtt0rA8wzKSlBdBJKfVcRn3Za5DoggVG2e+mJY/yzOnK7K/YgLo9PrG1LIclLkF4bnJ1fL0SODx1reP4elLTsaPVu7T0Pp54pdQmpdQ8pdTGBzE6Js4BYzGSyKdlMkYcWGGgHzDFFA0PcAPYopSKw0hQVuU+19s969dDv37XKekzlUURNWlSsQF1u39sa1kOy+4TYXz0yyv0a3KMYxe8+GGNA+avU8pIpJ1kn1vs7MLOK6UWKqV+B66mbhcRL6Ab8K5S6qZSajOwCHjcdMp2oJJpJ3UIcE+CMnsmNhbeGp5I06bvEB9XlNORz7O3YBF2dRxht0NkRyAhKZ7wY3/h5b6Uh0pE89ZPlYi65WABtNu3mwL51tpayX2xC8OTARWARKXU4VRtuzGNbJRSV4HfgA3Ap8B9NzqKyCARCRORsMuXL1tYctaIiIB6dW/z8Sf9gQ9p51OIbzqN5NKzM0nwyG9reQ5NjeA6+Hj6sWn/Sga1PsXlKDc+XuBgq4U1axrVF6dNs7WS+2LvhscbSJsyMBIjNQcASqnJSqnGSqlHlVL3XVlTSk0DRgM73dxsn0s4PBzq1oV/9/QB5vJq2foMevEXilZvhziks8G+cHF2oUGl5vx9eBPFC16hdtlIpq8KcqwVLjc3ePJJo1LFBftbb8nwUy4iziLSVUQeExHnVO3W2ix4E0hb99YHiL7PuRliL7Faq1YZO5FdXeHV2mX5sE4vmvSZoA2OlWlcpSVxCbFsP7yJ0PI3uBTpzqGz3raWZV4qVIDERCOA1M7I7NP+ExCK4T/ZKCLJ49HnLarqPw4DLiKSOhyjBrAvuzeyh5zLCxdC2zZz8EwYwZjHwmnUthdVW+tsIragUskaPFyxCZ7uXtQua1RlWhGeC8seZ0SxYkY52ClTjMz/dkRmhqe4UuptpdQooC8wXUSamluEiLiYNgM6A84i4iEiLqbKpAuBMSLiJSINgM7AzOz2YcsRz/Hjxqi3W7fdqKSnKcUX+HroDYG2xEmceLvHJ9Qp35DCfnGUKHib5TsDbS3L/HTubJQYmZntr4xFyczwuImIO4BS6hTQAXgdY7+NORmBsRFwONDf9PsI07EXgHzAJWAOxibBXDPi+fFHIy3pnDnH8PTsgm8Bf4YO/B5nN70h0B64eTuK89ciqF0ukg3/FuRWrAOlRwUj439oKIwfb1dL65kZnleAlOQlSqlooBPwP3OKUEqNUkpJmsco07FrSqnHlFJeSqmg5M2DOejD6iOeOXOMIM/gAtNxj6+IU+Jl3n7sfXwLOuC+kVyIUorXZgzgu1UTqFX2BrHxzqzbU9DWssyLiJE57vBhWL7c1mpSyNDwKKV2KKUuJj8XkUClVKJSapblpeVelILp06F//wQe9tzDqMsv8JCrG1/3nkDF4lVtLU9jQkSoV7Ex/xzbRunAc3i6J7A83AGnWx4eULq0ke7ATsjuUspci6iwAtaaal24AJ07JfLMM4MplTSQPxI7ULPZIN59bSV+JatbtG9N9mlYqQUJSQmEH9tA8xpXWBJW2J5mJObBxQUmToSmZnfP5pjsGp5cu53WGlOtRYugatVEli4dCEyjffF/2P7yLI7VfxzJhbXI8wLli1Um0LcoWw6splzRGE5d8mTLAX9byzI/IsbS+uHDmZ9rBbJreHLtFitLjngSEmDwgEt07jyGW9dqk6R+5IWgmrR86lsSdVS5XSMiNKjUjF3Ht1Mj+AweronMXFfC1rIsw+TJxiayONuXwMtWdLqIrFVKNbOgHotj7uj0O3egTx/4/ffbuOBLJRdnOpSrzyPdPzJbHxrLcjnyInEJsRQvGMSEP8qw87gvF35ciYebg825Dh2C4GB47DFj+mUBshqdnt3ec+1Uy5xERSbxvyf+JMG3Bbv3OLN7tzCo9SVaV1qEq7cDDtMdnEK+/2UibFLtCuv2BrBkR2G6NzhvQ1UWoKKpyvjy5TZJmZGa7E61eltERS7jkRKDmbGoHT//XJtTB5+nskcFWlXco41OLubExcN8tnAEZYucpaj/HcedbikF8+bBmCwVjrEY2TI8qZfWcxvm8vEsG7uT/TdDKOVag4L5L3MjdhoB7lG4ejvYdvs8RlxCHJv2ryLs6Cb6Nj7LsvBAbtx0sFQZYDiZz5+HsWPhwAGbyTBbZKKIpK2nbleYY1Xr+PrjPPteUSo5t2Tii5OY/Nw8Xun4Hs8O+FEHeeZyKhSrQiGfImzev5pCvrEkJDqxcpeD/jN56inw8oK33rKZBHN+Wxqb8V52xdy5UK3iUco27cgF9Q8vdTmDs1c+3F09aF6jPb56ipXrEREaVm7OruN/U8z/LPnzJbB0h4NVoUjGz8+ow7V4MZw5YxMJZjM8SimblLexJPHxRtngPn1Oc/hwMzw4xsim/1LiIZ2oyxFpVKUlCUkJbD+8nlplb7AsPJDERFurshCtWhn+nu++s0n3eWZ+kFUfj0pKYtu6KwwfDlWrwhdfLMKDynhwhokdhlGzQRPrCNZYnbJFHqJmmXo4OzkTWu4GV6Lc2XHEweqsJ1OkCNSqBd9+a/yHtTLZ9p5lUEc9FogAVtijE1optRhYHBoa+mxG573f92NGzvsKJ1lFpQoBeNGDMi7Ca93G41++gZXUamyBiDC67xcARN+OxNkpiaVhhXn4oRs2VmYh2rY1nMxLlkAX69bCzMmIpwLwJtAUo3Z6U9PzmhgJwo6LSBuzKbQi27ffZNT8qbiJDzOGRDK83h88FVSFMUOXaaOTh0hITOBO3FnqP3Td8aqNpqZ2bShZ0kiZYeU5ZU4MjxPQWynVSCnVVynVCOiJkZT9YYz8ObmuNsuuXdCixWiUOsOrXYfj7+eCb1AIrZ+YirOHg6XE1GTIyNkv8+nCd2hf5yL/HPfl7FUHDXtxdjYShW3eDBs3WrXrnGxUaA30SdO2hP+yAs4CvnoQUdYk7ugxtmw8ze0iScyf3xx3t0745NPOY4dEQXyU4uqOeJIyCFeqVfZhflw7mVt3TgGV+WxhWSY+m+3cc7mDli2hbFmrR67nxPAcw5hSpTYuz/Ff7fQAwG7zel66BG+8Aa1bG8/PHoyhQu1qBLjG4ON2Bx+/IoiTA24c06CUIirmBnCZy1vSd6g2qNSCH9dO5tTlZfh6tWDPybT1BhwIEcPwAGzdCiVKQFCQxbvNyVTrGeB1ETkjIttE5AxGGeGBpuMVgXfNJdDcREbC+HE7adFC0aLFagqVvk4JryRKFgnA17+ENjoOjIjg4+WHq0/GIYdFChSjfLHKbDmwmuqloth9wsexSt/cjzlzoHlzw9lsBbJteJRSO4HyGMnfJ2CUFS5vasdU0vhbs6o0A8nL6Z5Jp3BKqk1o/moUc+6Hi9NtvPwcdA6vuQcRyVKoc8NKzTl6/iDBgXu4dtONgxEO7ufz9oYRIwxHsxXI6b/3Jhh+nkClVAcRCRURH6WUfdZL5b/l9HJFKz3bv+hDTDiymUQgv5sX4mT5BN+RUZH0e7oXAHv27aF6leqUCgpm8vipOb7nmvWr+WT8h7i4ulDAz5+fp89l7Kej+XP1Cnx9jNCQaV/OYMbMb/lz9Qp88vtQvlwFJn7yJe++/zZNGjWjeZMWAPy+ZCHHThyjRZOWvPneMJydnEhSSSyatwxXV9cHfwPuw+69u3BycqJaFfvLzNikWlsqFK+Cf35fZq6H1bsDqFTypq1lWZbq1Q0DFBtrPHd3t1hXOdnH8xJGEvjvMOqag1EV4gugvvmkWQYRaNxrHMHb53P82F+4WsmR7Ovjy5IFKwBo06Vlyu9g+B4MbdnLOjLuy09ZOGcR3l7e3LhxPaX9o1Gf0LD+3REsyW0v/G8wh48eokPbTsxdMDvF8CxdsYTXX3mDMR+PYtqX0wkqEURkVCQuFsrbAobhcXFxsUvDU8C7IAW8CwIJFPG7w5rdhXipw0lby7I8s2cbTtBu3WDSJIt1k5NP1VCguVLqpIi8aWo7iOHbyTUE1e1JUN2ed7V5d898+1F8i7bEPvdKyvlxPfsT17M/cu0Kyj8gWxrGfjqaixcvEnHuDE/2G8DBQwcY/to7/DT7B9zdPejVrTcffvY+W//egrOzM1+Nn0rJ4v9VqEhKUmzdtoVmjzbHz69AlvqMijJ2btetXY/h7w1DKUVCQgInTh6nYvmHyOeRj/Ub19Kza++UUVNqFi39nS+nTsLDw4O3h71LpQqVePalgdyMiaZm9Vp8OOoTxn46miYNm9KwfmMGvTSQkW+NYc36VWzcsoHIqEicnZ2Z8/18fvz5e6JvRrFxywae7DeAd8e8jaenJ3169KNPj37Zei8tweXIi/y+bRblir3Kur21SEgUXJwd3Nnj7Q1OTkbOHgsanpw4l/MDyZFlyX8FV8Dq+RRFJFhELovIetMj14UTly9Xgd/mLMbP917Dseff3Vy7cY0lC1bw4ahPmPT13fPvLz+fzLyFcwhtHMJnE//bOvXWqDfp0L0NHbq3IdG0MeytUW9SrV4lfHx8qVCuIiJCrZDahO3cwcYtG2jU4FEAxr73EWE7d/BI8zq89vZQUmeoTEhIYOLX41n8y3IW/7KceqEPM2PWdHp1683yhau4EXmDXXv+Sfe1FilchF9mLqSgf0EOHj7Ak/0GMHTIa0ydZEwF33/vQxbNX0bv7n0f6D01FyKweMd8YC5Rt1wJP2rb8tdWQQTatYMjR2DHDot1k5MRz0aMwnsfpGp7GVhnFkXZZ4NSqrs5bnQz1fQnu+dnd7STTEj1EMD4e6fcy/RlP3TkEBs3r6eDaSRWrGjxu64tX7YC0yf/QHx8PH0G9OT4CWNHQ3pTrerVQnhqcH9iY2Nxd3enQ5tOLP1zCZFRN+jf6wnAMA5ffD4ZpRRD33iJDZvX06SRscfj0uVLBJcqjYeH4Yx3cnLi5MkTdGzbyfRaanL85LG7poypDVelipUBKFqkKJFRd8fMPfvUYMZ9+Sk/zJrB888MIaR6zey+lWYnwKcwlUuGcOLCImA8Xy8Lpl7FXbaWZXmaNjUqUU6ZAnXqWKSLnIx4XgK6iMhJIL+IHAJ6ALYqAt5ARDaJyIeSXSeJHeAkxp/AJ78vFy5dAGD/QWOzWvmy5WnRtBVLFqxgyYIV9ziijx0/CoCrqyu+vr4kqYxzBPvk96Ftq/bM+3UOAI3qN2brts3s3ruLWiG177qniODvX5CkVLVeAgsFcvr0KWJNzsekpCRKlQpm1x7jy7hrzz+ULlUm5bUkJSVx6MjBlOvTGiRXV9eUEZl/AX/GfzSJEW+O5OPxH2b5/bM0jau24uy14wT67uCAo69sJePpCU2aGEvsP/8M33wDERFm7SLbIx6l1HkRqQPUAUphTLu2K5XJpz4DRORF4CmM0shzlFJPpTrmD0wHWgFXgLdSVRM9jxEvdgv4FugK/JpTHbaketUanIk4TY/Hu1LA5K8JqV6TlWtW0KF7G5ycnOjRpReP93ky5ZoJk8dx+MghnJydebjOI5QrUx4wplXJ/pnPP5xwVz89u/Si11M9eKLvU7i4uFC2TDm8vbxTjMKcBT+zYdN63NzcCC5VmkcbNkm51sXFhZeee4X23Vrj6enJW6+P4On+A3n2pYHMmPkt1avUoGaNWhT0L8jjz/Zl6YrFKa9Dj7BuAAAeKUlEQVTlftSpXZeXXn+BfQf+pUTxkixbuYSYmBheffF1s7yn5qDBQ82YtmIcnu4/czBiHElJhgvE4WnbFv78E/r3N543bAibNpnt9tmqMmEpRKQrkIQRjpEvjeGZgzEyGwiEAEuB+mnrp4tIO+BhpdR7GfVVvlglNX7gjynPi7Vzo2xQeTO9Ek1u4NjpI5xblnWX5GcL3+FOXAV2HJ3Mvq/WUTnIwZfVk4mIMOqt//GHUY/r2DEjvisDzFplIoNUGHeR2Zc+g+sWmvoJBVKybIuIF8aSfVWl1E1gs4gsAh4HhotIflM9d4BGgO2SyGoclmFdP+DcNXd2HIUtB/zzjuEpYfoqDhoEbm6ZGp3skNVBY8lUj/IYzuXmGNOcZqbnlhg2VMCIek9d/nA3UMX0e0MRCReRTUBxYHbaGwCIyCARCRORsMgYB82torEoRfzukD9fBJv358E0t+7uxupHQoLZbpklw6OUGpD8wNhw3kcp1cCUFqMhlit74w2kTRkYibGkj1JquVKqtilFxxNKqfu+M0qpaUqpUKVUqK+Xg2aU01iUKcs/Jjb+YTbvz6Ofnx07jKyFZnIy58RN1hb4PU3bH4Alci7fBNKGBvsA0fc5N0OSY7Vi7uSRYbLGrFQOCiEu4SzHL/7LheuWCyWwW4oXhzZtjHCKiRPhiSceKHlYTgzPUWBImrYX+C8thjk5DLiISOppXA0g28lRksvbeOmkXpocULd8I5ydXIFf2HIgD063ihWDXr1g/37YuxdmzoTPP8/x7XKaFuNVEYkQkb9FJAJ4zdSeI0TERUQ8AGfAWUQ8RMRFKRUDLATGiIiXiDQAOvNf0rHs9GHTEU9kVGTKbuKgSsXo0L0NQ1597oHuuWb9alp1aka7bq3oN9CY7TZr/9/GwdBGIWz9ewsAoz8ayep1qxj00kBOnznFT7N/oO6jNenYoy2PP9OHeFPC71oNqqfofGPEaw+kLzNmzf3Jovc3J14e3tQq+wjwC5v354EdzBnRuTPUrw/vvmuk7swBOUmL8Q+GI7kPMB4jPUZKWowcMgIj0HQ40N/0+wjTsReAfMAlYA7wfNql9CzqtumIJzlIdMmCFVR+qMpdGwKVUnft8M0qyUGiy35dyeRxUwAoWqQY586fIzLyBoGBgfyz2/iz7N77D7VCat11/dAhr7H4l+VUqVSVDZvXAxAYWDhF56djxz3AK86cWfOy/f/DpjSs3AyI4Pdtp2wtxbaIwAsvQEAA9OsH0dn2fOQsLYZSKh4w224ipdQoYFQ6x64Bjz1oHyLSEehYtMD9a2K/NdKNvfsebGdYtSpJfDQ66/tDLBEkWiukNjt3heHtnZ/e3fqyPfxvlFLciLyBf4GC99WRNnwhPaJvRjPkf4O5dv0aZYLL8sXnk5k9fxbfz5yOs4sL4z6cQMXyD/FY7w4sWbCC4yeOMXHyeL74fDItOzalWtXq7AjfzpgRHxATc5P9B/cZI6v/DWfW3JmcPReBk5MTi39ZnuX30JrUq9CIWmWnse9UXRITN5hzdTn34eMDzz0Ho0dDx46wbJmx4zmLZOmbJiLvZ/G80Vnu2crYesSTHuYOEq0dEsrO3TvZuTucOrXrEhcXx/ETxygdXOae+0+cPI46jWty8PCB/+KxLl1MmWp9/sWnd53//czptGrehiULVjDx0y+Jj4/n+5nTWf7bKqZM+IYPPkv/Y3L12lXeGfYuP8+Yx48/z6BD204pI786tepx5epllv76J3/MW5rt99BaeLp707hKM2ITfDh01r4+RzahZk343/+MRPE9epCdNI1ZHfEMFZEZZJ677WVgZJZ7tyKZjXiyM1IxJ+YOEq1ZvSZfTJlAYKHCDH3hVfx8/Vi7cQ21atS+p++hQ16jW+cePDGoHzcirxNQsFDKVOt+HDt+lDaDXwKMANHzF84TFFQKFxcXSgeX4fqNa+kGiBYqFEhB/wB8fRLuGWHly5ePbp17MOilgQQHBTP8tXdwstO4hFKB14GP+G6lJ+OfudeY5zkefRTu3AEXF6M+VxbJ6l/XC2M1K7OH3a4z2uuIx9xBor6+fkRHR5OYmIiTkxM1qoXw/czpKUGgaXF3d+fpJ55h6vSvM9Varmx5wnYaqRKSkpIoFFCIU6dOpuTz8fMtgLOzMzG3Yu56HXBvgKjRZjxPSEigZ9feTPtyOucunGPPv7sz1WIrggolIHzH/C1zbC3Ffmjd2sjXnA2yuoHQSSnlbPqZ0SPrkzzNXaQOEr1uyiYYUr0mfr5+dOjehk492zF3wd0bsydMHkerTs1o06UlJYsHpQSJlildlsBCgQDUqFaTw0cPUaNaSLp9t2jSkrUb1hAXF3fXVOu5V+4uuvpU/6dZ9ucSOnRvw//efBk3Nzee6v80bbu05Lmhg3h7mJHjv9mjzWnbtSXbdvyV4WsOqVaT/gN7E/bPDjr3bk/rx1pw6fIlKpZ/KHtvnhVxdXaikG9/zl79myPnztlaTq7FLoJErUGqqdaz3wz5L4BdB4nmPbIbJJqWyUvd+fOf2rze5TE+G/Bk5hfkIaRTpywFidrnRNoC2OtUS5P7qFrKG+jAD2vWkWDl0r+OQp4xPBqNuShXNAZ4hmL+VbgSFWVrObkSXb1Oo8kmRf3v4OPZmgaVqlKkwF5by8mVZNnwiMg6/kvufj+UUip7rm0rktlyukaTVZwEggJus2ZPAKcvX8bX0xNfLy9by8pVZGeqNQv4+T6P9UB14BFzizMn2sejMSdBgbc5dekCwc88w49r7baOpd2SZcOjlJqe+oGRGqMSRoDoQoykXZp02Lx1I9XqVaJD9zb0fboXd+7cyfb1H48zCnukF7x5+swpNm5Zn+37JbN33x7adm1J+26tadOlJbGxsXTs0ZabMUZgbY/Hu/KzKb5q5pwfmfb9VD4e9wHrN61LeX0de7SlS5+OXLt+FYAO3dvQvltrOnRvw+PP9MnWa84us+fPuis5vSUJKnSb2PhyVAkqy/erV1ulT0ci285lEfExhVAcBQoDtZRSg5RS5k1D74D06tabJQtWULd2Pf5Y+l9Ko+wGiaYXvGkYng051vf5pE/4evw3LP31T+b/9Cuurq6EVK/JblOtLDc3t5S6Wf/s2XnPbuhe3Xqz+Jfl9O7ejwW//5LS/se8pSxZsIKZ31l2051VDU/AbQAaV+7MrhMn2HnMEllhHJfs+HjyYVQRfQ1jetUwJ1HitiIrPp4O96kk+ljHrjzz5CBu3b5Fz8e73nO8b8/+9O3Zn6vXrlAwi7W1qlWpzp5/d/PxuA84c/YM5y+cY9qX05nx03ds3LIBJycnvho3haCSpXjxteeJOHuGksWDKF7MCJlo06UlK35bxbYdfzHyg3dxc3VlwOMDWfbnUv4O28aO8O38MW8pn074KEv3SyZfPk/WbVpLkcJF8clv5F+rFVKb8F3hlCwRRKWKlTgTYdRy/Hf/v1SrUp3V61be8/oio7KWXjYpKYmhb7zEsRNHyZfPkwWzfmP9pnV88KmR4vudN96jSaOmKa8XjL/RkgVG5Y06teuxYfM6BvQfSOWHqrB3/14692rP432e4vDRQ2zZthk3V1emfTmDokWKZklTVgkqZBieov5d8XD7kumrVlGrbFmz9uHIZGdV6wRGvpxPgTCgsIgUTn2CUspuJ7tKqcXA4vLFKj2b6ckWZuvfm6lZvRb/7t9L2TLlmDx+KvsO/Mu5C+dYsmAFh44cZPxXn/N47ydxdnLm97lLGPflZ8TH3b3pbfRH7zF7xlwK+geQlJREYEAgwaWCGfHGyJzd752xfPT5+3w5dRKPNmzC+I8mUTsklEVLfyeoZBC1atTm4sWLRN+MRkRwd787Qmber3NZsWo5SUlJLDcZCoDOvdojIlSs8BDjPpyY0r7szyUEBBTii88np4xUPh7/Ib/O/gOA7v27pASv3o8uHbsy/NW36dK3I8t+XUm1ytX4fe4SXFxcaN+tNcsXrsTJySlHKUcyw8czAV+veE5cLEH3+vWZv3kzE595BlcL1pp3JLLzLt3BWNV6Pp3jCsjVUXPpBUcCeObzzPB4VkY7836dy987tlGxwkO0bdWef/fvJaSaUTHz8NFDbPlrc8qoq3DhIpw8fYJqVasDRnjBjvC/0+03bVBlTu5XOLAwEz/9CqUUr771Cms3rKFF05aciTjDP7v/YfDTz3PqzCnm/TqHqpWr3qOlV7fevPXaCF4eNoSIs2dSanv9MW8pLvf5Qh49fpS6ofXu0i9CymjL2fnu15TWgFSqWBlXV9eUeLfUvPzC/3h+6CD8C/gz4s2ReHmaf9UpKOA2+07nZ+6wfnw+YIA2OtkgO87lYKVU6QweudroWINkH8jnH0zA2ZTMJfkLV65sBZo2bpYSEDp14reUKhnMv/uNfSL3C5wUkRQnblJSEi6uriQmJuX4fqmriAYUDEgZhRQsGMDefbspVrQYNaqFMOOn76gdcv9d8c7Ozgwd8irjv8o8LWb5suUJ27k9Rb/xUxEVHUVUdFTKa1FKERsby74Dd8/s0xaOTV2ZtHH9R/nmi+8ICCjEn6stk9+nZKHb7D/jTanAwhQukH7hQs296J3LdkK1ytUIDCxMh+5t6NijLT/Pm0lorTrExcXRuVd7jp04es817w0fTe+netKxR1t+X7KQShUr83fYXzz9/BM5ut8vv82jRYcmtO/WmjMRZ2jepAUAtWrUwt3dqJdevWoNDh05mG60OxjpOq5evcLFSxcBY6rVoXsbOvdqf9d5bVu15+LFi7Tr1opeT3YH4I2hw+nWtzNd+3Tizf+9BUDfHv1o27UlfyxdmOF72Kp5G/oN7M2ipb/T/5netO3aktXrVtLg4UYZXpdTggrdJvq2KxFXPNh3+jRN33mHw2fPWqQvRyPPBIkmoyuJah40SDSZf0/l5+2ZlVg+chs1Sh+i5NNP83qXLnz8ZN4NHNVBommwdbJ3jeORvLK1/0x+ivr70y40lB/XrtWBo1kgzxgevXNZY258PBMI9I1l3+n8AAxs0YIL16+zPDzcxsrsnzxjeNJFQVx8rEWWXDX2hVKKuPjYjCMOs0nloOgUw9MuNJTCfn5MX7Uqk6s0eX797+quOJJiT+PkYWslGmuQdAeuHzBffm1vjwT+OliAuHjBzdWFET173rPaprmXPG94Ys/DhfO2SfSuyf1ULRXNkh1F+PtwARpVucaLHTrYWlKuQE+1NJoHoHpwFE6iWL37vw2kN2/f5uf1660WN5YbydWGR0ReFJEwEYkVkR9srUeT9/D2SKRcsRhW7SqU0rZo+3b6jx/P2j17bKjMvsnV+3hEpCuQBLQG8imlnsrsmoCAABUcHGxhZRpN3iQ8PFwppTId0ORqH49SaiGAiIQCWUot6OPjQ7t27SyqS/PgODk54eTkhLe3N76+vpQvX5569erdE5iqsS9EZGdWzsvVhieriMggYBCAp6cnY8eOtbEiTUakNwr38PCgXbt2TJo0iRIldArb3EyeMDxKqWnANIDQ0FAVFhZmY0WazFBKkZiYSExMDNevX2fPnj2sXbuWb7/9lmrVqjF58mT69u1ra5maHJKrncsax0VEcHFxwdfXl+DgYDp16sTEiRPZvXs3lStXpl+/fvz222+2lqnJIdrwaHIV5cqVY926ddSpU4cBAwZw/PhxW0vS5IBcbXhExEVEPDAyIzqLiIeI5InpY17Gzc2N+fPnIyL07NmT2NhYW0vSZJNcbXiAEcBtYDjQ3/T7CJsq0liF4OBgZsyYQXh4OFOnTrW1HE02ydX7eHKCdi47Fs2aNePAgQMcO3YMT09PW8vJ84iIzsejcXxGjx7NhQsX9Kgnl6ENjyZX06hRI1q0aMEnn3xCTEyMreVosog2PJpcz+jRo7l06RLffvutraVosog2PJpcT/369WnQoAGTJ0/WEeG5BG14NA7Biy++yNGjR1m58t7Kphr7QxsejUPQtWtXihQpwpdffmlrKZosoA2PxiFwc3Nj8ODBLF++nKNH760ZprEvtOHROAyDBg3C2dmZKVOm2FqKJhO04dE4DMWKFaNr1658//333Lp1y9ZyNBmgDY/GoRgyZAjXr19n7ty5tpaiyQCrGh4RqSUiL2Rw/AURCbGmJo1j0ahRI6pUqcLkyZN1rTQ7xtojnjEY+ZHToyUw2kpaNA6IiDBkyBB27tzJ9u3bbS1Hkw7WNjyhwKYMjm8E6lhJi8ZB6d+/P/nz52fy5Mm2lqJJB2sbHj8go4CaO4C/lbRoHJT8+fMzYMAA5s6dy5kzZ2wtR3MfrG14TgH1MzjeEIiwkhaNA/Pqq6+SlJTEhAkTbC1Fcx+sbXh+AfqKyLNpD5gqQfQGFlhZk8YBKVWqFH379mXatGlcu3bN1nI0abC24fkQCAemishxEVksIotE5DgwBfgHeN/KmjQOyhtvvEFMTIz29dghVjU8SqlbQGNgJBANNAdamH5/D2iolNJJVTRmoWrVqnTo0IFJkyZx48YNW8vRpMLqGwiVUneUUmOVUjWUUp6mRw2l1AdKqTvW1qNxbN5//32uXbvG6NF6l4Y9oXcuaxyakJAQnn32Wb766isOHDhgazkaE9rwaByesWPH4u3tzdChQ/VuZjtBGx6Nw1OoUCHGjBnDypUrdVJ4O0EbHk2eYMiQIbRr145XXnmFv/76y9Zy8jza8GjyBE5OTsyaNYuSJUvSvXt3zp49a2tJeRpteDR5hgIFCrBw4UKioqKoX7++djbbEG14NHmKGjVqsGHDBmJjY2nYsCF//vmnrSXlSXK14RERfxH5TURiROSUiPS1tSaN/VOrVi22bt1KYGAgbdq0oUePHpw4ccLWsvIUudrwAJOBOKAw0A+YIiJVbCtJkxsoU6YMu3bt4v3332fJkiWULVuWdu3aMW/ePC5fvmxreQ6P5NZ9DSLiBVwHqiqlDpvaZgJnlVLD07suNDRUhYWFWUmlJjcQERHBtGnTmD59OufOnQOgcuXKVKtWjcqVKxMUFETx4sUJCAjA398fHx8fvLy8cHNzs7Fy+0NEwpVSoZmel4sNT01gq1IqX6q214FHlVId07tOGx5NeiQkJBAeHs6aNWvYunUr+/fvz3AK5uzsjLu7O+7u7ri6uuLq6oqLiwvOzs44Ozvj5OSEiKT8vN8DSPdnatK23e+cB8Ec9/vqq6+oW7dulgyPywP3Zju8gcg0bZFA/rQnmlJuDAIICgqyvDJNrsTFxYV69epRr169lLY7d+5w7tw5zp49y9WrV7l27RrR0dFER0dz+/ZtYmNjiY2NJT4+nvj4eBITE0lISCApKYnExESUUiQlJaGUuucBpPszNWnbzD1YMNf9XF1ds3xubjY8NwGfNG0+GJHud6GUmgZMA2PEY3lpGkfBw8ODMmXKUKZMGVtLcShys+E5DLiISHml1BFTWw1gX0YXhYeHXxGRUxbSFABcsdC9LYXWbD1yo+7sai6VlZNyrY8HQETmAgp4BggBlgH1lVIZGh8L6gnLyvzWntCarUdu1G0pzbl9Of0FIB9wCZgDPG8ro6PRaLJObp5qoZS6Bjxmax0ajSZ75PYRj70xzdYCcoDWbD1yo26LaM7VPh6NRpM70SMejUZjdbTh0Wg0VkcbnmwiIi+KSJiIxIrID2mO9RSRAyISLSL7RcQuHN8i4i4i000R/NEi8o+ItE11vLmIHBSRWyKyTkSytBfD0mSkW0QeFpFVInJNRC6LyC8iUtSeNac5b6SIKBFpYQudacnCZ8RTRL4WkSsiEikiGx+kP214ss85YCwwI3WjiBQHZgGvYuygHgbMFpFAqyu8FxfgDPAo4Au8C8wXkWARCQAWmtr8gTBgnq2EpiFd3UABDMdnMMamtWjge1uITENGmgEQkbJAd+C8DfSlR2a6p2F8PiqZfv7vgXq7XwyJfmT+wDA+P6R6Xg+4lOacy8Ajttaajv49QDeMGLatqdq9gNvAQ7bWmJHu+7TXAqJtrS8rmoHlQDvgJNDC1vqy8BmpCEQBPua6tx7xmI8w4ICIdBIRZ9M0Kxbjj2dXiEhhoAJGeEkVYHfyMWVUcj1marcr0uhOS+N02m1KWs0i0gOIU0ots6mwTEijux5wChhtmmrtFZFuD3L/XL2B0J5QSiWKyE/AbMADI0FZD2VnJZlFxBX4GfhRKXVQRLwxRmapuW+Uvy1JqzvNseoYJbA720JbeqTzXn8ItLKtsoy5j+6uQFXgV6AY8AiwVET2K6VylLhaj3jMhMlJ+CnQBHDDmCt/JyIhttSVGhFxAmZiGMUXTc1ZjvK3FenoTj5WDmPq8opSapMN5N2XdDSPBmYqpew2z2o6um8D8cBYpVScUmoDsI4HMKDa8JiPEGCjUipMKZWklNoB/A3Yy6qFANMx0sR2U0rFmw7tw4jqTz7PCyiLnUxbMtCNafVtNfC+UmqmjSTeQwaamwMvi8gFEbkAlMRw4L5pI6l3kYFu87sLbO3Aym0PjOmpB/ARxn8GD1PboxjpA0JM59UErgKtbK3ZpGcqsA3wTtNeCGNq1c30Wj4BttlabxZ0F8fwRQ2ztcZsaC4IFEn1OAP0SHueHep2BY5irHS5AA0wRsQ5XoCw+YvNbQ9gFEYqjtSPUaZjL5r+QNHAceA1W+s16Spl0nkHY2qV/OhnOt4COIgxpF4PBNtac2a6gZGmY6nbb9qz5vucexI7WdXKwmekCvAXEAPsB7o8SH86Vkuj0Vgd7ePRaDRWRxsejUZjdbTh0Wg0VkcbHo1GY3W04dFoNFZHGx6NRmN1tOHRWAQR2SciTcx4v5M5yV0jIk+JSKKI3BSRSubSk43+R4tIjCn3jo6NNKENjwNi+pLGmXLtpG7fZfoCBFtag1KqilJqvanfUSIyy9J9ZsBfSilvlcOAxgdBKTUSO4z0tzXa8DguJ4A+yU9EpBpGDTKNxuZow+O4zASeSPX8SeCn1CeISHtTissoETkjIqPSHH/ClArzqoi8m3q6YxrFzBeRn0ypMveJSGiqa0+KSAsRaQO8DfQyTXd2pz6e6vy7RkUi8niqvt9Jo8tJRIaLyDHT8fki4p/VN8bU1y8iMsukfa+IVBCRt0Tkkum9aJXq/AHyX0rb4yIyOM393hCR8yJyTkSeMY0qy2VVT15EGx7HZRvgIyKVRMQZ6IWRmjU1MRjGyQ9oDzxvSmCGiFQGvsaIiyqKkQ6zeJrrOwFzTdcvAr5KK0IptQIjB80803SnRtpz0mLqewrwOEb+l4JAiVSnvIxRyPFR0/HrwOTM7puGjhjGuQDwD/AnxvehODAG+CbVuZeADhjpQgYAE0SklklrG4x0ty2AciZNmkzQhsexSR71tMQIAj2b+qBSar1Saq8y0njswSgDnfzF6Q4sVkptVkrFYSTaShvYt1kptUwplWjqK1OjkkW6A0uUUhuVUrEYUdFJqY4PBt5RSkWYjo8CumfTebtJKfWnUioB+AUjSv9jZaSCmAsEi4gfgFJqqVLqmDLYAKwEGpnu0xP4Xim1Tyl1CyPnjiYTtJfdsZkJbARKk2aaBSAi9YCPMbLLuQHuGF9CMEYSZ5LPVUrdEpGraW5xIdXvtwAPEXExfZkfhLR9x6TpuxTwm4ikNkaJGHlk7jKuGXAx1e+3gSsmA5r8HMAbuCFGtYWRGKlAnQBPYG8qrWGp7nUGTaboEY8Do5Q6heFkbodRSSItszGmSCWVUr4Y+VjEdOw8qaY3IpIPY8qTIyn3aYvB+AInUyTV7+cxkmQl9+2Zpu8zQFullF+qh4dSKqtGJ8uIiDtGys/PgcJKKT9gGem8T6l1a9JHGx7HZyDQTN0/93N+4JpS6o6I1AX6pjq2AOgoIvVFxA1jCiH3uUdWuIgxdUn9edsF9BYRV5NTunuavjuISENT32O4+7M6FfjAlIEQESkkIpbKt5w8ErwMJJhGP6lTfs4HBph8aZ4YU1JNJmjD4+CYfBNh6Rx+ARgjItEYX5j5qa7bB7yE4e84j5Hc7BJG5Yzskjx9uyoiO02/v4uRYvU6hlGbnabvIaa286ZzIlLdbxLGSG2lSfs2jEoIZkcpFY3hzJ5v0tHX1Hfy8eXAFxg5iI9iJMuCnL1PeQadCEyTJcSokHADKK/sOFl5WkTkcYwVqjiMGmcW3URo2h39L+CulEoQkZEYq17ugFcqP1KeRhseTbqISEdgDcYUaxzGqKKW0h+auxCRLsBSjGKIPwJJSim7KF9tr+ipliYjOmOUbD4HlAd6a6NzXwZj+ICOYayuPW9bOfaPHvFoNBqro0c8Go3G6mjDo9ForI42PBqNxupow6PRaKyONjwajcbqaMOj0Wiszv8B2gPbJbCi8PQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11086eb38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m_1 = 17\n",
    "m_2 = 26.5\n",
    "bins = np.linspace(m_1,m_2,100)\n",
    "bin_cent = get_center(bins)\n",
    "\n",
    "\n",
    "#fig = plt.figure()\n",
    "#fig.figsize=[5,3]\n",
    "#ax1 = fig.add_axes([0.1, 0.6, 1.5, 1.2])\n",
    "#ax2 = fig.add_axes([0.1, 0.1, 1.5, 0.5])\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1]) \n",
    "ax1 = plt.subplot(gs[0])\n",
    "ax2 = plt.subplot(gs[1])\n",
    "fig.subplots_adjust(hspace=0, wspace=0)\n",
    "\n",
    "\n",
    "mask = ~(np.isnan(SERVS_elais['f_servs_irac_i1']))\n",
    "plotting_servs = np.histogram(flux_to_mag(SERVS_elais['f_servs_irac_i1'][mask]*1.e-6)[0],bins=bins)\n",
    "ax1.plot(bin_cent,plotting_servs[0]/elais_servs_area,\n",
    "         c='red', linestyle='-.', label='True SERVS counts')\n",
    "weights = np.zeros(len(SERVS_elais[mask]))\n",
    "weights = weights + 1/elais_servs_area\n",
    "ax1.hist(flux_to_mag(SERVS_elais['f_servs_irac_i1'][mask]*1.e-6)[0],\n",
    "         bins=bins,color='red',alpha=0.3,weights=weights,log=True, label=None)\n",
    "\n",
    "mask = ~(np.isnan(SWIRE_elais['f_swire_irac_i1']))\n",
    "plotting_swire = np.histogram(flux_to_mag(SWIRE_elais['f_swire_irac_i1'][mask]*1.e-6)[0],bins=bins)\n",
    "ax1.plot(bin_cent,plotting_swire[0]/elais_swire_area,\n",
    "         c='blue', label='True SWIRE counts')\n",
    "weights = np.zeros(len(SWIRE_elais[mask]))\n",
    "weights = weights + 1/elais_swire_area\n",
    "ax1.hist(flux_to_mag(SWIRE_elais['f_swire_irac_i1'][mask]*1.e-6)[0],\n",
    "         bins=bins,color='blue',alpha=0.3,weights=weights,log=True, label=None)\n",
    "\n",
    "errors = np.arange(0.5,20.5,0.1)\n",
    "mask = ~np.isnan(depth_elais['ferr_irac_i1_mean'])\n",
    "irac_error_elais = depth_elais['ferr_irac_i1_mean'][mask]\n",
    "irac_error_elais = np.histogram(irac_error_elais,errors)\n",
    "probg_swire = gaus_prob(errors,5,irac_error_elais[0])\n",
    "cdf = gaus_cdf(mag_to_flux(bin_cent)[0]*1.e6,4,1)\n",
    "\n",
    "probg_swire_interp = np.interp(bin_cent,np.arange(0,19.9,0.1),probg_swire)\n",
    "#ax1.plot(np.arange(0,19.9,0.1),probg_swire_interp*plotting_servs[0]/elais_servs_area,c='red')\n",
    "ax1.plot(bin_cent,cdf*plotting_servs[0]/elais_servs_area,\n",
    "         c='black', linestyle='--', label='Predicted SWIRE counts')\n",
    "\n",
    "ax2.plot(bin_cent,cdf,c='black')\n",
    "#ax2.plot(np.arange(0,19.9,0.1),probg_swire,c='red')\n",
    "ax2.set_xlim([m_1,m_2])\n",
    "ax1.set_xlim([m_1,m_2])\n",
    "ax1.set_ylim([10**4.1,10**7.5])\n",
    "ax2.set_ylabel('c',fontsize='xx-large')\n",
    "ax2.set_xlabel('Magnitude [mag]')\n",
    "ax1.set_ylabel('N [deg.$^{-2}$ dex$^{-1}$]')#,fontsize='xx-large')\n",
    "\n",
    "\n",
    "ax1.legend(loc='lower left', prop={'size': 8})\n",
    "\n",
    "plt.rc('font', family='serif', serif='Times')\n",
    "plt.rc('text') #, usetex=True)\n",
    "plt.rc('xtick', labelsize=12)\n",
    "plt.rc('ytick', labelsize=12)\n",
    "plt.rc('axes', labelsize=12)\n",
    "\n",
    "column_width_cm = 6.9\n",
    "width_cm = 1.4 * column_width_cm\n",
    "hieght_cm = width_cm / 1.3 # 1.67\n",
    "width_inches = width_cm/2.5\n",
    "hieght_inches = hieght_cm/2.5\n",
    "fig.set_size_inches(width_inches, hieght_inches)\n",
    "\n",
    "\n",
    "plt.savefig('./figs/completeness_mags.pdf', bbox_inches='tight')\n",
    "plt.savefig('./figs/completeness_mags.png', bbox_inches='tight')\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:herschelhelp_internal]",
   "language": "python",
   "name": "conda-env-herschelhelp_internal-py"
  },
  "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
}