skymap_class.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import healpy as hp\n",
    "import pandas as pd\n",
    "import scipy.constants as const\n",
    "from astropy.io import fits\n",
    "from astropy.cosmology import WMAP5\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from mpmath import *\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "from scipy.interpolate import interp1d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoYAAAFzCAYAAACjPchOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztvXmUrdlZ3vfsc2q+U9/b9+r27YnuRkIDQpMRMgaM5ECCLYztZJlRxPYi8cJ2cAJokYWBBDCsxApYsRNjh+AABovBJrYxgwEJEEFoAk1IaGikVneru3V7uvf2HWo8tfPHObtq11ffvIf3/c55fmvd7qo601d1qk796nn2uz9jrQUhhBBCCCEj6QMghBBCCCE6oBgSQgghhBAAFENCCCGEEDKDYkgIIYQQQgBQDAkhhBBCyAyKISGEEEIIAUAxJISQA4wx9xpjbhhjxj1vf8MY80CkY/l+Y8zPxrgvQghpC8WQEBIdY4w1xjy/8DE1omOM+V1jzH9T/Li19hFr7Ulr7aTP/c5u+6nZY/yUMeaHGo7jrxhjPmCMec4Y87Qx5reNMff3eWxCCInBkvQBEEJISowx476il5KZOP8rAP8lgN8GcBLAfw5A3bESQhYHJoaEkOwYY15rjPmMMeYfzJKyTxtjvsm7/KeMMf/CGPNbxpjrxpi3G2M+x7v8RbPLnjXGfNwY87WF2/5zY8yvGWNuAnhdh+O6b5Z2Ls3e/11jzA8ZY/5gVhP/R2PM7caYfz1L+d5rjLnPu701xjzfGPO3AXwTgO9ytyt5uFcAeMha+zY75bq19pestY+UHNeyMebnjDG/NKu7bxljbvcuf5Ux5qnZ9f6mMeYdxpg3G2OuGmM+ZYz5c7OPP2qMedIY8ze8277eGPP+2efzqDHm+0u+Hn/bGPO4MeYJY8wb2349CSHDg2JICJHiDgDnAdwF4G8A+HFjzAu9y78JwD+cXecDAP41ABhjTgD4LQBvAfA8AF8P4MeMMS/xbvuNAH4YwCkAvx94nF8P4Jtnx/m5AN4J4CcBnAPwUQD/c/EG1tofnx3vm2b18l8uud/3AXjRTOBeZ4w5Wfbgxph1AP8ewDaAr52J4+8C+Frvat8M4Oettbuz918D4EMAbsf06/TzAF4N4PkA3gDg//Qe7yaA/xrAbQBeD+DvGGP+auEwXgfgBZgmmv+jMeYryo6VEDJ8KIaEEEm+z1q7ba19O4BfxVHZ+VVr7e9Za7cBfA+ALzbG3APgqwF82lr7k9baPWvt+wH8EoC/7t32P1hr32Gt3bfWbgUe409aaz9prb0G4NcBfNJa+1Zr7R6AfwPglX3udLYW8bWYCucvAnh6lnb6gngawH8C8EkAf8urxH8aU8HDbFDmGwD8jHe7h2ZfnwmAXwBwD4AfnH2tfxPADqaSCGvt71pr/3j2tfoQgJ8D8OWFw/0Ba+1Na+0fYyrF39DncyaE6IdiSAhJwQTAcuFjywB2vfevWGtveu8/DOBO7/1H3RvW2hsAnp1d/jkAXjOrSa8aY65imi7eUXbbCFz23t4seb806WuDtfZd1tqvtdZeAPBlAP48phLs+LMAXgbgf7XWWu/j/wHAS2aDKl8J4Jq19j01xwxrbelxG2NeY4z5nVkVfQ3At2Ka0vr4X8/i80QImSMohoSQFDwC4L7Cx+7HVCocZ2e1sONeAI9779/j3pilaOdmlz8K4O3W2tu8fyettX/Hu60vUVJ0OgZr7XsB/L8AXup9+DcB/C8A3maMuehddwvTlPENmNbIflrYlbcA+GUA91hrzwD4FwBM4Tr3eG8XnydCyBxBMSSEpOAXAHyvMeZuY8xotibtLwP4t4Xr/YAxZsUY82WYVsT/xrvsLxljvtQYs4LpWsN3WWsfBfArAD7PGPPNs2GLZWPMq40xL+54jEvGmDXvXzHhDOUygMo9DWef239rjHne7P0XAfgaAO/yr2etfROm8vY2Y4yf5P0rAH9zdpsQMTwF4Flr7ZYx5oswXZ9Z5PuMMRvGmM8H8LcwfX4JIXMIxZAQkoIfBPAHmA5+XAHwJgDfZK39sHedz84uexzTQY1vtdZ+zLv8LZgOdjwL4M9gtqbOWnsd0yGIr5/d9rMA/hGA1Y7H+M8xrVTdv5/sePsm/iWmde9VY8y/L7n8KqZS98fGmBuYriX8d5h+rY5grf2HmA6gvNUYc272sXcA2AfwPmvtw8XbdODvAvhBY8x1AP8TpklkkbcD+FMAbwPwI7N1ioSQOcQcXbZCCCHpMca8FsDPWmvvrrj8pwB8xlr7vTmPa2gYY34bwFustT+R6P7vA/AQgOXZsA0hZM7hBteEEDJAjDGvBvAqAH9F+lgIIfMDq2RCCBkYxpifBvBWAP/DrFonhJAosEomhBBCCCEAmBgSQgghhJAZFENCCCGEEAJg2MMn7MAJIYQQQo5T3KS+NUMWQ0LIgvDyv//mqPdnJnH/rjz9sP6dXH7vV75L+hAIIQNgyMMngz1wQsiU2MJXRmwJLGMIYtgWCiQhc0HvxJBiSAhJRg7xK5JDBIvMkxg2QXEkZBBQDAkh6ZEQvSokBLCKRRLDrlAkCRGBYkgIiYcmAfTRJIM+FMPuUBgJSQrFkBDSD60S6KNVCB0UwzhQFgmJBsWQEHKUIQhfGdolsAyKYV4okIQ0QjEkZNEZqgg6hiiEDoqhLBRFQo5BMSRk0Ri6CDqGLIQOiqEuKIqEUAwJmTte8t1vxvJN6aOIx62LFuOdw9eqtafn40d484LBicemn8vJJ+ZDEFef3T5420wsti6sCR5NXLZvG+PdP/ud0odBSGoohoQMnZd89/EEcMhieOvi0R9RXwqBYYvh5oWjn4sTQ58hS6IvhsDxVHfIorh92/jYxyiKZA6hGBIyFMoEsIh2ISxKXxlFEaxCuyAWJbCMMjEsQ7ssFoWwijb1v3Z5LBPEIhRGMmAohoRopY0I+miUwjYi6NNWCh0a5bCNEDraiqFDoyC2lUJH17WhGkWxjRz6UBTJgKAYEqKJrjLo0CKFXUXQ0VUIHZrEsIsQOrqKoUOTIHYVQ0ff4SEtothVDh2URKIciiEhEvQVwDJyS2Ff+SvSVwbLyC2IfSSwjL5iWEZuWewrhGXEmjDPLY195bAMCiNRAsWQkFzElEFHDimMJYI+MaUQyCeGsYTQEVMMHbkEMaYYAmm2H8ohijHl0EFJJIJQDAlJSQoZBNILYQoZBOILoSO1GMYWQkcKMXSkFsTYYuhItT9laklMIYgAJZFkh2JISCxSSWCRmFKYSgCLpBLCIrEEMZUIFkkphkViiWIqISySawPzmMKYSg6LUBZJQiiGhISQSwYdoVKYSwQduYTQJ1QOc0khkFcMgXA5zCWFPrnPcBMqirnk0EFJJJGhGBLSh9xCCPSXwtwy6JCQQqC/GOYUQkduMXT0FUQJMQTkTn/YRxJzi6GDgkgiQTEkpAkJCSzSVgqlJNAhJYM+XcRQQgZ9pMTQp4skSomhj/Q5stvKopQg+lAWSQ8ohoSUoUEGHU1SKC2DPhrEEGiWQ2khdGgQQ0eTIGqQQkBeDH2aJFGDHDooiaQlFENCAF0i6CgKoSYB9NEig0WKcqhFBn00iaFPURK1SGERTZLoKJNFTYLooCiSCiiGZLHRKITAUSnUKoSAXikEjoqhRikE9IohcFQOtYohoFMOHb4kapRDgIJIjkExJIuHVhn02T2p/9tUsxQ6rPJD1CyGjts/ouR8izVolkPHtReckD6ERiiJBBRDsgi85B94Iqjw2b911+Tg7eVrI8EjqWbn0u7hO+PDL+L6p1YFjqaezUuHSdfaZ5cEj6Sa7fPT5/zEI4cp0uoVhd+cAE595vC5X7m2I3gkFewfft2Mzi/hEbZvP0wRt87pSxH3Vg+94I/+5XcIHgkRgmJI5pcjQuhQ9Oz7Qgjok8IjMugYH/0CahNDXwoBfWLohNDhi6FDmyD6YggolMP9o1+vockhoEsQfTF0UBAXCoohmT9KhRBQ88wXhdChQQxLZdAxLv8CapDDohD6aJDDohA6ysTQoUEQi1Loo0IQ98u/RtrlsCiGDi2CWCaHAAVxQaAYkuFTKYJFBJ75Kgn0kRLCWgn0qRBCHyk5rBNCh5QYVsmgT50Y+khJYp0YOsQEsUIKfYYqiD4SslglhkUoinMJxZAMl9ZCCGR/1tsIoSO3GLYWQqCVFDpyy2EbKXTklsM2Ugi0F0Mgvxy2kUJHdjlsIYUOzXLYRgwduQWxrRwCFMQ5g2JIhkMnEfTJ8Ix3EUGf1FLYSQJ9OgihTw457CKEjhxi2FYGfbqIoU8OSewiho4sgthBCn20CmIXOfTJIYpd5NCHojhoKIZkGPSWQiDpM95XCIG0UthbCIHeUgikFcM+QuhIKYZ9hNDRVwwdKQWxjxg6kgpiTzEE5k8OgbSC2FcMAcrhgKEYEr0EyaAj0bMdIoRAOikMEkIgSAqBdGIYIoVAOjEMkUIgXAyBdHIYIoZAQjkMEENgPuUQSCeIIXLooCQOCooh0UUUGfSJ9GyHiqBPLCkMlkCfQCEsEkMQQ2WwSCw5DJVBnxhi6BNLEkOlsEgUSQwUwiIaBTFUDn1iiWIMMfShJKqHYkh0EF0IgSjPdEwhBOJIYVQhBKJLITC/YhhTCoH4YgjEkcNFEENg/uUQiCOIseUQoCAqhmJI5EkihUDvZzq2DDpCpTC6EDoSiCHQXw5jC6FPiBzGlkIgjRgCYXIYWwp9egtiAikEdIohEF8OHX0lMYUYApRDpVAMSX6SiaBPx2c5lQz6dBHDZBLok0gIHV3FMKUQOrqKYQoZ9Eklhj5dJTGlGDo6C2IiMXRoE8RUYujTVRJTyaEPRVEFvZ9o+VM0kEGiTQpv3TVRJYU7l3bzSGEGNh/Ybn/dDFLYldRSmIvts+l/oXdl58xK+ysnlkKNrD6zlfwx1p6dYO3Z9t/jS9vpn4c/8y3/OPljkHQwMSSdyCKEjhbPcA4ZdLSRwqwymDgp9GmbGuYUw6bUMKcQ5kgMfZrSwxxpoaN1aphRDBcxOXS0SRBzpIYOpodisEom6cgqgz4Vz3BOGXRUSaFYKphRCn3KBFEqJawSQ4mEMLcYOqoEMacY+pRKolBSuMhy6KiSxJxi6ENJzArFkMRHTAiB0mdXQggdRTEUrYmFpBDQL4ZStbGUGALlckgxnKJJDiXE0FEmiFJyCFAQM0ExJHEQlUHH7JmVFEGHE0IV6wUFhdDhxFDDWkJfDKXXEUqKocMXRCkx9DmQRAVrCymIR3GiKCmHDkpiMjh8QsJRIYUzNEihQ4UUKqHLIEoupKVQCxqHUzRIITlOl2GV1HBQRR9MDAle9H1vxkjJ68StO5UcCIDli5vYud5h6jIhr/y8h/H+T94rfRhTNuXTMQBYuX0Lu5fXpQ8DgI7E0HH64QlWntPxc7RyRckfEiMDM9HxK2P31Ar2V/RkMjcupTv/eBcmq8AHfozpYURYJZPuvOj7pgmhtBTevP8wkTPbci+Yyxc3D96WFsI77roCALh04rmDj0mL4Z2Xrhy8/finzosdx8rth1uASIuhPTH94Vl9bPngY8vPySZ3px8+/IGWFMSlW4ePPdoWXnowOnxO9lenEj++JXtMu6cOX2MkRXH7zOEfNbsbst+7k9kSZgpiFFglk244KZTk5v27R6RQEl8KJbnjrisHUljklZ/7SOajOcSXQkl8KZTEnpgcSGGR3dMWu6d1/N26c1pHkrm/KphKjcp/P042ljDZ0JGWjXb2pQ8BALB8y2L5lvz37iv+LutlSZgYLhhlQiiRGJYJoURaWCaEEmlhmQz6aaFP7uSwTAolEsMqKcydGpYJoZ8Y+uROD/200EciOfQTQ0f25LBCCl1q6CORIPqpoUMiPfRTQ4dEejgp2S6V6WFvWCWTaprSwVxi2JQO5hLDpnQwhxhWpYKOKil05JLDuqQwpxzWJYW5xLAqIQSqxdCRSxCrxBDIK4dlUujIJocVUugok0OfHKJYJoY+uSSxTAx9cklimRj6UBI7wSqZlDOUyjiHFC5f3ByEFGrgzktXVNTHK7dvqaiP66SwDRqq5Z3TYxXV8v7qkmy13JIcNfPy9fqzxox29rPUzKvX6r+/WTEvFhTDOUaLFGqgzRrC1FJYt37Qpykt1MKdDzyd9P41CCEQLoUODXIIpF93WJcWamK03XycOdYhNskhoGsNojSUw/SwSp4z+shg7Cq5jwzGTgy7DpOkksIu6WAfIYxdKfdJCWNXyn2EMEWd3EUIm6rkIimq5boauYrY9XIfKYxeLTdUyGU01co+qSrmplq5SOyaualOLiN2xdxUJZfBerkSVsmEUuiQnjB2yeAQKmNtaEgJ6yaOY6FlcllDtTw0XIooPdEcO0VsqpPLiJ0gjntse8kEMT5MDOeAkMo4hhiG1MWxpDBEBmOlhSEiGFIfx0gNQ9YTxkoMQ6QwRmIYIoNdE8MiMRLEPomhI1ZyGFIjR0kOe6SFji6pYZFYKWLX1NAnVoLYJzl0xEgQ+6SGDqaHR+BU8iISYw1hqBiGriEMFcPQdDBUCmOkgqFrCkPFMHTIJIYYhiaFoWIYmhCGiiEQLochYgjEkcPQ9YXBchgghkCYHDpCJTFEDoFwQQwRQyBcDkPE0EFBBMAqefHgYImOyjiUGIMmkhtfA+FDKNL1cerauC3S1XJorSw+dBIohUC7gZQm5q1i7goHVIYPE8OBEVMI+6SFMWWwT1oYWwa7JoYx1w3Gnj7umhzG3o6ma3IYUwj7JIYxhTBGYujTJz0MTQx9uqaHsaWwc3IYQQp9YiSHjq4JYmhiWKRPghiaGvr0SRBjpIaOBU4PmRguAtIpoWRC2GYPwq50kcJ5GybRsEehJFpSwiqk00NphrDPYVu6Dqq02b6mC7n2QqxCOkFketgdiuFAkJRC6XMap6iMu0phbIayV2EqYtfHXb5HtEuho4scxkwLAU4rx6iUi0jKISBbMUtvkE057AarZOWkEsI2NXJKGWxTI6daQ9hGClOmgymlsKlOTp0UNtXJKdcTNtXJKYUwdpVcpKlaji2GPk21cuq1hY21cuQa2SdmpVykTcUcu1Z2tKmXY9bJRdrUyzHrZJ8FqpY5lTxvpE4Im8QwdUJYJ4aph0qaxDB1ZZw6LaySwxz1cZ0Yph4yqRPD1ClhajEE6uUwpRgC9XKYY+ikUg4TSiGQVgyBZjlMJYaOOkFMKYZAsxymEkPHAghi7x+O+VnIQaKwyJPGOdYQLmqFLD15TMLYOT2OfpaUITDaniSVQ1cvpzqbShOjnf3oZ1Bpi6uWY589hYTDxFAROdcRFhPDnEJYTAtzymAxLcw5UJJTCoupYa5hk7LEMKcU+qlhzrWEORJDHz89TJ0W+pTJYa5tao6lhonTQp/UyaFPURJTp4Y+RUlMnRr6FAUxdWLoM6fpIaeSh470cEkuJKWwyDxNGdeRcwK5uKchk8L5ojiUknPvwnmaVK5Dch/E4oBKn9Pk9YXDKXpgYiiMlBCOJjK1sdkeicngzvUVERmUqo/f/8l7RbalcamhhBTuXl4XmTrOnRg6lp8zWRNDh0sOJTa1Hm3vZU0LfXImh47xrb2sqaHP/sooa2ro2N0wWRNDnzlKD5kYDpFFSQl9pKTw3OmbC5MQOqT2KrzzgaeZFM45O6fHYmc6WZTk0DHZWMJoIrPVjNQWN0wPZaEYCvEF3ykjhVt3TLB5r4wUrt8uewo7CV5/10dEHvdLz31S5HEB4ML6TbHHXr10S+Rxt58nI0lGsDfZPSknaGYi84mPb8oO50mwfFNGDteuWCwJ/X0p9ftZC6ySM+O+4XY38j7u1h3TX1x2Jf8PuS+Ek0m+v0Xuv/DMwdvXtteyPS5wVAjfd/WerI/tS+EvfeYVWR/bl8KPfvZitscdjQ5fDraey9xB7R5+T68+mbd2W7l22BadeiTjerAb+97beSdqx1uHj2fHeStls3f4eU9O5F0+sD+WyXH2Ng6/p3dP5D+GvfXpc7yX9yUck9nj/vGPfnveB44Hq+QhIJkSSrB++6YKKcyNVEoILG5SKMZuYYpTKDmUZBGTQ2Bx0kN/ycDyzX2xBJHpYT4ohgLkSgu37piISqEUlML8SEqhnxYCwNrpbaEjyYufFgLA9XvlTmNHOVwsJOUwlyCONwdZSkaBVXJiin9tpJbCJhFMXSXXCWHqxLBOCFNXyXVCmLpKbhLC1HVylRTmqJKLUuiTvFLerf5+Tl0pF6XQJ0el7FfJRz+etlb2a+QyUlfLfpVcJHW1LFUlA0fr5CKp62VXJVdenrhinhQef0DVMqtkIlcZA8dr4yKSUpgayZRQmrqk8MV3XE762HVSuMikTg6rpBCQTQ6lmef0sG4CXbJeBuQq5nmGYpiQnGlhGylMlRZKTxsvshQuan0sTk1aCCzmWkPHotbKwHzLYROp5HCpRaWbUg6LlfIirDlklZyAqm+c2GLYNSGMLYZdhDBFYthFCGNXyV2FMHad3EUIU1TJXaQwdqXcJSlMUic3SKFPikq5rkYuErtWrksLj183fq3cVCX7xK6V62rkMmJXy5JVMlBfJxeJXS831cnHrh+5Xi7WyQ7ltTKr5EVjSFIYm/svPMOUUBBNgyZNzNsgShcpBBZ3IAVgeihJ7PSwTWp45Pqsl4NgYhiRuog5RloYsoYwlhj2FcJYiWFfIYyRGIYIYYzEMEQIY6WGfaUwRmoYsqYwSnLYISksEiM57CqFPjGSwy5p4dHbxUkOu6SFRWKkh10TQ58Y6eGQEkOfWOlh19TwyG0jJIhVqSGgNjlkYihN6nUHkoMlQPNwSWoWPSXUwJCSwugESCEgv+ZwkZNDDSx6eig5nAKkTxDnbd0hxTACiyCFJIxX3fZo0O1D6+P/6u4PBN0+VApTTyiTtPRNCx2hchiSFgLhtXJIWkimUA6HA6vkQNp8M/SpkWPKYJ8aObYM9qmSYyaEfavkmElhnzo55nrCPnVy7JSwT6UcMy3sVSkHpoU+fSrlkAq5SJ9KOVQKj99fd8kLFUOfPrVyTDHsWysPtUquok/FHFInH7uvHr8S6upkh6JamVWyBPP0F4I2JGtjh3R9vMhDJo6hV8hFFrlSdkhXyxxKISmZBy+gGPYk1ZMveRo7QH4tISAvha+/6yOUwjmVwqFNKcdMCx2UQx1yuOiCKL32MOXp9YYuh6ySe9DlSW9TI6cWwTZVcmoZbFMlpxTCtlVySiFsUyWnFsI2dXJqKWxTJ6dMClvXyZHTQp+2lXIKMXS0qZVj18jH77+5Io5ZI5fRplpOucawbbU8b1VykTbVcswqufT+W/yaaFMnO4RrZVbJuYj9l4D0YAmgY7gkdUp4ZrX+T8McKWHTAIp0SgjMb1Lo0yo1TCiFQLtKOaUUAkwOHU3pYerBkzbpobQU5kB6OAWInyAONTmc/++2iFAK0yBdHS8SodPJMaibUM61prBWDhNLoaNODlNLoaNODlOnhW1InRYSXVAOdcAquSV9n9yyKjm3EJZVybmFsKxKzi2EZXWyxFrCYqUskRQWK2WJpLBYKeceNCmtlDNJoU9ZrZxLDB3FWllCCou1soQUFmtliW1qyqplDYlh6iq5SFm1nLpKPvZ4JdVylyrZR6BWZpWckljGLz1Y4ljUlLBYJ0sPmACsjx0S08daB1FySyHAWtkhPZQCHJ9a1iCFEmhJD2MliENKDpkYNhD6ZO5uyFfGdmVfhQxOJiPx2vja9pq4EL7v6j3iQugSQ2kp/OhnL4pvSbP13KpIUlhk9cmxiBT6nHpkIl4hu+RQvEZW8LtxcmJZhRjmTgvL2D0xyp4YFtlb658YOjImh70PlGJYQwzDv/6AfEK4dkk+EQKAe89dkT4EfOl5+YRuY7QjfQgAgN975gXSh4B9GDz45AXpw8DWM+vSh4BTD8onZsvXLU49Kr+ub+3pxKepaIOS3427pyOc5zsQDWIIAJvn5Y9j+2y4nGaSQ1bJGtk5LX0EwNJ5+RfYe89dUSGFf/WOD+D80nXpw8B/cVK+wv6claelDwH7/V+3orKzLS9kAHDrTvnqDACu3yP/9dg72e/sIDFRcRo8Y7B8Xf4PyfHOPsY78l+P5Vvysr5yTf4YUsPEsISQpLAog9vnZRJDXwiXlmWOQYMMAlMh9Hl675TIcfy10+8/ePs3bny+yDEUhfBnnvhikeMoSqFEalgUwv0bcjIyvnn0b/SNx/P/zb58/ehLqlRyuLR19PVq6YbMRtBm9/A47JJQhmKO/pzsnloROQy7dPQ4JisyX4+9jcPH3d2Q+cPSFj71nTP9jyNxcsjEUAMaEkLgeEq4tysfv0tRlEIpfCkEZFLDspTwmy+9M/txaEkKi4xOyghIUQq1IJEcFqUQkEkPfSkElKSHgEh6WJRCKXwpBHSkh8B8Jog6X5EE6ZsWapDCpfNbKqpjQEdaqFUKyXFe8LynpA8BQH451CKFxbSQHGeR5bCIhloZmMqhBkHsK4daJ5VZJXv0eZKahDBHldxGBnPVyUMRwlx1cpMU5qiU26wnzFUpN6WFOSrlNmsKc9TKbaQwR6XcJIW5KuWytPDYdTLUysW0sEi2Wtk0J3U5quWmxDBXrVxMDMvIUS8Xq+QifarlRJUyq+RQukrhzmkdKaEWtAyYaIJJ4VG0VsikHg3DKJrQkhxqQctgihZWrtnOCaK25JCJIbo9KV1lMFVi2KcyTpUaahDCPrVxqtSwqxCmSg37TB6nSA67CmGq1LDP9HGK5LBPfZwqOexaIadID9skhcdukyg5bEoLiyRLD1ukhUVSpYdd1himSg/bpIVFUqWHTYlhkS4JYuTkkIlhDpgQHmXoKWGKrWu0pIQatqPpS4r1hlq2pOlLiq1shryuUMN2NkCi9LCHFKai6+AJ08PjDHE4ZeHFsG1aqEUKNQ2YaGHoQyYa9jVMRd/6eF6HUYY+bDLPtXLXtPDgdkqqZQ1DKZrQMpgCtJdDLZXyQlfJbZ6EGEIYo06OJYOhdbKmhDCGEMaok2OlhKGVcqyUMFadHLqmMFalHCMtjFEpx5DCWJVyaFoYq1LuUyOX3k9gtdxXCotEqZYjJIYxauVY29SE1st9auQyYlTLXWvkMtpUy5HCy0spAAAgAElEQVQqZVbJKdCUEpKjDD0l1EyMvQ1jDJrESA2HXiEXiVEpx6iQYySHsaRQE0wP9TK09FCShRXDurRQ08SxJinUkhbGlEINp8iLgaY1hTGnj0PkMKYUhlbKMSvkEDmMua5QU62sZc1hMBHXF4bIYcxNrTWtOdQkh3WCKF0pL2SVXPVFTymDXevkVELYp0qeRyEs0rVSTpUU9qmTU0lhn0o51ZY0XWvllElhl1o55ZrCrrVyqmGTPrVyqrSwT6Ucq0Yu0rlWTjh00rVaTnW2kz61cqwquUjXajlGjVxFVb0cWCmzSg5FS0IIMCUsQ0t1DKStj7sOosxrUhjCvNXHVaSYVO7DkJPDVFII6KmVAT3VMtPDcrTVywsnhmVp4SJJYdvzJg99K5pU5FhT2FYOU0thl7WGqaVQy5TyUEm9NU0XOUy9trCtHKaUwoPHUCSHWuiypU2qtNChXQ6lKuWFqpKLX+TcQlhXJ+dOCesqZS1CmDslrKuTcw+ZNFXKOZPCuko5d0pYVynnTgrrKuXc29LU1cq59yusq5ZzDpw01co5xNCntlrOvHdhXbWcqkYuo65aTi2FReqq5ZQ1chnFarlnpcwquQlpKaxC276EWqRQgqEMomipj7VUx4BMfVw1kCKxV6GWWhnQUy0PZiBlwBtah6KtWtaSIBbTw9zJ4cIkhu4LKymExcRQUgiLiaEmIZReT+gnh9Lb0fjJobQQ+smhtBT6yaH0mkI/OZTewNpPDqXPbOInh9Jb0/jpYe6ksMiR5FBYCv3kMLcUFimmh7kTQ59iepg7MfRx6WGP1JCJYROatqABOGBShbQU+khLoY+0FPpISyFwuOZQWgqBw+RQWgqBw+RQWgoBPckhoCs91LTuUMtQCqAvPdRC09Y2KZB/JcuA9J5AjtWnp4MfGqTQDaFQCsuhFB7nmy+9U4UUOjRIoUa01crSaaGDcliORjmUTAsdTg4l00KfnB6j5FNOx8u/7c1qksKtS3sqpBAAzp25idMbOo4FAF534RPSh3DA3zv7PulDOOBbb/uY9CEcsGvH+MY73i19GAdcOKdoTajRkzCYXT3yvnbVYvNC+OnZYhHjVHGxMFtxTisYg9GODnkHgMmqgZno+HlavmWxtKXjWFauWbz82/LI4dyLIQCsKQhcti5NXwR2r60KH8lUCrXwugufOJDCa5MN4aMBvu7UgwCAO8fyLwb3L01T3b+4If8NvGvbbXOUiydvngSgQw73b06Ty8lJ+V+uTgqfe0D4QApokEM7W8+3c25d+EgAzBJDsxl2Rp0YTFanP9vjW3sY39IjqxrkcH9WSmiQw731fH/wLYQYSrJ1ae9ACjVQlMKrW3IvkppSwq879eCBFGrASaFDUg6LUiidGjopdEjKoZNCh6QcFpNCaTlcu3r0l6mkHNrCkIcKOZxhNndVCKJDUg53Th19rdEgh46lLT3pYWrmWgxzxa5VVAmhRGp47sxNdUmhFqqEUCo1LEqhJFVJoZQcFqVQkqIUOjQkhw4pOSxKoUNDcugQk8OK9YWUw3LMxIoI4n7FEmZpOczhNXMthlJoTwkl8avjMnLXyZpSQqBeCnOnhk31cW45rJPC3KlhlRQ6csth3brC3HJYJYWO3HJYTAt9NCWHQH45dDVyGayWq5n39HBuxbBo1RrWGWpFsk4uQ8NaQ0fO1LBNUphLDrWuKawjlxw2SaEjlxy2GTaRrpWL5JLDOil0ZJVDRdPIdVIoQbFGJocU1xemTg3ncoPrui/a1vkkxzK9744p4fKZ7URH0i8lvG1tM8GRTOlTHZ8Z30pwJFO6JoWPT9It/O1THf/6rXTfyF2l8C2ffU2iI5nStT5+6tnqUxvGoK0YAsD4Rtpfdl0nkE9/KtGBzGhKC4usP5Vuq5Q2Ulhk5dl0r4FdpdCup9tep48UTjbSbQ/VRwrtON1rclWNXMXeWrpjqRs6+eD/UbvpNTe4lqZPdaxhQtknVXKoaT0hoKs+1rSeEOiXFKaslPusKUyZHHaRQiBtathnW5qUyWFXKQR0rTkEdFXLmtYcAunWHWpLCrtKISC/7jA2C7dD7NrTcVNDriWsJlQIr002oqaGIULoKuWYyWGIFLpKOWZyGFIff+Md746eHIYMmlw4dz1qcthVCH2cHMZMD0P2KnzugfjJYR8pdGxeWImaHPZJCn12zq3HTQ4D6mMnhzHTw5AK2clhyvSwLf6aw5TpYVucHMZMD3NuUeMzd4lhzknkGFKoLTWcVzSlhMB8JIVFYiaHmqaP55GYyWGIFDqYHNajLT2cV/qkhUVypoepfGfuxDAXmpJCIF5aGKtOjlUfaxpEAeIMo8SUwhjDKEMcNGlDrEo5JC30iVUrxzqzyTwOpISmhdFRNGwCxBs4iVUrx6qRNU0sA8Ovludq+KSLPfetk1MJYd9BlFT1cd9BlFTrCftWyqmSwr6VcqqksG+lnEoK+9bKqZLCvrVyLCn06VsppzrVXd9aOUZSWEbfWjmVFPaulRNJYd9aOdUUct9aOdXawr61coy0sIy+1XKXGrliCIXDJ13ps31NypSwT6U8b2sKY8P6uJ6USWGfWllbfZxCCgFdG2AD/ZLDVFIIsFZuQlutrGmvQ6BfephKCoF+6aHU2kLHworh0EkthV0r5dRS2LVSTi2FGs6l7NO1Up7X+riKrrVyKil0dJXDVGmhY+i18qJVyF3lMPWehV3lMPUksrZqeWjMTZXcdxFmU6Wcey1hm0o5Z1LYplLOmRS2qZRzJoVNlXLulLBNpZxTCttUyjmTwjaVcmop9GlTK6eWQp82tXLKtLBIm1o5txQ2VssZ1xW2qZVzbmTdplbOuT1Nm1o5ZVpYRlO13DctLKmTWSXPC3WVssT5jpuSw9z1cVNyqK0+zk1Tcpg7KWyqlHPXx03JYU4pBJqTw5xSCDQnhzmlEGhODtUlhZlpSg5zn92kKTnMvWdh03mWc0vhUFh4Maxba6ht8lgbXFNYXylLrSmsksNFq4+7klsKHfOw5jAlVXIoJYW1aw4FppC55nBY1K05lF5b6KAv4/im19JCuHtt9UilLD1kcnVr/UilLC2EZRtfSyaFxc2vtQ2ZALJS6FJDv1aWlEKXGvq1spQUOiYnJ0dq5dxJYZHiJti5k8IisTfBDsXJ4UGtLLwtjdncPVYrS54LuWwjbMkznLjU0K+WJdPCss2wtUghMCeJYcxNHqWlsIi0FDpcpSwthQ6/Ul70+rgMPzXUkhQ6QWRSWI7W5FBaCh1+crjoFXIZfnIoKYU+ThC1nfZOCzH3O4zpQXMhhjEY7eqSQjOyMCMdL8iOjSU9f7EDwKvWP61KCj9/ZV1VWvgXN57Gq1Yflz6MI2iSwpTnVO7L5OREPC300VYr722MVUnhzrl1TE7o2V7HbO5idGNL+jCOsH1Gz2uimVjsL+laW7i0ZTHR8y0EgGJ4hI1HFH23KOPOk9ekD+EIrzv5JwCAB/finUM0BhsjPT/hlye6RP69N+7Hl1yMfJLeAFyVPDqh5w9CANg7q+d4Vp4zuPp8Xb8mJut6jme8rSvlNZPp8Sw/o6Npclg9T5k6tk/r+UPHMfinK0Z8unPm8G0NcrhyVtdffL4UvvvKfXIHMsNJoSbOjnWdus+XwtesPSp4JFPee+P+g7c1ySGpZuW5w19YGuRw+eahhGmSQwCYnFmTPoRjaJDDW5cOd9nQIIeTVV0S5kthjL0XY9XJCp4qOXbOHJVCaVbObh2TwivXN3Dlupx0lCWFknJYJoWSqeHZ8cYxKdwYrYglh5cnO6VJoaQc+lLokJTDp549dWw/w9GJPVXJ4d7ZPdHk0JdCh6Qc+lLomKyPxARxvD05lhZOzqyJCqJLC30k5dCXQoekHGqTwjKattfJxcKKYZ0QakgNpbnz5LXa+lhCDuuSQm2VMpC/VtZWHQPlUujQmBxqkkMpyqTQoSE5LJJbDrXVx0C5FDok5LBMCh12pCM9lKauQpaWw0E/PX1j0zYpYW45bKqPJVNDDbSpj3PLobb6uIncqWGdFDpyy2GbM59oIndqWCeFjtxyWJYWaiZ3algnhQ4NtbIk2tLCNusK+8phjDp50GLYhy7VcQ45LKuPq8ghh01Joc+7r9yXJTnssqYwhxyW1cdV5EgNq+rjMl6z9mgWQWwjhY4cclhWH1ehLTXMVSu3kUJHDjlcvjlpLYW5auW2aWGuWrmNFDpyyOGtS6u1aaFPruRwiFLokEoOF0oMNa0nBPoNmSxacqhx0KQrmiaVc9BFCh3aamVtcpiaLlLo0Fgrp0RjhdyVRUsOtUlhHyTkcGF+svtKYarUMGTyOJUcat2SpispU8O+9XEqOey7rjBVathHClPTtz7WJoeatrFxpJLDvvVxqtSwrxSmTA27pIU+2uQwVWqoUQr7bk2TWw7nXgxjTB7HlkON29GESGGKSjk0KYwth13q4ypiymGX+riK2HIYKoWxU8Mu9XEVGuUwtiD2SQt9Ysth6JrC2LVyaFKYolbuK4WO2HLYpUIuYxEGUkL3K8w5sTznT0U8YslhLCnUWCnHksNY9TEnlZuJJYexkkJtlTKgTw5jEiqFjlhyqG3QRGN9HCqFjlhyGCKEqdCWFmrcxLqOuRXDFHsUhsph7KQwdI/D0KSwjFA5jL2m8MG95SBBjJEUxiRGUlgkVA5j18dfcvFTQYIYIyksok0OQ5PDledMNCl0hMphbCkMTQ5jS2Focmgmk2hS6AiVw9hSGJocTlbN3EthjuRwLsUw5ZBJXznUWB+noq8czsOgSRv6poYp9ynsK4ca1xSmQpscaqSvHKZMCrWdJUUbfeUwZVI4L7VyyqQwpRzOyZdfN6mlUGOt3BWNUpgyKdRWKQPd5TC1FPZJDVPvU6hNDvukhrGTwiJd5VBbfQykrZB5+rw0aEsKh8xciWGuU9x1SQ1zJYVt5TBFfVxGl9QwhxR2qZRz1cdt5TBFfVxFWznMlRS2rZVT1MdVaJTDNoKYoj6uoq0c5pLCLrVyjnWFXWrlFBVyGV3kMMe6wi61skYpzLGuMFWtPDdimHuPwjZyqK0+zk0bOcydFGobSBlicrhI9XEV2uRQI01yuGhJYR9yCKFPGznUOGyijdzDJrHlcG7EUII6OZSQwqbUUGKfwjo5XLT6uIo6OdR2/mMpKaxLDaVOc6dNDutSw1xJYRFtm2DXpYYSUji0WllCCptSQ21p4dAmkMvQ9VPbg1z1cRVlciiZFJZNKueqj6sok0NJKSyrlaWnj8vkUFIKc506rwtFOcxZH1ehUQ6LgiglhVV0OdVdCsrkUDIpLKuVc9XHVZTJoWRSWFYrL8IEchdi1sqDFkMtp7jz5XDR6+MqfDnUkhRqrpW1JIW+HGqokLnPYTc0SKGfGmqpj305XPT6uApfDlkfN6MlKYwhh4MWQ01sPLKkSgpdaqjpNHfvvnKfGin00bRPocY1h4AOKXRolEONaJBCh7ZK2aFJCrXWypqk0KWGTArTovOndYBsnd/H3mN6BGM0srh2cx0fvXyH9KEccP+JZ/ALz75G+jCO8PjE4CM7m9KHccB7tnfx4V1dLzJXJxt4wfpl6cM44B2XH5A+hOMoeyXdOz3Brbv1JJkbn7XYOanr+3rp1i7MZF/6MA5YfvwKRjf1vBYBwNZdp7F6Rc/30fbpEbZuG2F3Xc/30mTFYElPJhQFZS9n3RgreTK2zh++uGiSQ03cf+KZg7e1yOEJo+cFD5hKoeOpfR1/pf/GzRcevK1BDjVK4f7mbCmJwldTDXK48dnDauvmxbjnne/L6jPbB29rkkMAMFs6lpFs3XX64G1NcqiJycqhoGqRw72NcGlW+FLWjfGWrCD6UuiQlMPRyGI0OrrG4KOX7xBNDn0pdEjK4Qmzd0wKP7KzKZoc+lLokJTD37j5wiNSqIGiFF44dx0Xzl0XOpopB1LokH5FnZjpPw8Ncuhz8+KSqCD6UuiQlMPlx69g+fErRz5mtnZEBdGXQoekHG6fHmH79NEfrt11I5oc+lLokJTDvQ0TRQoB+ZexQVMmhQ4JOSwKYREJOSyTQoeW5NBHQg7LpNChJTl0SKWGdUmhlBwek0KH1KvqpPqXgpQc+mmhBsqk0CEhh0Uh1ECZFDok5LAohEUk5LBMCh1aksMQ5kYMc6eGdVLoWPRauU4KHbnlUHN9rIWmpFBDpSxNpRQ65uaVtT9NUqilVvZZ9Fq5Tgodi14r10mhI7ccxkoKHcZaXX/RdcC+9I1vLr1gkni4q40U+izddSvRkRzSlBYWefHFzyY6kiltpNDn6869O9GRTOkjhJ+/sp7gSA7pKoUXRtVpRyy61McPbl5MeCRTuq4pzLGvYaMUHrlyuuM4oCYpLGPjM+mFrEtSeOJyetGoSwrLsOO0Zt8nKbRraXcsaCOFPttn034fNSWFRZY307tMGyn02UvsIlVC+OEf+XYA6G2L/Ls2A6mTw65SqJFFq5WHmBQWSZ0c9hk0SV0rd5JCIP0rbEcpBNLXyl3r40VLDodWHw+F1JVyVykcMnMphilr5a5poUNbrZxyvWHXtNChUQ5T0FcKU6437Dtoski1cmcpdMzlq2w5fdcUppTDrmnhotFXClNWyl3TwtT0lcKUlXLs+thH11c/IimmlftKoSO2HJZNIHchxbRyXyl0xJbDsgnkLkhPKxd5an81uiCGTh/HlsN3XH4gaFuaFNPKvaXQEfuVtmT6uAu37t6LnhyGDprElsPVZ7aDpDB2alg2fdyFFJPKoUlhbDksmz7uQopJ5dCkMLYcxpw+rmJuxTA2oVLoiCWHGuvjUCl0xJJDbYMmgL4KOdaWNLHkUPU+haEofLWNJYexpo+11crahlGAeAMpsepjjcMoseQwVn08tEllhS9VcYmRHMaSQofGWjkkObz/xDPRpNARKoexpTBGchhTCmMkh7H3KQyVw9hSGCM5jCaFjtBX3MCksIxQOYy9JU2oHIYmhUVC5TA0KSwjVA5jrylcvbIXJIihSWEKYq8pDJXDHEmhQ9czkZC+chhbCh0hcpgqLdR0+jygvxymTAr7ymGqpLCvHKbavFrjmsO+chhdCh19X3UjC2EMUu1T2FcOU60n7CuHKQdN+sqhtkGTVEIYkhqmGjTpK4e5hNCxMGII6DmFnqOPHGqskFMyDwMp2urj1PSRQ20VcjIpdCh75e2TGmrbvDo1GmvlrmiTwtT0kUNt08e5pRCY030Mm2i7z2GqtLCMNnsd5pTCNvscxq6P62i7z2HOdYVt9jnMKYVt9jnMeZq7tvsc5pTCNvscJpfCIw/W4joZk8K2exzmlMI2+xzmnDxus8dhzi1p2u5vmFMK2+xxmLM6brPHYU4hbLu/YV8pXOh9DGeffGfaJIc5pRBoTg8XLSks0iY5zD1somlauQ25z33cJjnMnRRKn1u5M5nr4zbJYe6kUNtAShO59ylsUykvWlJYpCk5zJ0StqmUA6UwiEGLYQh1cphbCjXSNJCSMy101MkhJ5D1nVcZqJdDbfUxkDktBFS+AtfJocb6OPc+hRor5To5lJDCukEUjYMmEtTJoUR97LPQz06KvQ5D2Hts41hyGLpXYShlcighhY4yOZSUwrJp5fds74qtK6yaVs6dFvqUyaGkFJZNK+9vLuWXQscIx1+JE0wfd6FMDiWlsCw1jD193AUz2T8miCmmj7tQJoeSSWGZHEoKYVlqOFkxqtYU5pw8rmPwYhgjNvXlUENaqHE7G4ekFDo0DqRoq5V9OZSUQocvhxqTQnIcXw41JIW+HGo7mwlPc1eOL4caUkJfDjUIoZ8axhDCGD4EzIEYxkSDFGpF61Y22ipkTRPI2qrlF6xfViWFLjUUSwqLuFdjRdvS3Lp7T4UU+miSQm21sksNNUihjwYp1IrGza+VvCLqYO3pkRo53L08nXhdvdQ8rZyLTz57XkVi6Pjla68CAHzDbe8RPpIp/+zqPQCAV69/SvhIpvzU018GAHjZiUeFj+SQ5524gSdvnpQ+DADtJpSzM7bTfzs6fpGuPbkEQI8YLt2aYLK+hPGmjj8Ixw/Ndm/YaN6hIAu7ezC7ewD0iOHOSR3fywCwtGVh9fzdBQDYPaHsgDAniWGs+BSYyqEmtp/QUSuvLE9fiN/6qHwtCQDr48Nk7ueufpHgkejESSEAfOjmPYJHcsivPfUFAKZyqInRug7JgNEjYICTQuDGvTp+cZ156DApnKwz0zjG7uH38frHmrcby8H1e6Zb6RgdeYs6YkphTA8a9D6G/jt99zQEqvc1lEoP7VL5cyKVHjopLPIV93w885FM8aXQRyo5/N1bn1f6ccnk0BdDh1Ry6ISwiFRyWJUUilbKVVIolBw6KSxy8hGZ3xe+FDokU8ODpLCAkUoOd8u/FpsvklkC5ISwiBXKXZa2yr9v95dl/uipEsJJwOqfEjFczH0MY1C32bW29FAbWtJDh0RyWCWFAPDeTZn1dGVSCOhJDh1MDmcoTQq1UCaFgFxqWCWFAGBvCQyhVUihFFVSSKbUpYRjJctnaT7KqEoLgWmtnLtarkoLHbnlsCotdCx6rVwlhY7ccliVFjpyy2HTusLsctgkhSt5W4smKcxdK1dJoWPhK+UGKcxdKTdJYe5KeWnLVqaFADDa1fVHmRbmpkoGutfJbU+N50hdLddJYZEctXKTFPrkqJWbpNAnR61clxaWkbpabpJCnxy1cpMU+uSolbsMm2SplbskhRkq5S5JYY5KuUkKfVLXynUpYRlZKuUOSWGOSrlLUpijUq4TwiKpK+Wuawm7VsoV6wtZJS8aqZPDLlIIpE8Ou0ghkD457CqFqekihcDi1cpdJ5CTJ4dd6+PEyeFQ6uMqtCWHySvljvVx6uSwa32cOjnsIoXkOAsrhl3TQiDtmsMuaaFDy8Syg2sO65Fac1hFSjnskhY6FmbN4RysKUxZKXeVwtR0TQuTwzWFtfSRwpSVcp/JY+m1hnNVJQPt6uQ+UlgkZq3cRwqLxKyWu6aFZcSslrumhWXErJZD08LYlXLXtLBIzFq5jxAWiVkrx9irMGqtHEMKI9bKMZLCmLVyqBTGrJRjCGHUSjmCEMaslGMIYcxKOUZKGLNSjrEVTZtKuWabGlbJudE2sRwrPYwhhUC89DCGFMYkRoUcMzkMlUJgfmvlWBtYR0sOYyWFkWrlodfHZSxcpdyRWJVyrJQwVqWsrTrWuGl1F3TZTQRibvLYRAw5jJEWkmZYK9cTQw5jpIUObbXyvBFTCmPUyqyPG1BWH88rMSrlnFKYynfmrkoG6uvkGDVyGX2r5VRi2LdajpUY+vStlVOlhX1q5VTDJn1r5RhJYRl9a+WYUujTt1ZOdbq7XrVyqjWFPSvlVElh30o5lRT2rZVTSWGvWjmREPatlFOsJ+xbJ6dKCfvWyamEsK5ObhBDVsnS9EkPtaWFKaQQ4FBKHfOQHKaSQqBfcpjyHMida+WUgyY9KuWU9XGf5JBJYQ3KUsJUQybaTpfXJzUcenVcZC4TQ6A8NUyVFhZpkx7mlMI26WEqKSzSNj3MtbawTXqYa2uatulhqrTQp21ymFIKfdomhyml0NEpNcwxgdwyOcy1prBtcphDCtumhrmEsFVqmEkI26aGuaaO2ySHOdcStkkOcwlhWWrYokZmYlik+EXLJYVDJJcUtkXTwEnO/QrbpIc5pBAY5kBKDikEOqSGubalaZEczuOgSRvaDKPkTAk1DaPkPitKE9qSQ00Ut69JPUsxt2IoydrTI1VTy5r2O3zroy+srZZzS6GmWhmol8NcUuhoksNcaaGjTg5zSaGjUQ5z71VYI4e5pbCpUs5dH9fJoarqGMheH9fJ4fV7VlTtUahp8nj3hJm7+thHj70kIOeEchlVciixtlDiPMt1aFp3WCWHUmc30bTusEoOc0uhQ9O0cqUcKtrAWioprJJDTWsKpaSwMjVUtKZQmxBKSGHVWkNpIczhNXMthhrQnh5K1sja5FBTeliUw9xpoU9RDqWk0FGUw9xpoc8xOZSUwkJqKF0fF+VQUgqLqaGqpHB3j1I4o1gnMyWUQY+xJETD+kInhxomkbUlh04QNawtdHKo4VzITg4lpdChdc2hpBQ6DuRQQ1I4k0NpKSyiISl0cqhBCg9SQwVC6OpkVsfHcamhBiHMeZq8uRfDD//It2O8JX0UU9aeHmHlqo4v+fYTG9h+Uo8gvuOynvr0Bz7z1dKHcMC3f/zrpA/hgA/dvEc8LfTRIIUHaJBCR6Qzo8Tgxr1GhRQ6lh5+UvoQDth76GHpQzhgvKvne2Z5U6Y6LsNYHVIITCeTcy2P02EpiZFea+iYrE2/2VXI4Xj6v+uPnZY9Do+3Pq6nWn77s/KJoeODz9wpfQgAgN9/7H58+spZ6cM4YGVVPm3RyPianrTw1EMWz93b4oSvGVj92OPSh3DA5LIeQd353OdJH8IBS9vT35Ea/s5yx7B2RYc05/QYBYaSh/EW1CSHgBI5nEE5LIdyWI4GOXzwyQsAdMjh/tZY+hAOcFK4e0r+l9mphxT8dp9BKSxHoxRqQIOYOiar9Wc/ScHcbnBdpLjhtcS6Q5cY+uzcJvQCXvG77NRdz+U9DgDrK+VrC7/izn6n0gvhEzfKXyi//Nwnsh7HWx55deVlL7897y+533/s/tKP33f2StbjAA6FsMjOtkxKplEKiyxfz/9HaJUUnn4kb61cJ4R2P+9rb5UQLt19V9bjAKqF8OoD+X8xVgmhFWpwy6Rw62z+n6GiDPZIDLnBdRPFL6qW9HDlqp51h0D+9LBKCiWokkJtaEkPNSSHDonkcAhSKIGWpFBTSljH3mcekz6EA277VN5fjNpSQi1JYQQpDGJhEkOH1KnyytLCMrIkiC1+n+VIDrtIYer0sK0U5koO6xJDR47ksCotLJIjPaxKC4ukTg+HKIQ5UsO2QpgjNWwjhbkSw7bVcerksEttnCezFGkAACAASURBVCM5bCOFOVLDtjKYIzXseeq7KpgYtqXsi6wlPQQyrD1s+Tvt+mOnVa091EKOdYdtpBBInxy2lUIgfXrYVgpJOanXHHZJCVMPo7RNCs0owy96ResJtbC0bVsnhakTPC0JIRBdCoNYODGsYqHksANa5HBRhlLaSqFDS60M6KmWNQyk5EBLfaylOga618cp5VCTFGoZMtFWHWsh93BJEwtXJTvKKmVH7Gq5bY1cRvRqOaAFi10v911fGLtWDllbGLNa7iqFPrFr5S5pYZHYtXLftDB2pTzECrlI7Eo5RApjV8ohawpj1sohQhi7Tu4rhLGr5BAhjF0nhwhh7Dq5TggjpIWskrtS90XXtLXNvKaHIUMnMdPDoQycNDGvyWFIhRwzOZwHKQTiVspDTgq1EnMQJSQljDmEwpTwOE1b0EjvvazHOhRCOdSLlmpZ016HsQhJCx0x5DDGusJ5q5Vj1Mcx5DCGFGrZ+DomWupjLdWxJjRJoXYWtkp21FXKPn3r5ZAauYygajly6NG3Wo69RU1ItRwzMQyplUNq5DL6VssxpNAnpFaOOXASUitrSQtjrikMqZRjJ4V9K+XYKWFInRxTCEPq5NhC2LdSjp0ShtTJsYWwb53cVggjpoWskvvS9knQlB4yQTzK0NPD2FKoib7JoZYpZC1SGJu+qaGm+jg2fQdRmBLqRPO+hFVIV8gOGkYHtMihFvpsaTPvG1prqZY1rTnsSgopHHqlnGICWcNp8/qiZU1hCikc8obXXbaiWRSGUB0XWfgq2dG2Una0rZZjV8lltKqXEwcfbWvlHGLYplrOMXTSplrOkRa2qZVjV8hltKmVcySFbWplLUlhji1p2tTKOZLCNpVyDiFsUynnSAnbVMo5UsK2dXJqIWxTJ+dICNtWyV2FMEFayCo5lK5PiqbJZQ1wQ+zjaEkPtaBln0PSjXmuj/ugpTrWAlPCozRNHJehpUJ2MDEs0DU5dFQliDkSQ0dtcpg5AKlKEHNWyVXJYe4taqqSw9xrC6uSwxxpoU9Vcph7XWFVcqghLcy9eXVVaphbCqtSw9zVcVVqmFsKq1LDnOsJqxLD3DJYlxjmXEtYlRj2rYwTSiETw1j0fZI0pIeaBlM0pIdDH0pZJDhsIsuQ1xvmQktSqGHIREtCqGXARKEUBqHDIuYEDXIIlOx7qOR3ncTgiVY5lJhELhtIyZ0WAnor5UWVwiokKuSyvQ0lBk2KE8qLLIUxN7sOoSiAGoQQGOZwSROskmvoWys7Jmt5q+Qydm7bFxdDVytLTyS7aln6bCdffu4T4lvUvPz2x0Wk0MdVytJp4c72khoplD7/sauUpdcVnn5kW3zy2O7viwuhq5KlU0JXJ0snhdbIC+HW2VGwDGZKClkla0RDgqihWtYymPLWx18oLoWAjn0LNWxnoyU51LKVjbQUOqSlENCxHY20FALTrWukpRDQM2AiLYXAfCaEReStQTExrH68FfkM4D0Y3dLxND/55BnpQ8CHP3tJ+hDw9LWT0oeAxx++Hdcvyx/HH3/8HulDwNbVnqc1io2CX3prT0kfAXD24zex9znyMqSB8UX5r8NkVUeavvKc/DrYyUr473Ot6wp9WCW3IKRS9qeVpWrlvfWj7+9vCPyAnTiayjzvedfyHwOA524dPiEvveMJkWP44GOHk4bnz9wQOYbHH7794O1TF2WO4fqTh2K6dlYmXj8mhQKRxPhqISkU+lvy1KcO/4BcviHzWnX24zcP3l56WCaxmzz97MHbdndH5Bh8Idx9Uf9T5IVQFMLN8/kT7aIM7pyWCTl8Idw61/8HNLMUskpOSawnc7xlmCDOkEgPfSkEmB46NCSHW1fyp3ZqksIiAk7mS6EUvhRK4UuhFMWUcPlj+c+EoiEl1JIQxkgJgWEkhQ4mhh3okxzWnSElV4JYTAx9sqWHJ6rXceVKD4ti6JMrPfTTwiK50kM/LSySKz3000KfnMlhpRhmTA2PpYUHx5DtECqlMFdqWCeEuVLDOiHMlRrW1ca5UsM6IcyVGNYJYa7EsE4G+ySGQlLIxDAHsZ9cpodTNKw9XBTqpDAXVVII5EsOa9PCNufeikClFALZ/uytSwp3T8q/PnGtYT6YEk6JlRA6hpQUOuStYGBQDtOQWg7r0kJgWiunrpbr0kKAtbIjtRy2qpATy2GtFB4cQ9JDYH08o6k+NssrSR9/fPF5jUMmqetkSuEUSuEU+VeGAdLlyW6zZY1beygpiaNbI3FBfPLJM+LpofS6w6evnUwmiG3TwpRyWJcW+kisOTxGIjlsJYWJaSuFKVPDtlKYKjWcPP2s+JpC6anjyepYXApXntsXlUK3jrCNFHapkYcqhQDFsDepnnQKony1nCM9bEI6PZzX5FDtwEkVCVJDJoXyQyZtUsLUUAjjDZYUGbIUAhTDIIb+5GtGWg4BHemhJPMqh52InBpKp4WUQh1SKI20FM4z8+AFnEqOQJtp5brp5Obbhn+qdZPJTUSZXK6ZSm5D6ORy0xrDJmJMLTetMawjxsRyyOBJjGnltjVyGaHTysFJYYRJ5WApDHTUECmMMaEcIoQxppNDhDDGZHKoEIZOJofKYIyp5JCEMMZEckhC2KZGViaFnEqWJPU3g/SAinS1PA+ESCEQnh6GTiOHpochUkjkCV1rGJoShq4zDE0JUw+gLALSwyWpamOHMikMgr/xI7EIcigpiNKDKdK1sgYkq+WhV8qLXCFLV8cakKyPpQdMpNcSApTCrrBKTkBVtRxSJx+/r26ffkiVXEbnejmwSvbpUyuHVsk+fWrl0MTQp2utHHvvwq61csy0sE+lHH3gpEOtHF0Ie/x+iymFXSvl2FLYtVKOuZ6wT50cUwj7VMkxhbBrlRxbBrtWyTFlsKpGVi6ErJIXja4J4tJm3MdneijHIg+ldE0OBzeF3ETHP4cXOSlc5CETDSlhTCSlcBFhYpiQYnIYMzEs0iZBjJ0aOlqlhxETwyJtEsSYiaFPm/QwZlpYpE16mOpsJ22Sw5RrC5vSw6RS2CI1TFofN/zeSymEbVLDVFLYJjFMKYRtUsNUQtgmMUwpg20Sw1SVcRsxTCmDxcRQeVLoYGKokeI3T5vNrvvSJkGMnRo6pIdTmB7KoGE7GzEa1hxKrymURDIpXPSUUBJK4fxAMZwjVG+MfTPtL8pFlkNJ6uQw9SSy+ECKUlLXx3UTytL1sSTzLIV1aaH0cAlr4/iwSs6Iq5ZTVso+ZfVyqjq5SGm9nLBO9imrllNVyUXKquWUVbJPWa2cqkYuUlYr59qipqxSzra2sKRSzpYWlvw+zLWmsKxOziWFZXVyrqSwrErOJYRlVXKulLBMDHPJYFlimEsGXVo40JSQVfKQSFkpH30cudPrSQ+nSCGZHi5qrVxMDudu4KSKgpst6qCJZH08zylhHVIJYcrT2BXpcl7keYNimBGpvzok5fCIICauk30oh/lZyDWHhbWGi7i2kFKYn5xS6KeFkrWxVGU80LQwCFbJQrz0jW/OVin7TNZstjq5yP7GfrY6ucjayW2Rx33pHU9kq5J9zp+5ka1GLnLq4g2xM52snd2SSQyNlZFCI5cULt+wYlK49PCTIlJod3fEhHDr5feKPO7m+SUxGdw5PRIRwq1zZh6EkFXy0JBMD1evLF69LIVUerioyeH2ExsijyuVFK5cZX2cE8mUUAopKdw8P2ZKKMTi/aZWxEd/SO6bL9c6Ry1s3VjFzvay9GHkRWiJzOjmGKMbsltnkLSceWgX2+dWRR576dOXRR5XCrOxLpYWjm/JNDy7J+TU5A//n+8Qe2wtsEpWwou/t/w0eqnYWz/88uWstPcvdD+tVExGy4d//a6s7mZ97Mkk34vdzjXvl3bmn5TRzUMp3D85yfrYZuvwa2zX8yUd42t5RXjl2uHnuZo5PDvz0OHPzeqzeZdo+FI4eeZKtscdnbst22MBUxl0bL7gQtbHLsrgzm0r2R7bF8K99bx/2c6hELJKHjqLkh6Onsr3ItPEwiSIGV9ffSkEwORwzvClEEDW1HBRkkJfCnOjRQpzM4dSGAQTQ6WkThD9xNAnV3oomRz6qaFPjgQxR2p4JC30yfATUxRDR47k0E8Li6ROD3Mkhn5KWCRHaliUwsPHTp8aVklhjtQwR1pYJYM50sKqujiXFFYJYerEcAFksPcXcPH2ViC1uPRQYmJaGpcgphTE8Xg/a6WshdGNcfZamcSjSgrnndRSKJkOAnJrCAHZhJDUw8RQOamSw6rE0CelHGpMDH1SymFKMaxMCx2Jf2qqEkNHSjmsSwyBtKlh6sSwLi0E0iaGTVKYOjFsqpBTpoYpxbCNFKZKDNsIYcrEsI0UpkoMFyApdPT+AlIMB0IKQWwjh0A6QdQuh44UkphCDhul0CfBT0+TFPqkEMQmMXSkEMRUYtgkhI4UYtglJUwhh13WFKaQwxRS2CUhTCGFbRPCVFLYNiVMIYULJIQODp/MO9LDKfO2vc3+bvtv/RRDKuOx3EnnNRB7KKWtFAKA2RzGy15bKQSA7XNxH3tRq2PHIkthCnZPjDhcMiCG8QpJAMjKIUA5nKsp5sU9DSgA/XLYRQrnkXmaQDYb66JSOL61Jy6FklAKu8MqeaDEqpbb1sllxKiYNe1r2JUYFXOsSrlTjVwkwk9Slxq5SIxauUtiWCRGtRyrSg4Rwhh1ckhSGKNODhHCWHVyjLQwZKgkhhiGiGCMGjlEBmPVyBRCVskLh3R6CMxfgtiVGAlijEo5SAojECKFgPxeh1rSQ+mUMLQ+ljoTimN8+9ng+5CWwhhIpoOAfEIIUApDYWI4B4SkhyGJoSMkOZRODIGw1NAnJEEMSQ6jiGHAt0GoGDpCksOQxNARkhyGJoYxpDAkMYy1pjAkNYxRH4ekhiFSGEsGQ9LCWEIYkhjGkMKQxJBCeAROJZN+ghhDDH36SOI8ySHQXxD7yGH0tLDHt0MsMQT6y2EMMXT0EcS+YhgzJewrhjEHTfqIYez1hH3ksK8UxkwH+0hh7HSwjxTGTgj7iCGFsBRWyaRfvby0GXcKoU+9rOk0eTFwgypzNaySCelaGZhWy6nr5ZVrI/HqGOD0cR/cMMm8VcaUQuJgYjjHtE0QY6eGPm0TxHlLDctomyS2SQ6Tryts+S0RMy0s0jY9jJkWltEmQWybGKaWwbapYUohbJsapp48bpMctk0KU0tgm7Qw9drBtmKYcg1hWzGkDLaCiSE5DgdUutFl+5o+tE0Rh7LHYUopBHSkh0C84RQNCSGQPiWUHkJpSxspzJEMapDCtnCwZDGQf5ZJUrTIYZMgaqmTU8shkGbD7OgYqNjrcHRjrEIQQ+VQgxSeeWhXRXW89OnLg9inMEdV3CSFufYgbEoLpTeodlAK88AqecEoq5dTVslVlFXMGupkIH2lXEZZzVxVKWffnqbk2yN1WlhGVbWcukouUlYtV1XJuYWwqkqWEMKySjm3EFbVyWVpocSawTIxlEgHy8RQQgTLqmTKYG84lUy64QuihBgClMMyioJYlEOxPQsL3yISYggcl8PcUugoymGZGEqlhEU5lEoJi2IolRIW5bAohVJDJEUplKqLtUghcFQMKYTBcI0h6cZQKuZFo24douhG1gpqZUDXusOqelnLxDGgZ+pYsjqu2vhaw2SxQ9MaQlbGhIkhwYu/981iqWGR3Xt0JIYOqeSwyHhpIn6GE4cWOds/ORFLDIusPaXjawJME0MtQnjifY9IHwKAaWo4vvuS9GEAmCaFWkQQAG7eFeHcphHYWzcUwriwSibhvOCH/7H0IQAA9k5a2LM6frGZpX0YJWnZ5OaS9CEAAMzOCGZH/otiJtNj2F+Vl3ctYrh8ffr/s5+Q//lZubaL5U8+IX0YU9Z1yM/k3ElM1nUMn23esYr9sfzPMQC86+e+U/oQ5hFWySScB79Hz19r5oqOF08AGO7fTvExO/peMkbb+o5JAieFmtj9XB0pnQYm505KH8IBm3foaB8ASqFGmBiSSiQTxL2Tx59eyRTRLB1NpSRTRMnksEwMpdJDlxj6SKWHkolhmRBKJYYr144/rmhqKJgUFkVQMiksE0HJtJAymAVWySQdUoJYJoeAjCAWxfDg40KvrVJyqF0MARk5lBLDqpSQYggxKaxKBiXEsCoZlJJCCmFWKIYkDzklsUoMfXJKYpUcHrlOxtfbnHLYtkLOKYlVYujIKYi5xbCpNs4phmUyWEZWQcwohW0q4pxS2KYmzimGlEExKIYkLzkEsY0YOnIJYhs5BOZPELWJYZMU+uQQxFxi2GUdYS45VCWGyoQQyCeFXdYN5hBDCqE4HD4heckxqLJ0o/33tbmynGVgxe61+5Gxdn6GVroMnNgVfZ/0vAynaBwuaSuFwPwMokzOnVQlhZt3rFIKSVSYGJJopEoRuySHPqlSxLapYeXtE70up0oO+04ip0wPuySGRVIkiKkSwxAZTJkYdhFCn2SpYaKkMHSSOJUY9p0qTiWFFEGVsEomOkghh33F0Ce2JIbKIRBfEFOIYYztaWILYogUOmLLYQoxjJEQxpbDvkLok0QOI4thjK1lYkthjC1mUoghpVAtFEOii9iCGEMOgbiCGEMOgbiCGFsO51UMgbhyGFsMY9XGCyGGEaUw1l6DMaUw1p6DsaWQQqgeiiHRSwxJjCWGRUJFMZYcHrvfwNfwUEFMtZF1qCTGksIioZIYQwxTrSEMlcMYMlhGsCAGCmGqDadDpTDV5tMxxJAyOCgohmQYhEhiKjkEwgQxlRwCYYIYIocpz3ASIoepxBAIk8NQMUw5WBIihqmkEAgUwwApTHkGkhApTHk2khAppAwOFk4lk2Gg6bR7PiFTzW0nlXvdd8B08/jEXtyDiYRdsWqnl3NPMC9f1zttnFIKg+gphV2miXvdf08p7DpVnBNK4WLCxJCI0yVFTJkaltElSUyZHJY+Xoe/B7ukh7nPh9wlQUyZGJbRNkXskhjmFsEuiWFuGeyUGnYQwtznJe4ihbklsEtaSBGcK1glk+HTVhBzyyHQThBzi+HB47b88W8jh7ml8PBxmz+J3FLoaCOHbcVQKh1sI4dSCWErOWwphbmF8OBxW4ihRCrYVgophHMJxZDMH1WiKCGGZVTJopQgHjmGmpeEKkGUksKjx1B94FJi6FMliXViqKEqrhNDDZVxpRzWCKGUBB45hgoh1FINV4khRXAhoBiS+aVMELXIIXBcEDWIoU+ZJBblUIMU+hQFUYMU+hQFsUwMNQihT1EONQihzzE5LJFCDTLoUxRDLUIIlEshhXChoBiS+ccXRE1i6PAFUZscOnxJ9OVQmxgCR+VQmxg6nCD6YqhNCB2+GGqTQqAghp4UapNBhy+FmoTQ4YshhXAhoRiSxeP+f/Kj0ofQzIVt6SNoxF5dkT6ERkab+sTV58RjOsXV5+J7b0kfQiOjLX3CWuTG/TpF1ecPfvGN0odA5KEYksVGqyQaLzi0F/VKomY5XHruUAr305x6NhjNYrh88/Cl8tyfbAoeST2apfDaC08dvD3e0fmrhzJIClAMCQF0CqIptMpaBVGrHPpiCOiUQ61i6EshoFcMNUqhL4MOjVJIISQVUAwJKUOTKBYFsYgmYdQiiUUpLKJFEjWJYVEGi2iRQ00yWCaBPpqEkCJIWkIxJKQJDZLYJIcOaUnUIIZNUugjLYgaxLBJCH00yKG0GDbJoEODFFIGSQ8ohoT0QUIW28phESlZlJLELmLoIyGJUmLYRQZ9pMRQSgbbSqCPlBBSAkkkKIaEhJBbEPvKoSO3JErIYV8xdOQUxNxi2FcIHRJimFsK+8igT24xpBCSyFAMCYlFLkkMlUOfXKKYSxBDpdCRSw5zimGoFDpyyWEuIQwVQZ9cUkgZJAmhGBKSklSyGFMOy0gljKkEMZYQVpFKFFOKYSwRrCKVIKYSwpgCWEYqKaQEksxQDAnJRQpJTC2IQBpJjC2IqcXQEVsQU4hhaiF0xBbDFEKYWgaBNEJIGSSCUAwJkSCmJOaQQ59YohhTDnOJoSOWIMYUw1xC6IgphrGkMIcI+sSUQsogUQLFkBBN9BXG3HJYRV9p7CuJuYWwjj6y2FcMc0tgHX0Fsa8M5pa/KvpKIQWQKIdiSIhWukqiFjn06SqKXQVRkxg6ughiVzHUJISOrmLYVQi1iKBPVymkDJIBQTEkZCi0EUWNcujTRhTbyqFGKfRpI4htxVCjEPq0lcM2UqhRBH3aSCFFkAwYiiEhQ6dMGLULYh1FeSyKonYhrKMoi2ViqF0C6ygKYlEEtUtfHWVCSAEkcwjFkJB54/5/8qODFsMi9uL2ETkcshj67C8fiuGQZdDHF8PR1u6gRbDIeMdSBMkiQDEkZNF44M3y536OwWhb/jzDoZx+SPoI4jAvYvvun/1O6UMgRBqKISGLztBFcciCOHQxHLoQUgQJOQbFkBBylKGK4hAFcahiOFQhpAgS0gjFkBDSjyEIpHZZHIIYDkECKXyERINiSAiJh1ZZ1CqImsVQqxBSAglJCsWQEJIeTcKoSRK1iaEmGaQAEiICxZAQog8JkZQQRikxlBBAih4hg4BiSAgZHjnEMYco5hLDHCJI8SNkLqAYEkLml9gCGVsWU4hhbAmk8BGyUFAMCSGEEEIIgAAxXIp5FJnRs/KcEEIIIWQOmI+TlRJCCCGEkGAohoQQQgghBADFkBBCCCGEzKAYEkIIIYQQABRDQgghhBAyg2JICCGEEEIAUAwJIYQQQsgMiiEhhBBCCAFAMSSEEEIIITMohoQQQgghBADFkBBCCCGEzKAYEkIIIYQQABRDQgghhBAyg2JICCGEEEIAUAwJIYQQQsgMiiEhhBBCCAFAMSSEEEIIITMohoQQQgghBACwJH0AhEhy3lyyO9ievmMMjLvAHPznEOO9X3p54WOmcGHt+033O8VW3a7x8QBbctfltwVgzPT6ddc5dmw116u6bov79en0ORQ+3urzqXvsro9bepmtvk3j/dtjHy+9iXH/O3p9U3Is/ve7KVxWvI1/f/63qfu4f/vp5WUfL96vrb38+H0fXt8cuf+j1zt6u+P3Wbyef93y+y9cduRt//i96xnAe0U5dv3jl00/8kcf2v4Na+1XgRAhKIZkodnBNl4z+kqYkQHM6OD/GJmpKI5moboxwGg0+7+B8a5z7HL378jl3vVLLrfGTPN773Lrbut93L+ePfxNN33/4OPu+jhy33b2MXf5gegYACN3u/LLj3zMYPpYVZfN/j997KMfO3IMB8dedr2Ky3F432W3a3PfxWM78nbh9ocft9W39y9H8bbWuw4OZKl4ufGve/A5usuO/h84vL4xdvat4V9mZ98Sh5ePZm8fXO59bITD993lB5cdXG4PLh+5j88u898/uPzg4/sH7x/edh/jg+tMLx8fXD5731iM4C7bP7jvsdnHCHb6f7OP8ew+/NuPvdu56x/edvrY0+vsz64zu+3s+mMcPsYY0+Nwtx/PPq8x/Ldn9w+LkZn9H8DYAGOYg1pubAxGMLOPmdn7o9ljG4zNCCOMML704HkQIgirZEIIIYQQAoBiSAghhBBCZlAMCSGEEEIIAIohIYQQQgiZQTEkhBBCCCEAKIaEEEIIIWQGxZAQQgghhACgGBJCCCGEkBkUQ0IIIYQQAoBiSAghhBBCZhhrj51BlJCFwRjzYQBb0sdBCCEz1qy1L5U+CLK48FzJZNHZstZ+ofRBEEIIABhj/lD6GMhiwyqZEEIIIYQAoBgSQgghhJAZFEOy6Py49AEQQogHX5OIKBw+IYQQQgghAJgYEkIIIYSQGRRDQgghhBACgGJIBogx5kXGmHcaY7aNMW8suXxsjHm/MeZXKm7/HcaYPzHGfMgY8zZjzOd4l73JGPMRY8xHjTH/1BhjZh//YWPMo8aYG4X7utcY8zuzx/uQMeYvxf58CSH6qXpdMsasGWPeY4z54Oy15Qcqbl/3uvSPjDEfnv37upLb/lP/tanuvghpgmJIhsizAP4+gB+puPy/B/DRmtu/H8AXWmtfBuDfAngTABhj/hyALwHwMgAvBfBqAF8+u81/BPBFJff1vQB+0Vr7SgBfD+DHOn0mhJB5oep1aRvAX7DWvhzAKwB8lTHmz5bcvup16fUAXjW77WsAvNEYc9rdyBjzhQDOtrkvQtpAMSSDw1r7pLX2vQB2i5cZY+4G8HoAP1Fz+9+x1t6avfsuAHe7iwCsAVgBsApgGcDl2W3eZa19ouzuALgX6TMAHu/8CRFCBk/V65Kd4tK85dm/Y1OfNa9LLwHwe9baPWvtTQAfAvBVwLQdAfC/AfiulvdFSCMUQzJv/O+Yvkjut7z+twD4dQCw1r4TwO8AeGL27zestXXJIwB8P4A3GGM+A+DXAHxbj2MmhMwxs+UtHwDwJIDfsta+u+EmB69LAD6Iacq4YYw5D+B1AO6ZXfbfAfjlij9ay+6LkEZ4SjwyNxhjvhrAk9baPzLGvLbF9d8A4Asxq4uNMc8H8GIc/nX9W8aYL7PW/n81d/MNAH7KWvujxpgvBvAzxpiXWmvbiikhZM6x1k4AvMIYcxuAfzd7jfhw2XWLr0vW2t80xrwawB8AeArAOwFMjDF3AvjrAF5b9bjF+yKkDUwMySAwxvw9Y8wHZv/urLjalwD4GmPMpwH8PIC/YIz52Yr7+woA3wPga6y127MP/zUA77LW3phVP78O4IsbDu1bAPwicJA4rgE43+FTI4QMlJavSwdYa69i2kp8VcX9lb0uwVr7w9baV1hrvxKAAfAJAK8E8HwAfzp7zdswxvxp030R0gTFkAwCa+0/m70wvsJaW7qOz1r73dbau62192E6CPLb1to3FK9njHklgP8L0xfMJ72LHgHw5caYJWPMMqZ/ZTdVyY8A+M9m9/tiTMXwqY6fHiFkgLR5XTLGXJglhTDGrAP4SgAfK7le6evSrIa+ffb2yzAduCraMwAAAQhJREFUjvtNa+2vWmvvsNbeN3vNu2WtfX7dfRHSBp75hAwOY8wdAP4Q06GPfQA3ALzEWvucd53XAnijtfarZ+//IIA/tNb+sjHmrQC+ANN1hADwiLX2a2YLuX8MwJ/HdHH4f7LWfsfs9m8C8I0A7sR0wOQnrLXfb4x5CYD/G8DJ2W2+y1r7m0m/AIQQdVS9LgG4D8BPAxhjGsb8orX2B2e3afO6tAbgfbOPPQfgW621Hyh5/BvW2pOzt0vvK/KnTOYUiiEhhBBCCAHAKpkQQgghhMygGBJCCCGEEAAUQ0IIIYQQMoNiSAghhBBCAFAMCSGEEELIDIohIYQQQggBQDEkhBBCCCEz/n8EvnpCnSZ/pAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efdcc499320>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "class skymap:\n",
    "    def __init__(self, ulfile, pixelfile):\n",
    "        assert(type(ulfile) == str and type(pixelfile) == str)\n",
    "        assert(os.path.isfile(ulfile) and os.path.isfile(pixelfile))\n",
    "        \n",
    "        ul = np.loadtxt(ulfile)\n",
    "        pixel = np.genfromtxt(pixelfile, names = True)\n",
    "        \n",
    "        self.theta = pixel['theta']\n",
    "        self.phi = pixel['phi']\n",
    "        \n",
    "        self.map_dictionary = {}\n",
    "        self.map_dictionary[95] = self.create_map(ul[:,1])\n",
    "        self.map_dictionary[75] = self.create_map(ul[:,2])\n",
    "        self.map_dictionary[50] = self.create_map(ul[:,3])\n",
    "        self.map_dictionary[25] = self.create_map(ul[:,4])\n",
    "        self.map_dictionary[1] = self.create_map(ul[:,5])\n",
    "    \n",
    "    def create_map(self, ulcolumn):\n",
    "        npix = len(self.theta)\n",
    "        nside = hp.npix2nside(npix)\n",
    "        indices = hp.ang2pix(nside, self.theta, self.phi)\n",
    "        skymap = np.zeros(npix, dtype=np.float)\n",
    "        skymap[indices] += ulcolumn[indices]\n",
    "        return skymap\n",
    "    \n",
    "    def interpolate(self, conf, theta, phi):\n",
    "        if conf == 'all':\n",
    "            maps = []\n",
    "            for key in sorted(self.map_dictionary.keys()):\n",
    "                maps.append(hp.pixelfunc.get_interp_val(self.map_dictionary[key], theta, phi))\n",
    "            return maps\n",
    "        else:\n",
    "            return hp.pixelfunc.get_interp_val(self.map_dictionary[conf], theta, phi)\n",
    "    \n",
    "    def plot(self, confidence):\n",
    "        hp.mollview(np.log10(self.map_dictionary[confidence]), title= 'Upper Limit Skymap')\n",
    "        \n",
    "skymap3nHz = skymap('uls_smoothed_3nHz.txt', 'skymap_pixels.txt')\n",
    "skymap10nHz = skymap('uls_smoothed_10nHz.txt', 'skymap_pixels.txt')\n",
    "skymap31nHz = skymap('uls_smoothed_31nHz.txt','skymap_pixels.txt')\n",
    "skymap100nHz = skymap('uls_smoothed_100nHz.txt','skymap_pixels.txt')\n",
    "skymap318nHz = skymap('uls_smoothed_318nHz.txt','skymap_pixels.txt')\n",
    "\n",
    "skymap3nHz.plot(95)\n",
    "\n",
    "input_theta = .3\n",
    "input_phi = .4\n",
    "input_freq = 5\n",
    "input_strain = 5e-15\n",
    "\n",
    "#print(skymap3nHz.interpolate('all', input_theta, input_phi))\n",
    "#print(skymap10nHz.interpolate('all', input_theta, input_phi))\n",
    "#print(skymap31nHz.interpolate('all', input_theta, input_phi))\n",
    "#print(skymap100nHz.interpolate('all', input_theta, input_phi))\n",
    "#print(skymap318nHz.interpolate('all', input_theta, input_phi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "69.3674195194486\n"
     ]
    }
   ],
   "source": [
    "interp3nHz= np.array(skymap3nHz.interpolate('all', input_theta, input_phi))\n",
    "interp10nHz = np.array(skymap10nHz.interpolate('all', input_theta, input_phi))\n",
    "interp31nHz = np.array(skymap31nHz.interpolate('all', input_theta, input_phi))\n",
    "interp100nHz = np.array(skymap100nHz.interpolate('all', input_theta, input_phi))\n",
    "interp318nHz = np.array(skymap318nHz.interpolate('all', input_theta, input_phi))\n",
    "\n",
    "combined_columns = np.array([interp3nHz, interp10nHz, interp31nHz, interp100nHz, interp318nHz])\n",
    "\n",
    "#print(interp3nHz)\n",
    "#print(interp10nHz)\n",
    "\n",
    "#f = np.concatenate((interp1, interp2), axis=1)\n",
    "\n",
    "freq_array = np.array([3, 10, 31, 100, 318])\n",
    "#print(interp1)\n",
    "#print(freq_array)\n",
    "\n",
    "row1 = np.array(combined_columns[:,0])\n",
    "row25 = np.array(combined_columns[:,1])\n",
    "row50 = np.array(combined_columns[:,2])\n",
    "row75 = np.array(combined_columns[:,3])\n",
    "row95 = np.array(combined_columns[:,4])\n",
    "\n",
    "combined_rows = np.array([row1, row25, row50, row75, row95])\n",
    "\n",
    "#freq_interp = interp1d(freq_array, combined)\n",
    "freq_interp1 = interp1d(freq_array, combined_rows[0], kind = 'cubic')\n",
    "freq_interp2 = interp1d(freq_array, combined_rows[1], kind = 'cubic')\n",
    "freq_interp3 = interp1d(freq_array, combined_rows[2], kind = 'cubic')\n",
    "freq_interp4 = interp1d(freq_array, combined_rows[3], kind = 'cubic')\n",
    "freq_interp5 = interp1d(freq_array, combined_rows[4], kind = 'cubic')\n",
    "\n",
    "# = np.array([freq_interp1(input_freq), freq_interp2(input_freq), freq_interp3(input_freq), freq_interp4(input_freq), freq_interp5(input_freq)])\n",
    "\n",
    "\n",
    "confidence_array = [1, 25, 50, 75, 95]\n",
    "\n",
    "new_freq_rows = np.array([freq_interp1(input_freq), freq_interp2(input_freq), freq_interp3(input_freq), freq_interp4(input_freq), freq_interp5(input_freq)])\n",
    "\n",
    "#print(new_freq_rows)\n",
    "\n",
    "conf_interp = interp1d(new_freq_rows, confidence_array, kind='cubic')\n",
    "\n",
    "print(float(conf_interp(input_strain)))\n",
    "\n",
    "#confidence_interp = interp1d(new_freq_array,)\n",
    "#print(freq_interp2(input_freq))\n",
    "#print(freq_interp(input_freq))\n",
    "#freq_interp2(input_freq)\n",
    "#freq_interp3(input_freq)\n",
    "\n",
    "#print(full)\n",
    "#print(full[0][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}