{
 "cells": [
  {
   "cell_type": "raw",
   "id": "2622d571-d984-44f7-ac27-6ffb7fa61ee2",
   "metadata": {
    "vscode": {
     "languageId": "raw"
    }
   },
   "source": [
    "---\n",
    "date: 2024-07-29\n",
    "categories:\n",
    "  - cosmology\n",
    "  - pysm\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f229c7d7-453a-4810-bdb2-2c90b87b7bf4",
   "metadata": {},
   "source": [
    "# Generate point source maps with pixell and convert to HEALPix\n",
    "\n",
    "Testing the pixell `sim_objects` functionality to create maps of point sources pre-smoothed with a gaussian beam.\n",
    "The purpose is to include this functionality in PySM to be able to generate on the fly maps of source starting from a catalog.\n",
    "\n",
    "Compared to the [previous notebook](./2024-05-02-pysm-point-source-pixell.ipynb), here I also use `reproject` to convert the map to HEALPix, plot the result and check the flux in HEALPix also agrees with the input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6bc32901-4ebd-49d0-9b3d-788648643510",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from pixell import enmap, utils, resample, curvedsky as cs, reproject, pointsrcs\n",
    "import numpy as np\n",
    "import healpy as hp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "455c650f-5f4e-4cff-b1b1-84310cc2e2ba",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fwhm = 5 * utils.degree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a8ce7814-70d2-427d-b319-b692f2519ac3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "shape, wcs = enmap.fullsky_geometry(res=fwhm / 3, proj=\"car\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f984cc3e-441e-491a-a6d3-d9b1d6d26af2",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(109, 216)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "284621fa-b279-4ec6-8f25-f819063400a1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def fwhm2sigma(fwhm):\n",
    "    return fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e8789614-4b2a-44f4-98c5-4adfe8aa4b56",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def flux2amp(flux, fwhm):\n",
    "    sigma = fwhm2sigma(fwhm)\n",
    "    return flux / (2 * np.pi * sigma**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "bee0e3cf-96b5-4257-93bf-d72c89c5e3ac",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "assert flux2amp((2 * np.pi * fwhm2sigma(5) ** 2), 5) == 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3d492f49-3236-4f47-beae-7326425253fe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "n_sources = 1\n",
    "flux_sources = np.arange(n_sources) + 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a590ec0b-6028-4033-874b-7f383bb74c04",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "amplitude_sources = flux2amp(flux_sources, fwhm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "a0c43d61-b439-4c28-8be5-05076bc0315f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "source_pos = np.array([[np.pi/4], [np.pi/3]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "c66c3b78-ad8a-4049-86c5-6c3f0b2aa596",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "r, p = pointsrcs.expand_beam(fwhm2sigma(fwhm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "57f02ead-2d33-4cbb-9262-9755aae395e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "source_map = pointsrcs.sim_objects(shape, wcs, source_pos, amplitude_sources, ((r, p)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "8871cee2-0d12-488d-9513-ed9444e59063",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "d3967f9b-87ee-40e0-9d02-300756459462",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.78539816],\n",
       "       [1.04719755]])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "source_pos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "4957a65f-e376-41d1-b93f-ed4553fb25b0",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f81203692d0>"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAEqCAYAAAA/LasTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhZklEQVR4nO3de3DU9b3/8dd3L9lcTFa57WYlYPQXfyqhHAkeKrWS2pJzqBccOl4KpwdPexwsYE2xB2SoFZ02KZw51JlyRO1YpcdS/EfU32hb4qlEGeo0cmk17UEdU4hKjHrCbshlr5/fH5GtSyAJuGE/G56Pme80+/l+dnl/+tnv7MvP97v7dYwxRgAAABZx5boAAACA4xFQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1chpQHnroIVVWVqqwsFA1NTV65ZVXclkOAACwRM4CylNPPaX6+nqtXbtW+/bt0xe/+EXNnz9fhw4dylVJAADAEk6ubhY4e/ZszZw5U5s3b063XXrppbrxxhvV2Ng45HNTqZTef/99lZaWynGc0S4VAABkgTFG3d3dCoVCcrmGXiPxnKGaMsRiMe3Zs0f33HNPRntdXZ127949qH80GlU0Gk0/fu+993TZZZeNep0AACD72tvbNXny5CH75CSgfPTRR0omkwoEAhntgUBAHR0dg/o3Njbq/vvvH9R+lb4qj7yjVicAAMiehOLapRdUWlo6bN+cBJRjjj89Y4w54SmbNWvWaOXKlenHkUhEFRUV8sgrj0NAAQAgL3xyUclILs/ISUCZMGGC3G73oNWSzs7OQasqkuTz+eTz+c5UeQAAIMdy8i2egoIC1dTUqKmpKaO9qalJc+bMyUVJAADAIjk7xbNy5Up94xvf0KxZs3TllVfq0Ucf1aFDh3THHXfkqiQAAGCJnAWUW265RR9//LEeeOABHT58WNXV1XrhhRc0derUXJUEAAAskbPfQfksIpGI/H6/arWAi2QBAMgTCRPXTj2rcDissrKyIftyLx4AAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWCfrAaWxsVFXXHGFSktLNWnSJN144406cOBARh9jjNatW6dQKKSioiLV1taqtbU126UAAIA8lfWA0tzcrOXLl+vVV19VU1OTEomE6urq1NPTk+6zYcMGbdy4UZs2bVJLS4uCwaDmzZun7u7ubJcDAADykGOMMaP5D3z44YeaNGmSmpubdfXVV8sYo1AopPr6eq1evVqSFI1GFQgEtH79ei1dunTY14xEIvL7/arVAnkc72iWDwAAsiRh4tqpZxUOh1VWVjZk31G/BiUcDkuSxo0bJ0lqa2tTR0eH6urq0n18Pp/mzp2r3bt3n/A1otGoIpFIxgYAAMauUQ0oxhitXLlSV111laqrqyVJHR0dkqRAIJDRNxAIpPcdr7GxUX6/P71VVFSMZtkAACDHRjWgrFixQn/605/0q1/9atA+x3EyHhtjBrUds2bNGoXD4fTW3t4+KvUCAAA7eEbrhe+8804999xzevnllzV58uR0ezAYlDSwklJeXp5u7+zsHLSqcozP55PP5xutUgEAgGWyvoJijNGKFSv09NNP63e/+50qKysz9ldWVioYDKqpqSndFovF1NzcrDlz5mS7HAAAkIeyvoKyfPlybd26Vc8++6xKS0vT15X4/X4VFRXJcRzV19eroaFBVVVVqqqqUkNDg4qLi7Vo0aJslwMAAPJQ1gPK5s2bJUm1tbUZ7Y8//rhuu+02SdKqVavU19enZcuWqaurS7Nnz9aOHTtUWlqa7XIAAEAeGvXfQRkN/A4KAAD5x6rfQQEAADhVBBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsM6oB5TGxkY5jqP6+vp0mzFG69atUygUUlFRkWpra9Xa2jrapQAAgDwxqgGlpaVFjz76qD73uc9ltG/YsEEbN27Upk2b1NLSomAwqHnz5qm7u3s0ywEAAHli1ALK0aNHtXjxYv3sZz/Teeedl243xujBBx/U2rVrtXDhQlVXV2vLli3q7e3V1q1bR6scAACQR0YtoCxfvlzXXnutvvKVr2S0t7W1qaOjQ3V1dek2n8+nuXPnavfu3Sd8rWg0qkgkkrEBAICxyzMaL7pt2zbt3btXLS0tg/Z1dHRIkgKBQEZ7IBDQwYMHT/h6jY2Nuv/++7NfKAAAsFLWV1Da29t111136cknn1RhYeFJ+zmOk/HYGDOo7Zg1a9YoHA6nt/b29qzWDAAA7JL1FZQ9e/aos7NTNTU16bZkMqmXX35ZmzZt0oEDByQNrKSUl5en+3R2dg5aVTnG5/PJ5/Nlu1QAAGCprK+gfPnLX9brr7+u/fv3p7dZs2Zp8eLF2r9/vy688EIFg0E1NTWlnxOLxdTc3Kw5c+ZkuxwAAJCHsr6CUlpaqurq6oy2kpISjR8/Pt1eX1+vhoYGVVVVqaqqSg0NDSouLtaiRYuyXQ4AAMhDo3KR7HBWrVqlvr4+LVu2TF1dXZo9e7Z27Nih0tLSXJQDAAAs4xhjTK6LOFWRSER+v1+1WiCP4811OQAAYAQSJq6delbhcFhlZWVD9uVePAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwzqgElPfee0//9E//pPHjx6u4uFh/93d/pz179qT3G2O0bt06hUIhFRUVqba2Vq2traNRCgAAyENZDyhdXV36whe+IK/Xq1//+tf685//rP/4j//Queeem+6zYcMGbdy4UZs2bVJLS4uCwaDmzZun7u7ubJcDAADykCfbL7h+/XpVVFTo8ccfT7ddcMEF6b+NMXrwwQe1du1aLVy4UJK0ZcsWBQIBbd26VUuXLs12SQAAIM9kfQXlueee06xZs3TTTTdp0qRJuvzyy/Wzn/0svb+trU0dHR2qq6tLt/l8Ps2dO1e7d+8+4WtGo1FFIpGMDQAAjF1ZDyjvvPOONm/erKqqKv32t7/VHXfcoe985zv6xS9+IUnq6OiQJAUCgYznBQKB9L7jNTY2yu/3p7eKiopslw0AACyS9YCSSqU0c+ZMNTQ06PLLL9fSpUt1++23a/PmzRn9HMfJeGyMGdR2zJo1axQOh9Nbe3t7tssGAAAWyXpAKS8v12WXXZbRdumll+rQoUOSpGAwKEmDVks6OzsHraoc4/P5VFZWlrEBAICxK+sB5Qtf+IIOHDiQ0fbmm29q6tSpkqTKykoFg0E1NTWl98diMTU3N2vOnDnZLgcAAOShrH+L57vf/a7mzJmjhoYG3XzzzfrDH/6gRx99VI8++qikgVM79fX1amhoUFVVlaqqqtTQ0KDi4mItWrQo2+UAAIA8lPWAcsUVV2j79u1as2aNHnjgAVVWVurBBx/U4sWL031WrVqlvr4+LVu2TF1dXZo9e7Z27Nih0tLSbJcDAADykGOMMbku4lRFIhH5/X7VaoE8jjfX5QAAgBFImLh26lmFw+FhryflXjwAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANbx5LoAYExwnMFtxpz5OgBgjCCgAJ/FiYLJ8fsIKgBwyggowKk4FjoclxyXM/C/7k/OlLo+dcY0lZIkmWRKMimZlJHMQBuBBQCGR0ABRspx/hZM3G45brfkOHIKCiSXk94vkxoIIcmkFE8M/B1PSMYZCCpKEVIAYBgEFGAkjoUTt3tgxcTrleMrkOPxSIU+yeWS8QwEFhkjJ5EcWEXp6x8IKtGYFI9LyZRMUiKkAMDQsv4tnkQioe9///uqrKxUUVGRLrzwQj3wwANKfbLkLUnGGK1bt06hUEhFRUWqra1Va2trtksBPjtnYGXEcbvleD1yFRXKKS2V61y/THCikudPUP9FE9Xzfyfq6KXjFaker6OXjlfvxRPVf+FEpUITpcAEuc71DzyvqFCO15NefRnyGhYAOItlfQVl/fr1evjhh7VlyxZNmzZNr732mv7lX/5Ffr9fd911lyRpw4YN2rhxo5544gldfPHF+uEPf6h58+bpwIEDKi0tzXZJwGfjuAZO6TiO5PXIKfDKFPmUOqdAKZ9b0XM9ShY4SvocpbySKy65o0aeqEuuZIFM1C13LC4nlZJJJeXEYjKOS1Iy1yMDAGtlPaD8/ve/14IFC3TttddKki644AL96le/0muvvSZpYPXkwQcf1Nq1a7Vw4UJJ0pYtWxQIBLR161YtXbo02yUBp+dT15w4BQOncxx/mVJlxYqNK1J3RYHiJY56phglilMy5yTkKUwo0eeR0+OWp8elcw4VytsjlbZ75O3qlyvcI6WMlEhIfSmuSQGAk8j6KZ6rrrpK//3f/60333xTkvTHP/5Ru3bt0le/+lVJUltbmzo6OlRXV5d+js/n09y5c7V79+4TvmY0GlUkEsnYgDMh45s6Ho9MYYGSJQWK+T2KjnPUP0FKnt+vkopu/Z+pH6hm6iFdOLVTxZOPKlEeU/94R/3jHMXKPEqWFMgUFkgez8DrHbvgFgAwSNZXUFavXq1wOKxLLrlEbrdbyWRSP/rRj/T1r39dktTR0SFJCgQCGc8LBAI6ePDgCV+zsbFR999/f7ZLBYb3qXDiFHiV9HmVKPEodo5L0fOMYn6jKYEuTT7niKaXvqepBR/p7WhA/1McVFvxOHV0TZLxuBQ7xyVPn0fuHq9cHreM8Qx88yfFqR4AOJGsB5SnnnpKTz75pLZu3app06Zp//79qq+vVygU0pIlS9L9nOMuDjTGDGo7Zs2aNVq5cmX6cSQSUUVFRbZLBzI5f/t9E8fjkTxupYo8ShS7FCt1FDsvKde4mGZP/Ks+V9yuzxce1EXec/QX33tqKZiivb4L9P8+9CvqLlD8HJcSPS55irxyFXilVEqO48jIfPLVZEIKAHxa1gPKv/3bv+mee+7RrbfeKkmaPn26Dh48qMbGRi1ZskTBYFDSwEpKeXl5+nmdnZ2DVlWO8fl88vl82S4VGJlUaiBAp4ycREquuJE7ZuTucynR41Fbz3h5naQKnbi6Ux/qzXi5/qcvpEM95ynV45WnzyVXTHLFjVyJlJRMDVyHcoxJnfzfBoCzVNYDSm9vr1yuzEtb3G53+mvGlZWVCgaDampq0uWXXy5JisViam5u1vr167NdDnD6TEqSe+DvZFJKJOSKJuTpTcoXccv3vy654h7tf/d8tZ0zXm+UhjSx8Kg+6C/V4e4yHYkUq+BDtwq6Hfm6k/L0JuVE4wOvlUrKcGEsAJxU1gPK9ddfrx/96EeaMmWKpk2bpn379mnjxo365je/KWng1E59fb0aGhpUVVWlqqoqNTQ0qLi4WIsWLcp2OcBnY1IyxpGTTMm4U3LiSbliKbmjKXl6XDKOo96wT/+bcCuWcOujwhJF+n3qOVqo1FGvSnoceXold7+RKzbwfJP61AoKqycAcEJZDyg//elPde+992rZsmXq7OxUKBTS0qVL9YMf/CDdZ9WqVerr69OyZcvU1dWl2bNna8eOHfwGCqxjUkaOkjLRqJxkUk74qDzxhIp7i+TuL1Ky0KWiDz1KFniVLCrSEa/kikll/ZKn36j4w4Tc0ZR8nb1yevrlHO2V6euXSSSkZPKTrxkDAI7nmDxcZ45EIvL7/arVAnkcb67LwVh27LdQPvn1V6e4SI7PJ1NcqJS/WKkCt2LnFihV4ChZ4FLKI7kSkjuWkitmVHAkJlcsKdeRHjl9UZn+/oGAkkzKxBN/u28PAJwFEiaunXpW4XBYZWVlQ/blXjzAcEzqb6sdff0DP7IWi8sdi8vtdssdKZQ8LplPNieRkpNISYmUXL0D9+Ixvf0y8ZhMLP6pYMLpHQA4GQIKMJRPVjdMMikpOXBNSiwuxx2V6e1N36NHjiPHcX3y2ybmk1CTUiqRGLhpYDIpkxwIJQOvJVZOAGAIBBRgJIyRHGfgmhTXwB2JHX0SXJLJgd9KcZx0QDHGpIOJpL+Fk/TFsYQTABgKAQUYKWMkkxw4M+M4MolP2p0h7hjx6dM4hBIAGDECCnA6MsLGENeSEEoA4LQQUIDPihACAFmX9bsZAwAAfFYEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArHPKAeXll1/W9ddfr1AoJMdx9Mwzz2TsN8Zo3bp1CoVCKioqUm1trVpbWzP6RKNR3XnnnZowYYJKSkp0ww036N133/1MAwEAAGPHKQeUnp4ezZgxQ5s2bTrh/g0bNmjjxo3atGmTWlpaFAwGNW/ePHV3d6f71NfXa/v27dq2bZt27dqlo0eP6rrrrlMymTz9kQAAgDHDMcaY036y42j79u268cYbJQ2snoRCIdXX12v16tWSBlZLAoGA1q9fr6VLlyocDmvixIn6r//6L91yyy2SpPfff18VFRV64YUX9A//8A/D/ruRSER+v1+1WiCP4z3d8gEAwBmUMHHt1LMKh8MqKysbsm9Wr0Fpa2tTR0eH6urq0m0+n09z587V7t27JUl79uxRPB7P6BMKhVRdXZ3uc7xoNKpIJJKxAQCAsSurAaWjo0OSFAgEMtoDgUB6X0dHhwoKCnTeeeedtM/xGhsb5ff701tFRUU2ywYAAJYZlW/xOI6T8dgYM6jteEP1WbNmjcLhcHprb2/PWq0AAMA+WQ0owWBQkgathHR2dqZXVYLBoGKxmLq6uk7a53g+n09lZWUZGwAAGLuyGlAqKysVDAbV1NSUbovFYmpubtacOXMkSTU1NfJ6vRl9Dh8+rDfeeCPdBwAAnN08p/qEo0eP6u23304/bmtr0/79+zVu3DhNmTJF9fX1amhoUFVVlaqqqtTQ0KDi4mItWrRIkuT3+/Wtb31Ld999t8aPH69x48bpe9/7nqZPn66vfOUr2RsZAADIW6ccUF577TV96UtfSj9euXKlJGnJkiV64okntGrVKvX19WnZsmXq6urS7NmztWPHDpWWlqaf85Of/EQej0c333yz+vr69OUvf1lPPPGE3G53FoYEAADy3Wf6HZRc4XdQAADIPzn7HRQAAIBsIKAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHVO+WaBNjh2+6CE4lLe3UkIAICzU0JxSX/7HB9KXgaU7u5uSdIuvZDjSgAAwKnq7u6W3+8fsk9e3s04lUrpwIEDuuyyy9Te3j7sHRHzWSQSUUVFxZge59kwRunsGOfZMEbp7BgnYxw7bBqnMUbd3d0KhUJyuYa+yiQvV1BcLpfOP/98SVJZWVnO/w8/E86GcZ4NY5TOjnGeDWOUzo5xMsaxw5ZxDrdycgwXyQIAAOsQUAAAgHXyNqD4fD7dd9998vl8uS5lVJ0N4zwbxiidHeM8G8YonR3jZIxjR76OMy8vkgUAAGNb3q6gAACAsYuAAgAArENAAQAA1iGgAAAA6+RtQHnooYdUWVmpwsJC1dTU6JVXXsl1SaetsbFRV1xxhUpLSzVp0iTdeOONOnDgQEaf2267TY7jZGyf//znc1TxqVu3bt2g+oPBYHq/MUbr1q1TKBRSUVGRamtr1dramsOKT88FF1wwaJyO42j58uWS8nMeX375ZV1//fUKhUJyHEfPPPNMxv6RzF00GtWdd96pCRMmqKSkRDfccIPefffdMziK4Q01zng8rtWrV2v69OkqKSlRKBTSP//zP+v999/PeI3a2tpB83vrrbee4ZGc3HBzOZL3Z77PpaQTHqOO4+jf//3f031snsuRfGaMheMyLwPKU089pfr6eq1du1b79u3TF7/4Rc2fP1+HDh3KdWmnpbm5WcuXL9err76qpqYmJRIJ1dXVqaenJ6PfP/7jP+rw4cPp7YUX8uteRNOmTcuo//XXX0/v27BhgzZu3KhNmzappaVFwWBQ8+bNS993KV+0tLRkjLGpqUmSdNNNN6X75Ns89vT0aMaMGdq0adMJ949k7urr67V9+3Zt27ZNu3bt0tGjR3XdddcpmUyeqWEMa6hx9vb2au/evbr33nu1d+9ePf3003rzzTd1ww03DOp7++23Z8zvI488cibKH5Hh5lIa/v2Z73MpKWN8hw8f1s9//nM5jqOvfe1rGf1sncuRfGaMiePS5KG///u/N3fccUdG2yWXXGLuueeeHFWUXZ2dnUaSaW5uTrctWbLELFiwIHdFfUb33XefmTFjxgn3pVIpEwwGzY9//ON0W39/v/H7/ebhhx8+QxWOjrvuustcdNFFJpVKGWPyfx4lme3bt6cfj2Tujhw5Yrxer9m2bVu6z3vvvWdcLpf5zW9+c8ZqPxXHj/NE/vCHPxhJ5uDBg+m2uXPnmrvuumt0i8uSE41xuPfnWJ3LBQsWmGuuuSajLZ/m8vjPjLFyXObdCkosFtOePXtUV1eX0V5XV6fdu3fnqKrsCofDkqRx48ZltO/cuVOTJk3SxRdfrNtvv12dnZ25KO+0vfXWWwqFQqqsrNStt96qd955R5LU1tamjo6OjDn1+XyaO3duXs9pLBbTk08+qW9+85tyHCfdnu/z+Gkjmbs9e/YoHo9n9AmFQqqurs7r+Q2Hw3IcR+eee25G+y9/+UtNmDBB06ZN0/e+9728WwUc6v05Fufygw8+0PPPP69vfetbg/bly1we/5kxVo7LvLtZ4EcffaRkMqlAIJDRHggE1NHRkaOqsscYo5UrV+qqq65SdXV1un3+/Pm66aabNHXqVLW1tenee+/VNddcoz179uTFrwPOnj1bv/jFL3TxxRfrgw8+0A9/+EPNmTNHra2t6Xk70ZwePHgwF+VmxTPPPKMjR47otttuS7fl+zwebyRz19HRoYKCAp133nmD+uTrMdvf36977rlHixYtyrj52uLFi1VZWalgMKg33nhDa9as0R//+Mf0qT7bDff+HItzuWXLFpWWlmrhwoUZ7fkylyf6zBgrx2XeBZRjPv1fpNLAJB3flo9WrFihP/3pT9q1a1dG+y233JL+u7q6WrNmzdLUqVP1/PPPDzqwbDR//vz039OnT9eVV16piy66SFu2bElfhDfW5vSxxx7T/PnzFQqF0m35Po8nczpzl6/zG4/HdeuttyqVSumhhx7K2Hf77ben/66urlZVVZVmzZqlvXv3aubMmWe61FN2uu/PfJ1LSfr5z3+uxYsXq7CwMKM9X+byZJ8ZUv4fl3l3imfChAlyu92DEl5nZ+egtJhv7rzzTj333HN66aWXNHny5CH7lpeXa+rUqXrrrbfOUHXZVVJSounTp+utt95Kf5tnLM3pwYMH9eKLL+pf//Vfh+yX7/M4krkLBoOKxWLq6uo6aZ98EY/HdfPNN6utrU1NTU3D3rp+5syZ8nq9eTu/x78/x9JcStIrr7yiAwcODHucSnbO5ck+M8bKcZl3AaWgoEA1NTWDltmampo0Z86cHFX12RhjtGLFCj399NP63e9+p8rKymGf8/HHH6u9vV3l5eVnoMLsi0aj+stf/qLy8vL0Muqn5zQWi6m5uTlv5/Txxx/XpEmTdO211w7ZL9/ncSRzV1NTI6/Xm9Hn8OHDeuONN/Jqfo+Fk7feeksvvviixo8fP+xzWltbFY/H83Z+j39/jpW5POaxxx5TTU2NZsyYMWxfm+ZyuM+MMXNc5uji3M9k27Ztxuv1mscee8z8+c9/NvX19aakpMT89a9/zXVpp+Xb3/628fv9ZufOnebw4cPprbe31xhjTHd3t7n77rvN7t27TVtbm3nppZfMlVdeac4//3wTiURyXP3I3H333Wbnzp3mnXfeMa+++qq57rrrTGlpaXrOfvzjHxu/32+efvpp8/rrr5uvf/3rpry8PG/G92nJZNJMmTLFrF69OqM9X+exu7vb7Nu3z+zbt89IMhs3bjT79u1Lf3tlJHN3xx13mMmTJ5sXX3zR7N2711xzzTVmxowZJpFI5GpYgww1zng8bm644QYzefJks3///ozjNBqNGmOMefvtt839999vWlpaTFtbm3n++efNJZdcYi6//HJrxjnUGEf6/sz3uTwmHA6b4uJis3nz5kHPt30uh/vMMGZsHJd5GVCMMeY///M/zdSpU01BQYGZOXNmxldy842kE26PP/64McaY3t5eU1dXZyZOnGi8Xq+ZMmWKWbJkiTl06FBuCz8Ft9xyiykvLzder9eEQiGzcOFC09ramt6fSqXMfffdZ4LBoPH5fObqq682r7/+eg4rPn2//e1vjSRz4MCBjPZ8nceXXnrphO/PJUuWGGNGNnd9fX1mxYoVZty4caaoqMhcd9111o17qHG2tbWd9Dh96aWXjDHGHDp0yFx99dVm3LhxpqCgwFx00UXmO9/5jvn4449zO7BPGWqMI31/5vtcHvPII4+YoqIic+TIkUHPt30uh/vMMGZsHJeOMcaM0uIMAADAacm7a1AAAMDYR0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHX+PxbmeJONZlc4AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(source_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "52b78304-d1af-43d9-bab6-138d32baa181",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def aperture_photometry(\n",
    "    thumbs, aperture_radius, annulus_width=None, modrmap=None, pixsizemap=None\n",
    "):\n",
    "    \"\"\"\n",
    "    Flux from aperture photometry.\n",
    "\n",
    "    from https://github.com/msyriac/orphics/blob/master/orphics/maps.py\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    thumb : ndmap\n",
    "        An (...,Ny,Nx) ndmap (i.e. a pixell enmap) containing the thumbnails.\n",
    "    aperture_radius : float\n",
    "        Aperture inner radius in radians\n",
    "    annulus_width : float\n",
    "        Annulus width for mean subtraction in radians.\n",
    "        Defaults to sqrt(2)-1 times the aperture inner radius.\n",
    "    modrmap : ndmap, optional\n",
    "        An (Ny,Nx) ndmap containing distances of each pixel from the center in radians.\n",
    "    modrmap : ndmap, optional\n",
    "        An (Ny,Nx) ndmap containing pixel areas in steradians.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    flux : ndarray\n",
    "        (...,) array of aperture photometry fluxes.\n",
    "\n",
    "    \"\"\"\n",
    "    if modrmap is None:\n",
    "        modrmap = thumbs.modrmap()\n",
    "    if annulus_width is None:\n",
    "        annulus_width = (np.sqrt(2.0) - 1.0) * aperture_radius\n",
    "    # Get the mean background level from the annulus\n",
    "    mean = thumbs[\n",
    "        ...,\n",
    "        np.logical_and(\n",
    "            modrmap > aperture_radius, modrmap < (aperture_radius + annulus_width)\n",
    "        ),\n",
    "    ].mean()\n",
    "    if pixsizemap is None:\n",
    "        pixsizemap = thumbs.pixsizemap()\n",
    "    # Subtract the mean, multiply by pixel areas and sum\n",
    "    return (((thumbs - mean) * pixsizemap)[..., modrmap <= aperture_radius]).sum(\n",
    "        axis=-1\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "e5c8025a-8dbc-49fe-b2c9-2ba88e3a2110",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from astropy import units as u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "e7939667-c6b9-46ff-a577-068a4aafc5de",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "box_half_size_rad = 2 * fwhm\n",
    "box_center = [source_pos[0, -1], source_pos[1, -1]]\n",
    "box = np.array(\n",
    "    [\n",
    "        [box_center[0] - box_half_size_rad, box_center[1] - box_half_size_rad],\n",
    "        [box_center[0] + box_half_size_rad, box_center[1] + box_half_size_rad],\n",
    "    ]\n",
    ")  # in radians"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "33152991-7bf1-4aad-bbb8-7ed5591109c5",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.7853981633974483, 1.0471975511965976]"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "box_center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "e2851ab0-8af5-44e7-9bf0-49d2856ab87d",
   "metadata": {},
   "outputs": [],
   "source": [
    "cutout = source_map.submap(box)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "55e7302f-2bbd-4afa-b3bf-7bf10ccecd51",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f810b46c210>"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAGdCAYAAAC7EMwUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAX7ElEQVR4nO3df2zV9b3H8dehhdOWe3oQSH/NgmW3CUj9gS1bgCoYtQki0Zi4qaBEtgViQWoTBww3lYWegRtpYmdJ+YNhSJE/Nn4s0WnjBpUgsRRQwgwEZbQBu1522Tml6CltP/cPZ+8qFWF+v333nD4fyfmj337l/T6Bnqff9vScgHPOCQAAAyOsFwAADF9ECABghggBAMwQIQCAGSIEADBDhAAAZogQAMAMEQIAmEm1XuCrent7de7cOYVCIQUCAet1AADXyTmnjo4O5eXlacSIq1/rDLkInTt3Tvn5+dZrAAC+pdbWVt14441XPWfIRSgUCkmSSnW/UjXSeBsAwPXq1mXt1xt9j+dXM+Qi9OW34FI1UqkBIgQACedfr0h6LT9S4YkJAAAzRAgAYIYIAQDMECEAgBkiBAAwQ4QAAGaIEADADBECAJjxLUKvvvqqCgoKlJaWpuLiYr377rt+jQIAJChfIrRjxw5VVFRozZo1OnLkiO68807NnTtXLS0tfowDACQoXyK0ceNG/ehHP9KPf/xjTZkyRdXV1crPz1dtba0f4wAACcrzCHV1dam5uVllZWX9jpeVlenAgQNXnB+PxxWLxfrdAADDg+cROn/+vHp6epSdnd3veHZ2ttra2q44PxKJKBwO9914GwcAGD58e2LCV1891Tk34Cuqrl69WtFotO/W2trq10oAgCHG87dyGD9+vFJSUq646mlvb7/i6kiSgsGggsGg12sAABKA51dCo0aNUnFxsRoaGvodb2ho0MyZM70eBwBIYL68qV1lZaWeeOIJlZSUaMaMGaqrq1NLS4uWLl3qxzgAQILyJUI//OEP9Y9//ENr167Vp59+qqKiIr3xxhuaOHGiH+MAAAkq4Jxz1kv8u1gspnA4rDl6kLf3BoAE1O0ua692KxqNKjMz86rn8tpxAAAzRAgAYIYIAQDMECEAgBkiBAAwQ4QAAGaIEADADBECAJghQgAAM0QIAGCGCAEAzBAhAIAZIgQAMEOEAABmiBAAwAwRAgCYIUIAADNECABghggBAMwQIQCAGSIEADBDhAAAZogQAMAMEQIAmCFCAAAzRAgAYIYIAQDMECEAgBkiBAAwQ4QAAGaIEADADBECAJghQgAAM0QIAGCGCAEAzBAhAIAZIgQAMEOEAABmiBAAwEyq9QKAmUDAeoPhxznrDTDEcCUEADBDhAAAZogQAMAMEQIAmCFCAAAzRAgAYIYIAQDMECEAgBkiBAAwQ4QAAGY8j1AkEtH06dMVCoWUlZWlhx56SCdOnPB6DAAgCXgeoX379qm8vFwHDx5UQ0ODuru7VVZWps7OTq9HAQASnOcvYPqnP/2p38dbtmxRVlaWmpubddddd3k9DgCQwHz/mVA0GpUkjR071u9RAIAE4+tbOTjnVFlZqdLSUhUVFQ14TjweVzwe7/s4Fov5uRIAYAjx9Upo2bJl+vDDD7V9+/avPScSiSgcDvfd8vPz/VwJADCEBJzz512mli9frl27dqmxsVEFBQVfe95AV0L5+fmaoweVGhjpx2rAF3hTu8HHm9oNC93usvZqt6LRqDIzM696ruffjnPOafny5dq5c6f27t171QBJUjAYVDAY9HoNAEAC8DxC5eXlqq+v1+7duxUKhdTW1iZJCofDSk9P93ocACCBef4zodraWkWjUc2ZM0e5ubl9tx07dng9CgCQ4Hz5dhwAANeC144DAJghQgAAM0QIAGCGCAEAzBAhAIAZIgQAMEOEAABmiBAAwIyvb+WA62T0gpqBlBSTuZIUMHzdwECa4WsWphp+6XV3m412n8e/+SQ/5sZt5kqS6+kxm50ILxjLlRAAwAwRAgCYIUIAADNECABghggBAMwQIQCAGSIEADBDhAAAZogQAMAMEQIAmCFCAAAzRAgAYIYIAQDMECEAgBkiBAAwQ4QAAGaIEADADBECAJghQgAAM0QIAGCGCAEAzBAhAIAZIgQAMEOEAABmiBAAwAwRAgCYIUIAADNECABghggBAMwQIQCAmVTrBYacQMBu9KhRJnNTbhhjMleSer4z3mx2Z/5os9ldo+3+/29UZ6/Z7IzWTpO5KWfPm8yVpJ4L/zSb7bq6jCYHJHdtZ3IlBAAwQ4QAAGaIEADADBECAJghQgAAM0QIAGCGCAEAzBAhAIAZIgQAMEOEAABmiBAAwIzvEYpEIgoEAqqoqPB7FAAgwfgaoaamJtXV1enWW2/1cwwAIEH5FqGLFy9qwYIF2rx5s2644Qa/xgAAEphvESovL9e8efN07733XvW8eDyuWCzW7wYAGB58eT+h119/XYcPH1ZTU9M3nhuJRPTSSy/5sQYAYIjz/EqotbVVK1as0LZt25SWlvaN569evVrRaLTv1tra6vVKAIAhyvMroebmZrW3t6u4uLjvWE9PjxobG1VTU6N4PK6UlJS+zwWDQQWDQa/XAAAkAM8jdM899+jYsWP9jj311FOaPHmyVq5c2S9AAIDhzfMIhUIhFRUV9Ts2evRojRs37orjAIDhjVdMAACY8eXZcV+1d+/ewRgDAEgwXAkBAMwQIQCAGSIEADBDhAAAZogQAMAMEQIAmCFCAAAzRAgAYGZQflk1kQQMX9tuROi/TOZe/m6uyVxJ+nRGhtns3hlRs9m35HxqNvtYm93f94X3wiZzc98baTJXklJPdJvN7v2nzb/xgOuVrvFucyUEADBDhAAAZogQAMAMEQIAmCFCAAAzRAgAYIYIAQDMECEAgBkiBAAwQ4QAAGaIEADADBECAJghQgAAM0QIAGCGCAEAzBAhAIAZIgQAMEOEAABmiBAAwAwRAgCYIUIAADNECABghggBAMwQIQCAGSIEADBDhAAAZogQAMAMEQIAmCFCAAAzRAgAYCbVeoGvFQh8cRtsKSmDP/NfAunpJnM/ywqazJWki//dbTb7palvmc1+MvO82ezXxo03m/3C/zxsMvezj+3+jWe22HxdS5I6LtrMdb3SNX5pcyUEADBDhAAAZogQAMAMEQIAmCFCAAAzRAgAYIYIAQDMECEAgBkiBAAwQ4QAAGaIEADAjC8ROnv2rBYuXKhx48YpIyNDt99+u5qbm/0YBQBIYJ6/gOmFCxc0a9Ys3X333XrzzTeVlZWljz/+WGPGjPF6FAAgwXkeofXr1ys/P19btmzpO3bTTTd5PQYAkAQ8/3bcnj17VFJSokceeURZWVmaNm2aNm/e/LXnx+NxxWKxfjcAwPDgeYQ++eQT1dbWqrCwUG+99ZaWLl2qZ555Rq+99tqA50ciEYXD4b5bfn6+1ysBAIYozyPU29urO+64Q1VVVZo2bZqWLFmin/zkJ6qtrR3w/NWrVysajfbdWltbvV4JADBEeR6h3Nxc3Xzzzf2OTZkyRS0tLQOeHwwGlZmZ2e8GABgePI/QrFmzdOLEiX7HTp48qYkTJ3o9CgCQ4DyP0LPPPquDBw+qqqpKp06dUn19verq6lReXu71KABAgvM8QtOnT9fOnTu1fft2FRUV6Ze//KWqq6u1YMECr0cBABKc578nJEkPPPCAHnjgAT/+aABAEuG14wAAZogQAMAMEQIAmCFCAAAzRAgAYIYIAQDMECEAgBkiBAAw48svq3rCOUlu8Of2Gsz80uXLJmODF2zmStLov6WZzX75ozKz2W9ktZnNPt6eYzZ79N9sHnKCFz43mSvJ7Otakt3jmbv2uVwJAQDMECEAgBkiBAAwQ4QAAGaIEADADBECAJghQgAAM0QIAGCGCAEAzBAhAIAZIgQAMEOEAABmiBAAwAwRAgCYIUIAADNECABghggBAMwQIQCAGSIEADBDhAAAZogQAMAMEQIAmCFCAAAzRAgAYIYIAQDMECEAgBkiBAAwQ4QAAGaIEADADBECAJhJtV5gqHE9PWaze2MdJnODp/5uMleSvnNprNnsSydCZrP/lpFpNnvMJWc2O+Oszb/x1HP/azJXsvu6luwez5y79rlcCQEAzBAhAIAZIgQAMEOEAABmiBAAwAwRAgCYIUIAADNECABghggBAMwQIQCAGc8j1N3dreeff14FBQVKT0/XpEmTtHbtWvX29no9CgCQ4Dx/7bj169dr06ZN2rp1q6ZOnapDhw7pqaeeUjgc1ooVK7weBwBIYJ5H6L333tODDz6oefPmSZJuuukmbd++XYcOHfJ6FAAgwXn+7bjS0lK98847OnnypCTpgw8+0P79+3X//fcPeH48HlcsFut3AwAMD55fCa1cuVLRaFSTJ09WSkqKenp6tG7dOj322GMDnh+JRPTSSy95vQYAIAF4fiW0Y8cObdu2TfX19Tp8+LC2bt2qX//619q6deuA569evVrRaLTv1tra6vVKAIAhyvMroeeee06rVq3So48+Kkm65ZZbdObMGUUiES1atOiK84PBoILBoNdrAAASgOdXQpcuXdKIEf3/2JSUFJ6iDQC4gudXQvPnz9e6des0YcIETZ06VUeOHNHGjRu1ePFir0cBABKc5xF65ZVX9POf/1xPP/202tvblZeXpyVLlugXv/iF16MAAAnO8wiFQiFVV1erurra6z8aAJBkeO04AIAZIgQAMEOEAABmiBAAwAwRAgCYIUIAADNECABghggBAMx4/suqCa+3x270Z5+ZzHV/7zaZK0kjonbvHxX6eKTZbKUaful12/19u67LJnN7Po+bzJUk121zn78Y7ozmXvvjKFdCAAAzRAgAYIYIAQDMECEAgBkiBAAwQ4QAAGaIEADADBECAJghQgAAM0QIAGCGCAEAzBAhAIAZIgQAMEOEAABmiBAAwAwRAgCYIUIAADNECABghggBAMwQIQCAGSIEADBDhAAAZogQAMAMEQIAmCFCAAAzRAgAYIYIAQDMECEAgBkiBAAwQ4QAAGZSrRfAv3HOZuzlLpO5kuS6L5vNxjBi9LWFb8aVEADADBECAJghQgAAM0QIAGCGCAEAzBAhAIAZIgQAMEOEAABmiBAAwAwRAgCYIUIAADPXHaHGxkbNnz9feXl5CgQC2rVrV7/PO+f04osvKi8vT+np6ZozZ46OHz/u1b4AgCRy3RHq7OzUbbfdppqamgE/v2HDBm3cuFE1NTVqampSTk6O7rvvPnV0dHzrZQEAyeW6X0V77ty5mjt37oCfc86purpaa9as0cMPPyxJ2rp1q7Kzs1VfX68lS5Z8u20BAEnF058JnT59Wm1tbSorK+s7FgwGNXv2bB04cGDA/yYejysWi/W7AQCGB08j1NbWJknKzs7udzw7O7vvc18ViUQUDof7bvn5+V6uBAAYwnx5dlwgEOj3sXPuimNfWr16taLRaN+ttbXVj5UAAEOQp++smpOTI+mLK6Lc3Ny+4+3t7VdcHX0pGAwqGAx6uQYAIEF4eiVUUFCgnJwcNTQ09B3r6urSvn37NHPmTC9HAQCSwHVfCV28eFGnTp3q+/j06dM6evSoxo4dqwkTJqiiokJVVVUqLCxUYWGhqqqqlJGRoccff9zTxQEAie+6I3To0CHdfffdfR9XVlZKkhYtWqTf/e53+ulPf6rPPvtMTz/9tC5cuKDvf//7evvttxUKhbzbGgCQFALOOWe9xL+LxWIKh8OaoweVGhhpvQ789jVPWAE8NbQe5pJet7usvdqtaDSqzMzMq57La8cBAMwQIQCAGSIEADBDhAAAZogQAMAMEQIAmCFCAAAzRAgAYMbTFzAFrhu/RAgMa1wJAQDMECEAgBkiBAAwQ4QAAGaIEADADBECAJghQgAAM0QIAGCGCAEAzBAhAIAZIgQAMEOEAABmiBAAwAwRAgCYIUIAADNECABghggBAMwQIQCAGSIEADBDhAAAZogQAMAMEQIAmCFCAAAzRAgAYIYIAQDMECEAgBkiBAAwQ4QAAGaIEADATKr1Al/lnJMkdeuy5IyXAQBct25dlvT/j+dXM+Qi1NHRIUnarzeMNwEAfBsdHR0Kh8NXPSfgriVVg6i3t1fnzp1TKBRSIBC47v8+FospPz9fra2tyszM9GHDoWc43meJ+z2c7vdwvM9S4t5v55w6OjqUl5enESOu/lOfIXclNGLECN14443f+s/JzMxMqL80LwzH+yxxv4eT4XifpcS83990BfQlnpgAADBDhAAAZpIuQsFgUC+88IKCwaD1KoNmON5nifs9nO73cLzP0vC430PuiQkAgOEj6a6EAACJgwgBAMwQIQCAGSIEADCTVBF69dVXVVBQoLS0NBUXF+vdd9+1XslXkUhE06dPVygUUlZWlh566CGdOHHCeq1BFYlEFAgEVFFRYb2K786ePauFCxdq3LhxysjI0O23367m5mbrtXzV3d2t559/XgUFBUpPT9ekSZO0du1a9fb2Wq/mqcbGRs2fP195eXkKBALatWtXv8875/Tiiy8qLy9P6enpmjNnjo4fP26zrMeSJkI7duxQRUWF1qxZoyNHjujOO+/U3Llz1dLSYr2ab/bt26fy8nIdPHhQDQ0N6u7uVllZmTo7O61XGxRNTU2qq6vTrbfear2K7y5cuKBZs2Zp5MiRevPNN/XXv/5Vv/nNbzRmzBjr1Xy1fv16bdq0STU1Nfroo4+0YcMGvfzyy3rllVesV/NUZ2enbrvtNtXU1Az4+Q0bNmjjxo2qqalRU1OTcnJydN999/W91mZCc0nie9/7nlu6dGm/Y5MnT3arVq0y2mjwtbe3O0lu37591qv4rqOjwxUWFrqGhgY3e/Zst2LFCuuVfLVy5UpXWlpqvcagmzdvnlu8eHG/Yw8//LBbuHCh0Ub+k+R27tzZ93Fvb6/Lyclxv/rVr/qOff755y4cDrtNmzYZbOitpLgS6urqUnNzs8rKyvodLysr04EDB4y2GnzRaFSSNHbsWONN/FdeXq558+bp3nvvtV5lUOzZs0clJSV65JFHlJWVpWnTpmnz5s3Wa/mutLRU77zzjk6ePClJ+uCDD7R//37df//9xpsNntOnT6utra3f41swGNTs2bOT4vFtyL2A6X/i/Pnz6unpUXZ2dr/j2dnZamtrM9pqcDnnVFlZqdLSUhUVFVmv46vXX39dhw8fVlNTk/Uqg+aTTz5RbW2tKisr9bOf/Uzvv/++nnnmGQWDQT355JPW6/lm5cqVikajmjx5slJSUtTT06N169bpscces15t0Hz5GDbQ49uZM2csVvJUUkToS1996wfn3H/0dhCJaNmyZfrwww+1f/9+61V81draqhUrVujtt99WWlqa9TqDpre3VyUlJaqqqpIkTZs2TcePH1dtbW1SR2jHjh3atm2b6uvrNXXqVB09elQVFRXKy8vTokWLrNcbVMn6+JYUERo/frxSUlKuuOppb2+/4v8ektHy5cu1Z88eNTY2evI2GENZc3Oz2tvbVVxc3Hesp6dHjY2NqqmpUTweV0pKiuGG/sjNzdXNN9/c79iUKVP0+9//3mijwfHcc89p1apVevTRRyVJt9xyi86cOaNIJDJsIpSTkyPpiyui3NzcvuPJ8viWFD8TGjVqlIqLi9XQ0NDveENDg2bOnGm0lf+cc1q2bJn+8Ic/6M9//rMKCgqsV/LdPffco2PHjuno0aN9t5KSEi1YsEBHjx5NygBJ0qxZs654+v3Jkyc1ceJEo40Gx6VLl654U7SUlJSke4r21RQUFCgnJ6ff41tXV5f27duXFI9vSXElJEmVlZV64oknVFJSohkzZqiurk4tLS1aunSp9Wq+KS8vV319vXbv3q1QKNR3JRgOh5Wenm68nT9CodAVP/MaPXq0xo0bl9Q/C3v22Wc1c+ZMVVVV6Qc/+IHef/991dXVqa6uzno1X82fP1/r1q3ThAkTNHXqVB05ckQbN27U4sWLrVfz1MWLF3Xq1Km+j0+fPq2jR49q7NixmjBhgioqKlRVVaXCwkIVFhaqqqpKGRkZevzxxw239ojtk/O89dvf/tZNnDjRjRo1yt1xxx1J/1RlSQPetmzZYr3aoBoOT9F2zrk//vGPrqioyAWDQTd58mRXV1dnvZLvYrGYW7FihZswYYJLS0tzkyZNcmvWrHHxeNx6NU/95S9/GfBredGiRc65L56m/cILL7icnBwXDAbdXXfd5Y4dO2a7tEd4KwcAgJmk+JkQACAxESEAgBkiBAAwQ4QAAGaIEADADBECAJghQgAAM0QIAGCGCAEAzBAhAIAZIgQAMEOEAABm/g889PTJz+A+1AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(cutout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "d6bce5a1-1307-4808-a5d1-b2f88890493b",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.982631206252375"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "aperture_photometry(cutout, 2 * fwhm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "850ce4f8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.02908882086657216"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fwhm/3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "ea721fbf-b428-4738-9dd5-11be003120d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pixell import reproject\n",
    "\n",
    "source_map_healpix = reproject.map2healpix(source_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "ac0823bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "64"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hp.npix2nside(source_map_healpix.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "cdadf0e8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.015989479811663883"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hp.nside2resol(_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "1094cab4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0QAAAICCAYAAADvbw3rAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAy6klEQVR4nO3deZRcdZ3473dVd7rT2SAkkCbBkLA6YRFEVBYJIQkwLjiisii7w+A4R3BDz4xKXM4AKjrqceag54tBlEEjERgcFJIRGXEUGAFlGf2xBCJCCEmAJJCt0/f3R3XXkq7u9FJd232ec3Jyc7e6XVXdfV/5VN3KJEmSBAAAQApla30AAAAAtSKIAACA1BJEAABAagkiAAAgtQQRAACQWoIIAABILUEEAACkliACAABSSxABAACpJYgA6sC1114bmUwmMplM/PKXv+yzPEmS2G+//SKTycTxxx8/rNuYNWtWnHfeefl/P/XUU5HJZOLaa68d1v6G43Of+1xkMplBrbvj8VZLrW4XgNporfUBAFAwceLEuOaaa/pEz1133RVPPPFETJw4sTYHViF/+7d/GyeffHKtD2NAN910U0yaNKnWhwFAlRghAqgjp59+eixdujTWr19fMv+aa66Jo446KmbOnFmjI6uMvfbaK9785jfX+jAGdPjhh8e+++5b68MAoEoEEUAdOfPMMyMi4oYbbsjPe/nll2Pp0qVxwQUXlN1m3bp18aEPfShmzJgRbW1tsc8++8SnP/3p2LJly5Bu+5FHHolMJhM//vGP8/N+97vfRSaTiYMOOqhk3VNOOSWOOOKIknk/+tGP4qijjorx48fHhAkT4qSTTooHHnigZJ1yL5nbtm1bfPKTn4zOzs4YN25cHHvssXHvvfeWPcZVq1bFRRddFHvttVe0tbXF7Nmz4/Of/3x0dXUN+LX9zd/8Tey9997R3d3dZ9mb3vSmeP3rX5//d7mXzK1fvz4+8YlPxOzZs6OtrS1mzJgRH/nIR+KVV17Jr/Pe9763z/30jne8o899ev/990cmk4lbb711wGMGoDoEEUAdmTRpUrznPe+J7373u/l5N9xwQ2Sz2Tj99NP7rL958+aYN29eXHfddfGxj30s/vM//zPOOuus+PKXvxynnnrqkG77oIMOij333DOWL1+en7d8+fLo6OiIRx99NJ599tmIiOjq6oq77rorFixYkF/v8ssvjzPPPDPmzJkTS5Ysie9///uxYcOGeMtb3hKPPvrogLd74YUXxlVXXRXnnHNO3HLLLfHud787Tj311HjxxRdL1lu1alW88Y1vjNtvvz0uu+yy+NnPfhYf+MAH4oorrogLL7xwwNu44IILYuXKlfGLX/yiZP4f//jHuPfee+P888/vd9tXX3015s6dG9/73vfi4osvjp/97GfxqU99Kq699to45ZRTIkmSiIhYsGBBPProo/Hcc8+V3E8dHR2xbNmykvu0tbV12O8FA6DCEgBqbvHixUlEJPfdd19y5513JhGRPPzww0mSJMmRRx6ZnHfeeUmSJMlBBx2UzJ07N7/d1VdfnUREsmTJkpL9felLX0oiIrnjjjvy8/bee+/k3HPPzf97xYoVSUQkixcvzs8766yzkn322Sf/7wULFiQXXnhhMnny5OR73/tekiRJ8utf/7pk3ytXrkxaW1uTD3/4wyXHsGHDhqSzszM57bTT8vMWLVqUFP/q+b//+78kIpKPfvSjJdtef/31SUSUHO9FF12UTJgwIXn66adL1r3qqquSiEgeeeSRpD/btm1Lpk2blrzvfe8rmf/JT34yaWtrS9asWdPv/XTFFVck2Ww2ue+++0q2vfHGG5OISG677bYkSZLk8ccfTyIiue6665IkSZK77747iYjkk5/8ZDJ79uz8dgsXLkyOPvrofo8VgOoyQgRQZ+bOnRv77rtvfPe7342HHnoo7rvvvn5fLveLX/wixo8fH+95z3tK5ve+5Ou//uu/hnTb8+fPjyeffDJWrFgRmzdvjrvvvjtOPvnkmDdvXn6UY/ny5dHe3h7HHntsRETcfvvt0dXVFeecc050dXXl/4wdOzbmzp1b9qp5ve68886IiHj/+99fMv+0006L1tbS6/789Kc/jXnz5sX06dNLbuev//qvIyJ34Yn+tLa2xllnnRU/+clP4uWXX46IiO3bt8f3v//9eOc73xlTpkzpd9uf/vSncfDBB8dhhx1WcrsnnXRSyVUB991335g1a1Z+hG3ZsmVxyCGHxFlnnRUrVqyIJ554IrZs2RJ33313yegaALXlKnMAdSaTycT5558f3/zmN2Pz5s1xwAEHxFve8pay665duzY6Ozv7vC9njz32iNbW1li7du2Qbrv3RH358uUxe/bs2LZtW5xwwgnx/PPPxxe/+MX8smOOOSY6OjoiIuL555+PiIgjjzyy7D6z2f7/7633+Do7O0vmt7a29omU559/Pm699dYYM2ZM2X2tWbNmwK/tggsuiK9+9avxwx/+MC666KK4/fbb47nnnhvw5XK9t/v4448P6nbnz58fP//5zyMidz8tXLgwDjnkkJg2bVosX7489t9//9i0aZMgAqgjggigDp133nlx2WWXxdVXXx3//M//3O96U6ZMiXvuuSeSJCmJotWrV0dXV1dMnTp1SLe71157xQEHHBDLly+PWbNmxRve8IbYddddY/78+fGhD30o7rnnnvjtb38bn//85/Pb9N7GjTfeGHvvvfeQbq83elatWhUzZszIz+/q6uoTc1OnTo1DDz203/tj+vTpA97WnDlz4o1vfGMsXrw4Lrrooli8eHFMnz49TjzxxAG3mzp1anR0dJS8r2vH5b3mz58f11xzTdx7771xzz33xGc+85mIiDjhhBNi2bJl8fTTT8eECRPq/kp7AGkiiADq0IwZM+LSSy+NP/7xj3Huuef2u978+fNjyZIlcfPNN8e73vWu/Pzrrrsuv3yoFixYEEuWLInXvOY18ba3vS0iIg444ICYOXNmXHbZZbFt27aSEY6TTjopWltb44knnoh3v/vdQ7qt3gsLXH/99SVXrVuyZEmfK8e9/e1vj9tuuy323XffmDx58pC/roiI888/P/7+7/8+7r777rj11lvjYx/7WLS0tAy4zdvf/va4/PLLY8qUKTF79uwB150/f35kMpn47Gc/G9lsNo477riIyN2nl156aTz99NNx3HHH9TvaBED1CSKAOnXllVfudJ1zzjkn/vVf/zXOPffceOqpp+KQQw6Ju+++Oy6//PJ461vfOqyXZs2fPz/+7d/+LdasWRNf//rXS+YvXrw4Jk+eXBIvs2bNii984Qvx6U9/Op588sk4+eSTY/LkyfH888/HvffeG+PHjy8ZUSr2V3/1V3HWWWfF17/+9RgzZkwsWLAgHn744bjqqqv6fDjqF77whVi2bFkcffTRcfHFF8eBBx4Ymzdvjqeeeipuu+22uPrqq2OvvfYa8Gs788wz42Mf+1iceeaZsWXLlj6X1y7nIx/5SCxdujSOO+64+OhHPxqHHnpodHd3x8qVK+OOO+6Ij3/84/GmN70pInIvVTz44IPjjjvuiHnz5sW4ceMiIhdE69ati3Xr1sXXvva1nd4mANUjiAAa2NixY+POO++MT3/60/GVr3wlXnjhhZgxY0Z84hOfiEWLFg1rnyeccEJks9no6OiIo446Kj9/wYIFsXjx4pg3b16f9wX94z/+Y8yZMye+8Y1vxA033BBbtmyJzs7OOPLII+ODH/zggLd3zTXXxLRp0+Laa6+Nb37zm3HYYYfF0qVL44wzzihZb88994z//d//jS9+8Yvxla98JZ555pmYOHFizJ49Ox9hO7PLLrvEu971rvj3f//3OOaYY+KAAw7Y6Tbjx4+PX/3qV3HllVfGd77znVixYkV0dHTEzJkzY8GCBTFr1qyS9RcsWBAPPfRQSYzOnDkz9t9//3jssce8fwigzmSSpOcDFAAAAFLGZbcBAIDUEkQAAEBqCSIAACC1BBEAAJBagggAAEgtQQQAAKSWIAIAAFJLEAEAAKkliAAAgNQSRAAAQGq11voAABg9C7PvrfUhNIVl3T+u9SEAMEoySZIktT4IAMoTNM1BUAHUL0EEUCXihqEQUQDVIYgARkDkUA/EE8DwCSKAfogdmoloAihPEAGpJXigQDABaSWIgKYkdqDyRBPQjAQR0LBED9QPsQQ0KkEE1DXRA41PLAH1TBABdUH4QPoIJaAeCCKgqoQPsDNCCagmQQSMGvEDVIpIAkaLIAJGTPgAtSKUgJESRMCQiB+g3okkYCgEEdAv8QM0C5EE9EcQAREhfoD0EUlAhCCC1BJAAKUEEqSTIIIUED8AwyOSoPkJImgy4gdgdIkkaC6CCBqcAAKoLYEEjU0QQYMRQAD1TSBBYxFEUOcEEEBjE0hQ3wQR1BkBBNDcBBLUF0EENSaAANJNIEFtCSKoMgEEwEAEElSXIIIqEEEADIc4gtEniGAUCCAARoNAgsoTRFAhIgiAahJHUBmCCEZABAFQD8QRDJ8ggiEQQAA0AoEEgyeIYCdEEACNTBzBwAQRlCGCAGhG4gj6EkTQQwQBkCbiCHKytT4AAACAWjFCRKoZFQIAo0WkmyAidUQQAPRPHJE2gohUEEEAMHTiiDQQRDQtEQQAlSOOaFaCiKYhgACgegQSzUIQ0fCEEADUjjCi0QkiGpIIAoD6I45oRIKIhiGCAKBxiCMahSCi7gkhAGhcwoh6J4ioSyIIAJqPOKIeCSLqihACgOYnjKgngoiaE0EAkF7iiFoTRNSMEAIAegkjakUQUVUiCADYGXFENQkiqkIIAQBDJYyoBkHEqBJCAMBICSNGkyCi4kQQADBaxBGVJoioGCEEAFSLMKJSBBEjJoQAgFoRRoyUIGLYhBAAUC+EEcMliBgyIQQA1CthxFAJIgZNCAEAjUIYMViCiAGJIACg0YkjBiKIKEsIAQDNRhhRjiCihBACAJqdMKKYICIihBAAkD7CiAhBlHpCCABIO2GUboIopYQQAEApYZRO2VofANUnhgAA+nKOlE5GiFLENzkAwOAYLUoPQZQCQgh6ZDK5v3t/7GUyhWkAKEMYNT9B1MSEEBTpjaFyBBIAOyGMmpcgakJCiNQqHgHKtkR0b8/Nbm2NZPv2yLS05FdNtm+PyPS8jbJ7e9/RIwAoQxg1H0HURIQQqVY8AtQTOtmx7ZFs3Zqb1dERyaZNhVXa26N70+bc9JjW/HqCCIDBEEbNQxA1ASFE6mV7Rn6S7si0jomIXORERCRbt0amoyM33RNEmfb2iIjo3rS5ZD0jRgAMhShqDoKogQkhUqv4vT49MZQdPy6SzVty07vtGsmGjbnpKbtF90sv56Z33SUiIrrXvZjbzYTx+WXJtq7cvJaWSLq2FW7Lj0gABkEcNS5B1ICEEKlWNHKTGdOWmzW2vbB4XEdhui23vPull/MxFJELosyE8bndvFp4GV1vUEWEKAJgWIRR4xFEDUQIkXpF7xPKtrdHsr07srvuEsnmnvcC7dUZsS434tO1z54x5pm1ERGxddbu0fbUC/ltt0/bNbJ/Xh2Z8eOi+4XcOpnW1uje+Ep+HUEEwEgIo8YhiBqEGCL1doihiIjMxImF5VN3zU927TouP5205t4X1LKpKHAiomXdxsI6L76cny4bRX5MAjAMoqgxCKI6J4Qgcu8TSrpzk+PGRWzfHtlpu0dERPLq5th24IwYszYXMi8ePiUmrsi9DG7NYeNiysOb46X9xkZExJSHNkRExJapuZfVjfv/XohkTGvE6rX5m8oHUdKduzS3H5EAjJAwqm+CqE4JIeiRLXx2ULYjFzbZqbvl5219zZT89MaZY/PTW3bJjSi1bC7sasJzXSW77ljxYuEfPVGUdHXlr0aXdJWuDwAjIYzqU7bWB0BfYgh6FMfQ+NzL4DIzOiNpz10sYcNhe8aWKbnpZ99SeJncquNyo0kvvm57rHlT7sNZX+3MxOrDx8Tqw3OX5d7eXvS5RT3/L1QcQEl3UnL7ADBSzvHqkxGiOuKbBIqUi6Hp0/LxsvGvCiNDa+e05qdf3TsXNace+b/xk/veEBER454uLI+ImPJoIXwm/F/PRRU2bclfgrv7lVcLK3dvH/GXAgA7MlpUP4wQ1QkxBEX6i6EeLx2+e3SNzf34+su8bGyelhsR2vfkJyMiF0O9pu71Uow7Zk1ERCStuT9rDm2N1k3d0bopt11mU+Fy27FtW2TbcqNIkXSXXMwBACrFuV/9EEQ1tjD7Xt8QUKxcDHXunp/38uum5qf/Mq/wI6w3hv7nrV/Lz/vG/B/kp5OiQaKxawoD490TC+87im3Fl9ru7vnbIDoAo8N5YH0QRDXkGwB2kMkUribXMTaiu7s0hg4tvEzuubnd0d2WW/f8E34Zb5z8VHx79o355XMn/TE//eqWMZE57OUYuyYpiaG2tUUfyrrxlch05K4+l3R1RaalRQwBUBXOCWurdeerUGme9FBG8ecM9b5krcir+0zOT79weGHd80/4ZUREfGbqH6P3InJf3fP++I9XcqNLr27J7WvrY5OivWh/3S3lXwrXe2GFZLv3DgFQPb3nh95bVH1GiKpMDEEZZT50NSIiu9vkyGzeGq/uVxgZWn1E4cfWO990f6zZNiE+M7UwGrRn64T89C/Wz4mFs/4UWx+bFBERG/bO/Zn0VOGiCi2rcxdSSDZtimTTpsi0toohAGrGuWL1GSGqEk9u2LlMS+H9Q5mezxza3lkYGXqlsyUici9jmzg99yGrC3d5uOy+Ht70moiIuOW+1+d/0E39felL4LIbcxdT6L26XERE99ZtAQC1ZLSouowQVYEYgp1IktIYmjA+IgoxNHbN5nhpv8LL6Dpe+1JEFGLoP18dG891bcwvv3zNgRGRi6GIiK5dt/eJoXFP5j6UNR9DY8aUxpD3DwFQY84hq0MQjSJXDoFBKv4A1KIwalmzPiIiNu9RuBJcy5bCy+s+/6d3RETE7S8dkp93ymMn99n9Hr8u/YDVjr/kRpdi9dr8vGRT4QILYgiAeuF8cvT5YNZR4okLg1R8me2el8llOsZGZmLuvUCvzNkjv3zVG3OjRG2vezHGjsm9D+jN057KL3/qldx7jd48eUX8v/8+PiIidr+n8P8+u/3+xfx0ZktXPoi6N+QCKelOfBArAHXLS+hGhxGiClPxMEQ9AZLJFl31bXvuctrd4zui4+lcrHR1FH5cdd2XeyndCysL7y+69b7D89M/+NH8iIiY/svCLnd7oBBDEdEnhoqPBQDqkfPM0WGEqII8QWF4Mq2F67tkJ07Mzdt1UnSPz30u0IbX7pJfvu7A3IjSpr26YkeHzFkZj9+xT0RE7PanQtxM/GPP+4R6mirz5+cjojSGei+3DQCNwGhR5RghqhAxBMOXv8x10fuHko2v9lmv9ZXu/PTuv82tu9sDhW1e+PasiIiY/uvNMXZN7gIJkx5cnV+eWflcZFY+FxE7xJDLbAPQYJx7Vo4RohHyZIQR6vkMokxbW35WdlLuc4OSaYXPH9r0mon56Vf3yI0obS/6pNW2jbkfZS/Pzsb0X2/OzXum6GVya9blJ7s3vpKfzseQH4UANCijRSNjhGgExBCMoq6uyPwl99K27CuFK8B1PJeLmV0fK4wgTXkoNz359y/GrJtz7w1qf2xVfnny/JqIyIWQGAKg2TgnHRkjRMPkiQcVkilcTKF3lCg7blxh3q650aLuXcbn520fV/hMop7PaY3W9ZsL27zUc9W4jvZ8DJVcVjvEEADNx0jR8AiiIRJCUGE7BtH27fkLKyRbt0Z2j6mRbMh96Gqy17TIvvBSbNt794iIGPPntbFtrykx5qncSFKy2y4Rz66OzLjcxRi6167LfeDqK6+WXMUu6erK3a4ffwA0IWE0NIJoCMQQjILiICq6qEKmPfcGoUzH2D7zIiIi2/OK36KrwyWbt+T+LhoN6t66rbCNy2oDkBKiaPC8h2iQxBCMkt7/k+n5O+lOch+QGhFJUcwkr+Yip/ull6P7pdxltLvXFS6a0P3iS7m/N2zIX0K7e8uWwu2IIQBSxLnr4BkhGgRPKKiSbNEIUfFo0ZiezykqGk0q0fN+oOLPEiq5lLYfcwCklJGinRNEAxBCUAPZloikOzKtuQsnJNu3R3Zse3T3vBwu2zE2ki1b8i+f6960ObJtY/KjQZmWlsJ7hCLEEACEMBqIIOqHGIIaKQ6ZohGj/OKeiyMkXV2F5UnhA1sFEACUJ4rK8x6iMsQQ1Fhx1HRvL7z/pyd8il8aF93b+7wPCQDoyzlueUaIduCJAnXKS+AAoCKMFJUyQgQAAKSWICpidAjqnNEhABgx57ylBFEPTwyoc2IIACrGuW9B6t9D5MkAAECapf09RakeIRJDAACkXdrPiVMbRGl/4AEAoFeaz41TGURpfsABAKCctJ4jpy6I0vpAAwDAzqTxXDlVQZTGBxgAAIYibefMqQmitD2wAAAwXGk6d05FEKXpAQUAgEpIyzl00wdRWh5IAACotDScSzd1EKXhAQQAgNHU7OfUTRtEzf7AAQBAtTTzuXVTBlEzP2AAAFALzXqO3XRB1KwPFAAA1Foznms3VRA14wMEAAD1pNnOuZsmiJrtgQEAgHrVTOfeTRFEzfSAAABAI2iWc/CGD6JmeSAAAKDRNMO5eEMHUTM8AAAA0Mga/Zy8oYMIAABgJDJJkiS1PoihavQKBQCAZrSs+8e1PoQha7gRIjEEAAD1qRHP1RsqiBrxDgYAgDRptHP2hgmiRrtjAQAgrRrp3L0hgqiR7lAAAKBxzuEbIogAAABGQ90HUaOUJQAAUKoRzuXrOoga4Q4EAAD6V+/n9HUbRPV+xwEAAINTz+f2dRtEAAAAo60ug6ieCxIAABi6ej3Hr7sgqtc7CgAAGJl6PNevqyCqxzsIAAConHo756+rIAIAAKimugmieitFAABgdNTTuX9dBFE93SEAAMDoq5cGqIsgAgAAqIWaB1G9lCEAAFBd9dACNQ2iergDAACA2ql1E9R8hAgAAKBWahZEtS5BAACgPtSyDWoSRGIIAAAoVqtG8JI5AAAgtaoeREaHAACAcmrRCkaIAACA1KpqEBkdAgAABlLtZjBCBAAApFbVgsjoEAAAMBjVbAcjRAAAQGpVJYiMDgEAAENRrYYY9SASQwAAwHBUoyW8ZA4AAEgtQQQAAKTWqAaRl8sBAAAjMdpNYYQIAABIrVELIqNDAABAJYxmWxghAgAAUmtUgsjoEAAAUEmj1RhGiAAAgNSqeBAZHQIAAEbDaLSGESIAACC1BBEAAJBaFQ0iL5cDAABGU6WbwwgRAACQWoIIAABIrYoFkZfLAQAA1VDJ9jBCBAAApJYgAgAAUksQAQAAqVWRIPL+IQAAoJoq1SBGiAAAgNQSRAAAQGqNOIi8XA4AAKiFSrSIESIAACC1BBEAAJBagggAAEitEQWR9w8BAAC1NNImMUIEAACkliACAABSSxABAACpNewg8v4hAACgHoykTYwQAQAAqSWIAACA1BJEAABAagkiAAAgtQQRAACQWsMKIleYAwAA6slwG8UIEQAAkFqCCAAASC1BBAAApJYgAgAAUksQAQAAqSWIAACA1BJEAABAagkiAAAgtQQRAACQWoIIAABILUEEAACkliACAABSSxABAACpJYgAAIDUEkQAAEBqCSIAACC1BBEAAJBagggAAEgtQQQAAKSWIAIAAFJLEAEAAKkliAAAgNQSRAAAQGoJIgAAILUEEQAAkFqCCAAASC1BBAAApNawgmhZ948rfRwAAADDNtxGMUIEAACkliACAABSSxABAACpJYgAAIDUEkQAAEBqDTuIXGkOAACoByNpEyNEAABAagkiAAAgtQQRAACQWiMKIu8jAgAAammkTWKECAAASC1BBAAApJYgAgAAUmvEQeR9RAAAQC1UokWMEAEAAKkliAAAgNSqSBB52RwAAFBNlWoQI0QAAEBqCSIAACC1BBEAAJBaFQsi7yMCAACqoZLtYYQIAABILUEEAACkVkWDyMvmAACA0VTp5jBCBAAApJYgAgAAUqviQeRlcwAAwGgYjdYwQgQAAKTWqASRUSIAAKCSRqsxjBABAACpNWpBZJQIAACohNFsCyNEAABAao1qEBklAgAARmK0m8IIEQAAkFqCCAAASK1RDyIvmwMAAIajGi1RlREiUQQAAAxFtRrCS+YAAIDUqloQGSUCAAAGo5rtYIQIAABIraoGkVEiAABgINVuBiNEAABAalU9iIwSAQAA5dSiFYwQAQAAqVWTIDJKBAAAFKtVI9RshEgUAQAAEbVtAy+ZAwAAUqumQWSUCAAA0q3WTVDzEaJa3wEAAEBt1EML1DyIAAAAaqUugqgeyhAAAKieemmAugiiiPq5QwAAgNFVT+f+dRNEAAAA1VZXQVRPpQgAAFRevZ3z11UQRdTfHQQAAFRGPZ7r110QRdTnHQUAAAxfvZ7j12UQAQAAVEPdBlG9FiQAADA09XxuX7dBFFHfdxwAALBz9X5OX9dBFFH/dyAAAFBeI5zL130QAQAAjJaGCKJGKEsAAKCgUc7hGyKIIhrnDgUAgLRrpHP3hgmiiMa6YwEAII0a7Zy9oYIoovHuYAAASItGPFfPJEmS1Poghmth9r21PgQAAEi9RgyhXg03QgQAAFApDR1EjVyiAADQDBr9nLyhgyii8R8AAABoVM1wLt7wQRTRHA8EAAA0kmY5B2+KIIpongcEAADqXTOdezdNEEU01wMDAAD1qNnOuZsqiCKa7wECAIB60Yzn2k0XRBHN+UABAEAtNes5dlMGUUTzPmAAAFBtzXxu3bRBFNHcDxwAAFRDs59TN3UQRTT/AwgAAKMlDefSTR9EEel4IAEAoJLScg6diiCKSM8DCgAAI5Wmc+fUBFFEuh5YAAAYjrSdM6cqiCLS9wADAMBgpfFcOXVBFJHOBxoAAAaS1nPkVAZRRHofcAAA2FGaz41TG0QR6X7gAQAgwjlxJkmSpNYHUQ8WZt9b60MAAICqSXsI9Ur1CFExTwgAANLCuW+BIAIAAFJLEBVRygAANDvnvKUE0Q48QQAAaFbOdfsSRGV4ogAA0Gyc45YniPrhCQMAQLNwbts/l90eBJfkBgCgEQmhnTNCNAieSAAANBrnsIMjiAbJEwoAgEbh3HXwBNEQeGIBAFDvnLMOjfcQDZP3FQEAUE+E0PAYIRomTzgAAOqFc9PhE0Qj4IkHAECtOScdGS+ZqxAvoQMAoJqEUGUYIaoQT0gAAKrFuWflCKIK8sQEAGC0OeesLC+ZGyVeQgcAQCUJodFhhGiUeMICAFApzi1HjxGiKjBaBADAcAih0WeEqAo8kQEAGCrnkNVhhKjKjBYBADAQIVRdRoiqzBMcAID+OFesPiNENWS0CACACCFUS0aIasgTHwAA54S1ZYSoThgtAgBIFyFUH4wQ1QnfEAAA6eHcr34YIapDRosAAJqTEKo/RojqkG8UAIDm4xyvPhkhqnNGiwAAGpsQqm+CqEEIIwCAxiKEGoOXzDUI31AAAI3DuVvjMELUgIwWAQDUJyHUeARRAxNGAAD1QQg1LkHUBIQRAED1iaDmIIiaiDACAKgOMdQ8BFETEkYAAKNDCDUfQdTEhBEAQGUIoeYliFJAGAEADI8Qan6CKEWEEQDA4Aih9BBEKSSMAADKE0Lpk631AVB9vtEBAPpyjpRORohSzmgRAJB2QijdBBERIYwAgPQRQkQIInYgjACAZieEKCaIKEsYAQDNRghRjiBiQMIIAGh0QoiBCCIGTRwBAI1CBDFYgoghE0YAQL0SQgyVIGLYhBEAUC+EEMMliBgxYQQA1IoQYqQEERUjjACAahFCVIogouKEEQAwWoQQlSaIGFXiCAAYKRHEaBJEVIUwAgCGSghRDYKIqhJGAMDOCCGqSRBRM+IIAOglgqgVQUTNCSMASC8hRK0JIuqKOAKA5ieCqCeCiLokjACg+Qgh6pEgou6JIwBoXCKIeieIaBjCCAAahxCiUQgiGpI4AoD6I4JoRIKIhieOAKB2RBCNThDRNIQRAFSPEKJZCCKalkACgMoRQDQrQUQqiCMAGDoRRBoIIlJHHAFA/0QQaSOISDVxBAAiiHQTRNBDHAGQJiIIcrK1PgAAAIBaMUIEZRgtAqAZGRWCvgQR7IQ4AqCRiSAYmCCCIRBHADQCEQSDJ4hgBAQSAPVAAMHwCSKoEHEEQDWJIKgMQQSjQBwBMBpEEFSeIIIqEEgADIcAgtEniKDKxBEAAxFBUF2CCGpMIAGkmwCC2hJEUGcEEkBzE0BQXwQR1DmBBNDYBBDUN0EEDUYgAdQ3AQSNRRBBgxNIALUlgKCxCSJoMgIJYHQJIGgugghSQCQBDI/4geYniCClRBJAKfED6SSIgIgQSED6CCAgQhABAxBJQLMQP0B/BBEwJCIJqHfiBxgKQQSMmEgCakX8ACMliIBRI5SAShE+wGgRREBViSRgZ8QPUE2CCKgLQgnSR/gA9UAQAXVNKEHjEz5APRNEQMMSS1A/RA/QqAQR0JTEElSe6AGakSACUks0QYHYAdJKEAH0QzDRTAQPQHmCCGAERBP1QOwADJ8gAqgS8cRQiByA6hBEAHVMRDUHcQNQvwQRQBMTVJUhaACalyACAABSK1vrAwAAAKgVQQQAAKSWIAIAAFJLEAEAAKkliAAAgNQSRAAAQGoJIgAAILUEEQAAkFqCCAAASC1BBAAApJYgAgAAUksQAQAAqSWIAACA1BJEAABAagkiAAAgtQQRAACQWoIIAABILUEEAACkliACAABSSxABAACpJYgAAIDUEkQAAEBqCSIAACC1BBEAAJBagggAAEit1lofADSjzZs3x9atW2t9GAA0mba2thg7dmytDwOaiiCCCtu8eXPs0jE5tsbmWh8KAE2ms7MzVqxYIYqgggQRVNjWrVtja2yOY+Ot0Zppj0w2k1uQyRZN9/ydzeSnM9ls0fxsYb3e5Zlsbv3i7TOZHdaNovk7rrvj9v0fS5LJFF5QW3xbAy7vmd8zXbLujvOyxcuL9tMzXbz/pGR+Yd0keqcj/3UVlu+wbtH8/HrFt987P1t++7yS7fuZLnNflNx+2XXLTEc/y3c4ln6372/eTm6/V7/zyuyr3H0RmWTQxxKZ3FeV/9r6rJuUv83i+fnbL8zLDLR9JEXfAoVbzvSzfe/8kn0WbZ/ZYfuib5fIlmxfWJ4tmpft+eqL99O7bnaH6YiIbJTOy5aZ7t1Xf8t7bzM3r7twW7Hj8u5oKdqmsG5u/y2RRKZo+8K6RfOKp3vW7b2dlkx3fp8tPbfXu9/8bZXZV0umO3+MLUXr9X4bt0Txfnu3SfL7ym3fu11hPy1FX39L0bH0Pi75fUXhvuxdnpsXhfuq91gyES09j0hhXiay+XmF6ZZM77xs0bzc9PoN3bH3EU/F1q1bBRFUkCCCUdIaY6I1MyYyRUFSPJ37u3DWlMkUBVHx8mzR8p0GUabPdv0GUWaH5SMOokwhAkqCpxAGFQ2i4hPjUQ6i8tv3M53fPlM0XbT9ToJkZxEx6kFUZnmvHYNooPtlWEFUdFtVD6Li6TLb7xhEO25TeLoUTpyHFURlthluEJUGz+CDqHh+7u/+gqg4AoYeRNmyQZT0Mz34IGrJH1cmsj13aG+M5IKodzpTFClJ0bykaF+RP5bC7RfmDRRELUMIopZBBZG3fcNo8d0FAACkliACAABSSxABAACpJYgAAIDUEkQAAEBqCSIAACC1BBEAAJBagggAAEgtQQQAAKSWIAIAAFJLEAEAAKkliAAAgNRqrfUBQLPqim0RSTYySaZnTvF0z99JJj+dSbJF87OF9bp7lmeKlmeK/84WpvO7z5RZd8ftd1ieZPLTSSYTkcQO2+9sec/8iIjuKF2390vunZctXl60n/yhZPJfVlIyv7Bu781HJiKyO+5rh3XL3u1F25S520uOK8pt38/0jvdFpnR++XXLTEc/y3c4ln6372/eTm6/V7/zyuyr3H0RmWTQxxKZ3FeV/9r6rJuUv83i+fnbL8zLDLR9JEXfAoVbzvSzfe/8kn0WbZ/ZYfuib5dISrYvLE+K5iU9X33vfrqLtsnuMB0RkY3Sedky0/mfPP0sz0bxvO7CbcWOy7ujpWibwrq5/bdEEpmi7QvrFs0rnu5Zt/d2WjLd+X229Nxe737zt1VmXy2Z7vwxthSt1/tt3BLF++3dJsnvK7d973aF/bQUff0tRcfS+1jl9xWF+7J3eW5eFO6r3mPJRLT0PCKFeZnI5ucVplsyhfUK83L7W7+hO4DKE0RQYUmSxIQJE+LujbflzvC21/qIAGgWEyZMiCRJdr4iMGiCCCosk8nExo0b489//nNMmjSp1ocDQJNYv359vOY1r4lM8cg1MGKCCEbJpEmTBBEAQJ1zUQUAACC1BBEAAJBagggqrL29PRYtWhTt7e21PhQAmojfLzA6MolLlQAAACllhAgAAEgtQQQAAKSWIAIAAFJLEAEAAKkliAAAgNQSRDS1JEnic5/7XEyfPj06Ojri+OOPj0ceeWSn2y1dujTmzJkT7e3tMWfOnLjpppv6XfeKK66ITCYTH/nIR/pd56KLLopMJhNf//rXS+Yff/zxkclkSv6cccYZJevcf//9sXDhwth1111jypQp8Xd/93excePGknUuueSSOOKII6K9vT0OO+ywnX59AFTOf//3f8c73vGOmD59emQymbj55ptLlv/kJz+Jk046KaZOnRqZTCYefPDBkuXr1q2LD3/4w3HggQfGuHHjYubMmXHxxRfHyy+/XLLeKaecEjNnzoyxY8fGnnvuGWeffXY8++yzAx7bYH4PDuZ3ETQzQURT+/KXvxxf+9rX4lvf+lbcd9990dnZGQsXLowNGzb0u81vfvObOP300+Pss8+O3//+93H22WfHaaedFvfcc0+fde+77774zne+E4ceemi/+7v55pvjnnvuienTp5ddfuGFF8Zzzz2X//Ptb387v+zZZ5+NBQsWxH777Rf33HNP/PznP49HHnkkzjvvvJJ9JEkSF1xwQZx++uk7uUcAqLRXXnklXve618W3vvWtfpcfc8wxceWVV5Zd/uyzz8azzz4bV111VTz00ENx7bXXxs9//vP4wAc+ULLevHnzYsmSJfGnP/0pli5dGk888US85z3vGfDYBvt7cKDfRdD0EmhS3d3dSWdnZ3LllVfm523evDnZZZddkquvvrrf7U477bTk5JNPLpl30kknJWeccUbJvA0bNiT7779/smzZsmTu3LnJJZdc0mdfzzzzTDJjxozk4YcfTvbee+/kX/7lX0qW97ddr29/+9vJHnvskWzfvj0/74EHHkgiInnsscf6rL9o0aLkda97Xb/7A2B0RURy0003lV22YsWKJCKSBx54YKf7WbJkSdLW1pZs27at33VuueWWJJPJJFu3bi27fLC/B3f2uwianREimtaKFSti1apVceKJJ+bntbe3x9y5c+N//ud/+t3uN7/5Tck2EREnnXRSn23+4R/+Id72trfFggULyu6nu7s7zj777Lj00kvjoIMO6vf2rr/++pg6dWocdNBB8YlPfKLkf+22bNkSbW1tkc0WvlU7OjoiIuLuu+/ud58ANLaXX345Jk2aFK2trWWXr1u3Lq6//vo4+uijY8yYMWXXGcrvwYF+F0GzE0Q0rVWrVkVExLRp00rmT5s2Lb+sv+12ts0Pf/jDuP/+++OKK67odz9f+tKXorW1NS6++OJ+13n/+98fN9xwQ/zyl7+Mz372s7F06dI49dRT88tPOOGEWLVqVXzlK1+JrVu3xosvvhj/9E//FBERzz33XL/7BaBxrV27Nr74xS/GRRdd1GfZpz71qRg/fnxMmTIlVq5cGbfccku/+xns78Gd/S6CZieIaBrXX399TJgwIf9n27ZtERGRyWRK1kuSpM+8HQ20zZ///Oe45JJL4gc/+EGMHTu27Pa/+93v4hvf+EZce+21A97WhRdeGAsWLIiDDz44zjjjjLjxxhtj+fLlcf/990dExEEHHRTf+9734qtf/WqMGzcuOjs7Y5999olp06ZFS0vLwHcIAA1n/fr18ba3vS3mzJkTixYt6rP80ksvjQceeCDuuOOOaGlpiXPOOSeSJBlwnzv7Pbiz30XQ7AQRTeOUU06JBx98MP9n6tSpERF9RoNWr17d53/LinV2dg64ze9+97tYvXp1HHHEEdHa2hqtra1x1113xTe/+c1obW2N7du3x69+9atYvXp1zJw5M7/O008/HR//+Mdj1qxZ/d7261//+hgzZkw89thj+Xnve9/7YtWqVfGXv/wl1q5dG5/73OfihRdeiNmzZw/1LgKgjm3YsCFOPvnkmDBhQtx0001lXwo3derUOOCAA2LhwoXxwx/+MG677bb47W9/W3Z/nZ2dETH034PlfhdBMxNENI2JEyfGfvvtl/8zZ86c6OzsjGXLluXX2bp1a9x1111x9NFH97ufo446qmSbiIg77rgjv838+fPjoYceKomvN7zhDfH+978/HnzwwWhpaYmzzz47/vCHP5SsM3369Lj00kvj9ttv7/e2H3nkkdi2bVvsueeefZZNmzYtJkyYED/60Y9i7NixsXDhwqHeRQDUqfXr18eJJ54YbW1t8R//8R/9vgKhWO/I0JYtW8ounz179rB+Dw70uwiaUfl36kET6P1soMsvvzz233//2H///ePyyy+PcePGxfve9778euecc07MmDEj/36gSy65JI477rj40pe+FO985zvjlltuieXLl+cvYjBx4sQ4+OCDS26r9/XcvfOnTJkSU6ZMKVlnzJgx0dnZGQceeGBERDzxxBNx/fXXx1vf+taYOnVqPProo/Hxj388Dj/88DjmmGPy233rW9+Ko48+OiZMmBDLli2LSy+9NK688srYdddd8+s8/vjjsXHjxli1alVs2rQp/xkXc+bMiba2tsrcoQCUtXHjxnj88cfz/16xYkU8+OCDsdtuu8XMmTNj3bp1sXLlyvxnBv3pT3+KiNwITmdnZ2zYsCFOPPHEePXVV+MHP/hBrF+/PtavXx8REbvvvnu0tLTEvffeG/fee28ce+yxMXny5HjyySfjsssui3333TeOOuqo/G2/9rWvjSuuuCLe9a53Der34GB/F0FTq+k17mCUdXd3J4sWLUo6OzuT9vb25LjjjkseeuihknXmzp2bnHvuuSXzfvzjHycHHnhgMmbMmOS1r31tsnTp0gFvZzCXLN3xstsrV65MjjvuuGS33XZL2trakn333Te5+OKLk7Vr15Zsd/bZZ+fXOfTQQ5Prrruu7O1HRJ8/K1asGPCYABi5O++8s+zP4N7fLYsXLy67fNGiRQNuX/xz/A9/+EMyb968ZLfddkva29uTWbNmJR/84AeTZ555puRYIiJZvHhx/t87+z042N9F0MwySbKTd+IBAAA0Ke8hAgAAUksQAQAAqSWIAACA1BJEAABAagkiAAAgtQQRAACQWoIIAABILUEEAACkliACAABSSxABAACpJYgAAIDU+v8Beg8IB4C9K94AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 850x540 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hp.mollview(source_map_healpix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "5361935b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAJQCAYAAAANCFe9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAACCP0lEQVR4nO29eZgU1b3//66eGWYGBpFhBxHBIIhKRFGCG4vgQjR+UTQGEzHml5h4feI1aoz7jsQluZpco/cm0XhjjKK4a2SXJKBCAgZQjCiIssqgItvATNfvj6Gb6p5eTlWdrarer+fhYaanp8+Zme6qV78/n3PKcV3XBSGEEEJIGVKmJ0AIIYSQaEBpIIQQQogQlAZCCCGECEFpIIQQQogQlAZCCCGECEFpIIQQQogQlAZCCCGECEFpIIQQQogQlAZCCCGECEFpIEQD//rXv/C9730PBx98MGpra1FbW4v+/fvjkksuwaJFi0xPzwiO4+CWW26J/ZiExIlK0xMgJO48/PDDuOyyyzBgwABcfvnlOOyww+A4Dt5991088cQTOOaYY7By5UocfPDBpqeqlQULFuCAAw6I/ZiExAmH154gRB1///vfcdJJJ+HrX/86nn76abRp06bVfaZOnYrjjz8ePXv2NDBDQggRh+UJQhQyefJkVFRU4OGHHy4oDABw7rnn5gjDRRddhLq6OqxcuRLjxo1DXV0devfujSuvvBKNjY0537tlyxZceuml6NWrF9q0aYN+/frh+uuvb3U/x3Fw2WWX4ZFHHsGAAQNQW1uLoUOH4o033oDrurjnnnvQt29f1NXVYfTo0Vi5cmWref7+97/HV7/6VdTU1KC+vh7jx4/Hu+++m3MfP3MvVCpYu3YtfvCDH6B3795o06YNevbsiQkTJmDjxo1Ff8dDhgzBiSee2Or25uZm9OrVC2effXbJMTds2IBLLrkEBxxwANq0aYO+ffvi1ltvRVNTU/Y+xxxzDL7+9a/nfN8RRxwBx3GwcOHC7G3Tpk2D4zhYunRp0fkSEmlcQogSmpqa3NraWnf48OG+vm/SpElumzZt3EMPPdS999573ZkzZ7o33XST6ziOe+utt2bvt3PnTnfw4MFuu3bt3HvvvdedPn26e+ONN7qVlZXuuHHjch4TgNunTx/3uOOOc6dNm+Y+++yz7iGHHOLW19e7V1xxhXvWWWe5L730kvv444+73bp1cwcPHuym0+ns90+ePNkF4H7rW99yX375Zfexxx5z+/Xr53bo0MH997//7XvumTndfPPN2c8/+eQTt0ePHm7nzp3dX/ziF+7MmTPdJ5980r344ovdd999t+jv6/7773cB5MzDdV33lVdecQG4L7zwQtEx169f7/bu3dvt06eP+/DDD7szZ850b7/9dre6utq96KKLsvf72c9+5tbV1bm7d+92Xdd1N2zY4AJwa2tr3TvvvDN7vx/96Edut27dis6VkKhDaSBEEZkTy/nnn9/qa01NTe6ePXuy/7wn6EmTJrkA3Keeeirne8aNG+cOGDAg+/lDDz1U8H4///nPXQDu9OnTs7cBcLt37+5u27Yte9tzzz3nAnCPPPLInPH/67/+ywXg/utf/3Jd13U/++wzt7a2tpWIrFmzxq2urnYnTpzoe+6ZOXlP4BdffLFbVVXlvvPOO61+X6XYvHmz26ZNG/e6667Luf28885zu3Xr5u7Zs6fomJdccolbV1fnfvTRRznfe++997oA3OXLl7uu67ozZ850Abjz5s1zXdd1//jHP7rt27d3L730UnfUqFHZ7+vfv3/O74OQuMHyBCEGOProo1FVVZX9d9999+V83XEcnHnmmTm3DR48GB999FH289mzZ6Ndu3aYMGFCzv0uuugiAMCsWbNybh81ahTatWuX/fzQQw8FAJx++ulwHKfV7ZmxFixYgJ07d2YfN0Pv3r0xevToVuOIzL0Qr776KkaNGpUdX5ROnTrhzDPPxB/+8Aek02kAwGeffYbnn38eF154ISori/d7v/TSSxg1ahR69uyJpqam7L/TTz8dAPD6668DAI4//njU1NRg5syZAIAZM2Zg5MiROO200zB//nzs2LEDH3/8Md5//32MGTPG1/wJiRKUBkIU0blzZ9TW1hY8Wf7pT3/CwoUL8cILLxT83rZt26KmpibnturqauzatSv7eUNDA7p3755zwgeArl27orKyEg0NDTm319fX53ye6bEodntmrMzj9OjRo9U8e/bs2WockbkX4tNPPw28suHiiy/G2rVrMWPGDADAE088gcbGxlaik8/GjRvx4osv5ghcVVUVDjvsMADA5s2bAQA1NTU4/vjjs9Iwa9YsjB07FiNHjkRzczP++te/ZsemNJA4wyWXhCiioqICo0ePxvTp07F+/fqck+6gQYMAAKtXrw78+J06dcKbb74J13VzxGHTpk1oampC586dAz92/jgAsH79+lZfW7dunbRxunTpgk8++STQ95566qno2bMnHnnkEZx66ql45JFHMGzYsOzvuRidO3fG4MGDceeddxb8urdB9eSTT8ZNN92Et956C5988gnGjh2L9u3b45hjjsGMGTOwbt06HHLIIejdu3egn4GQKMCkgRCFXHvttWhubsYPf/hD7NmzR+pjn3zyydi2bRuee+65nNsfe+yx7NdlMHz4cNTW1uKPf/xjzu2ffPIJZs+eLW2c008/HXPmzMF7773n+3srKirwne98B8899xz++te/YtGiRbj44ovLft8ZZ5yBZcuW4eCDD8bQoUNb/fNKw5gxY9DU1IQbb7wRBxxwAAYOHJi9febMmZg9ezZTBhJ7KA2EKOT444/Hf//3f+OVV17BUUcdhV/96leYPXs25s6diyeeeAJXXXUVAGC//fbz/dgXXnghBg8ejEmTJuGXv/wlZs6ciVtuuQXXXXcdxo0bJ+0Etv/+++PGG2/ECy+8gAsvvBCvvvoq/vjHP2LUqFGoqanBzTffLGWc2267DZ07d8ZJJ52E+++/H7Nnz8a0adPwgx/8ACtWrCj7/RdffDEaGxsxceJE1NbW4pvf/KbQmFVVVTjuuOPwm9/8BrNnz8Yrr7yCBx98EGeccUZO8nH00UejY8eOmD59OsaOHZu9fcyYMXj77bexceNGSgOJPSxPEKKYH/7whxg+fDjuv/9+/PKXv8S6devgOA4OOOAAHHfccZg1axZGjx7t+3FramowZ84cXH/99bjnnnvw6aefolevXrjqqqukncgzXHvttejatSseeOABPPnkk6itrcXIkSMxefJk9O/fX8oYvXr1wltvvYWbb74ZU6ZMQUNDA7p06YITTjihVd9FIQ455BAcd9xxmD9/Pi644AJ06NCh7Pf06NEDixYtwu2334577rkHn3zyCdq3b4++ffvitNNOQ8eOHbP3TaVSGDlyJJ599tkcORg+fDjatWuHnTt3YtSoUcF+eEIiAneEJIQQQogQLE8QQgghRAhKAyGEEEKEoDQQQgghRAhKAyGEEEKEoDQQQgghRAhKAyGEEEKEoDQQQgghRIhIbe40NnWu6SkQQgghsWRGemrZ+zBpICTp5F0lM7FzIISUhdJACLHjpG3DHAghJaE0EJJkvCdqUydtG+ZACBGC0kAI2YcNJ20b5kAIKQilgZCkUuzkrPOkbcMcCCHCUBoIIa2x4aRtwxwIITlQGghJIiInZNUnbRvmQAjxBaWBEFIcG07aNsyBEAKA0kAIKYeKk7bfx6Q4EGIFlAZCkkaQEzBP2oQQUBoIIaLIEoegj0NxIcQ4lAZCkkTYE6/pE7fp8QlJOJQGQog+ZJz0KQ6EGIPSQEhSMF1ekIkNcyAkgVAaCCH+seGkbcMcCEkYlAZCkkBcl01SHAjRCqWBEBIcG07aNsyBkIRAaSAk7iRhO2iKAyFaoDQQQsJjw0nbhjkQEnMoDYTEmaRd5priQIhSKA0kfvDEYQ4bfvc2zIGQmEJpIPGEJw5zeH/3pv4O/PsTogRKA4kXNpywACBVYW7sDCZ/fhtO2pwDIdKhNJB4Y+KgnREGG8TBJDacMG34G9jweyBEEpQGEh9saMTLx9BJy6mshFNh9oRpevwsFAdCpEFpIEQ1Bk9apk/cRsd3PIc3E3+DfFGgOJAYQGkg8aDcAVnXAduCd7VOZWXu5wZO3N4xTYuLVVAcSMShNJDkkMAyRQbTJ27t4zsFDm06/walnmsUBxJhKA0k+vg5CKs8YJc7KWk4aeWnDDlf03TiLjaOaXEBYFzeslAcSEShNJDkwcQh3uMXShm8qP4biD6/KA4kglAaSLQJeuCVfcD2cyJSdNIqlTLk3E/hiVvksU2LCwDj8paF4kAiBqWBJBcmDvEbv1zK4EXF3yDIc4riQCIEpYEkGxkH7KAnH4knLdGUQSV+RcC0uAAwLm9ZKA4kIlAaSHThgTYUNpy0pc7BT8rgRZY4hH0+8vlMIgClgZAwB+uwJxwJJ6wwKYOsk3aYx7FBXpg4ECIGpYFEE9kHV/Y3RHsOQVMGLxQHQspCaSAkg9+DtcyTDMXBijkE/jvESWIJKQGlgUQPlQfUiCUOMhsgg560jZ/sZaQMXpg4EFIUSgMh+URMHGRiWgBMj5/Fz98hrhJLSAEoDSRa6DqIlhtH5cld8LFVLbP0c+JWcZKPpDiohOJALILSQEgxmDjYP77s0kQ+5f4OtkgsIZqgNJDoYOLAWWhMXSf0EuPo2Myp3IlbtViYFpcsTBwIyUJpIMRmmDiU+KLGw1ehv4MtEmsCx7FnLkQrlAYSDVIVek8SXrwHRxMn8bwxdW8ZXejErVMmTItLFiYOrbFpLkQLlAZCREj4wdH0ibvV+KYE0hZxsGUeQOJfG0mD0kCihamTBQCnqo2xsTMniShemCrq47YiVWHHidKUOBT62W34fRAtUBqI/eQfHE2Kg8kTV8L7G7JzMPj3twLvz29T4kASQcJffSSyaD5xOJVV+z42dPJM1VTDaWMw7QCQqqszOj4ApGprjI7vVFUa/zvkoFMcSiUKTBsSAaWB2E2pA2JCEwfTJyynttbc2NXVAMyLA2Do72B7ykJxiD2WPwMJKYOGg6g3Zci5neJgFBPi4FTl9pSY/jtk0ZE2iAoBxSHWUBqIvYgeCBOQOKRqqluPrfmElS8KusUhkzJ4SVTiUO55blN/A8UhtlAaSDxIgDgUHJuJgzZxyE8Zcr6WpMRBFIpDLKE0EDux5OBXrDShk0IpgxcdJ6xScpAkcSiF0r+DHylW8doJKgAUh9hBaSDxIaFpA2D+na5qcShUmtBJqZQh535MHFpDcYgVlAZiH2EOeBLFwW/KoEIcyqUMOeMrOmGJCoHpxMGGtAGIoTjIOOlTHGIDpYHEjwQnDqZRIQ5+UgYV4iCaMuR8j0xxCPN8ZuJAJENpIHYh7d1RuKd2mF4GWeLgJ2XIji35XW4QCWDi0EIsEgfZJ3qKQ+ShNJD4ktDEwYaTlSxxCNrLIEscgqQMOd8f9m8h6zlsW+JAeYgslAZiD0q6vv0/xWWtmAgjDkFShpyxJYiD6cQgLEwc8vD7+lJ9Yqc4RBJKA4k/TBzMjB9SOmSsmAgjDmFThpzHCvK3UPG8tSlxACgOEYTSQJKB4AFYxb4MfsUhbMqQM3ZAcZBWXrAgrWDikIeIOPBkTopAaSB2oGXvfCYORsYP0kwpeV+GyImD6ueqTYkDBSVSUBpIsihxMFa9+6OIOMhMGXLG9iEOSpZNRixxkFmaaPXYticOJk7iFIfIQGlIMrZ0Mdv0rifGmD5ZCW8SpXD3x0gkDjoTMZteezYci0hZKA0keS/WAgdlXdeYKJU2qEoZcsYvIw7Kt4OOQOKgMmXIGcfGxMH0scD0+KQslIakkv/iNPRidaoMHTgT2t8AmD9Zlbz4laZrTFibOJh6XjJxIIJQGsg+TImDqZPo3gO0iStZUhzsTBx0pQw5Y9qYOJiG4mAtlIYkUuoFmbQXqyWJg47SRKvx805Wuk/k+eOZuJKlLYkDAKPPRaDl+WjDpeCzJO1YFBEoDcQY3tKE0SWJBk7Y2bF5gSvTU8iKg4mUIYM1aQPMJG9FoThYB6WBtCZBZYpUu7YtYxsUh4ounYyNnTlZmTx52yQORuew97logvzXHsWBFIPSkDREX4CKX6jFGiATmzi0rzM3tgXvcit6dDc6vlPXDqn9O5gbf+/fwKQ45ENxIIWgNJDixDxxsOEAnarf3/QUkOpUb/SEaRM2/B50Py9Lvd4oDiQfSkOSCPKiU/BCFVlmaSpxSGraAJg7YWbGTdV3NDK+U9cu53Pdv4dCSY8NQpuB4kC8UBpIeWKYOJQ6KOsSh0Ipg25xSHWqz/3c8DttU+KQj+nfA6BHHERfYxQHkoHSkBTCvtAkvVCNbebkEyYO5sayRRx0UK6fhIlDESgOxqA0EHEMvFBVpA2iB2KV4lCul0GHOOSnDDlfS0jikF+ayJmDBWkDoE4cgry2KA6E0pAELHlxBU0ZuKLCDKpPmuUe34bEQeXvwM+qFSYORbDk2JYkKA3EHxHvbwhy8I2jOJRKGXLuF+PEoVTKkDOHGCYOYV9PFIfkQmmIOypeUBEXh0BjSxQHv8ss45g4+HnMOCYOQffGYOJQBIqDNigNJBg+X6SyGiDjIg6+x46hOPgaX7I4iKYMOXOIYeIQFopD8qA0kOBELHGw4WAbZjMnWeIgWppQRdCTb1wSBxk7cIZ5LssWb4pDsqA0xBkdLyCBMVQss0zi5k+mseGdtgxxCJIy5MzBgt8DYIcEWwnFQSmUBhKeCLxIpTaRBRQHGVtGh00bwqYMNpww45I4yMDv81qVbFuVNgCROCZFFUpDXLHkRaNyMyf2N5ghzAlT1sk2qDiETRly5hDgZ1FxcTBbEgeKQzKgNBA5WNzfoGxznIiJg8xeBhveaTNx2IfIc1yHZFMc4g+lIY6YeqHkjatry+goJA4qrmYZtcTB9AlWZsrgRfTnUn0JciYORaA4SIXSQORiWeKg5aI/EUgcVK2YMC0CNqQNgPnfQ4Ziz3fdYk1xiC+Uhrhhw4vDcYxcmMrWxEFFypAzdgQSB5UnVVvEoRSqUwYvTByKYMOxMQZQGogSnAozTy2vOOg+eNqaOOjYl8H0O+1y4qCqNJEzB0vSBiD3uW9UpikOsYPSECcse0HYIA6mUZ0yeLE1cdB1MrUhcSj0s+pMGbwwcSiCZcfJqEFpINJJVZvfAKnC0K6Hpjd/yhcH3bs/mn63XUgcdKQMOXOwNHEwCcUhPlAaiFJMpQ0myYiDzpTBJrwnTRMnUJsSB1MpQwanstLo+F4oDvEgeUf0uGLJC6BQymCsTFFTY2TclrHt7G9IChlx0J0y5MzBksQhZfB3kA/FIfpQGogWdIqD92BtShycA7oD9eZOGk77OqMXpkrt38H4SdN04uC0a4tUl07mxvekDBSHElAcfEFpiAOWPOnL9TIkMXEwKQ57DjB3wgKA3Qd1MTp+c7f9ke7d1egcABgVBy8UhxJYcgyNApQGohXV4lDs3a1OcXAO6J57gwFxaOrXAwDFAUAixaFYLwPFoQQUByEoDVHHkie6nxUTTBz0YkIcvLJgQhyau+2f87lucXAKrFpg4tAaikP0oDSQ2CBSQ0+COGRSBi9MHJKTOIismKA4lIDiUBJKQ5Sx5MkdZF8Gk0sxVYpDq9KERegSh2KCoEsc8lMGLzrEoVDK4IWJQ2soDtGB0kCMIVMcTHfqC6M4bSiUMnhh4hDvxMHvvgwUhxJQHApCaYgqMXlCx6m/QThlMNjfAFAcVFIuZfDCxKE1FAf7oTSQUMjYMjqsOARNGeLW31AuZfCiShxEhUCVOJQqTXixIW0A5IpDmN0fKQ4loDjkQGkgVhD1xCFQL4PhxCHpyBYHPymDFyYOraE42AulIYpY8gS24cJUMohb4iCK7LTBb3ogO20QTRm8xDFxCAPFoQSWHHdNQ2kg1hAkbZDVABlGHEKvmJAgDn5KE17Y3xAPcZB5YSqKQwkoDpSGyGHJk1ZVyhDXpZhliXjiEObkL0McgqQMXsKKQ9DSRD5MHFpDcbALSgOxDlFxULHM0q842LAvQ9CUwQsTh+gmDqouf01xIIWgNESJVIXpGQDQ08tgMnEwRkSXYso64Qd9nLApg5cg4iArZbARikMBXNf0DIySwCMziQqlxEHlZk5RKlPISBm8MHGwI3EQTRtUpQxeKA7EC6UhKmRSBsfsn0z3igmbl2IqK01EKHFQcZL385gyUwYvouKgMmWwpb8BoDiQfVAaoohhcdBNvjjo2jLa9sRBdsrghYmD/YmDjpTBi1XiUGGoVJvw0gRAaYgupsTB1IvVEMXEQUsDpOWJg+oTe7nHV5UyeCklDrp6GZg45OGmARgUh4RDaYgCxRogNYtDqq25hq9M2mDiwlQ2Jg4qUwYvTBzsTBx0pwxerBCHvWgVB6YMACgN0cdE4mDI8G3Zw0H7MsuEbzdNcWjBpsTBGHtTBi9MHPRCabAdkWWWGsShVcpgShzamnvXb0vioCtlyFAobdB9Is8fT0dpIh+vOJhaZmmLONiUNgAaxIEpQxZKAwmOZnFIdWs5cZgUh6avHmxsbJOYLlPYgg2Jg9Ozm+kpADAgDgVSBi9MHPRAabAZP5s5KUwbSvYyJDBx2NPJ0Lus+g747OR+ZsbGPnEwVS7YfVAX7D6oi5GUwcvOIX2Mje1W7e1l6GqHxKXq2lmVOigRB6YMOVAa4kSMV1RkUgYvusVhz4Be+z42JQ4Avuxba2xsGxKHhiPam54Cdhxivs/CFnEANKQOZVIGL0wc1EJpsJWgW0ZLFgfhFRMJTBySyOYj26LhcHO/88+/Yvbv3dh5n7BRHOyF4qAOSkMciVniUChl8KJDHLwpQ/Y2zWnDZ0P2nSBMpg0AjIoDYEfaAOgVh2xpIh9LxEFZ2uAjZfAiRRxYmmgFpSGuSBCHQPsyJMzwk1Sm2Hxk7vNBtzjkpwy6xcGbMlhH3MUhIEwc5ENpsBFZV7OMwXbT5VKGDCrThkIpQ87XNYiDN2XwwsTBfOKgI20omjJ4iaM4BEwZvAQWB6YMBYn+WYUoIdTujwnsb4h74pCfMnjRIQ6lehl0iEO5lMGK/gYgnuIgASYO8qA02IaslCFDhPsbRFMGL3EUh2IpgxcmDvFNHIRSBi8Uh4L4EgemDEWhNCSBCItDEGSKQ7nSRKv7xzBxKJUyeFElDqIrJlSJg59eBiYOuYQSBwmliXyYOISH0mATslMGLz7EQeqFqQK+SIOkDF7ikjiIpAxemDjEK3HwnTJ4iYM4KKCsODBlKAmlIUkwcfCF35Qh53tjkjiIpgxe4iQOQVdMMHHIxbc4KEgZvDBxCA6lwRZUpgxeyoiDsstf+3iRhk0ZvER58ye/KYMX04mDLIJu5hSHxCFUyuAlquKgGIpDMCgNSSQGSzFVEyZlyD6GwbTBNKbTBiAe4iCNKImD4pTBSytxYGmiLDx7kCzKUoYMES1ThCHKZYogpQkvYcVBxpbRYcRB1mZOFIdcmDhEG0qDDegqTXixtL9BZmnCix9xkJEy5DxeAHEIU5rwYrpMwcShBb/iIK00kY/t4qAxZfDiVFQwZRCE0pBkPOKgPGXwwsRBK0HEIWzK4CWIOMi+MJVfcVCxZTQTh1xsSxyMvHmLIJQG05h+olqUOKhKGbyUEwfZKUPOYwuKg6yUwQsTh+igLGWwkBxxMJQytAy9N2UwfTyOAJQGglRdnZmBE5g4mERUHGSmDF5ExUHV5a9F0waVF6Zi2tAaJg7RgtJgEpuenCmziYOOlMGLKXEolzaoSBm8JD1xsL2/QWvKQHEA4EkZiBCUBrIPU+JgiHxxUFmayBnH8FLMUuKgKmXwUkocVKUMOeOXEAddl79m4pCL29QEp9aivUVsekNnGck6S9iEJU/KVDuNDZBFcA7Uc7K2iULioDpl8GI6cTCNjYmDsV4GS8QBgHZxKJkyWHKMtg1KA8nFUNrgVrcxMm4mbdCVMnixLXHQkTJkKJQ26EgZcuaQJw66UgYv1iQOXeqNDe02NeV8zsTBbigNJrDkiVg0ZUioOJggIw46UwZbMN3fANiTOFixYsKgOORDcbAXSgMpjCZxcHp1z/nchDh8eWQPNHYyIyxA9PZwkAnFoYVtAw2esL2bGmkWh/yUwYtqcfDVAElxyEJpSChCvQwJSxxMisPWPuZeil/2rdVamsin4fAa7aWJfD4eu5/R8QHD4uCFiUNhKA4AKA36idoTT6E45KcMJvjyyB45n5sQh3UntpywTYrDZ19tNjY2AGweZnZ8ANg0pMrIuM3VTvZj7eJQbOtkDeJQKmXwokIcAi+zjNrxWwGUhgTie8WEgcTBVNoAmE0cTLDhpJad+EyJQ2ZcU+Kwo/u+k7YpcfDCxKE1TBzsgdKgk4Q/2byIpAyqxSE/ZfCiSxwyKUMGk2kDwMQB0CsO3pTBS9zFQTRl8CJLHKRs5pTgYzmlIWEE3pchYf0NJtEpDpmUwYtOcSg0lk5x8KYMXhKROIhe1ZGJQ2ESKg6UBl3E4QmWMHFQnTbkpwxemDjEP3EoljJ4iXviEIQw4iB9y+g4HNd9Qmkg/pAgDkEaIGWLQ6nShJc4r6golDJ4US0O5R5ftTgUSxlsQ4k4iKYMXiSJQ5DSRD5MHMxBadCBJU8qaVtGM3EITamUwQsTB7Pjq0obRFIGL0wcWuNXHJRemMqSY7wOKA0kGAHFwcZlliLELXEolzKoxrSM+EkZbOhvACSKQ5CUwUsIcZCRMlhLQsSB0pAQlFyYiksxAyGaMnjhHg5mkSkOflMGL9YkDpZgVZkCSIQ4UBpUk4AnkSiyUoYw4hAkZfASt8RBFIqDPYlDKMKmDBkCpA2qUgYRcVBamsgn5sd8SkMCUHr564T1NwDRF4egpQlZ4hD0cWSJQ5gGSBvEwZq0IcL9DcqJsThQGlQS4ydODgLioKKXwa84hE0ZvAQVhyClCZtg4hBOHMKUJrwEEgdZKYMXQXHQ0ctQTBy0pgxeYnr8pzTEHKUpQ85AyUscTBEmbZDRABlGHGRIRxhxkLXMkomDByYOhXHTgBONZb1+oDSoIqaWWRKLxUFmypDBb9ogM2XgUszoJQ6yUgYvURAH3SsmvOJgLGXwEjNxoDQQuRQQBx3LLOO0h4MofsVB9jJLv+IgWzT8ioOKzZwikzioKE3kw8RhH27eay1G4kBpUIElKYO20kSrgZP1tBIRB1W9DEwcopE4qEgZvEQhcdCNU1lpegq5xEQcknV0J9rRuZlTsbRBRWkiH9sTB5WbOYmIg0q5EBEH1VtGW5046EgZvHjEwfRmTqk2Bv4u+SmDlxiIA6UhphhLGbITsLe/QRW2i4NKmDgUR3XK4IWJQ2uMiEMpIi4OlAbZWFKasAHnAPXv8AvhFQcdKYOXQuKga5llMXFIypbRxcRB14WpbEgbgDxx0J0yeHDq2hkbG3v25HyqTRxKpQxeIiwOlIYYYjxl8GLoxcHEQT+m0wYbyBcHnSmDF1sSh9T+HUxPIQsTBzlQGmTClCGL07Ob5xMzL47PvtbLyLheor6Zk19Mi4MNZQpbEoddffY3NrazszH7sXZxyEsZvFAcwkNpiBlWpQxeDL04mmrMPMVtSRtMlCa84mBCIrzioKs0kY814nBQR9NTAJCQxEG0NJFPxMSB0iALpgxZclKGnC/oe3F8PqRL9mNT4vDhhBoj4wJsjLQhcdh4jDlxqNy57wSmWxy8KYMXLeJQImXwwsQhOJSGGGFtypBgdnUz14TY5ZLVxsYGgFH/35tGx8c5DWbHB7B5sB17BTBxaI1UcQiaMniJiDhQGmTAlEGchJQp1o7aN55JcTji0DXGxgaAs49ZZHT8tsdvNjKu63GFJIlDsZTBizJxEEwZvDBx8A+lgUilaGki505qXxje0oQXU2UKQL84HHzah9mPkygOnQ/4PPuxKXHwolMcvKWJfJg4tCa0OMhIGbxYLg6UhpgQudJEQhIHL0lKHEwnDPnoFAe3iB/YkjjYglRxCJAyeGHiIA6lISwsTWQRShlyvkH+C6NYyuBFtTh4SxMm8KYMXkwmDjolwpsyeElC4lAqZcigKm0QKU3kE/nEQXbK4MVScaA0xIDIpQxeEpY4mEwbAD3iUEwQbEgfVItDsZTBiw2Jgy1lCiAG4qASC8WB0hAGpgxZfKcMChBJGbyoEAeRlEGlOBRLGbzEOXEoljJ4iWviIJIyeJEpDkFSBi+hxCFkaSIfikNpKA0RJ9IpQwbLXhQ6MJ04qEJECuKaOIikDLbBxKEwQuKgsjRhMZSGoDBlyCIlZQgpDn5Thgwy0wa/vQxJaozMJ67i4AeZaYPflMFLWHEImzKEQnLK4MWqxMGiN1aUhggTi5TBS8L6GwC54iBSmvAiWxz8ioBscRApTdiGDf0NgD2Jg01pA1BCHEykDJaIA6UhCEwZ1BHghRE0ZfASVhzCrJhg4mAOWWlDmNJEWHEIkzJ4CSIOKlIGYXFQmDKQ4lAagmBBLcuWlEFJAyQTB9/4TRm8yBCHMCd/GeIQJmVoe/zmWJUqwmBT4mBL6pBqU5WbOFhw/DcJpcEvmRNawp84yhEUBxkpg5cg4iBrXwYmDtFEVgNkEHGQlTJ4ERUHHb0MRcXBQMpgvMfBdc2OvxdKAwmM8mWWTByECJMyeAkqDrJO+EEfR1Yvg+m0AWDiUAhbEgcAcFJ29BWYhNIQBkNpgy2liTgjKg6md3+UTdITB7/ioGKZZRTEIa4rJkRwKgz0tFmSMgCUBn8UeudrQhzS5ksj2jZzKpE2yC5N5BOFXSNlpQxeoiQOKlZMRCVxUFGayMeWxMGGtMFtajI9BSugNMhAozikamtaPrBAHLRhcKlRFMRBBaLioCoZiELioHozpygkDjpJ7d/BeMqQQWvaYFHKAFAaxCl34kpQ4mBky+i837/qlEEE1aWJcuKgImXwYjpxKIfqfRlsThx0pAxevOJgsjTh1NYaGbdQymCkTGEBlAaZKBaHbMpgGlMpBxsjtVNKHFSnAabTBsBucdCNLYmDKXEoRBLFgdIggiU7cRVE8wnc6d7FyLj7JuAYSRnyxUFnA2QhcVCdMthCMXHQuftjIXHQfZ0JrzjoThm8pNube+Pibtue/Tgx4mBZaQKgNMhHUdpQMmUwdQI3VR5Jm3khJTVxKJQ26EwBmDi0YEvi0NzVfFMioE8cRBogk5Q4UBrKESRliGl/QzZl0Dyuly8Gd9I6Xj5NNSljyyxtEwed2CQOJq9mua2nHScn3eLgTRm8xDpxsDBlACgN6oipONiCqbQBANJtzP2ed3VLGytNZMTB1Ak8M67JC1PZkDh8frCZnQnbNOw0Mm45VIqD32WWSUgcKA2lCNvLIEkcfDVAmhCHBJUp1o9o+VlNioNJkp447GisgnPkF0bGrtm87/luShy86EobiqUMXmKXOFiaMgCUBvXEJHEoWJpQPGY+hUoTSUscvjt6Lo7tuFr7uBke7vs0ruoyz9j4I/ZbgRsHvGRs/AymxMGLTnEoljLY0t8AxFAcLIXSYDnWLLMUIeaJQyZlsAGT4gDAqDgAMCIOOxrNv7vPJ+6Jg0jK4EWmOITdATKwOFicMgCUhuLIXGYZ8bShbMqgaFwv5RogTSUOOtOG746em/N50sRhxH4rcj43nTjoTBu8pQndiPQyMHEoTBwTB0qDLgKIQ+iUgUsxtWCyv0GnODzc9+lWtyUxcfCStDJFKeIkDjKvM+FLHCxPGQBKQ2FUbeYUwcTBV8ogcVwvppdZAuVLE6rFIT9l8JKExCE/ZdBNqdKEanEQSRlUiYPfFRMyxcFvaSIfJg5qoDToRlAcpPYycCmmFuKcOBRKGWzBdNoAMHHwEvXEQdXVLOMiDpSGfHRsGR2RxCFwyhBizHyCpAyyxcFPA2ScxaEUKtMGkZRBpTiINkCqEAe/vQxxEYewKYOXyCQOEShNAJQGopqE9TcA8sWhVGkiHxXiIJoyJL2/AYhX4hB2M6coJg6qUgYvUU8cKA2mKJE2KFtmGbHGyLC9DDLEIegyy6Ru/hQ3cQiyzFKWOIRZMRHlxEFmyuDF6sQhIikDQGnIRffVLC0uU4QuTQQcVzYmEwcZ+EkZMshMG4L0MsgUhyANkEwcWggjDjK3jI5i4qCDqCYOlAbT5ImDls2cIpY4hCWoOITdzCmp/Q1APBIHGzdz8kvUEgdVKYOXYuKgozSRj1NREamUAaA07EN3yuDFssRBesogOK4X2csso7iHQ5CUwUvUxSHsMkvTiUOYtEHmZk5+xUHVhamYOBTB5LknAJQGW7BMHIg8opo42LzMUpSg4iArZbChTAFEL3HQgVccTKQMAOA2N++dTHTEgdIAWPMHM3KdiTxxUJoyFBkzH1WbOflJG2RfZ8KvOIRNGbxEcSmmzM2copY4qNoyWkQcdFz+upg46ChNWI0l56FyUBpsImOdSSBi/Q0yiFLiIDNlMN3fAERPHFTBxCEXp7bWWMpQkAiIA6XBkj9Sqrq65QMT4rD3BK4lZSgwrhcdW0aXEwfTV7OUmTJ4iUriYPOW0WEQEQcdF6YqJg46UgYvXnEwmTI4lZVGxnWLHestOScVg9JgIwbFISnjRrExUgYi4qCql8F04mA6bQCYONiGu7NFlEyJQ1EsFgdKgwVkUwbDOLt2Gx1f94WpComDjpTBtDiYxHZx0LHMspg46L78tQ3i0Ny1gzW9DDrFoWjK4MVScUi2NFj6RwGgPW1I1XcEYEgcDK7isC1xUFWa8FIqbdCxYqKUOOgoTTBxaI3u0oQXp66dsbHzYeJQnmRLgwWUTBkMNUYmTRwy6O5liFJjpGxsTBx0b+bkFQfdKUMG02lDxaaW34EJcciUJvJRLQ5CKYMXy8QhudJg2R+iKBrEIZMymGTHVzqhaluy+huAXHHQkTJ4yRcH3fsy5IuD7gZIJg4tpCvsOBYycSiBReer5EpDlDCQOJjqbzAlDp8OMfeiZOJgjow4mNwy2gZx2NW1rfYxMymDF13iUCxl8KJCHHynDF4sEYdkSoMtv3w/DZCKxKFUyqBLHHZ8JbcB0pQ4VG41dwGZM09cZGxs0+IAmF1maTpx2P3+fviyj5mx91u9b48CE+JQCCYOJbDg3JVMaSDCJC1xMCkOm/fUGRv7hf5/MTa26bQBAMYe9J7pKRgTBy9JEAeRlMGLdeJgmORJgwWmBgRcZik5bRDtZTC9FFM1m47OfRnoFoezhv1T63jF6FFpRlp6VNbhG+12GBkbAGZvHQTAjDjsfn8/7WNm8KYMXnSIQ6HSRD5xSxxClSa8GD6HJU8aok7MVlTklya8mEobAHOJg4m04YbO+0oDpsQBgFFxyGA6cbAhbQDimzj4TRm8WJU4GBQHSoMBQm/mlKBrVKgWh/yUwUuSxMEU+ZKiWxwyKYMXXeJQLGWIuziIpAxe4pA4SEsZvBgSh2RJgyWlCSmEfBIGWWYpO20olTJ4MZk4qKZYaUKXOHhThgwm0wbAjsTBNKrFoVhpIp+4Jg5hSHrikCxpsACpW0ZzKWYoSqUMGZLaGJkEcSiUMmRQnTaI9DLEMXHwmzJ4CSsOYUoT+fgRByUpgxfN4pAcaYhTyhCSsJs5yRAH0ZTBS9xWVIg0QKoUh0IpgxfV4lDu8U0nDqb7GwA14iCaMnhh4tCapCYOyZGGuBKzxshyyBIHkZTBCxMHM6gSh1Ipg5e4ikMQwopDmJTBSxBxkJkyeCknDspTBi+axCEZ0mBJyqDsapY+npgyt4wOKg5BUgYbkCkOfpdZxk0c/Dxm3BKHIMssZYlDkJTBCxMHy9FwrkuGNCSBBCUOSVyKKZtypYl84pQ4iKYMXpg47CNq4qAqZcjOw6YyBaBcHJIhDa65CxJlUJYyeCkjDjZcmEoWYcTBb2kin7DiEHQzJ9NLMWWJg2kBCYoMcQi7mVMYcQibMoRBVmkiH1sSh0LioLU0kTOw2vNdMqQBsEIctGD5igqZpYkkJg4yxMFvyuCFmz8xcQDsSRuA0uKgOmXImYdtiYMikiMNgDFx0JIyGCYqjZFhU4awyNgyOsqJQ1jpCCsOQUoT+QQVB5lbRkdJHFSlDF5sSxyMpQwaSJY0JIUCT1gdpYly4qCqATJuSzFFCCoOYVIGW2Di0IIfcVBVmrApcbAFo4mDhjfGyZMGzWmDU1EBt8lALTFBjZGAmDioSBmiKA4yMN2XEEQcZKQMXvyIg8kLU6mmlDjoSBky5KcNOksTXtK79wBOfE+t8f3JSmGgTGFSHHQ3QBYSh6gusxRBVBxUXM3SjzjIThn8ioNs0WDiIJY26GiAtCVxsKVMAUC/OGg6ryVTGoDEiYNubFuKqbqXgYmDGUTFQXbKYBM29DcArcVBZ8rgxalrZzZlyJlM/E6x8fuJ/KBYHJyK1icS3eJg8gWUQWfKkMQVFSKo7GUQEQeVcmE6cSiXNugoTRQTB93LLG1JHFBVZXoG+9AhDhrfBCdbGoDELMU0IQ429zeoopg4qChNeInyigoZlBIHHSmD6TIFYG/ioJv053sTjqSJgybi85NYRqGUIYOutMGG+p4N4mB6maUubBUHXUJhY+KguwHSKw4mN3NKbWs0NnYOGsWhVWkin5iIQzx+irAkpL/BRNrQ3L0jajbv0j4uYM9STNUpg5dC4qBzmaXpxME0TBz24bat0T5mNmXwEvfEQfP5i9KQQeIvvlTKkDNkQsQBgDFx2N7dTJ9BUhsjgVxx0C0R+WmDiQbIjDiYXGZZ/bm5smvbDz/LfmxCHAqiWBzKpgxeIp44RHv2solR4mDD1qrN3XOXepoSh3ZrzVzltHJrhdaUwUtGHOKwmZNfTJcpADsSh8b97W3MNUIcEwcD5yxKQz4h/wiiKUPOkAlKHHTy+VfMHyRmf9zf2NhJX4ppcpnl8wuPQtP+5rcS1i0O3pQhg660oWBpIh8F4uArZfAS0cQhmrNWTcQTBxsuGZufMmSo2bzLSOJgIm2oHfg5AHPiMLbDMry8w5J42ABdq740PQUj4tD57dzjlw2Jg9u2JjGlCl+EEQdT11IyMmpMCZIymCauiUOhlMFUmQIwmziYFIf1TduMjDt58wAAZsTh+YVH5XyelMShUMqQjypxEEoZvEgSh8Apg5eIJQ7Rmq1OIpo22LDMsljK4CXu/Q2ZlMGLTnEY22FZzue6xcErC6bEIUOSEof8lMGLDYmDVUQ5cTC4vxCloRQRFQffY8ZsRUW5XgaTiYNJkpI4ZFIGLzaIgw2oEgeRlCGD7LTBd8rgJYQ4SEkZvEQkcYjGLE0iKA4ySxNBxSFMyiBLHERSBi9xTBwKpQwZdKQN+SmDbooJQhISh/zShBfVaUOplMGLDYmDNf0NQPQSB8O7GFMaRGDioJQ4ikMpktrfAKgXh0IpgxfTiYMN/Q1AfMQhVMrgJWriYBC7Z2cTJcTBhgZIWb0MUW+M9LvMkuKgn7gmDqVSBi8qxEE0ZfAiSxz8lCbyiWriIL00kY/F4mDvzIiZS2kjuDj4LU14MZU2yKZUaSIfFeIgWppQJQ6mhUAUJg72EFQcpKUMXmxPHCy4wCKlwQ8F/mCqUwZT4mCCsOIQdDMnLsXUjwq5KFeayCcu4hAkZcgQNm0IkzJ4iVLioDxl8GJh4mDfjGzHwv4GFcss/aYNYVIGL1Hub/CTMniRJQ5BGiBlioNfEYhKKiGCaGkiHxsSBxv6G4BoiYNWMuJgQcoAUBqCYaE4KBkzQo2RMraMZuKgH1ni4DdlyGA6bQCiKw6yUgYvouKgpDSRTwFx0JoyeLEocbBnJlHDdbU3QBYSB9WbOYmIg6yUwUvUEoegKYMswi6zDCsOYU7+phOHsOIQNGXwElQcwpQm8mHiUABbEoe0ebHMQGkgZbF9RYXsC1MlcUWFacKIQ9CUwQsThxZExUFFyuCllDhoSRlIUSgNQUlVwE2bLVPo3DK6mDioSBkymFxR4UccZKYMQcRB1mZOQdMGWUlBFBMHGSmDFz/iIDNl8MLEIY+9aYOx0oRFKQOgQBqaEtTtD8C4OGgd10DiEBVxkElS+xsA/+IgI2XwwsShBVvFwVjKYEuZwgJ8S8OkSZOws8jJY/Xq1TjhhBNCT8p6UrkvKBPiAAs2lNJFKXGQXZqwBYqDOUTFQXbK4AdVKYOXYuKgujSRjw2Jg7tzJ5xUMq9Zk49vaXjmmWdwzDHH4J133sm5/bnnnsNRRx2FjRs3SptclDCSOOzU/y7cmzaoLE3kY2tjpMoGSBFxUHWdCVFxUHWCj4o4qMKGtAGwK3GwoZdBuzhYVpoAAkjDW2+9Bdd1ceyxx+LRRx9FU1MTLr/8cpx99tk48cQT8c9//lPFPO0hVfxFlERx0Em+OOhKGbgU005klybysVUcdKQMXrzioDtl8OJUVhoZN/94l/TEwbc0DBo0CIsWLcKECRPwve99D3369MHDDz+Me++9F88//zw6dtT37jOpONXVpqcAVJp5B2JT4qBrmWUxcdBxNctS4qA6DTCdNgDFxUFXaYKJg51oEQcLUwYgYCNkbW0txo8fj+rqaqxfvx6HHXYYvv3tb8ueWyRJStoAABWbtxoZFzDTy2AycTCJjZs/qU4ZvNiaOOjGqDhsagBgLm0oRFITB9/S0NzcjKuuugpnn302xowZg6lTp2Lt2rUYMmQI5s6dq2CKFlGiNOFFpTgUSxl0ioPTvi77sQlxsGFFhe7NnPLTBh0pg5d8cdCZAtiQOHgx0QCZEQfdpQkv+y//HG61+ZO2TnEoV4pVJg6WpgxAAGk46aST8MADD2DKlCl44YUXcM4552Dx4sXo378/xo4di9tvv13FPCNHkhIH3ezqWoP9PzS0ZhrJXIppGq846EwZMphOGwB7Egft7E0ZvDBxMIdvafj4448xd+5cXH311dnbevTogdmzZ+Oaa67BbbfdJnWC1iCYMniRLQ4ivQyqxcGbMmQwVaYwJQ67OhsZFkCLOOhOGTJk0gZT7/xNJw5dq740usyy698rkK40c4Laf/nn2Y9tSBsA9eLgp+FbqjhYnDIAAaRh8eLFOO6441o/UCqFO+64A6+++qqUiZHgmEgckiYOu9821/B763tnGhvb9IqKb7x/mtHxnSbz7ypNiYMXbeJQIGXwwsRBP76loVOnTiW/PmbMmMCTsZYAKUMGIxs/GUKHOOzq2vqkRXHQy12bRhkbGwDe+KyvkXF/O28kADPi0PXvuccgneLgTRm8xDlxCLqsPAniIPTbnjdvHo466ijU1dVh3rx5Ze9/0kknhZ5YnHDTbugnk99llu7OXXBq5b4rLFSayKdi81Y0d95P6rgi7P/hHnzeT/2KivzSxO63O6LNV/WtXa+p2reF+K3vnYmbB7yobWwAeO3zIwC0iMO1XedoHdvLG5/1xdc6rjI2vtPkwK00+4YgXekg1WR2Dm51JZxGRdval0kZvDiVlca218/HSTnB3yxaXpoABKVh5MiReOONN3Dsscdi5MiRcJzCJ0DXdeE4Dpqb7f/BhQmRMniRIQ6+x1QgDiKoEodCKYMN6BYHLybEIYNucbhk1QRtY4mgSxzyUwYvqsWhWMrgRak4+MAmcYgzQtIwZ84cDBo0KPsxCUZQcQizmZMscRBJGUyjOm0o1QCpQxy8KYMJMimDF5OJg860IVOayMeGxMEG4iQOMna8DZQ2RCBlAADHdd3IPOPHps7VP6ikpMGLX3GQsQNkWHEIIg0y0wY/KYMqcRBZNaFSHEpJg460oZA0ZFAtDqVSBh3iUEwaMqgUh1JJQwZVaYNI0uBFmjj4KE0UIow4yNwm35c4WCANM9JTy94n1KWx161bh6VLl2LdunVhHsZeFAiDX2RtGR1mRUXQlCFOKypMLrMEyqcMqhsjSwkDYLY5UnVjZDlhANQ1R4oIA6CmMdKvMADRb46UfV0d4TeIFgiDKIGkYdq0aRgwYAB69+6NI488Er1798YhhxyCp59+Wvb8YglXVPgjSC8DV1TEB5FeBlMrKmzChqWYgARxCJkyZLBlOWbcVlT4loYnn3wSEyZMQEVFBW666SY8+OCDuPHGG1FRUYFvfvObePLJJ1XMUz+KUwYRcZB9YaogaYOMXoaoJw5+UwbZ4uCnl0GFOJRLGTLEcSmmSMqQQXbaIJoyeJElDkFSBhvxIw4qr94bJ3Hw3dNw2GGH4aCDDsKLL76IVGqfc6TTaXz961/HmjVrsHz5cukTBTT3NGgqTZR6Mqm6mqWf/gaZDZBBexzCrpoI2+MQtDQhq78hSAOkzB4HUWnIILO/IciKCZk9Dn6kIYOs/oYg0pAhbI+DDGkI1N8gKWXIR6THQaU0ZMco9GbRotKEkp6GDz74AJdeemmOMAAtO0Jeeuml+OCDD/w+pH1o7GVI0jUqgiQOti6zFCGppYq4JA5BhAGQkziEEQYgXOIgK2Wwpb8BYKlCJr6loU+fPtixY0fBr+3YsQO9e/cOPamkUUgcVKUM2TEFxCEKyyxFCFOmCNsAGVYcwiyzlCEOflOGDHERh6Bwu+kWfImDopRBBB0pQ4YccbAoZRDFtzRceeWVuO2227B58+ac2zdt2oQ77rgDV111lbTJEbXYfo0KmSmDyatimkwcTBJWHMJu5hRVcQibMnjxKw4qehlsSRxsSRuAaCcOvn+Ly5Ytw9atW3HQQQfh5JNPRvfu3bFhwwbMmjULnTt3xvLly/HjH/8YAOA4Du6//37pk1ZOuln7ckvvxk+qUwYRVKYMUdlqWuYyyyCbP8nYzCnMjpFBU4Y4ELQ0kY8Nmz/FfrtpHxTa/ElnypBDBFMGIEAjZH4vQ8kHl7yltPbNnQzs0+CkHO3SUKgxUkdpopQ4qOxlEBUH2XszmJCGDEHEQZY0BGmMlLlldJDGSFnSkMGPOMhMGryIiIPqVRNFxUFzacIrDqakwcYtr5U0QqbTaeF/kb8GhQETtKExUlcvg81LMVVs5uSnTCF7y2i//Q0yU4ao9TfIFgY/qBIGoHypQscyS9tKFRQG/4TaETIRaBYHJ+XAbWzUOiZg14oKHSsmuPmTXvyIg4oLU0W1vyGOtBIHQw2QNvU4RAlKA8liShySRjlxUHlhKhFxUNXLEIXEQWXKUE4cVKYMGYqlDbo3c7IhcUh/+aXpKUQSSoMImtIGb0etibQBANCc1j6kbWUKHdeZYOJQGNWXv2biYMdSTGCvOBhcZmmKKJcmAEqDOCb6G0yJgwEy4qB7MyeTSzFNYrM4qKaYOOjqZSgkDjpSBi9ecYjLltFED5QGPygUh2LrdnWKg9OmTcuYn5t5529D4qDzapaF0gaVpQkRdC2ztFUcdMHEwSymShNRTxkASoN/EpI4mBCHdLta1H5k5sXMxkj95IuD6tJEPhQHs+LgfLzR2NgkOIGkYffu3di4cSM2bdqE3bt3y54TMUAmZfBiKnEwRd16M0uEM+JgImXwikMSN3PKiIOpZZZOk6O9NOGlfvFngP42JqMwZQiHsDQ0NDTg2muvxcCBA9G2bVv07NkTPXr0QNu2bTFw4EBcf/31aGhISFOL5LRBZEtRU/0NusQh3a42+7GptAEAOi81Kw4msKG/QXfKYBNVOy04a2ueAlOG6CK0I+SqVatw4okn4tNPP8WoUaMwePBg1NfXAwC2bNmCpUuXYs6cOejatStef/119O2rJvbTviNkOSTtGOlnH3IVu0UWShla3Wd/tds+e6Uhw84+7ZWO6aWpdp8/bz5C/zu/1N7Aru64T7WPDQCfrumIM49ZbGRsAFi9vZOxsVdO74ddPczIYs+5+z7eU6u3Wly/uMAOpZqmYEoaTC6zjELSILIjpNBi2auuugodO3bE/PnzceCBBxa8z5o1a3DGGWfg6quvxtNPP+1vplFFwjUq/F64xG1stOLaFDqo/ehLreKQofPSZiPiAADb5ncxJg4vLhxiRBxeXDgEAHDEoDXax85Qs77CmDhkqNqZ1i4OrUhDuTgkMWWIgjCIIvT0mD17Nm6//faiwgAABx54IG699VbMmjVL2uQiARsjQ1MoZchgqlRhqkwBtIiDTj5ds680kjmBm2DpO8WPL6pYOb1f9uOa9XpF0Zsy6KZgypDBgmqJCriZkxyEpKGpqQm1tcUP7Blqa2vRFCOjEiagOIS5PKoscRApTWTHNNQYqVocmoq8u9MlDinLeol1ikP+WCbEwYtuccjHiv4GQJk4MGWIPkLSMGzYMEyZMgXbt28vep/t27djypQpGD58uLTJEfuQLQ6lUgYvSUscdKcNNqFLHLwpgxcd4lAqZVAtDiVTBkLKINTTcPfdd2PUqFHo168fJkyYgCOOOAL19fVwHAcNDQ1YunQppk2bhh07dmDu3LmKp2wpPvsbwqQMGcL2N/hJGXLG/Xyr8sZIXRRLGbyo7HEolTLo6G/wlia8mOpvyLD0nQMT3eOQlP4GXXCZpTyEVk8AwHvvvYcbb7wRL7/8MnbmXU60trYWZ555Jm699VYMGDBAyUQBC1dPFEJQHGRIQ/axAopDUGnIfr8EcRBNGjLIbowUkYYMKsRBpDShUhyKSUMGleIgUgZRKQ7FkgYvqsRBtJ9BhTj4ThokTcFkaYLSIIbI6glhacjQ3NyMDz74ILsnQ6dOnXDwwQejokJ9pBcJaQCExEGmNAD+xSGsMADhpcGvMGSQJQ5+hAGQLw1+ehlUiEM5YcigShxMSoOIMGSQLQ5+GyBlikPg0oSEKSRtmWXUhAEQkwbfT4WKigoccsghGD58OIYPH45DDjlEizBEijKNkbKFAYjfiopSJK2/ATDb46CiMVL0MU03RtqAFc2RIaeQxAbIuOLrouau6+LNN9/EsmXL0NDQAMdxUF9fj8MPPxzDhg2D45jfS53oJWh/Q9CUIUPYPRz8pgwZZPU3BFkxIbPHQTRlyGCyx0F2f4OflAGQ298QdJmljB6H0A2QEexx4DJL+QhLw5///GdcffXVWLduHfIrGo7joGfPnrjnnntw/vnnS59kJCnSGKkiZcgg2hgpozSRM66hxsgkbv5kElniECS5SHpjJBDd5sgkpgxRLE2IIvTnf/LJJzFx4kQMGjQIjz/+OJYtW4Z169Zh3bp1WLZsGR5//HEcfvjhuOCCCzB1avmaSGJIyMZPgL9SRdiUwQaivBTTb8rgJeqbP/lNGbyEXYopYzOnoKUKqcssLaiWEHMINUIOGTIExx57LB5++OGS9/vBD36AhQsXYvFiNTFmZBoh89mbOKhMGfIpljjIThlaPb5A4qBCGvwkDkFLE4UIkjjI2MwpTJkijDRkCJo4yJCOMIlDGGnIEDRxkLkDpN/EQcneDAJT4IqJaCGtEXLFihWYOHFi2ftNnDgRK1asEHnIZMHEIUscUgYvUUscZAgDEM3EQYYwAMESB9lbRvtJHJRt5sTEIZEISUN9fT3ef//9svdbuXJl9uqXJA9X/yssXxxUpwwmEV1RITNlyOBHHGRuGR21XSNliobpVRWmt5sG4rGqQhXGUoZms30vOhA6gp577rm45ppr8NRTTyGdbv0sSafTmDp1Kn72s5/hvPPOkz7JuGDiCcWlmPEmbksx/RAVcbD2wlSyKCIOSWyABADEfBWhUE/D9u3bMX78eMycORPt27fHoYcemrON9Lvvvott27ZhzJgxePbZZ9G2bVslk41sTwOQfSI5Bva0cKqrjaQM+f0NukoTxfobVKQMXsr1N6i8MJVIj4Os0kQ+ov0NqiRDpMdBVmmiEOV6HFRLQ6n+Bm3XmSgwhcRt5pT/ptDfvolWIK2noV27dpg+fTpeeuklnHfeeUilUvjggw+wcuVKpFIpnH/++Xj55Zfx2muvKROGSOMxz6SkDUDyEgdu/hTuPkGJSuKgChvLFIlNGRKA722kTRLZpKFAXKU1cdg7VsqQ0Dn772ekAdKbOKhOGbwUSxx0XAK7WOKgKmXwUipxUF3KKJU2qEwZvBRKHHSWJvITByNXs9w7haSlDECRN4TROb0CULSNNPFJkfqWicQhvWOH9jEBc4mDKQolDjqEwTTFxEBH74PptMEGbEkckpgyFD2ex7C/wdc20u+//z6effZZ/Otf/0JDQwNSqRR69eqFESNGYMKECagOcZnmJOI2N2vvcUjv2KE/caj09TSThqkdIwFzu0YW2mpaR8qQwbbtpnWlDEDrXSNNNEBmdo00kjJkcNOAw/ejWRwncolDKYTKE+l0Gj/5yU/w3//932j2GJXjONktpXv16oXHH38cJ510krLJRq48IWiZSsWhyGNrFYe90uB266RvTA9fDuxgZFygpVRhImXwioNOaciQEQcTKyy84qBTGjJkxMHkqon2K74wMq7zyXrPJ3rFwZoGyKJ3tF8cpJUn7rvvPvzqV7/C9773Pbzwwgt47bXXcMMNN2C//fbD73//e8yfPx/Dhg3DuHHj8M4774SeOIkRnpTB2dhgZAqmDqAmidoeDjLJlCpMCAPQkjiYFIb9lmyCs8tM83MOBvamsZqYlCqEkoYBAwbgG9/4Bu65556c25944glcfvnlWLduHSorKzF69Gh07twZTz31lJLJxjVpABSlDWUeU0vaUKA0kaTEoXJ7GpuOqtI+boadB5jb0vbMYxYb3cuh9hMzZTEA6Pn3XdjV2czffb8lm7IfuzV6S8Y5SUP2RvWJg/UpQ8432Zs4SEsaPvroI5x22mmtbh83bhw2b96c3S3yu9/9Ll5//XWf04wpPq0ySY2RSUscuv5zj5FxAaDLG+aWA/79f4YaGxsAOqwy+063ZrO5v3sGnYlDQWEAmDjkE/HEQUgaOnfuXHAb6ffeew+O46Bdu3YAgD59+mDr1mR1ystEqjgIJhdKxaFEAyTFQT11n7S8ozEpDvWLzYzd/a8thzYT4tDz77uyH+sWB2/KkCHupQqTyywDE2FxEJKGs88+GzfccAOef/55NDY2Ip1OY/78+bj44otx6KGH4sADW2qI69evR/fu3ZVOOBKEeEIwcYgPldtzD5RJShy8smBKHDIwcVBP0ZTBS8wSh9DH6oiKg1DRb/LkyXjrrbcwfvx4pFIppFIpNDc3o2PHjnjttdey91u6dCnGjRunbLJEHdKXYhpaZilC+xVfGF1RoYNMymAL9YsrsGWIHiHOpAwm8KYMuimUMmRwdjVq728oiOTlmJFMGbxEcDmm8I6Qzc3NmDp1KubPn4/du3dj0KBBuOCCC9Cpk76mtkg0Qkqyx1CNkQG/15Q0xLExMj9l8KKjMbKYNHz6NT0n7kLpgmlp+KKvepkoJQ2qGyNLSUMGFeIglDLkfEP0pUF6ImyJOIg0QnIbadlIjJwCi0MI4ZAiDgFShriJQylpANSKQ7mUQbU4lCpHqBaHcimDSnEQSRlUiYOIMGSQLQ6+pQGQJg6xkQbACnHgNtK6kVyjCvTEDLl0M2n9DSoaI8sJA5Cs/gYv7G8w398Ql8bIWAkDEJkeB6nS8Je//AX9+pnZUIXIw5Q4mCJJKyoyqBIHESlQJQ6ivQxxEwc/KUMGWeIQKGXIELPGSClEQBykSsP27dvx0UcfyXzIxOPLaiVuEBVYHEI0QMZhNYVIyuBFtjj4aYBk4iAPvw2QTBz2ElAcIt8AWQrLxUHoCD9v3jyhB1u+fHmoyUQahX9oExe2Asxc3MrZ2GCkv8Hkioqu/9xjbNfILm9USOtx8CsCOldUFKLDqrSW5shi1GzeY2zXyAxhVlWEShm8ROgCVyaWxNuGUCNkKpWCI3BSdF0XjuPkXNRKJlY3Qmqww7LioEAsfEmDxGWWUW2M9Js0eAkrDmGWWcoQh6DpgQxxCLPMMqw4hF1mGUYcgpQmChFEHKRJQ/YBxf4OJlMGrdJgoDFSpBFS6Cjftm1bjBkzBt/5zndK3m/BggX45S9/KTa7OKEpTiqZOChKIoxcShvRTBzCCEOSMZ04mCaKiYN0YQCsTxy0pwyW7uEgJA2DBw9Gc3MzzjnnnLL3TaQ0aMREqUJIHBRs5hRFcQhDmDJF2M2cwpYpTPYohN3MyXSZAggmDrJShgxWbABluThox0JxEPrrHHnkkVi8eLHQA0Zo2wc52NC0okEiuKKiNLJShjiuqBAhqo2RJneAjC0lmiNjt8wygghJw2WXXYbJkyeXvd+4ceOwatWq0JOKFAYkydQTuKg4KNwy2uSKCi7FFEPWCT/I48jcMjpKSzFlpwwZRFZUKClN5MPlmC1Y+CZc6BU3aNAgXHjhhWXvV1tbiz59+oSeFClPVhwMlCp0Y/tSTBW9DH7EQfZ1Jpg4iKEiZeBSzOIwZbAD6cWjpqYm2Q9pP4ZsMElP5jjtGClKFBIHFSf5KImDCmwWBy0pQ4akpw0WpgxAAGmYNGkSdu7cWfBrq1evxgknnBB6UsRuctIGjVezpDi0RuXVLG1PHFRezbKcOKjuZSglDqpKE/lYkTjsFYdYb+YUMXy/6p555hkcc8wxeOedd3Juf+6553DUUUdh48aN0iYXKUylDbt3GxmX16hoIcnLLFUnAkwc7E0ctGIwcTCW5lqaMgABpOGtt96C67o49thj8eijj6KpqQmXX345zj77bJx44on45z//qWKe0cDiP7QKKA56KJY2qEwZMphMG4Di4qAyZfBSSBx0rpjIFwddKYOXjDhoLU14SG/bbmRcUhjfr7xBgwZh0aJFmDBhAr73ve+hT58+ePjhh3Hvvffi+eefR8eOHVXMk5TAVNqQRLzioDNliEJ/A1EDEwczMGUoTCBdr62txfjx41FdXY3169fjsMMOw7e//W3Zc4smSStTbN1qZFxnrbkymA2Jg46UwUu+OOgsHeSPpStlyGC6TAG0iIOJlCGDu3GzkXETlzJYLgxAAGlobm7GVVddhbPPPhtjxozB1KlTsXbtWgwZMgRz585VMEVCCmNSHEz1MiQ1cbClvyHJmzm5zeblSRdJWpnmF9/ScNJJJ+GBBx7AlClT8MILL+Ccc87B4sWL0b9/f4wdOxa33367inlGC6YNarFgWW+7Febe9ZkWB1MncFvEwQTV72+As9NMicCbMugUB1MpA8sSpfEtDR9//DHmzp2Lq6++Ontbjx49MHv2bFxzzTW47bbbpE6QRANt4pAnDCbShtT2liXHpsShdr3ZyLbTMnNbivd/zFxtvePbn6HNp2Z/97rFoVBZIkmJA2mNb2lYvHgxjjvuuNYPlErhjjvuwKuvviplYpEnYWlDUsgIgykywtDpbTN9Ffu/b04YOi1tGbv9Kv1/g45vf6Z9zAzV72/I+dxU4qATpgz24lsaOnUqfdXBMWPGBJ4MiTbK04YiZQmTvQ0myxS6xcErDCbTBtOYTht0Uar5kWlDcvG1nd+KFSvw3HPPYdmyZWhoaIDjOKivr8cRRxyBs846CwMHDlQ1z2jiukaugunu3g2nTRvt4yqlTB+Ds3Yj3F7dlE7BlpQhiWRShgztV+3El31rtYxdKGVo8+l27O7STvnY+SlDBmdnI9xas5exdpvTcCrkr2RhymA3jitwLeumpiZcdtll+O1vf4t0Oo3OnTujvr4eALBlyxZs3rwZqVQK3//+9/HrX/8aFYouojQ2da6Sx1WKoUtnm5SG1H77yX1AwcZHldJQThi2D+yqbGygtDA0fLWD0rGB4mWJhsPbKh87Xxi86BCHUqUJleJQTBi8qBIHP0ssZYqDySWWlAZgRnpq2fsI/bUnT56MRx99FLfffjvWrVuHTZs2YcWKFVixYgU2bdqE9evX44477sCjjz4qdAntRMHeBm2oKlOIJAwqyxSmE4ZSfQxxL1OY7GUQQUV/g6k9GUxCYRBHSBoeffRR3Hnnnbj22mvRvXv3Vl/v1q0bfvazn+H222/H73//e+mTJNFCam+DBcsrbcdUU6QOSqUMgNqmSBFhUNXfIJIy2IKs/obEbeQUUYSkYd26dTjmmGPK3u/YY4/Fhg3RebJrg2lDMAIIg+y0wU8fg4q0QTRlUCUOIqslVKUN5YQhg4nVFF5MN0bKTBuCpgxRboxkyuAPIWno27cvZsyYUfZ+06dPx0EHHRR2TkQikd3wyYKEIemNj36WV8atTGHTEksRbFiGGUYcmDJEByFpuPTSSzF58mRceumleOutt7Bz576D6a5du7Bw4UL8x3/8B6ZMmYL/+I//UDbZSBNRqwxDlK9LEVQYkrQEUyWiKUMGmWlDEGGQlTaYLEuwl0HnwNE9HwitngCAO++8E3feeScaG1uMtrq6Go7jYNeuXdnPb7jhBlx33XXKJhvJ1RNeDK2kAMytpgi0kkJSyhBmNUXYlCHsaoowKYOM1RRBN3GSsZrCrzB4kbGaIkzKEHY1RVhpCLOaQpY0+F1NwRUT9iCyekJYGoCW5ZWvvvoqli9fjoaGBgAtmz0dccQROO2005RfFjvy0gBwCWY5JJYlgkqDjLJEGGkIW5YIKw1hd30MKw4mpSFsWSKMNMhKGYKIg+yUQVQcKAx2ISINvjZ3qq+vxwUXXBB4QsQckdjwSXIfQ5ANn2T1MbRbsUn53g3F6PT2F1r2blBBGGEA9G76VAhdmz6Vwu/GTyrKEqo2fpIFr2IZnEB/1d27d2Pjxo3YtGkTdke9Q183FlumKoR6GyxofLSBKDU/FsN0U2TQ/gZZzY9B+huitMRSFolsfozB8V9YGhoaGnDttddi4MCBaNu2LXr27IkePXqgbdu2GDhwIK6//vpsyYKUgUswteGnKVL2agm/TZEyhSFIU6TMi1EFEYewKUMYorZaohyiqylUNj/augyTKUM4hHoaVq1ahRNPPBGffvopRo0ahcGDB+dsI7106VLMmTMHXbt2xeuvv46+ffsqmWwsehoysLdhHxpShnJlCpXLK0XKFKoSBj9lCtlXsPTT26BCGPyUKVRIg2iZQmXKUK5MoWPFRKEyBXsZ7ERaT8NVV12Fjh07Yv78+TjwwAML3mfNmjU444wzcPXVV+Ppp5/2N9MkksCLWaW3bm0tDhaUJUzvx2ADKi553WnZDi3XpiiGaH+DqpRBpL9BdVmiVH+DriWWNvU3MGUIj9Bfcvbs2bj99tuLCgMAHHjggbj11lsxa9YsaZMjRBY2Xz5bZR+DSJlChTD4Ic5liVL9DUndk4G9DNFGSBqamppQW1ve2Gtra9FkwTvHyJDA3oacpkgLniu6UoZi4qCj8dHkpk/lehtUC4PpLaZtwKbdIo2WJZqazBxzYyQMgKA0DBs2DFOmTMH27cX/4Nu3b8eUKVMwfPhwaZMj6jDeFGlAGPLTBpYl9KQMtq6m0NX8WCht0J0yeMXBVMpgTWNkzE7iuhHqabj77rsxatQo9OvXDxMmTMARRxyB+vp6OI6DhoYGLF26FNOmTcOOHTswd+5cxVOOGYZ6G0yS/uwzpNq3NzJ2Zu8GE8KQv3eDzuWVhfZuSHJZQjdR3L9BNuntLX9vJ2Wglyv/TYqu424MBUVIGo466ii89dZbuPHGG/Hoo4/mXHsCaClLnHnmmbj11lsxYMAAJRMl8jHaFPnll8bEwQZM7MdgctOn/KZI3cKQ3xQZtyWWothwfQk37RoRh9YTSd4bNhn42kYaAJqbm/HBBx/kbCN98MEHo6KiQskEvcRqyWU+SVmC6eleNiUN7u7dSHXtbGRsAEh3MPeOMyMNplKGjDiYShm+7FtrVBicz780NjYApBu2AFVV+sfdnvv31ikNrVKGnIkonEcEUwbp154wTaylAYi/OBRY7qRbHLy9HKbEwf1yG9wDgl9MKyzNbfWfNHIweMSpWr0Rbr3BLbbXbYLT1sw21+mGLfs+0SgO+cKQQYc4lBQGLyqOvdE5tWaRfu0JANi4cSOWLl2KhoYGpFIp9OrVC0cffTSqq83Vykg4tJQpLFgfbbz5Ey3CYJLUp58jBWBPny5Gxq/6uAF7DuhkZGzjrGtZQePu2KldHHKEAQD27DGSOHixpkyhgggKgyjC0vDmm2/iqquuwvz581t9rV27dvjBD36A2267DW3bmtvMJfLEtcZWQhh09TYUEob0ps1a0wavMDifbDSaNpig6uOWkmbVJ2bEoWp1y+oZZ8sX+tOGdblLbk2IQys0iEOxlCGDSnEQThmA+B57FSAkDfPnz8fJJ5+M9u3bY/z48aiursabb76Jjz76CFdffTUaGhrwu9/9Dn/7298wZ84coT0diF1E4iqYJBSpTz/Pflz10afG0gbAnDhkMCIOhmiVMliGNYmDLHGIccoACPY0jB49Go2NjXjttddQV1cHoKUh8pJLLsGKFSvwt7/9DWvWrMHQoUNx6aWX4pZbblEy2dj3NGSIU2+DYFlCZdpQriyhI20oVpbQlTZ4hcGLLnHIpAw5Y2uUhkzKkI8WcVhXfEdQ1WmDkDAoShvKpQz5yBQHXylDq4mEnEeEpUGkp0Foc6eFCxfipz/9aVYYAKCiogI333wzFixYgE8++QQHHnggrrnmGjzxxBPBZ0yMIr3m76OPIf2lmq5ykZ8pvUntMrRSfQzOJ+a2t9ZFIWEAWtIGLeMXEQYbcHdYsMHYnj2mZyCVUMIAhDvpR1gYRBGSBsdxkEq1vmtFRQVc18XWvVsDH3nkkVizZo3cGSaRODzx2PhoDcVSBqClTGESXeJQDGeL4i22S6QMqvFVlpAsDn5TBqClTBFp4nDcFkBIGoYPH45f/OIXaGzM3cP8zjvvRPv27dG/f38AQGNjY04aQaKHlBNtQGFQlTYIja0obRBZLaEybSglDBlUikOxlEEXIimDMnEQFAYr0gaJBBGGDGHFIXTKkH2gZAhAEIQaIW+//XaMGDECffv2xciRI1FdXY033ngD//73vzFlyhRU7a2JLVq0CIcddpjSCScGg928cWiKtCFl8LO8UsVqChFhUImoMKhqirS5LJGP7NUUgZofLViGCUS0MTJBkiGUNBx77LGYN28eBg8ejJdeeglPPPEEamtr8dvf/hZXX3119n7nnHMO/vd//1fZZEkECFmWkJE2BBUGmWmD6f0Y/MIyheS0IUBZQlbiEGq1RMgyRZiUISzSUoacB02ODIjCHSFtx+DaYd9pg6Q+hjArKWQkDDJWUwSVBllpQ9CUQdZqiiBlCZlpQ9CUQcpqihB9DDLSBilLLAMmDjKlwW/aoEQagPLH4OicQssibfUEIWWR2PhosrdBBmFSBhn9DabLEkGRlTaEKUsob4wsQ9i0weSeDLJTBj/9DcqEAYiVFMhAqjS8++67uO2222Q+JDH4hBV+165gpUQQcZDVxxCmTBG1skQ+MsoUYZofTZcpQiNhtYQVjZE+yxSqyhIi4qBUGLKDFJlHAoVCqjS88847uPXWW2U+JCFCyG58DCIOsoQhTNogI2UIIw5RWC1RjsBpg8TllUHEQXrKYMn+DdYsxUygIBSC5YkoYHPaoHA/BtG0wYaVEjYQ1bJEPkHTBpmrJUyXKfyirCwhIA6xa34sOaBb+OMEIbTksl+/fkIPtmOHuSdP7LFxCaYFGzipxM8FrWSXJUxf0CrItSlkpgymr03hGwWbOFlxUSvAiqWY1izDBBIrCxmEpOGjjz7CwQcfjAEDBpS834YNG/Dpp2aXbhFNaBKGclfBtCFlUNXH4EccVKQMfsRBRVnCjzio2JPBhotaiYhDnJofS5EvDtpTBgJAUBr69++PYcOG4Q9/+EPJ+z399NP45je/KWVipAA2pg0aKCYOOoShXNpgQ+NjXMoSQVG5iZOQOCjeKrqUOGgTBgvSBsCyxCGhCPU0DBkyBEuWLCl7P4fXI08GFpQldCYMqi9oVYpyTZGqhUGkKVJl86MNqymi1t+gjLz+BlO9DG7aZcpgECFpGDduHLp1Kx+TDho0CDfddFPoSZESmG6KNCQMNu7doCtlMH0lzFLioGO1RClxML5VtKYLUhVaTWGkLLFXHEw2PyJt/k1LkuGOkFHEVKLjunAqhSpaSsiUKEz1MXjLFLrLEoV6G3SWJQr1NuhcXlmot0G3MLQqUxi4gqW3TGGsl6GqitIQU7gjZFwx4Xl7xzQZC6a//NJo46NNZQrdfQxJvzYFYEeZIpM4GG1+3Gawj4fCYBxKQ1QxWaYwJA5u2kV65y4jY+fMw1Dzo01lChObOHnFISlliUKYFIYsbtr0DIghKA2kPNGpYCmn+ZP1pqdgfLVEHHZ9DIqz5QujwgAYTvsaGz0T0SwOTBmsgNIQZXSczIuMofvA5d1K1lTa4O62YFvdf682OjzLFIC709y1ITINwVw9QExBaSDFsSRhsGbv+b2kvzCzkiNbz/7scyPjA0D68y+Q3vKZufG3fIb0VnMraTK/e5PikEG3OOSkDNlJaEobmDJYA6Uh6qg6sQs8ro6DVjFh0J025KcMpsTBJOnPzTYCemXFhDiYlDXA7LLjgsKQQbU4UBisgtJAWuNDRFSKQ7mEQZc42FCWyF+nb/wEZjBtsAWdaUMxYbCmTMHGyMRAaYgDlpQR4kgpYdCVNhS7TLJOcSiUMugUh0Jj6Uwbiv2uk1CmKJkyqIYpg3VQGuKCLHEI8DgqDlqifQyml2CqFodiwqATm8oSrb6mQRzKyZlqcRApS1iRODBtSASUBrIPSxILv42PqsTBhrJEOVimiDemt0/3nTLIFAemDFZCaYgTYU76IYVB1jsdW1ZK+BEGVWmDaMqgUhxEUgaV4iDy2CrTBtHfbRzLFIHLEjLEgcJgLZQGIi1hCHvQCiMMMtMGGxIGliX8yYjJZZiqCJIyWFGmAFiqiDGUhrhhSYkhaZhegik7bfArDHEsU/j9ncpMG8KUJWSIA5sfSTEoDUlHsmQEPWDJKEvISBvCpAyyxMGGlCEIMsUhyGPJTBuCSpgNZQprYNoQSygNcURUBCxJJWT2MYQRh6iXJWSlDVEqS7T6XgvKFGHFQUbzY5i0QWrK4FccmDJYD6UhrpQTAoXCYE1d1QeyhCFM2iAjYQgrDmGFwYYyRVhxML0iRRZBXodKyhJMHGIFpSGJaEgYRA9YKlZL+E0bZCcMpvsbTBMqKTAsHbKEIWjaIHuJpTUCLyIOTBkiAaUhzhSSA0tKEoDa5ZWmN33yi8w+hqAnPplliUA9CTJ7IiJYplC1J4OoOLD5kYggVRpSqRQOOOAAPPTQQ2iyxXCJMUodrKK4H4Mf/KQNKhof/YqD6T4GFfgVBxVliag0RmoRBpYpYoFUaTjppJPQo0cP/PjHP8Yhhxwi86FJULzJgiUpgy5hKJc2qG58FBEHG1ZKqBIGX/ssxKQsEXh8xTs/Wl2mYMoQKSplPtjcuXMBANu2bcO8efNkPjQJiyFhcJua4FRKfZpJwYaVEqpJf/Y5Uh33Nz2NsijdUXLrl0jt117Z44vg7twJp7a26Nd1bRVd7LWovSzhpgGHlfGoouQvV1dXh3Hjxql4aBIEwwmD912O7rKE6d6GUmmDjpSh3Dto1WUJ0wkCUL5MoSNlsKVMYV3iwJQhcviWhuXLl5f8+gsvvBB4MiTemOpjyBcH3SlDIXHQWZYodlLU1cdQ8iqVmqSimDjoLEsUEofIXZBK6uAUhijiWxqGDRuGRx99tNXtTU1NuOKKKzB+/HgZ8yIxw5Z3ODaUJWzoY7ABG1IIk5gShsxr0agwkMjiWxomTJiAiy++GBdddBF27jXn1atX4/jjj8eDDz6Iu+++W/okCQlLeucuo8Jgeu+G/HfUuldL5AuCCWHITxtMND+yTJGZgB1N2cQ/vqXh0Ucfxe9+9zs8/fTTOOaYY/DrX/8aQ4YMwcaNGzFv3jxceeWVKuZJ4oDhONJtMpsypL/40mjKkDlJmlpeaUOykBEHk6sl3J07jZcl3KYmnrhJIBzXDfbMWbJkCY477jg0Njbi6KOPxowZM9ChQwfZ88thbOpcpY9PNJGq0D/m3sYrp7JK/9iZKTQ3I1VTbWx8AIDjmB3fFprNCmy6sRFOhYHXAQqkDLqfE5QVa5mRnlr2PoFWT3z88cf40Y9+hKamJgwePBiLFy/G/fffH+ShCFGPZ2248bRhl9k6sunVJK4FdXTTv4NE9xJQGCKPb2l4+eWXMWTIEKxduxavv/46/vGPf+Caa67BbbfdhrFjx2LTpk0q5knihM4yhSW70LmG39kC+4TF9EnTpDhkfva0of4WrzCYeE4U7GXgiZz4wHd5oqKiAqeffjoee+wx1NfXZ2+fMWMGvv3tb6OyshJr166VPlGA5YlYoaNEUUIYdJYpCp0cdJcp8hOOVG2N1vGB1rLgVOsv1eQLU6qN3nJVoZRBZ5miZAOk6jIF5cR6lJQnJk+ejJdeeilHGABg7NixWLJkCbePJmIkpCmy2LvJpJUpWJYoXpbQlTiUXTHBkzoRIHAjZDHS6TRSKTVbhDJpiCGqEgeBsoSOtKHUCUFX2lBKUHQlDsWkQVfaUEoYdKQN5foYVKcNvpZYqkgcKCSRQFkjZMkHVCQMhAgj2MegOm0o9w5SR9pgOtEASqcMViQQNmz4ZUHPizIoDLFC6EpCo0ePxoMPPoiBAwdi9OjRJe/rOA5mzZolZXIkAaSb5aYNPhsf3aY9ShIH0ZNAelejssRBRBjSO3cpTRtEpMBtbFSaOIiUJdK79yhLHERXS7jNzUoSB98bObkul+aSogjFAt4KRjqdhuu6Rf+l03Z0qxNiiqi9a1RV67ciRbC0j6EYsp87gXd+lJUOMGWIHdJ7GlTCnoYYIyNtCLG8UmbaEOTALztt8FuWUJE2+JUG2WlDEGGQnTYE3ZNBVuIQervoMIlDdE4tZC9GehoICUTY1RQh92OQ1d8Q9J2izN6DII8l+x15kJTBimRCYn+D6U2cjF9fgsQSoZ6GfJqbm/HUU09hzpw5aGhoQKdOnTBq1Cice+65qKwM9JCEBIcbOGUJIx+q+xt0EkaCZPQ3hBUGVf0N/icSsL+BKUNs8V2e2Lx5M0477TT885//RGVlJTp16oSGhgY0NTVhyJAheO2119C5c2clk2V5IgH4LVNIFoagZQpZwhC2TCEjsQgrDmETg7BlClmpSRhxkJUyBBUH6SmDH3GgMEQWJeWJK664Au+99x4ef/xx7Ny5E+vXr8fOnTvxxz/+Ee+//z6uuOKKQJMlBIC/MoWChCFImUJmwhAqKbB8eaXOxzCJzLJEkOeWkrIERYDsxXct4cUXX8Qdd9yBb33rW9nbKioqMHHiRGzatAm33HKLzPkRQgSQ2hMRsEwh82QfdBmmzN6MIGUKFX0MfkoVSvsYREoVlIvY4ztpcF0Xhx12WMGvHX744YjQYgxiKyJpg8I+Bj9pg4o+Br8CoCJhML1UMQgq5mzDxk9WweN74vEtDWPGjMHMmTMLfm3GjBkYOXJk2DkRUhoNjY8i4qCy8dGGUoMfVJQUbClTiIqDytUSIs8146slKBSJwHd54sYbb8TZZ5+N5uZmTJw4Ed27d8eGDRvw+OOPY9q0aZg2bRq2bNmSvX/+ha0IEaLYTpFcKZGDSrkQLVOoPLmLlilMJyM6lldavaKCwpAYfK+e8F5bwvE8cTIP4+Q9mZolHly5eiKBeMXBgDAUW02hSxpKrabQlUaUEgddaUApcdAlDKX6G3TuyVBIHIykDN5jPaUhFoisnvCdNNx0002txICQuFLo2hQ6U4Zi16aIWvlCFToThmKNkYndxCmTOFAYEgW3kSb2k6owWpbwSoOJsoQN0lAobdDdc1AobTBRlvCKgylh8KYNxnsZSGyQtk/DkCFDcOedd+Ldd98NPSlCfBN2i+mQZJoiTfUx5AuCiZQh/+Rsokkxf0zTfQwmyTwXKQxEN0LSMGrUKPz2t7/F4YcfjkGDBuHGG2/EkiVLFE+NEHuQdW2KoGREwWRZwoaTdEYcTM4ls5oisWUJkmh8lScWLlyIZ555BtOmTcPKlSvRt29fnHPOOZgwYQKOPfZYlfMEwPJE4jHZS+O6cq7EGXFStTXGl0I61dXmBSbdbPb5CLCXgEhHpDwRuKfh7bffxjPPPINnnnkG7777Lnr16oWzzz4bEyZMwAknnKCkWZLSQIwcqL0vEZPiUGwZqk7ctPFlf25Tk/nfQ6ZkZkocKAxEAUqlwcu7776bFYi3334b3bp1w/r168M+bCsoDUT7QTr/5WHqZOXt6zA1B08zqilxyInkbfhbmJZYQiSiTRq8fPDBB5g2bRquvvpqmQ8LgNJA9qLrQF3spaH7ZFWoEdTECZPSUPhvYVpkCZGEEWlQCaWBADAvDYC+E1axlSO6T5gFlrzqFoeCjX82CBxgx3OSkJBI29xp9OjROZ/Pnj072IwIkYHI1fZkjGEzOvsbiuyRoXNb46IrBXT+Hkot/eVzkiQEIWno06eP6nkQ4g+VB2mRg7OOk1W5/SksaIzUIQ5llxba8LcgJCEIScMjjzyieh6E2IGfd3MqT1a2nKQsuUBYJDAtsoRowPelsQmxhrgeSP0Ig0q5EBQGlTtlCm9gpPL34OexVTwn4/o8J5GE0kCijcwDapDHsiERUDEHnwmDCnHwveOhit+DDX9fQiyC0kAIEE4+ZJ5YeJIKhw2/P9MiS4hCKA0k+oQ9sNpyYA5zwpN5sgzYxyAzbbDiugphfqcynlO2PC8J8UBpIPEg6AFW1oE57ElbxklfxmOEbHyUIQ6hhcGa3yVP+iR+UBoIkYUN0XgYJK2UMHUJ8Rxs+VuYlllCJENpIPHB74HWlgMzeyKySC1LBP1dmP4d2vK8JKQAQtKwadMmNOe9e1i6dCnGjx+PHj16oFevXpgwYQLeeecdJZMkRBjRA66qA7PfE44tHf+S92MIkjZEvo+hGJQAEiOEpKFHjx74xz/+kf18+fLlOO644zBr1iwMGTIEgwcPxvTp03H88cdj5cqVyiZLiBRUH8RFTzym39FmsGADJ2XCYM3v2LDMEiIJIWnIv6bVTTfdhI4dO2Lp0qV45ZVX8Oqrr+Ltt99Gu3btcMcddyiZKCHClDrw6joomz5ZiY6vUBis6G0A7JG4cs89CgOJAIF6GubOnYtrr70255oUffv2xdVXX41Zs2ZJmxwhsUWHVJgWF4iJg5ayhMh1PAghZQkkDV988QWOOOKIVrcPHjwYmzZtCj0pQkJT6F2b7ndyxU5EOk9QJa/MqKcsUUocYtvHUIxiz0GmDCQiCF2wCgDee+89VFa23L1Lly7YunVrq/ts3boVtbW18mZHSBi8FxAydVDOv6iVLe9oLehj0I4FVwUF0PrCVhQGEiGEpeGiiy7Kfuy6Lt544w2MGzcu5z7/+te/0Lt3b2mTI0QKST8o558sDQhDoUtoG0kZbJE4lVfEJEQhgS+N3b1791a3zZs3D6ecckr4WREiCxuEIXOiMpkyWPAu2ysORssSNvw9ADuem4T4xHHzl0ZYzNjUuaanQEh0SVXYUZaw4ZBjgzQQYhkz0lPL3oc7QhKSFHiS3Ad/F4QEQqo0fPzxx1izZo3MhySEyMT0u3zT4xNCQiHcCClCv3794LoummxYRkUIKYypJjwKAyGRR6o0XHjhhUinLaiZEkLsgsJASCyQKg2/+93vZD4cIUQVXPJHCAkAGyEJSSq63v0zZSAkNvhKGlzXxZtvvolly5ahoaEBjuOgvr4ehx9+OIYNGwaH71wIIV4oDITECmFp+POf/4yrr74a69ata3XVS8dx0LNnT9xzzz04//zzpU+SEKIIlikIIT4QKk88+eSTmDhxIgYNGoTHH38cy5Ytw7p167Bu3TosW7YMjz/+OA4//HBccMEFmDq1/OYQhBCLUJUGMGUgJHYI7Qg5ZMgQHHvssXj44YdL3u8HP/gBFi5ciMWLF0uboBfuCEmIQmQmDhQGQiKHtB0hV6xYgYkTJ5a938SJE7FixQqRhySEEEJIxBCShvr6erz//vtl77dy5UrU19eHnhQhxACy0gGmDITEFiFpOPfcc3HNNdfgqaeeKrh5UzqdxtSpU/Gzn/0M5513nvRJEkI0EfaET2EgJNYI9TRs374d48ePx8yZM9G+fXsceuihqK+vh+M4aGhowLvvvott27ZhzJgxePbZZ9G2bVslk2VPAyEaCNrbQGEgJNKI9DT4ujT2K6+8gmeffRbLly9HQ0MDAKBTp0444ogjMH78eJx22mnBZysApYEQTQQRB0oDIZFGujSYhtJAiEb8iEN0DiOEkCJIWz1BCEkgoiJAYSAkMVAaCCHBoTAQkigoDYSQ4lAKCCEeKA2EkNIUEwcKBSGJg9JACPEPhYGQREJpIISUxysJFAZCEgulgRAiBmWBkMRDaSCEiENxICTRRGpzJ0IIIYSYg0kDIYQQQoSgNBBCCCFECEoDIYQQQoSgNBBCCCFECEoDIYQQQoSgNBBCCCFECEoDIYQQQoSgNBCSYBoaGtC1a1esXr3a9FTK8tJLL2HIkCFIp9Omp0JIYqE0EJJg7rrrLpx55pk46KCDsrc9+uijGDx4MGpqatC9e3dcdtllOd+zdOlSjBgxArW1tejVqxduu+02iO4R19jYiCOPPBKO42DJkiU5X3Mcp9W/hx56KPv1M844A47j4E9/+lPgn5cQEo5K0xMghJhh586d+N3vfodXXnkle9svfvEL3HfffbjnnnswbNgw7Nq1Cx9++GH261u3bsXYsWMxatQoLFy4EP/+979x0UUXoV27drjyyivLjvnTn/4UPXv2xNtvv13w64888ghOO+207OcdOnTI+fp3v/td/OpXv8K3v/1tvz8uIUQGLiEkkTzzzDNu586ds59v2bLFra2tdWfOnFn0ex588EG3Q4cO7q5du7K33XXXXW7Pnj3ddDpdcrxXXnnFHThwoLt8+XIXgLt48eKcrwNwn3322ZKPsXr1aheA+8EHH5S8HyFEDSxPEJJQ5s2bh6FDh2Y/nzFjBtLpNNauXYtDDz0UBxxwAM477zx8/PHH2fssWLAAI0aMQHV1dfa2U089FevWrSvZF7Fx40Z8//vfx//93/+hbdu2Re932WWXoXPnzjjmmGPw0EMPtepf6NOnD7p27Yq//vWvAX5iQkhYKA2EJJTVq1ejZ8+e2c8//PBDpNNpTJ48Gf/1X/+Fp59+Glu2bMHYsWOxe/duAMCGDRvQrVu3nMfJfL5hw4aC47iui4suugg//OEPcyQln9tvvx1Tp07FzJkzcf755+PKK6/E5MmTW92vV69ekWjcJCSOsKeBkISyc+dO1NTUZD9Pp9PYs2cPHnjgAZxyyikAgCeeeALdu3fHnDlzcOqppwJoaVj04u5tgsy/PcOvfvUrbN26Fddee23J+dxwww3Zj4888kgAwG233ZZzOwDU1tZix44dAj8hIUQ2TBoISSidO3fGZ599lv28R48eAIBBgwZlb+vSpQs6d+6MNWvWAAC6d+/eKlHYtGkTALRKIDLMnj0bb7zxBqqrq1FZWYmvfOUrAIChQ4di0qRJRef3ta99DVu3bsXGjRtzbt+yZQu6dOki+mMSQiRCaSAkoQwZMgTvvPNO9vPjjz8eAPDee+9lb9uyZQs2b96MPn36AACGDx+OefPmZcsVADB9+nT07NkzZ9mmlwceeABvv/02lixZgiVLlmRXazz55JO48847i85v8eLFqKmpwf7775+9bdeuXfjggw8wZMgQ3z8vISQ8jusKLrAmhMSKpUuX4qijjsKmTZvQsWNHAMD/+3//DytXrsT//M//YL/99sO1116LDz/8EEuWLEFVVRW++OILDBgwAKNHj8Z1112H999/HxdddBFuuukmoSWXQEsvRd++fbF48eJsGeLFF1/Ehg0bMHz4cNTW1mLOnDm48sorcdFFF+H+++/Pfu/cuXNx5plnYuPGjSUbKgkhamDSQEhCOeKIIzB06FA89dRT2dsee+wxDBs2DF//+tcxYsQIVFVV4S9/+QuqqqoAtOybMGPGDHzyyScYOnQoLr30UvzkJz/BT37yk+xjrF69Go7jYO7cucJzqaqqwoMPPojhw4dj8ODBuP/++3Hbbbfhvvvuy7nfE088gQsuuIDCQIghmDQQkmBeeeUVXHXVVVi2bBlSKTnvIebOnYvx48fjww8/zCYYMvj0008xcOBALFq0CH379pX2uIQQcbh6gpAEM27cOLz//vtYu3YtevfuLeUx//KXv+C6666TKgwAsGrVKjz44IMUBkIMwqSBEEIIIUKwp4EQQgghQlAaCCGEECIEpYEQQgghQlAaCCGEECIEpYEQQgghQlAaCCGEECIEpYEQQgghQlAaCCGEECIEpYEQQgghQlAaCCGEECIEpYEQQgghQlAaCCGEECIEpYEQQgghQlAaCCGEECIEpYEQQgghQlAaCCGEECIEpYEQQgghQlAaCCGEECIEpYEQQgghQlAaCCGEECIEpYEQQgghQlAaCCGEECIEpYEQQgghQlSangBJFrt27cLu3btNT4MQIpE2bdqgpqbG9DSIBigNRBu7du1Ch9qO2I1dpqdCCJFI9+7dsWrVKopDAqA0EG3s3r0bu7ELJzhnoCrVBnBScFIO4DjA3v+dVGrv56m9n1dkv+a9HQ5yvg9Oat/HyH1M13FaCnHOvsfJvQ1wvffd+9ju3q+7jud+e7/PzX6+9z5oGRapzP09X/N83jIO9j72vvsA3s8L/e/kfF70vpnHSWHfnArct9xjAUU+zv+elFv0MeG4LXPI+dncvMdxW33Nyb8fXDiplkdy8u7jOG7O/R3H3fvnc/c9ZbL3c3M+T8F7Xxcpx0VFKt1y+97Pvf9a/rx5tyPzcXrvfTK3pVHh7Ps45QAVcOHsvV+V09wyHlo+r3DScNDyf8pJowLY+//e78/cD+mW+2TuizRSwN7HSe/9uovU3vtVohkVmcfZ+z0Ve3++CmQe30UF3L0fo2XOACocoALO3o8dpLDvX8vnKaTgYPs2oO/RH2H37t2UhgRAaSDaqUQVKp2qFmnInPD3HuGd7MflpCFXDHKkIe9r5aXBgVsRUhoyJ7iA0lBcFrBvXBFpyHwcVhqKCUT+fctJw97vlyINTnFp2Hd7rgiUlQbv52WkYZ8gFJcG7+e50uDuEwDHRZWTaiUNqRxpcPP+Ly4NFY7bcnJ3HFTA2Ttuy4m9Ci0n/BZpSO8VgYwooKA0VJSQhooC0lCR+RuTRMBGSEIIIYQIQWkghBBCiBCUBkIIIYQIQWkghBBCiBCUBkIIIYQIQWkghBBCiBBcckm004Q9cFrW6+39f99aPSe7FrHlfyddgewaO6fle3L2ach+LQXPWryW/13PMkoXrZdcZm8DXO99gyy5BEItuQRyfg3cp0HiPg1u3pLLzOfu3n0a0p4ll66mfRrSsdinAdl9GkhyoDQQbbiui7q6Ovxt20tAs+nZEEJkUVdXB9d1y9+RRB5KA9GG4zjYtm0bPv74Y+y3336mp0MIkcDWrVvRu3fvlo3aSOyhNBDt7LfffpQGQgiJIGyEJIQQQogQlAZCCCGECEFpINqorq7GzTffjOrqatNTIYRIgq/rZOG4bHklhBBCiABMGgghhBAiBKWBEEIIIUJQGgghhBAiBKWBEEIIIUJQGkhRXNfFLbfcgp49e6K2thYjR47E8uXLy37fM888g0GDBqG6uhqDBg3Cs88+m/P13/zmNxg8eHB2k6fhw4fj1VdfzbnPRRddBMdxcv597Wtfy7nPBx98gPHjx6NLly7Yb7/9cN5552Hjxo3hf3BCIsa8efNw5plnomfPnnAcB88991zJ+69fvx4TJ07EgAEDkEql8J//+Z+t7vO///u/OPHEE9GxY0d07NgRY8aMwVtvvaXmB/Dw+uuv4+ijj0ZNTQ369euHhx56KOfr06ZNw9ChQ7H//vujXbt2OPLII/F///d/yudFWqA0kKLcfffd+MUvfoFf//rXWLhwIbp3746xY8fiyy+/LPo9CxYswDe/+U185zvfwdtvv43vfOc7OO+88/Dmm29m73PAAQdgypQpWLRoERYtWoTRo0fjrLPOaiUkp512GtavX5/998orr2S/tn37dpxyyilwHAezZ8/G3//+d+zevRtnnnkm0um0/F8GIRazfft2fPWrX8Wvf/1rofs3NjaiS5cuuP766/HVr3614H3mzp2Lb33rW5gzZw4WLFiAAw88EKeccgrWrl0beJ5z587FQQcdVPTrq1atwrhx43DiiSdi8eLFuO666/DjH/8YzzzzTPY+9fX1uP7667FgwQL861//wne/+11897vfxWuvvRZ4XsQHLiEFSKfTbvfu3d0pU6Zkb9u1a5fboUMH96GHHir6feedd5572mmn5dx26qmnuueff37J8Tp27Oj+9re/zX4+adIk96yzzip6/9dee81NpVLuF198kb1ty5YtLgB3xowZJcciJM4AcJ999lnh+48YMcK9/PLLy96vqanJbd++vfuHP/whe1s6nXZ//vOfu3379nVramrcwYMHu1OnTi36GHPmzHH79OlT9Os//elP3YEDB+bcdskll7hf+9rXSs5tyJAh7g033FD2ZyDhYdJACrJq1Sps2LABp5xySva26upqjBgxAvPnzy/6fQsWLMj5HgA49dRTi35Pc3Mz/vznP2P79u0YPnx4ztfmzp2Lrl274pBDDsH3v/99bNq0Kfu1xsZGOI6Ts6FMTU0NUqkU/va3v/n6WQkh5dmxYwf27NmD+vr67G033HADHnnkEfzmN7/B8uXLccUVV+Db3/42Xn/99UBjFDt+LFq0CHv27Gl1f9d1MWvWLLz33ns46aSTAo1J/MELVpGCbNiwAQDQrVu3nNu7deuGjz76qOT3FfqezONlWLp0KYYPH45du3ahrq4Ozz77LAYNGpT9+umnn45zzz0Xffr0wapVq3DjjTdi9OjR+Mc//oHq6mp87WtfQ7t27XDNNddg8uTJcF0X11xzDdLpNNavXx/2xyeE5PGzn/0MvXr1wpgxYwC0lER+8YtfYPbs2Vnh79evH/72t7/h4YcfxogRI3yPUez40dTUhM2bN6NHjx4AgC+++AK9evVCY2MjKioq8OCDD2Ls2LEhf0IiAqWBAAAef/xxXHLJJdnPX375ZQBodblb13XLXgJX5HsGDBiAJUuW4PPPP8czzzyDSZMm4fXXX8+Kwze/+c3sfQ8//HAMHToUffr0wcsvv4yzzz4bXbp0wdSpU/GjH/0IDzzwAFKpFL71rW/hqKOOQkVFhf9fACGkKHfffTeeeOIJzJ07FzU1NQCAd955B7t27Wp1st69ezeGDBmS/byuri77cXNzMxobG3NuO/HEE3MaoQsdP/Jvb9++PZYsWYJt27Zh1qxZ+MlPfoJ+/fph5MiR4X9YUhJKAwEAfOMb38CwYcOynzc2NgJoMf+M3QPApk2bWr0T8NK9e/dWqUKh72nTpg2+8pWvAACGDh2KhQsX4v7778fDDz9c8HF79OiBPn364P3338/edsopp+CDDz7A5s2bUVlZif333x/du3dH3759BX9qQkg57r33XkyePBkzZ87E4MGDs7dnGo5ffvll9OrVK+d7vGXDJUuWZD9+8803cc0112Du3LnZ22pra7MfFzt+VFZWolOnTtnbUqlU9vhx5JFH4t1338Vdd91FadAApYEAaDH39u3bZz93XRfdu3fHjBkzsu8adu/ejddffx0///nPiz7O8OHDMWPGDFxxxRXZ26ZPn47jjjuu5Piu62ZFpRANDQ34+OOPcwQmQ+fOnQEAs2fPxqZNm/CNb3yj5FiEEDHuuece3HHHHXjttdcwdOjQnK9lllWvWbOmZCkic3IHgE8++QSVlZU5t3kZPnw4XnzxxZzbpk+fjqFDh6KqqqroGOWOH0QelAZSEMdx8J//+Z+YPHky+vfvj/79+2Py5Mlo27YtJk6cmL3fhRdeiF69euGuu+4CAFx++eU46aST8POf/xxnnXUWnn/+ecycOTOnOfG6667D6aefjt69e+PLL7/En//8Z8ydOxd/+ctfAADbtm3DLbfcgnPOOQc9evTA6tWrcd1116Fz584YP3589nEeeeQRHHrooejSpQsWLFiAyy+/HFdccQUGDBig6bdEiB1s27YNK1euzH6+atUqLFmyBPX19TjwwANx7bXXYu3atXjsscey98kkANu2bcOnn36KJUuWoE2bNtkS4d13340bb7wRf/rTn3DQQQdlE4C6ujrU1dWhffv2uOqqq3DFFVcgnU7jhBNOwNatWzF//nzU1dVh0qRJvn+OH/7wh/j1r3+Nn/zkJ/j+97+PBQsW4He/+x2eeOKJ7H3uuusuDB06FAcffDB2796NV155BY899hh+85vfBPnVEb+YW7hBbCedTrs333yz2717d7e6uto96aST3KVLl+bcZ8SIEe6kSZNybps6dao7YMAAt6qqyh04cKD7zDPP5Hz94osvdvv06eO2adPG7dKli3vyySe706dPz359x44d7imnnOJ26dLFraqqcg888EB30qRJ7po1a3Ie55prrnG7devmVlVVuf3793fvu+8+N51Oy/0lEBIB5syZ4wJo9S/z2pw0aZI7YsSInO8pdH/vcsg+ffoUvM/NN9+cvU86nXbvv//+7Ou9S5cu7qmnnuq+/vrrRedZasml67ru3Llz3SFDhrht2rRxDzroIPc3v/lNztevv/569ytf+YpbU1PjduzY0R0+fLj75z//WfRXRULCS2MTQgghRAju00AIIYQQISgNhBBCCBGC0kAIIYQQISgNhBBCCBGC0kAIIYQQISgNhBBCCBGC0kAIIYQQISgNhBBCCBGC0kAIIYQQISgNhBBCCBGC0kAIIYQQISgNhBBCCBHi/wcDq6YFLqc7mwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 580x640 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "hp.gnomview(source_map_healpix, rot=hp.vec2ang(hp.ang2vec(source_pos[0], source_pos[1]), lonlat=True), xsize=1900, reso=.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "cd1554a1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.995092243734167"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from pysm3.utils import healpix_aperture_photometry\n",
    "\n",
    "healpix_aperture_photometry(source_map_healpix, source_pos[0,0], source_pos[1,0], 2 * fwhm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7df1f3b3",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}