MaDCoWS VLASS Radio Contamination (1).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from astropy.io import fits\n",
    "from astropy.table import QTable\n",
    "from astroquery.esasky import ESASky\n",
    "from astropy.nddata import Cutout2D\n",
    "from astropy.convolution import Gaussian2DKernel, convolve\n",
    "from astropy import wcs\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle\n",
    "\n",
    "path = \"/home/alex/Downloads/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def background_emission(hmap, map_w, pixel_w, max_sep, get_background):\n",
    "    '''\n",
    "    Args:\n",
    "        hmap: data array\n",
    "        map_w: width of half the array/map\n",
    "        pixel_w: size of the pixels\n",
    "        max_sep: radius of map that will contain the central emission \n",
    "        get_background: determines whether to look for background values (boolean)\n",
    "    Returns:\n",
    "        hist: array of central pixel values\n",
    "        means: array of length 20 of mean values of random background areas\n",
    "    '''\n",
    "    # Returns the average surface brightness in Jy/pixel of \n",
    "    # pixels within max_sep of the center of hmap (assuming hmap has units MJy/pixel)    \n",
    "    center = map_w + 0.5\n",
    "    emission = 0\n",
    "    div = 0\n",
    "    hist = []\n",
    "    alpha = []\n",
    "    delta = []\n",
    "    for i in range(2*map_w):\n",
    "        for j in range(2*map_w):\n",
    "            if pixel_w*np.sqrt((center-i)**2+(center-j)**2) <= max_sep:\n",
    "                hist.append(hmap[i,j])\n",
    "                alpha.append(i)\n",
    "                delta.append(j)\n",
    "                div +=1\n",
    "    # Find background flux\n",
    "    # Randomly choose a point that is at least n RA pixels and m Dec pixels away from the center\n",
    "    # Shift the circle of pixels that we formed above so that it centers around the random point\n",
    "    # Calculates the mean of the pixel value in this new, displaced circle\n",
    "    if get_background:\n",
    "        means = []\n",
    "        for k in range(20): # Iterate 20 times\n",
    "            background_hist = []\n",
    "            sign_RA = np.random.choice([-1,1]) #Randomize direction of offset\n",
    "            sign_dec = np.random.choice([-1,1]) \n",
    "            offset_RA = sign_RA*(300.-10.)* np.random.random_sample() + 5. #Offset RA between 5 and 70 pixels away\n",
    "            offset_RA = int(np.round(offset_RA))\n",
    "            offset_dec = sign_dec * (300.-10.)* np.random.random_sample() + 5. #Offset dec 5 and 70 pixels away\n",
    "            offset_dec = int(np.round(offset_dec))\n",
    "            # Take the Circle of pixels around the center of the stamp from above, and shift it\n",
    "            # with respect to the offset you have randomly settled on\n",
    "            offset_alpha = []\n",
    "            offset_delta = []\n",
    "            for x in range(len(alpha)): \n",
    "                offset_alpha.append(alpha[x]+offset_RA)\n",
    "                offset_delta.append(delta[x]+offset_dec)\n",
    "            # Record pixel values of the new, displaced circle\n",
    "            center_RA = center + offset_RA\n",
    "            center_dec = center + offset_dec\n",
    "            for i in range(len(offset_alpha)):\n",
    "                background_hist.append(hmap[offset_alpha[i],offset_delta[i]])\n",
    "            means.append(np.mean(background_hist))\n",
    "    return  np.array(hist), np.array(means)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def central_emission(hmap, map_w, pixel_w, max_sep):\n",
    "    #Returns the average surface brightness in Jy/pixel of pixels within max_sep \n",
    "    # of the center of hmap (assuming hmap has units MJy/pixel)\n",
    "    center = map_w + 0.5\n",
    "    emission = 0\n",
    "    div = 0\n",
    "    hist = []\n",
    "    for i in range(2*map_w):\n",
    "        for j in range(2*map_w):\n",
    "            if pixel_w*np.sqrt((center-i)**2+(center-j)**2) <= max_sep:\n",
    "                hist.append(hmap[i,j])\n",
    "                div +=1\n",
    "    return  np.array(hist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "mdcw = fits.open('/home/alex/Desktop/Clusters/MADCOWSUnion.fits')\n",
    "\n",
    "mdcwRA = mdcw[1].data['RADeg']\n",
    "mdcwDec = mdcw[1].data['decDeg']\n",
    "mdcwName = mdcw[1].data['name']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: FITSFixedWarning: 'obsfix' made the change 'Set OBSGEO-L to  -107.618000 from OBSGEO-[XYZ].\n",
      "Set OBSGEO-B to    34.078827 from OBSGEO-[XYZ].\n",
      "Set OBSGEO-H to     2115.607 from OBSGEO-[XYZ]'. [astropy.wcs.wcs]\n"
     ]
    }
   ],
   "source": [
    "hdu = fits.open('/home/alex/Downloads/newMoo/'+mdcwName[0])\n",
    "w = wcs.WCS(hdu[0].header)\n",
    "header = hdu[0].header\n",
    "pixel_size = header.get('CDELT2')*u.deg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "stackArr = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# fill stackArr - 998 is weird for moo\n",
    "i = 0\n",
    "while i < 1:#len(mdcwRA):\n",
    "    cutout = fits.open(path + 'newMoo/' + mdcwName[i])\n",
    "    vmap = cutout[0].data[0][0]\n",
    "    #Convert from Jy/beam to Jy/pixel\n",
    "    conv_factor = (np.pi*(2.5/2)**2/np.log(2))\n",
    "    vmap = (vmap/conv_factor)\n",
    "    stackArr.append(vmap)\n",
    "    i+=1\n",
    "len(stackArr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cut excess rows/columns off evens\n",
    "i = 0\n",
    "while(i < len(stackArr)):\n",
    "    if(stackArr[i].shape[0] > 720 and stackArr[i].shape[0] % 2 == 0):\n",
    "        while(stackArr[i].shape[0] > 720):\n",
    "            stackArr[i] = np.delete(stackArr[i],stackArr[i].shape[0]-1,0)\n",
    "            stackArr[i] = np.delete(stackArr[i],0,0)\n",
    "    if(stackArr[i].shape[1] > 720 and stackArr[i].shape[1] % 2 == 0):\n",
    "        while(stackArr[i].shape[1] > 720):\n",
    "            stackArr[i] = np.delete(stackArr[i],stackArr[i].shape[1]-1,1)\n",
    "            stackArr[i] = np.delete(stackArr[i],0,1)\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cut excess rows/columns off odds\n",
    "i = 0\n",
    "while(i < len(stackArr)):\n",
    "    if(stackArr[i].size > 518400):\n",
    "        if(stackArr[i].shape[0] > 720 and stackArr[i].shape[0] % 2 != 0):\n",
    "            if(stackArr[i].shape[0] == 721):\n",
    "                stackArr[i] = stackArr[i] = np.delete(stackArr[i],stackArr[i].shape[0]-1,0) \n",
    "            else:\n",
    "                while(stackArr[i].shape[0] > 721):\n",
    "                    stackArr[i] = np.delete(stackArr[i],stackArr[i].shape[0]-1,0)\n",
    "                    stackArr[i] = np.delete(stackArr[i],0,0)     \n",
    "                stackArr[i] = stackArr[i] = np.delete(stackArr[i],stackArr[i].shape[0]-1,0) \n",
    "        if(stackArr[i].shape[1] > 720 and stackArr[i].shape[1] % 2 != 0):\n",
    "            if(stackArr[i].shape[1] == 721):\n",
    "                stackArr[i] = stackArr[i] = np.delete(stackArr[i],stackArr[i].shape[1]-1,1) \n",
    "            else:\n",
    "                while(stackArr[i].shape[1] > 721):\n",
    "                    stackArr[i] = np.delete(stackArr[i],stackArr[i].shape[1]-1,1)\n",
    "                    stackArr[i] = np.delete(stackArr[i],0,1)\n",
    "                stackArr[i] = stackArr[i] = np.delete(stackArr[i],stackArr[i].shape[1]-1,1) \n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# array for just 720x720\n",
    "sevenTwenty = []\n",
    "# array for indeces\n",
    "indexArr = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0\n",
    "while i < len(stackArr):\n",
    "    if(stackArr[i].shape == (720,720)):\n",
    "        sevenTwenty.append(stackArr[i])\n",
    "        indexArr.append(0+i)\n",
    "    i+=1\n",
    "len(sevenTwenty)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0\n",
    "while i < len(sevenTwenty):\n",
    "    cent_em = np.sum(central_emission(sevenTwenty[i], int(100), pixel_size.to(u.arcmin).value, 1.0))\n",
    "    back_em = np.median(background_emission(sevenTwenty[i], int(100), pixel_size.to(u.arcmin).value, 1.0, True)[1])\n",
    "    thisdict = dict('name': mdcwName[indexArr[i]],'source flux': cent_em,'background flux': back_em)\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "pickle.dump(dictionary, open(path+\"mooEmissionCatalog.pickle\",\"wb\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "done = pickle.load(open(path+\"mooEmissionCatalog.pickle\",\"rb\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'MOO J0001+1428': [-8.8941675e-08, 6.816236e-08]}\n"
     ]
    }
   ],
   "source": [
    "print(done[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ignore the rest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0\n",
    "dictionary = []\n",
    "while i < len(txtArr):\n",
    "    a = {txtArr[i][0:14]:float(txtArr[i][16:])}\n",
    "    dictionary.append(a)\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "textArr = np.array(dictionary)\n",
    "a_file = open(\"/home/alex/Downloads/madcowsCentralFlux.txt\",\"w\")\n",
    "np.savetxt(a_file, textArr,fmt='%s',newline = ',\\n')\n",
    "\n",
    "a_file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define stacker\n",
    "def s18dStack(ras, decs, width = 20., weight = False):\n",
    "    stack = 0\n",
    "    divisor = 0\n",
    "    count = 0\n",
    "\n",
    "    while(divisor < len(stackArr)):\n",
    "        if len(stackArr[divisor]) == 720 and len(stackArr[divisor][0]) == 720:\n",
    "            stack += stackArr[divisor]\n",
    "            count +=1\n",
    "        divisor += 1\n",
    "\n",
    "    print(\"Number in stack: {}\".format(count))\n",
    "    return stack, count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extent in arcmin for plot axes\n",
    "w1 = -720*pixel_size.to(u.arcmin)/2\n",
    "x1 = 720*pixel_size.to(u.arcmin)/2\n",
    "y1 = -720*pixel_size.to(u.arcmin)/2\n",
    "z1 = 720*pixel_size.to(u.arcmin)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number in stack: 1652\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUUAAAEXCAYAAAA6MlahAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9e7xtSVUe+o219j6n3zTQBrEbBIXI6wpXWgIakARuwDcazUWjKBftm+SCJhoNeo2iMV41RgXFa4giQRHloSYapMUoIiooLS8bMDzk0bxfTb/POXuvkT+qRq2vvhpzrd3nrHN6n+Ma+7d+e6056zGqZtWobzyqprk7trSlLW1pS4VmdzQDW9rSlrZ0mGgrFLe0pS1tiWgrFLe0pS1tiWgrFLe0pS1tiWgrFLe0pS1tiWgrFLe0pS1tiWgrFG8nmdk3m9mr72g+tvS3l8zMzew+dzQf5yqtFIpm9nIz+6Hk+lea2YfMbMfMnmdmP7ymnOeZ2Z6Z3V2uX2pmz61l3Whm/9PMni71vMHMbjCzj5nZH5jZvSfquMLMXlrTfcrM/srMvrneu1cdSDur+DzMZGbXmtlN9Nkzs9+u9/6umf1XM/uomX3CzK42s8+R/P+q9vMNtc+PHrDem+Szb2Y/U+893MxeUev8qJm9mJ+xFfoxM/t4/fyYmdlEPY+uz+g35fqD6/VXHpDfZ5jZiTqeYkz9bDL2LjGznzaz99Z2vbP+vmxN+Y+o5c7p2n+euPbz9fsDzez3aj9db2bXmNmXrKjj7mb2i2b2wVru28zsB83swoP0wUHoIPP2byutQ4r/BcA3JAP5GwG8wN331lVQH+Q/BvApAN8gt38KwEUA7g/gTgC+AsA7ar77AHg+gO+s9+4N4NkA9ieq+mUA7wPwmQDuWnn88Dr+zhZy9we6+0XufhGAi1Ha+uJ6+1IA/w3A5wC4G4A/B/BfI6+ZPQ7A0wE8BqV/PgvADx6w3ouo3k8HcCvVe2cAzwFwr1rujQB+ibJfBeAJAB4M4HMBfDmA/3tFdR8F8Agzuytd+yYA//MgvBL9urtfDOAuAL6q8n1NCEYzOwLgfwB4IIDHA7gEwCMAfBzAw9aU/TqUefN5dO2RAK6Ta48C8Kr6/bcBvKLy8XcAfBuAG7LCzewuAP4MwPkAHlHb8X+gPOPPXsPbGaOzGWCsJXef/KA8mE8BeBRduzOA2wA8uP5+HoAfXlHGk1Am8LcD+Cu591cAnjCR72sAvGEVf5L+JgAPmbj3XgBe09yEMgE+G8AfoEyEjwF4AYBLKc89APwGykT9OICfrde/GcCrKd1/APBqAHdK6j0K4KcBfKB+fhrA0Xrv0SgT6TsBfATABwE8+YBt/SIUAXThxP271Pbetf7+VQA/QvcfA+BD9fsX1Pbfo/5+MIBPArhfUu43AXgXAJuo9/MA3Ei//xTAVfT7KQBeM5E3+uPnAfw/9docwPsBfD+AV1LaZ9YxdQOAawA8ku49A8CvSNlzAG8E8BP197egLJgXrejj+wN4JYDrAVwL4Cvo3v8A8J31+9+pffKDcs0BXAHgsvr90qm6pN4fBvBmALMVaRzAfer3VwL4FrrXxicAQwEeH6l99WYAD0JZrE4AOI4yH367pv8MAC9FGfN/A+DbpF9fAuBXalnfgrKAvK7+/jCAnzzofD3Mn5VI0d1vBfAiFMEW9E8AvM3d37gqL9E3AXghgF8DcD8zeyjdew2Af29mTzaz+0q+v6zpf8rM/oGZXbSmntcAeLaZPdHM7in3HlX/X+oF9fwZyoD5/1AGwv1RhOAzAKCqQb8D4D0oKOjyyn8jM5uZ2X9GQUD/yN0/lfD0/wJ4OICHoAibhwH4Prr/6Sgo+HIUgfFsM7vzmnYCpU9f6u43T9x/FIrQ+3j9/UAUoRD0RgB3M7O7uvufAvhPAP6LmZ2PMuj/rbu/baLe53udJRP1Xku/s3ofuKJdQNEOYrw9DmXh/ICk+QuUPr0LisB/sZmdN1Wgu++jIOdH1kuPBfByd78pS29muyjo7vdQBNzTALyATBKvwnJMPQplUXy1XPsbd78OZUF9B4BfMbMnmNndVra+8PYb7r5Yk+4g9I8qL38XZZz9EwAfd/fnoICAH6/z4cvNbIbS5jeijMfHAPiXVcsI+koUwXhpzf9MAM9090tQQMaLNsDzHU8HWLn+PspqeV79/ScA/hXdfx4mkCKAewJYoCI4AFejdGLcPx/A96Ks9idQBs8X0/2Ho3T0R1HQ6fMwsbqjINgfRZmU+wDeAODz6717oayuOyva+QQAr6/fH1HrHNKjrMSvBfDrKKvqkRVlvhPAl9DvxwF4d/3+aBRVdIfufwTAw9c8jwtQVuZHT9y/AgVdfZ3w8Xj6vVv74170+xoUJPFyJEgQRT3eB3DviXo/F8An0KO2fRDiBHDfWm9W/qMBXFe/vx3FFPBrAP4pCip55Yo++SSWmsszIEixXv9nAN5ev78CwI+uKO+RAD4EQmsoC/sziNePoyyszwTwrShmoA/TtV+SZ/Kz9TksUITqfSfqfjuAf7ZmDBwUKf5DFNPDwyHIEzJvAfw9AO+VNN8T7aj9+iq5/yoUhHzZKn7Pts9a77O7vxpFvXqCmX02Ctr51XX5Kn0jgLe6+xvq7xcA+Pq6EsPdb3X3H3H3h6LYAV+Esurfpd5/jbv/E3f/NJSB+igU9JXx+Ul3f7q7PxDFrvYGAL+1wrB/NzP7NTN7v5ndgIKQwsh+DwDv8Wmb6X1QVs0fdPfjK9r/GShoM+g99VrQx6WOW1Am1yr6ahTh80d6w8w+DQXd/Jy7v5Bu3YRiNwuK7zcCgLufQJkkDwLwH72OeKFvRJlsf5PUex8Avwvg2939j9fUe9NE+Uy/DOCpAP4BgN/Um2b2r83srVYcatejoKCVDhIU9POJ+v3jAO6+Iu1nAHif92jtPbUMoGglF6H016MA/LEX1Pk+uhb2RLj7de7+VHf/bJTF5WYURJzROt4OTO7+ByjC+NkAPmJmzzGzSyaSfyaAz6iOoOtrv34vylwKep/keQoKCn2bmf2FmX3ZJvi+o+mgITmh0nwDgKvd/aAOjCcB+Kzq9fwQgJ9EGbyD583dbwDwIwAuRHGq6P2/QLHxPWhdpe7+MQA/gTK4w76m9CP1+v/mBf5/A8oqD5SHf88VxuS3AngygN818fIKfQBlsAXdE6MqeHspVWGr2v17AP6bu/97yXMtivoe9GAAH/aqXpvZ5QB+AMVJ8h8t90w/CcXx1pGZfSaA3wfw79z9lw9Q77VYT78M4F8AeJm73yL1PRLAd6Oognd290tR7N7p4lfzzFCcPCGwfx/A41Z4cz8A4B41X9A9URA43P02FBX+ywHc3Zemhj+u1z4XJBSZ3P19KEJqahz/PoCvkrpX0c0o2kPQp0t9z6qg4wEoAuy74paU8z4Ulf9S+lzs7jxXuzzu/nZ3/zoUE8OPAXjJJj3kdxTdHqH4WBQ1YZgYAOZmdh59jphZODMehmL/eQjKQPhVVJuRmf1bM/v8mv48FGfM9QD+2sz+vpl9q5n9nZr2fije6ddkDFoJ93iQlTChiwH8cwDvqBP/oyhqy2dRlotRkMynqlD4Lrr35yiOjx81swtrm76Q66tI7HsB/H5F0Bm9EMD3mdmn1VCP70dBpCdFZnYFCnr6L3L9EhTTxJ+4+9OTrM8H8BQze4CZXYpi13xezWv1+y+irPwfBPDvpPwvQEFJL5brl6M4q37W3X9+ot7vMLPLzewzUJxKz1vXzopGvwi5VnAxgD1U84aZfT96NMr87ZjZ/VGew6ejLMrAMlLhpWZ2v2ofvquZfW8NlXktCmr/bjPbNbNHowg7tiu/CmW8/ilde3W99kF3f2fl4c41nOY+tZ7LAPxfmBjHlcdLUGy8n1nLuNzMftLMPjdJ/wYAX21mF1TE/hRq/+eb2d+rmtnNKCaoQL8fRj8f/hzAjWb2b8zsfDOb1/n0+RN8wsy+wcw+rSLq6+vlTdhC71g6qJ6NYrv4JKr3lK4/D2UF4c+rUbyIL03KeRiAYygI7vtQDOk3oKg2rwTwBTXdg1AMvx9GEV7vRlmNdif4+xkUe8xNKBPmdwDcn+7/UL1+PYqN5YEodrSbUAbWd6LatGr6ewL4LSy9089ysdnU39+K6pBJeDoPwLNQBM0H6/ewzT6a66vX3g3gsSuewfegqGp6/Ztqv9+MpYf9JgD3pDTfUfvyBhREGF7wb0cxrh+pvz+j9hPbBv8TgF9O6v0B9F79m1DU47hvAH68PttP1O9TnuuhP+hesymieJKfW9vxQRTU2PoNxfZ1ovJycx0TPwfgcinzTijRAO+rad+JIpDCY/9AFBPFpwC8BcBXSf7H1bZ/B127W732Qrp2Icoi9u5az4dQhPTlWVvpGTy3pr0RwNtqX19Q77NN8TIUDeFGFHv/M7C0KT4GwJtqvRFhcVG9d1+UcX89gN+iel9Y6/0kiuDmflWv/q+g2MFvQtEA0kiSs+1jtXFb2tKWtrQlbLf5bWlLW9pSR1uhuKUtbWlLRFuhuKUtbWlLRFuhuKUtbWlLRFuheJaSrTn5x8qpOo8+w2zdLjKzV5rZt9zRfGxpS0xbobhhMrN7Wn/UlpvZzfT7kRP5NnpOo5dTdV55MnltxZFtVo7mOulYy02RmT3VzF5nZsfM7HnJ/QvM7OdseZTcq+heHC/Gz+mz6r21x7Bt6dymc/f4nzuI3P29oK16ZuYo+3LfccdxdXCy5ZFtX40SmH0RysECU0e23VH0AZQTZR6Hsode6Tko4/v+KDGSD5H7v+7uepQdsDyG7ckosX/fj3KYxP02w/aWDjttkeIZJDO7k5k9v6KQ95jZ99VdDvdHCXZ/REUt19f0X2pmr6+I7X1m9ozbUde7zeyx9fszzOxFte4bq2p95UTWh6Bs9/ofXuhGd3+pu7/XzB6Psovn/6x8vrGW/2Qre5FvNLN3mVl3ZqIgz3fWcpTfu5vZm8zsu/ReRu7+G+4ewfVaVux+usrdP+ru++5+zQHL/XN3/0V3/4SXPeE/BeBzrD/jcUvnMG2F4pmln0HZSfFZKNvYnoRyhuJbUU5x+TMvRzldWtPfXNNcCuBLAfxzM3vCSdb9FSjb1AIJ/exEuskj29z95Sh7xn+98hn7mj8C4MtQtqc9GcBPmdnnAYCZPQwFeX5XrftRKLs7GlXV/I9Qtgv+h3rt6Wb2OyfZ1oeh7DL6wao+v9nM/rGk+fKqHl9rZv98RVl6DNuWznHaCsUzRFbOaHwigO+p6OvdAP4jyukzKbn7K939ze6+cPc3oWzB+qKTZOHV7v4yL2cL/jL6gxq4znehbLm7HOXUoo9ZObp+8vQed//v7v7Oiiz/CGXbWdhOnwLgue7+itqO93t/VuMDAPwhgB/wcs5flPmj7n6yp65cgbJN9FMoW9eeirKX+P71/otQ1OpPQ9mm+f1m9nVaiJW95s9G2SK5pb8ltBWKZ44uQzm3UI8SuzxPDtTN/H9Y1e1PoaDJdUdkTdGH6PstAM6b8lz77TiyrfL5xWb2moq8rkc5BYmPYXvnCr7+KcrpMy85eFPW0q0o+59/2N2PV0H9hyi2Ubj7W9z9A1Wt/lOU8w+/Rto0dQzbls5x2grFM0cfQ5moepTY++v3bBP6r6Kouvdw9zuh2B0nj8g6HeTjkW16ZNlRlMN2fwLA3arq/zLi831Y/W6RZ6D0za8avfjpFOlNybVVm/wd1K+2+hi2LZ3jtBWKZ4iq2voilNcvXGzlWKjvwPIosQ8DuMLKS5WCLgbwCXe/rdrmvv5082nrj2z7MIB72fK8vyMo76L5KIA9M/tiVERW6RcBPNnMHlOdSpfXMoNOAPhalNNknm8HPEfQyrFg56GcmhNH1wXyfRXKe3m+p6b7QpQj166ueb/SypFeVvv121Bf9GXrj2Hb0jlOW6F4ZulpKM6Td6Ecr/arKEdEASX85VoAHzKzj9Vr/wLAD5lZhIaciXdgXI8iBN9sZjehvJ7gN1GO/QKWZyp+3Mz+0t1vRBEqL0I5burrUdAtgOLNRXW+oNj4/gg9WoaX08u/GuXoredW4fm9Zva7K/j8PhQ1+ekoBwTfWq+heo2/EkWN/xSA/wzgSWTLfCLKqy9uRHEC/Zi7xxmVXwXg81EEOccx6nt/tnSO0vbosC1taUtbItoixS1taUtbIjptQtHMLjWzl5jZ22pg7yNOV11b2tKWtrQpOp3b/J6J8m7dr6nOgwvWZdjSlra0pTuaTgtSNLM7ocS2/SJQDOnufv3qXFva0pb+NpGZPdfMPmJmf7Wh8r7JzN5eP99E1x9adzW9w8yeZZa/9rilPx2OFjN7CMqG/Leg7Jy4BuWdwDdTmqsAXAUAF1544UPvd7/74e2vf3dfkPI21ZZIF/c5n1n5Lfciha0qN9Lz/aQsrocD3iT4bbqOVXVymevKyHjNKCsn8nJfRt0T6Ye2an7lYd3zOwh/WremrencfXmf0qf5sjbT7674qXvc3qnxktSV1sv36R6Udy17qs567Qb/xMdqQP5J0+P+wYX+8U+sPxvkmjcdu9rdhz3uTGb2KJSXXj3f3de+upjyvRLAN9ddYXHtLgBeB+BKlMd0DYCHuvsnzezPUSIkXosSQ/ssd5+MbDhd6vMOgM8D8DR3f62ZPRMldOLfRoK6pes5AHDllVf66173Ojzu4m+Gzefw4/X98rMKZBeL8n1Bb0+cz5eDKu4D8P395SSeERDWMmazknaexAtTGiwWpa7CNLBYlAkXE8+spEsmUFeW8qvtC4o6te36fapfuAyijteMtAxtf237UP+aele2Q/8D4zMxA/b3+/K1noyneL7cZuV/Xdv5mdbnjoXDdneWZWXPPcbf/v76/tYxMPHsAPT9QuOwlRG8ZO0B8Hu3/sp70hu3gz7+iX38+dXro5Pmd3/72p1X7v4qM7sXX7PyuuBno2zBvAXAt8q20Cl6HIBXuPsnajmvAPD4KkAvcffX1OvPB/AEAJNC8XQ5Wq5DeV3la+vvl6AIydVUB/IwkObz8UHzw9/dbQOxDRIebBODpENAZuNg3N8vA1IGfQzONkgXi+XEjTK0Xr6fTQTmWe/rhOGya3kd4mchQ9c9hPr+fl+GCsvgJQRBTe8i1Ls6I/98vlxEmKhNKa+Vr/ZMgj/3JT9mpezgdz7vhYg8v+75LBbwfWl7fKIs5bUVtBxbNqdntAotMi0WPXJLBGKqsdV2W5a38hPPtCtfxxGPrVMkB7A4wB+Ay6ycdxmfqw5YxXNQwNRDAfxrlNfTHoQuR9k9FXRdvXZ5/a7XJ+m0IEV3/5CVo64+x93/GuX9s29ZmzEmgVIVTlaFY1ONeNAHTa28MYj5WtB83oSW1/9lAsyX94SHqKOt1tqODM3whM9QFfOsCGdCsEea4KGhinrPaj5Gz6ZlKS8kVNvEUx6BEWWzgKHFraE/5iXpgzWmnmWZIRz3RY1jxBTPi/qlCTROy2VzO6bQMQvD4EeFHKeN/LqAyNhoZghFt4gmifZT29+lXVXHhsjhOOEHOlrzY+4+dTxdSvXQkS8A8GIaC0frvSejvKMcAO4D4GVmdhzlmLuvuj31rKPT6X1+GoAXVM/zu1B2NRycWHWLiVwHWAishmKm8gO9EM2QDdCQEGazYXIwemmTGSR8YlAy6qh5lY90xWaVh9Pu7jYEqoI3VYVDiAnqGNKq4HUHTpwYkKWR0BqEVab68sIjyG1KCKaLxtQzEmSsC0WndsdLzYnv7jll9SsFHyEco+wQhLogZ2aNbNGGLFxM0ddcZ4zH+Ry+t7dsS9QZZbPmRCifx+ymqCLB00EzANe7ux4IDHf/JQC/BOQ2RZQzBB5Nv68A8Mp6/Qq5/n6soNOlPsPd3+DuV7r757r7E9z9kwfI1P/OVEoAdt5R4MjuMp2qD3EtUdWcB88Epba3yksb0Ds7A1Lq1MvgXexAVElvZ1SqqDnKbGXPZqXeTMUWIYHFovHo+/vNHND6QIVxLDgsQFSQ8/PI+jHMDaLSu3tD4ZNqXdzXhYbKaG0LtZr4HlAnq9jRR2QaWKlShpkh8mmdwtdgkuD2SF80LUTV4vm8LIZqey6V9O3iPlFzBoGHwTRziuRw7Pv6z0mV7X4DgL8xs68FgLo3PT3iLqGrAfyjuqf9zij776929w8CuMHMHl69zk9C3ec+RadNKJ4Uqa0nJqjYv+yCC4D9OlDVjmPWJlWHFgK17Oz0NkpFeMwLRKULIcPqUCf05FpMJlXleZLpPTbSo6ioLKicy1OkFTaoHVIAOnvYfBR6Uk5qJ8woQ6oi9GdHjy6R/ny+XEQ0r6JXbZMsbt1ixKYGUtcbP1kfTwlEWYCN1PCBeJGqfdoWHG5THY+dgIrxp4JU7YO8iEQZrXrvtJgyLhy+v8j7doO0gK/9HITM7IUA/gzlZPPrzOwpKEfJPcXKqe7XouxhX0vVwfLvAPxF/fxQOF1QzhD4BZT97u/ECicLcAjf0dLsduUH/PjxXgV0x+Ljn8g9i0A3OeO+ZZOAfytiy1Qh9TKTMF7yYoOACdteNxm5DDYT1PIm+wWCFNTuRTw17zj3adJHSp03nfmJOtWeRoKoQ2shwKOdLJzUOaFIjuuJZKqKZmaHen9Q9zP7Y1BmAuD2AqMNNuvHGdlqtX/ZFkr91o3hrG+4T8Lm7b6sZ+FwFFtpUbur/XfKxr0BcgD7BxR6a8tyHw72rbQylMfdHz1x/blYHrDC11+H5dF3a+lwIUVFE6wGVOqQEhGv0KxqDgMiJmpQoL4QpJm6m3kZqRzOo2pShwRCZWJBwXxFOSFoTuzlfDFvKwZ8IB3jclV1ZYGk6lrwqZMsVP/IF+pgEKO0KI4RF5cZdahpIZ6J2ACbuqiIeRWPxEPSSf3vDIVHm8SZNdgTVU1VIcnqOGsdnCZDyh27tkSn81kJDyIejMfmhgSh0qaQ4mGlQ4MUH3fBN04jmFA5aGUNh0IMoc74HM4JYEBgQ1AvC41EwLg7jNWYEF6CEhoPdH3wTGs52sbB5kZtyuxPiipYbZyym6kw5vozgZ3ZeaPccNDoIqOomvtB00cd3J/aD2zfm+JZ1WQWMhRR0JUbeWXcqbe8pWe7Huj5rkLgzDe3h8bu0N/83BjlZuaFoAxla3s3QA7gxEnaDM8WOjRCkSe57y/68AnU0JAsXjHykpHednb6QRFCMpwjdQB1Kh+ruJWX1MvLAjFb3UnADg6Lep3zDoKThEqHYXiyMOpQIn6bqsV1KxJmAaKTekrF4z5N1NehzdyHqpZzvvq/M6FUPozL1jzKm7YzE0Zqj6X/bdHMTDMI9st46haiCdtsFmdorDUkDq9hfKlg47AnXSB5gdC2niI5fGPq82GlwyMUK/Jo8WS6i4TV3mzixMAQYRhlD4IH6ARku05q32T8IQUzGwtzRpUawsF1k2DsamBEopOWDerMG4XutLyV0qDfzPHAfUjpoy7fX0zH0anA0T6LPtA6JlB545uFAi96jJIZKWblr6KZ7GhSlEX91NqVmTyShSS1a3KbZWHs1GQz+N7eWF+kYySfaCuZJjD5XE6GHNg/t2XiIbIp8goqAnHYITGjEIwMNalQqdd8b29pz+IVOAYM2xDVlhVqM6kkwUOq2sXE4oGr6h2A2V3vgtld7zLaz26PoNZQGymjeWNZJWM1V9XNmmapbi7rZNSocZupsGRioR/1U58M6m2Q8q+oqNbf5Yv+lz7haIRh25wsvEObswUne2b8O/qa7YeaJsZkFWjsaOzSxRiZsvVGHuGpxfRugBzA4gCfs5kOD1JkRKWrtPuIkHiAxaA9cWJ5j5EaC1xWVVjghGG+Cs8ORYYA0IFV04cK2al4K9SzxrM7/FM3LO8xwp2aQKqW17SNPxbAtY3NW6mqKweG7+6OCHiqTQk66bzWKlgC6WWqNZXb2eim+izKFcEx5JsSVmy20HT0uzPhZGq0qMpGi8nQz2rjm4pPnTIpJKi8owwxcxlVMG6GDPtn9t1pZ5wOj1BUm1OlsL1MoohqSB+i/SmEIa53e2spjap+lgze2D3TCSRWvRJ1bYiTJOTWjPVx+AWV09UVvLGwylRDnhDZfRWaLDDDJBBJo45o31RIC5XVhYpE+av2EwfPJBCxcPgsQcPaRlYjM0Gl13gRJEed2g1ZcHQcZA4a6pPOPhgLM7dbzQGYUGlVHdYFLuqOMrPFM1HlN2tTBE74ViieGeLBkBnm1WYiKzWAVGB0W700JIbQZSp0M3Vpd7cIWBZmnG+xGIWv2O0YoXS2UUa0rFby9czpQROk6wvmLTMpxD1FcBmSy4QaTdpVHvypycttt8UCmCcOhWifCP0QKt22xEhbBVZ7DuyUYMHKfPHzizIjCF4RmDo4dExp//DvOr6NFxlGgGzGyMYUzw+uW3mkMWLnHS3hXRsgB855pHgobIpxjmK3/SxsJhN2Hqa2W2ICrYWwaLa/GEg8uHlSq8CJ+3Ui+d7eOPjZBlXLmp1/3hB72O3CYaLJ2Nn6Ii8Lwiw/qdYr97tyuUqKfPl5RDu5zei9rC0vq4/ufbuUhwzJ6TNg1VN4590k7X595kbPjDL0vGq7qMw2RqYQXfAVpoHMbMDX457ukuE+4DbzwsgaRyDSxaKMxZq/xYKKucD3F6PQPAVauK39nM10OJBiZlBmDyOjSLYvMvKL+xlVJGM7O001LQXkq2ynCmr9tTwAq4/IAuDHSZVKbFwaBGwsTJI2MOLrHBuZh1vzriIWZNEu7hsVhlT+gMQS21ZDkeGs2pdDLjJ0rsg16Wvug8Z9CCh1Vkz0SceHovCpdi8kznBWPdlcr9o3BXGvstV25hLmWdThbgwsxOPNZW5YfT7XkeLhEIqVWhiCqgY02YaQEJ2MsXIDy8kSed1hR3bhx9GOA4t03YDNbDhBPJDZKyqrfRckrhM8eMvsPxPU2RkBZHa3zua4Skj2nd73YYYyuX8zHhWdTKnpysdUmyMd2zJpPAzoVBcwVpO1L3gBmc+BaovOBG8TOhoUP5/DeKzNZr3gJ7469T7rx2TxbSYFGcvaR6ZlZOk2TA7D/uFQME8bHQ6hmAk5GoB8eEK3AT6yx4DOUKOuvqQWDjsSWI0Lvlg4M7HtRzCDslwAACAASURBVAN4ZxQyxINUnTRTA5gnQ7s0qn0NnVTe0y2KtT5FZukJ5ar+yWLB/TiUsUzQCa0MvUzGzYUQWTgws7w9VK7zggmMAeW6oEEWF7JLRmRCN67M+jGiaNV9OFG9c+zJd6s8AsjPa0wWiZUxhjp+Ms0hGUunSme7eryODodQhHhkeacJD0YSdOlAqQKsUz8CWdUywksNoPckR/mrVHXmgdUzVh8rfykqYoSgdrbEFrQSsZl1dtRU5Wc1VwSMHpY71BN2KkGo3AfDawOoXSwUOkG1oDMDiZ8o1wA45HAJpvkcVid+8yCr8yTKTFRum89727JoE8vm9OOms7Emtr8sDrBD7RQaxv00iZ7reBmiF6ZU8ig7c+BsLBynIMXjvrnyDiMdDqEoK7bxRCcU0IhXWSZZFTskFMggtgDqCh35GfGpA4cnTjbAM2HI10Ut53jIbLtfQzDz+TLUI+xy7HwACXjuI1Xvqa2K8qIu39/vbGWq4ne/E5MAL1iDqs9IKkEvwdPK9+aw460w3ps+GN1xPuqPNDg6+JpCp5mtkZ5put2REKvGvQ59KN9XxuVynmxHDvdLpNmQo8UBLLbq8+knB9rAGl4HACztgVNe15omPMVGA6KpLSxE4l6dhGqwXjImtp9VK+7UYJ+ySdY2tkkV98nm1CaSTvApO1MmfLOFI+kjAPATe33AsuZVM0MQC5mk/kydDjXZ5rNBmHb1cqwkmzIm2tohNvLyWrYlj9vJ/9XOSjGvQxiQLqzRx9xnWZoM+QkNAnFqDLI5KBY2zrshgRi0dbScAWqeOOQ2KAA9YuOHLOruYJfi8mT1ZMHAA3cYVFm9Osl4gE/t4OBJoXafSBfqOBOrobowrEIe3C5SyeJ3lLc8ro0QsaphmRCKcjmN2EpVkAQVFdo6m2B6AIdO6ANMcF7oOL50IOY9s2tzfYSSh+D+DAFi4nlFfVPCOdowFdmgdky1jwd/2bjYALkb9v3cRoqHq3WJ86EZ0xP1r+XZ3YUdKZ+lvauP7Rr2Jyt6oe1ZMaiGPFJOxwOlYftgix0LG50Ky+CFBy3bnyidJf3DfXLg72FLY0FAKCvsrul+2akJ3TolEX4zec0BI6foL0o7IK0oWu1ynEbVfBHEtrMDO3qkL5MXJXpGWd2ZwBvshEI6rlpdNX+2Z757BYKo6639nCbay2NDg/wTh9Op0AK29nM206FAijBb7r0FOiFjMbDFO9jF6O3vwxdL9arbicACBehthvxbkZHREWSI4mTC8/0q6JpjYiK8o03WUN0jLTtepJ7BS57Zt4ICralAU0QdfUp2R33p1UqvZ4YYWchMeTyjfUne5k2fmMDqiBu834HmEkHlx4/3zg3uZ+4P/q0IktG19o0i6WzBI2TZxUYmaK/VT8851V4yoSyLoe3u9DGzp0DF0XI4xMbpokPTOj9+vEcuGpAcwjAReHyfd550NqYpw3UMRh3MycQczvgLT2bYPBe5V5zVme54LCSChydk2BuzgOCZnHXIaIdtb5yP+0uFU9TFyEyFcaRj/qQvVSB3kQD1/MF08hu9zGmK71p2dwAH8yWCqe1K4f+rzAJKMY7ot2Xjg8tfpbZm2x3VFq5tonEzeKIzdT0ZK37bsbFtJ0lbR8uZoilbGTAOCDJ6AxIyEWXNZ8V5s1iMtplkhW3CjR0v2eCXSdfFtJFDwKpnubtey+rQGAcnywRuiFiQEwvCoc900sS1DNmpPZUm2RBfJ5O9Oa9iwvKCw0Id6HeaUB9wn3Tq5T69XztTs2M8aLs1IiEEdPDEXn62TzMa5D7K0Ca3c4XaPOm4U4E4ZXYIPkKoKc+av7bROL+2YYPq8/42TvEMkaocQK/6Cspqg4TU1hisfvzEcqJoHUw0wbqXta+K4eNyEmdL895mam2kKY3onSqMaFngCc+ph5ZRYuThxSSb6CwIGVWwYyR4VXsd51H0lDm2Au1mTon4Te3t0CovVEDbhaJB1m2fMaNmrkc1Ba47o+CXvN4DktO+YX4m7nUqf6Yy1z5gxNwFj2e2z8gzNVYnTAonQ9sdLXcEsarFhmJILJ6GRgTF/UwFE2HbqYezWRr/mO7R5ZWcyw+BXNNqzN8gCIJPRjCMjiLt1L7fJN5OUZQGsk95RFuM4JRarwJEvre2ZkHsOiFDiOnupaTshoqXDVraFlnoswCbzcop7rIlVN842O0yIWSmfTogb76XqdE1T7ajRTWXbpOBepW5jkxVzoRpJmh1sTpFWpzj3ufDIxRpAtl83qtRmWEevbdS49CaYAV6VTSI7EWduhz3KnUqCQvDCSTX0lY1OmxluotiQHcLLyExPPmAPj6PA5AXcpyU2rEEwQU6HdBOqLmJ6twJ5ExocfmsUqrp4SACRfoudX5EuWFb1HyCirq+CR7ZiafxnvxMgoc69gZkmqnbK1TUVSr3EIY0VRaZJgbzzlSeDSHExgKwRYpnlHRwME3ZUyZUNx5ojDJ48Aw2M7LntW2B6kmO/JkqFfdFkHRea3UQtYVgwnapaiELbDYhaD2JkBx4jv4iVXPYVbO3Byx6GyLX1UwWYkfrECvzwgJPn10InsxzzchcBZfWF+kRzaNFb5VQzuITWcCqV38KOcfvbHukmDi6vNnzYRTLzyvaxekVJSbP61TJYTix3eZ3BiiEkwy47hQQVoHUuQIxggfx5KqoTVHVENxt/TYzvq/e2E6FVSHNNsLYmhcDe5VtSNEWl01oOM7ISwOVmTKbkpan16nPyjZD4S3pG7WjtaPaVJ3OduMcxD6WIR5xug0LgfA26UALWoeKE7ttJ6ykbe3ZyOI1tFM97vX/YB4SoZyhT46x7fabb2j/szu2wdunQmY2N7PXm9nvrE27s9MGlycTx/f3+8DWer85BGICx6rKAy4Gna7yag/KJpWqX8CyzqmBVnktnzK422GgU3wwLzw5WXgSCuNTZAY0krUv+lGRarQ9+8/9QpO0tb3WMaAQRcgirLkc/t79l2fZgvi5jpqmWxC4bu6DmbzAKfo02hZ7qtcEZLd6eYGNvCRAu0V6SgNKvOYaJ9r6mRdVSTe0OaH05WonResDt7fB26vp2wG8FcAla1OG3czy4GsAiK17TKuOlxoOTmB1RG0w2cRehSrUniS0VMlpsK8qc91EDP4DQbM6rWpo4jDq+FL+2bjPqrsIquaJlrYPjo+oV9S9zhYYdldC/U0dZNseoa4ORWdojj3Y2qesBgt1iJ/7Th1Z60wnxE8X9xpjUJ95gjwnnXpAdwp4Gg/LGs4UOj1FcmyR4kmTmV0B4EsB/MLaxKRq8tFe3SCfz2G7O/0rUIHVD57RJdusapndf1TBFXlo5e629rHNJupmZKfl1mtDbFz8F9TbUCUjXSUWZjyJMjUsyduhjAzFZSp8pBXEwvzoVjPTRYlQrytiBgk+RY9A2fGU2QCZ10Bsq9oPWqAYacez5DKjPBWu2nZdHDPhnSE1dvxkwlD7IFks2nVaGFjTalrNzubwzz5maz9nM51OpPjTAL4bwMUHSs0r+WyW255IsKWHB/Dq3LJYv1JHGQkKyIzinQCcz2G8IyNsN9lOhfm8f+WqtmUVimMeot0hRNRrTZO7Q3zRR9mWR/WektAa+lT5JRW+ix/UybxP7yGJ51nri0mbvoy+75Al79qXzKf0ZXr4bVesLXlhhxGX4Z4/jxCc3H8qwFRoZXwp4sycT/ydx2yki0VHd91Q0LpxH26AHGf/O1jW0WkRimb2ZQA+4u7XmNmjJ9JcBeAqADjPLixIoBqEU8O5DKKpQd8NPDZMs5BlWx1EwFZB2CYs8xCrLiGBFljO4THzeXntQdSr7aDf3RZDThPtz5ACp4lrLITrfWPUxWo12+1C2Ko5YJacrE0CsbGkW/NA9q5sVw4vWIySdLFTlKltDlLBRM+Ey0tj91q3WnFcsZ0WI00ehMyULSphkw5e41nJWBh+B/9xT5+9qvnLTDna3gA5gBPn+N7n04VzvxDAV5jZuwH8GoB/aGa/wgnc/TnufqW7X3kER5va5Ht7bXJriEnzZFY1LQtB6WyM4uENVVhPWzFVCTM0kgizbmJFGRVJ+K23peqvi3DuJiqXwXwEL/ydVU22lWYCNdKpCidmgcG+FWFJLFjDlKF22oU4F7hfKg9ehbft7HRpmskk+oM1BXlWrf+Y38gzIQA6gcgqfmgCs9nyVCU2x4hppgvzYtU3KMawqvU8hniMZYuFPL/BZs6fZQOXzzeZE6kJ5qTJsH+Az9lMp0Uouvv3uPsV7n4vAE8E8Afu/g1rMvW2PNreVTiddWmBWOH32/FcqWdVEEGmJgHo1aLw+rG3NXbYxCTSCaQDj21clYdUJeNJymp+pl6zmlx5jT3b7CXuVGqanM6LhNjJTIRFW0BU7eLJy4iTBRr1N6dvqnzkrYJFw3nCa9+hUhZGvLAkcY9OC0QT4KwhBHE/MVH/d/WrUJOFsbM/Twm9eAZiwki1I6VVvGo+eoYD4j4FcgALn639nM10qHAw21uaqkH2m25LnqgkU6oeMCGMok5VkUuG8p/q6kIaFv0Wri6P1h9lsNrNZZdC+omeCEY7cgR2wfnw6z+1TAdMtqvZoMibqrtfhp081N50p4z2i7ZdBaKq+JF2TXhIQWHjIpPakRPb8yBkyUQxBOMzJRpC82KrKUMRYKQdGzP5Iqk2hkALSCxO9b3OHW/uXR+kLx/TORS0IZsicO6fvH3aRbq7v9Ldv+zAGXSgZZ7dZFXsYhQ5jmzWx9R1KitPjswOxPd4Us8kLlI9yxnSiHuskiqaVHWZ+sKPH8ciBGLwF8J5incSYllcW0Pi0ucDCmQUmfVV8BufKVVWzBmNV05PqDr6OUWrasroqknCeXQ8RRlq8w1TgxKjeY2AUN44jzrhSIir8G7ls51QxtPUxoIox/f2xkVHF/VTIHfbGFI0s39lZtea2V+Z2QvN7LyD8mFmjzezvzazd5jZ0+n6U+s1N7PLTqKJh8x3HuhDDOndA52awHWl9QjrYVU6s8UkEzGd+FleRqtAP4nC1qSTMO7VskOgtgnME4QM5WwSCEExBFHH5M6IhEvav5FM8zMKpP7tysiQoAp27p9M+GhAfmgAVF9M/uHg3tgHzvbJzKamlF1XFKhjJXs+ybhITRTJuGu8h2NxlgTBK78Tmgar7elurw2ixOJoma/9rCMzuxzAtwG40t0fhLJv6okH4cHM5gCeDeCLATwAwNeZ2QPq7T8B8FgA77ndjat0aIRid9S6IIrhQfPETNTXDgkycmFkx+UtMy7zsFqSCc0oX8tJHBZNsMeHkECqftOkWhekW7KK51cFep+4a0vztKutlV8Sz+iJy2OUG2nCAx6CUupvtl+mCXvYcE8Qb/Deqc7VhtjVwwsYL1r8XHlhCuLFIPJmCDgzG0wtuAm6a0Kfx0IIygT1ARi3JC6mbZKZLf3kqbyjZd3ngLQD4Hwz2wFwAYAPmNlDzeyPzOwaM7vazO6e5HsYgHe4+7vc/TiKM/crAcDdX+/u7z6VFh4OoRhhDpkaBPRCRkk9o3TkWFMlQshJqEjbescUE52FCwtL5kNRXIYsWGXnCZJ5fnXS6bbFILFHDipVImS6/OzQySgTpNxODvmIPgvBMiPnT7SVVP1hgmZtFPTXaDbrPeLxO54vPSOL4HEVSswbX89MB7wYZJ7dzKwwYWZowk2ciUoBDuL/Kptxl2aWmBmAfD/3KVBxtNjaD4DLzOx19LmqK8f9/QB+AsB7AXwQwKcA/CGAnwHwNe7+UADPBfDvEzYuB/A++n1dvbYROjSOluH9yjVQGqiGcs2gRn+xI2E2609GUcEIWtXjehYwzvWpcBajduOTvcFBKtS47cEHC6sViCJVW+M6CZ92jdPFZGYVfIWaGbt8muBVG1zStjQAXBaCxgvzJ46yAcFWJ8XgONO87KWXBSR9fw0jtNo/zQlCi0H3RHjMZUSLRKmCTA6qAUX76Vm0eEhqaxZPu2xCIjyn7J6nSAfcsfIxd79y6qaZ3RkF3d0bwPUAXgzg3wB4EIBX1PbMUQTmGaVDIxRTFBYPPwKkswEStJC9szx5ePJR+mHiTQkuVr+YWHWsk169mxzInHrHma9MMCsfmbotNJwcLhO0e5+NeNUbsq2/tazBc6ttinImzBqpBzlIbLNDOhIcQ9A7awwsJLXtwR/zrAuePKt0W2Rmw83q4HJJSK/yDPPYZtPRJGrUscJ9orudTpE2uKPlsQD+xt0/CgBm9hsAngzgWnd/BCc0s3sA+O368+cBvBHAPSjJFQDevwmmgMMkFHmQxaDmgT6n9zTP5SU/uvLKqgtgKfxY5eGJN2V8D5vaOvUjqwvJCp7VpQOW0nShP5ka3L7Ke0FU6Or1rN7ajrSsVWkKg8t+YNQrwtlUWER79TubL4JiTNTn396rM5+PjoasP7NnOIX0umon+nWxmH4+KqjENgos1edOY4kFKWsD2WnTxYXr1AVwg7ShF1e9F8DDzewCALcCeAyA3wLwNDN7hLv/mZntAvi77n4tgIdExmqDvK+Z3RtFGD4RwNdvgikAh8SmmHkM2UgeCEpVIbJVAWKjiQGjgcWKGHUQs5OFr2eUTEB9UdUQ/BvfaWI0J9OinJPIQdODbSyEdPBa/2uwedeGuK72pcwpMJ+Pji0x7g8e3sViPJyB+qDbqcP9ORXaUm1qfmKv9VFrR9iKK2rl49M6e6rY5rrA6uAreONrfE9sgoPNbpaEZkW/qvahZTPPLMBUePJvup8djNLGEdu/uawNkDtwYjFb+1lfjr8WwEsA/CWAN6PIop8H8DUAfszM3gjgDQC+IMm7B+CpAK5GOYXrRVVwwsy+zcyuQ0GPbzKz9QfSCB0KpOioSCfU5P19YFHeypdupRJiVaNbfSlweRh4altiZKpOlvgeqmSmTjEftPJ3sWiJKqNxeDafjQKEbE5d+5lfsRuqTSy1gamamdn2gti2lqFf7g9VR/kZcH0TSLX1IQVwN7NE2Jozm6UgOJvPOxtdOLQ6M4u2YQJRds8pzA68yLLpg2gw4Uh97TBY3c7H+cmsMbSZeE7NCQfRcm4HFfV5M+W5+w8A+AG5/AYAjzpA3pcBeFly/VkAnnUqfB0KodhWW2D5cOc2TmQWChOD2HTSiarS0sUXVveA8QSUEC47O/2xZpGWeSkMLL+vEmRUh6kwygJ+uZyYlCzYVcDQ2ZPsJBjUPeafy+IFhfuNVcngI5uMMhEbD2qH5by64BGiGvoi0rNQYoEYeaiPmoAMJMWCLdozZXZQPpjfEHDq4AobrDrEKL+afJQGJ5g4hWw+h0+o9oMWsCE613e0HAqhmJIKE6BDAACmvYnAcHRXOA5a2ZE/rotq263QobIRHx2fGcrMJpY6A1SQ6HetK+qLvJmA4PZW4v5i4T5sVQxSNVP7jXnlaxkCF1vrUF8mfPXZT6G4CZXctD+zsnTh4YUxEyDMQ+ZsinGpSDl4oHTD82dBt0KADuh6f395EpMIvzZ+68K1qQBuB7ZHh50RYqgfpAKwXmuhDfyQdf+yCgy+znaaWf++6A5VsMCk/6mBO1NVdPIl6JPbwUeVdXt1a33p6zmjL3hXyMRigt3dJT+I6gXVZqqW2rhk8RhQSfYstb6s77ReNTWoGWLKxKCkbWIkmQkw5mVCuJdb8npVrUv5mjJJZPm5LasQHi+Q4u3PXn+7Gdqc+nxY6XAIxRiggjJM08T9EDQx2MhOGAKFbZTDxKXB24WJZBNEWR2M24ulHVAFBaUDkKuJtU0NZTDKW/RhGF2b4nemojPRhDEzOEQFBvoQHZ5EmVmAFhs+UELfsTz0W0QQ7O7Aj8s7mTMkp4JknZBgoT6lEcQ1DUPK2jtFddx12gqQL8KRPhsLLNhrus7EkJGaBphnbb+M+clY1JOgs/0dLOvocIh89+VHvb+cBugf+GIxhPKwsyO8l4r2VG1KB1l8p5W/29cKwM47Ctvd6cuu39ODGjRdtGuxGJEOCb6g9DCApE2tblWpyBPL5ft+LxyHnUQTtjQ+ELgLKwkKHmIBms+B3SNDW7KyW9+oiSOeOZsjWNBIXwwqZ3V+RVC6q7CN58saRZQti1o7loufteaNvfi0JbLtCpJ628KY8UB9vuweGo8BGmqZ3Uve+N4pkjtwYjFf+zmb6XAgRWCM2QL6wc5qaFxTVMPqy2JRBFYVNsOrURU9lspztTxR7coAnnXCxYgPRn6s+nKoTce7bosD+gmb7CSZEozDIRQkQMpEcmC3HPaK48fhM4x9zuTe8Q4kDgDKmwXWAyjCYW+v90Bn6Kaifjv/PGBvD37bsWU5apZgtJmp/qIKq13TVOgyX9oH4tBZuZ94NmvjVtMN41AdVonGke2wGfijcdAJz729zb3idPs6gjNE7H2empRBigrYDqP5WU2eCuplh4zaJiXuUc9e9L29fOALP2xLs6pywyYEJPMcv6dsXtTeLswkUz0DaVCoSzhcVgroIHZgED/Z2xSbSYLbw95emcCNfxFifuz4Mk2CxDvBrKp+zdO9QoDtgBNCXcvt8u/t9Yt25uTLFhdF0LV9Ld1MbMdJ7GS2q6vzsLP6HjzWcd1FGmyAznX1+XAIRbaH6IDIVF7J0wYSDUjf75FiZ1RfZZcC2e6CFM14fW1CpvIykgm+eTsd2QzTA24z5ACM73FpLMshC8EP2VvT9mSTJCYb28tCSLLDY3d3uZBMCGzTe5w/EbrDHuewHWtbCbGmTq8ohwVHNo7m8y6ioAWBr0C7A6rkcUXpBrSnzyb4I76Gd8Pw8yKBt4y3dGBeF594xhVhn07aep/PELWhwiuhrsaq6gLdQOkEV1Vvm81L3wRX6+r+k+pkITSyFZ6R6iqnBOflk1FosjVVTnkQW1ITaoy2Kn+DR1fUrpWkPFL/YTaDnVfVVxZO+/uDsEr5ZuL2ufd9x22K32HuULUzHFGiyg/9FOVwG2djWFaq/k6pxtnYAVLUmwpGpUyQ6rObMA8UjYCEqWoL2eK/SaS49T6ffuLdJwPC0IMdgKVqFCqHrurqMKA8HU3ZoHRwgwRYtoVK+OnCScR+1u12UaQxNXAX3iZB2h4RZt01nfzufX9yev6/vw8cO9Yb7HVRCkS7yraWkTrNEuGRmiPYpikIsRsDulDENfbUa9/peKrCpnMgZU6Zyn/bmcLfsUS3Q1uygz/02XFfMSXak+6M6mhFVMDtJXfD3lYoniEi9OH7Cxiv8BOrrQ5gftnVYEdRJCFCsptUKgAUnWToADSReXBmXnNW5bMYS0aj+/vLkB81Bahg1EkxEas2oKDMZrZYwI8f79EK9SOA3rKk/HBbBAV21zLhkS1W2jfES+yIsiO7JeKAUWjwNrngLKifa3+tO0SBBGSMuSmnS2orZ+E6ZcZQNTh5loqmTe22UW+m8ZwCbdXnM0WMbPi6rtRiJ7MQMiw8CSnyJJ/0NJrl+2mFv8kVn9FOlC9t6urjSbEKzU7tWNF+W3VfBaWiBhYs6qzKVDokQhU0QautrluUMmEnfTmcWzglGLumLW1rBe0lITqK7nj8RF+EQOf2aX0TaD6LzewW5gkbakatD+s70Fu9mXDX8YhkYec6N3RaztameKZIBm7EkDVVlCdMqIj1dBnjQVdJjebLanzc8zxlD4rvcU9sO90qTUJicGpw3uAz4bnjjctRflSwZIKGhZlsd8xsWSnCZSI0lfIrQr5zIDHCpHCYQd2NhUsRkiLNoBntGadnMuwTz9qhC6+YSaaC0FeGjcn9YQysMt3Qoh2Ok77iiYWTvnfOJ87LWsWG6FwXiofDOBArc/10g4rtSLyH1+TYpdlsucKqPWqF6uRO74teLLrv3dFWstp3ajgJhQ5lkK20OVsWizZIOcC7oYqKKjNDv1cHRByd1fHK/RC8HNnF7JKLazdWxEF93ZkDan90/7k8jROVdGm7k3QAlsHiU8TlsICfEhaZAIrnQuOq5dHxwEHxbJeNe7XfhjHHFCo8j5nWXtl7POW84/Q8XiStBt/3bIg9lW3CG6CIUzzA6wjOWjocQpFXMlZF438MVBV2UyqW+zB4gETty4QPsBRW7iOy4v8scGmlbvXywBZbWENFrF6rYyaES71mOzsjkoq8bGaodfiJPfjNNy/7idvDbWd1M8pkXrmvg5eMj7gf9QgaCqFi8z64OD2LMRaPbDJHv/MiM3Vfia9lqDSxRadRBrz4cnq+VsdAt8jHuIh00ufGwld5TgRkd57kVKxvZqM8BVrA1n7OZjoc6nNdiQdPXXxnlJaont2qSHYkVdOGAOGalx9hOFR8f7/s+JjYF9zFRlKZ3UodE2AdMU8sCBgtx+JAdUSwb1dv1B3q4G3HRvWbF5uwrylaZIG/wkHS7kcedgSpCs3fFdFkZowVJo7UMca7XLi+VSgpixUN4npVyKuQSeyG3RZINbGwQ0TUYNvZWZaZtN+YNy6T2xz/N+hgKfUAewc4RPZspsMhFIFii5IHqEdbsXe5XVe0KKpKp3JzkGumAtWBNtjYeMJRmm63gwrJyJcJBq6fA7WDWEAqoqhtzuL4BgEQ7UVBmSnimddDW8mj2uqhcjjcp+t/tfmxc2gVYlPHFgtDVWETm+LKbWtqO+UdH9zOLJyopusiGbI6xIOcRTx0tlUekxPqc6cK82IqC29nVqj3I7Z28OJHuil78UnQ2a4er6PDIRTZNtSpvlUlbWOBENIBHnKmGmfb0tIJWO+H0wdAt7e5K4vUGFWpBsfLlJGe07FQZuRFE5qdSU3dZxvmOoRAaLRDilgime4EIEKSiobTF7HzRFaBy5NXJ7wKNEK0LJSHhS3KBHqHVvSj1jNhp2ShwoJueI7ch6AFOhtTTIqEhTrH4dS4VB7o3gAUDqqtHJC2e5/PJMXgJ8GH3eXrBYA+5IBVJw3eHt5mB4xn8LFAasZoCQeZJV7KoCl7Zp04LWQo7hH6G1DeTHbRZEIhQbeBTrtXAIt+BQAAIABJREFUuUoZDb2owMVS+KmnNFLoIoAZ7dFlwRSoSsvifhHeV6WzCy+E33ZbH3OnKL0uVsOJ7dzfiraonm6xIyEfbWpmlEyDWeOFxmzWeeunzDZDu/seGevR8aECvn7X/d6bJj/HheLhMA7EQ89UhsW0Wh3p+AirzrkS6oMIGueVeDax3YvtQCywgtbYqZodKcoIuymHfGS2Pq5/2eCR94yPBXmUlT9GmurYCtWrqvHN+8m8SZndzp1aThojp21K+k2fnd9yC3DiRO+BZUTMPHBfaZum+ol4XbWdbwiMXiyWXmYWiPw9CQBvQjyzwzK6XpCjpyLyLu9EuwZv9KIcB5eakTZA57qj5bQIRTO7h5n9oZm9xcyuNbNvP0Cm1KBt81mU2XnaGtEA0nttUlG6QQjWesMr2tJGzJuigUww8uTfXx6eYPM57LyjS9QQzhvdKpipVFm4SbSzfqZ2UZTs/cRw4quznVKbJsNNojxFTaw2Z3ZaXUzi2UofGj2Drlzmpy5+6WIUQoGFLiH2VmcsnovFsLA2NR/JwhNpoi7ivcVGxnfqYw3F6Ras4EPf+8MLBC3a2V7vYTx3gn85H9Jg9JOkInPP7ZCc06U+7wH4Tnf/SzO7GMA1ZvYKd3/LZI5kNex2HQSFXYkN4RN2wsxL1xmpZSIPe6kzw/YUqlSVdzbDLY/8HNxwjx18+q+9pUyYRDXuDOOxG4RNArMZBg9zlEXe4+w06RDykyoUt4naOjt6tBcME06w7p0kq+xpqlrqMyHBlwZPZ9/D+VP5U7to53AjtbNbPFnlpHHUHeY7JUwSXqdEwSp75BBXWL8bjzVWraP/JBoirofZSefQ5tRow/457n0+La1z9w+6+1/W7zeivJv18hUZ0nfytgGnaiyqN3U2GwdbrPo6QXmlX4hzolJ4j7uVl5wdTfiwWh5CgVXgug93dmyBu/z1sXL8vnvZl0vqIECTUyftRMxZFwdJdsNQ6Yb8s9kyvpH7JFssANgFF8AuvqgPTldVOfgmHkKwdgHlSASyqpok9Lv+iPaSCt3SRTv5dwhosT1aMnZa3ZndLRPY2aK7WIxlS/ntuUT+bKcPL7zxn7UFNmUwT8JPhm479XmFVnF7yd3Wfs5mOu2OFjO7F4D/HcBr5fpVAK4CgPPswlEITKmsTBoQmyEKFlpcvtTVGcMzpKPCMhO8hPAcwNHX/s8+3IZfOTqlqk60dTi1JUufqeFZ+xUBc7jKYh/YT0wEmQeXiQT0YMtSFZqQW+TpitqXgz2irCowui2QGQJSoZ/9pj6ZUsvt6NFyKAYvhCr8dnf7sy/jf4LcU1OMaiTanhk54TCN+jOb4SbtiEEOnPXq8To6rTjYzC4C8FIA/9Ldb+B77v4cd7/S3a88gqORYRw0upLGPZ4wYoccBm5QZ6uc98huhT1tGLzMa0U6gTw7FTdskoQaar+MAiwxoEfdjDQmUY0KfUKf3e+axgMBcxnuBdXedltfj6DBlE/mVxESt7W2o5XH96OvtJ/4E9czO5q2fxUtEqcUl2UGu+Si5Tt4mBKkxsTjoaVnQZktEmGPVHsnUAQvjdOmDVQNoPXH1LjfJPlSeVr1OZvptCFFM9tFEYgvcPffOECG3l4CpKiCaXAeiDoMYPota1FnqMRaJ1PUz+iF0VBV5bvTV5QHLitBbcMZfKGySvsGoznb9DhPXA81MxBN8NAh2H6SuwjKyWfAfU4xpmkaFgRTz2TdbOKyuYwYM9zPLHCVBzVdaL9HdR/7xFi3lhPB3l7CudLnn5CGZbFteTLUJ+oHxhOb9uXcUeFzk6E5Z7t3eR2dFqFo5Wn/IoC3uvtPHigTr5jAaHPh7wdYDYfg4MLYoAKpc6VTbRNVcbBzBlVUmL7GIN65fOJEr4ZRms7gHgM8i7ObJc4IVoO5f1iQqSd+FQV/mT0wysqEmiBX3eUxxGcyAmVBDuTvP+F+0Drdx+B/YDkGWHAn42gQZurEiPRsX41nATEBKJJcgywx619ry9SNY+YpsaUb0D8jrn9DAdy+dbScNH0hgG8E8A/N7A318yUrc0ypPaEaxOSiLWC2uwM7sttNAt/bWwoZVV+CukHah5p0sW8q/HigBUqlkIx04MV1zjuFpqLJHE7C/RBfY28sC5Fabhjch7g4YHnqz+4u7MiR/hSgIPcivDNUyLqRoklgUHWHE2eyvgGW/Sfe49ZXio5JuDsLVqIulpC1AClj2I2TCV1+ZvTc+S2NsZBwPGNaXvCnvAF9CA+jxmWmsS+0T+IZ8bPYsEq9VZ9Pgtz91ZiOUDgYqZBM1DM+7KCzQzEyLAwt1cgoM/Ix2goERGpRx0vkXfShNB3K5DZMIZL5vEN7gTa6shV5KJKKdnH7CHEMuzsonc1ngFW1Pd5JrKianwOja11guPwMFZHgNeY1FjFFj3GPBZnajCNf7BiKdDxmJrz3XR/Oxl0kmnbQHLrbidOLTTIqEJPA7kaLxfjM+b86xRhR87Nm4U5xtptUeM927/I6Ojzb/IDxMAFabf14ed3lYK9hAahLFA8+HiyZnUZVVR7clZrQm/UxkEPwcUz8+E6qeLxxreVPVMthsmUoOmnHZJyllOH7C9h55aX08WKq1Cmi3zOUy2putsAsmVtOYFFxOzWbYx+ZB7XjZvvEuV9Yq2AVmvuiFZaMHTVrZIt0BG1zf+tiymOJ26vhOVMLC7eRF6+gTMDyApFpDKdApRlboXj6KVtxgaUDZLFo6GpIw/aWqQGldko2ytdr7fADLYMGdWfcnjKmK3qL+rg8aQcL2MkJwaSLR2ZIp3IGG+T+fnkp/az2TbwYiwVI9EvY5rhMFpLc3in+Ig8JxLZAUBtSYoGi5gIgR4TMGwu1TDVWJCYmmsG2mfA5mGoyLz2PwyiHnWJAh4iHZ6ZCn9uaOQBPow57rofkHA6hOEUi7AavXExe9cCqKsXqJzCeU2e0xY8HYYYoa/r2W9XTZMAyoijZk/yEPAY+MsTH9c56Q73GU2ZvIPTjxxu/vL2x22GTIaWsP+j+gIKYVAiqQIr2RloWSMz/lBNG+eRxwXY+XhRFpeWdMe1MTaG0jcyvRjPwf14MmHixtOWe587GnS02UnZKqxadk6Cz3Wa4jg6FUFTjchrWwGpmhMCoACqZc7tc3JPyugkp3tsuNCbusz0nESKr1N4W2JwNXpe9xWwf4vbFvVXqVpDa5RjhEcLUAOnOc9vY8+7Q3aZyaxu53czf1AIT/JXKR5upIsUVbR8WHVUtMzsfL0TR/n16Q5+OmZlsuWOUx+UrAgR67zejbHbEKWVqPZb93733hu2z2TjcADkMi633+QwRCaKB6AGboIS27Q3oPbQsEGWSNsHENidNy+qjutQy2xpP1KTOoa3KnxlsZ6e1x7lOQWyDZzn4UQEpNtEBoTH6jk94bRf9c9DdEWksHttIM3WeVENXVZKFBKuy7IXmMjnNYmnbzYRC16dRX1yf2l6qvGV2yFJIL/TqIpohzEHwabkisJudNcYLzwO2Py985JvaOSzWp0h+gM9ByMwuNbOXmNnbzOytZvaIg/JgZg81szeb2TvM7Fk1DBBm9rVWDqFZmNmVt6thlQ6FUORH2XZ+JA+x21cbEyF2QOjgCqGXqH+pwORQllndj7uYeLwcFqSeycw2WvnpDgaItkRZdXJFezqnzYkTWBw71t7KNwx+HvQklJ2RyxRaiD4Su5bNZ0vPMbcpc1hx3B6jNTZvlIuTKrcKuKDJHTzrdDgS/kO4FQm/ycUTKPGlmQNKBHTbn02bCXxvb+RxSo3VMZP9ZicLje0yplaUy892E1QdLRva+/xMAC939/sBeDDKGQkHpf8fwLcCuG/9PL5e/ysAXw3gVbejrI4OhVBMDdNANxA7VXYqrQYos0pD6mKbaBxzxpO5Cdx+InQHMSRIsDuqjCdAzZ8KM7rfBnDi8ezyRh+oQEaC0FRwcX3RB7wziInVVv4di0Zi3wr12/epTUJ6HiHHZQ4IOYullL7R/92hGYz+VLizNiCCqB3+4T79Qq/6v4VW8fjgfk60h5WHUDB/utjzc4tFVhbfVqZ4vTdGG4CKZnYnAI9C2eQBdz/u7teb2Web2cvN7Boz+2Mzu1+S9+4ALnH313jpyOcDeEIt563u/ten0rzDIRSDplDNTMJdNC0PvDph26RMytFdFhzKUqqg0BaOs6M02WRsQlsGcjoB4n9MyECLIaRm44kzjXQ7G09ubmvwnnn3E4HW0GVMcg1KDz5UTVVvulnZM0wIuLWVeRCBoc8mbTuiqIm+AWkPUyaNxaJ7F3ZTdYOnqq1w+12RGreD6u3qifva1yRIu/GhC7qUz/3f9R8/XwYOam7YEB0QKV5mZq+jz1VSzL0BfBTAL5nZ683sF8zsQgDPAfA0d38ogH8N4OcSFi4HcB39vg6rTuG6nXQoHC0AcvUz0BlNuu5Y92wlRR2cc4w2NqIQfDy5w0nQnUhTEvdpmMTW170cqyVZDuLBGRAqNE+CzNak/aNOAK1r2dBesGk/8K6MyF/zBI9de7hP4zvzNoX6p4jVa24ft1+jBRR1T23fy9TVTMgyGhZU2sYJx0WSlxgq0Pj5aD+s4KUzrWj/alvYhqlN2dvrNwcE3Z5nsoIcwGJxIAH7MXdfZdPbAfB5KALwtWb2TAA/DOALALyYeD96CuyeFB0eoTjx0LpdH5kata5MnrjzeRt8w2NVew3QHdLQpWEK9cp93E2ieWb13R/7C2Bu4yAXvlXdjgWis2NODXYyObCHeSrkRnduMOpcdbZj6uFPFqruOSj6IfNGqdBy4Tjr32LYCRIWUEnIDreveWwzFZxRGredKdKLw2/oC1pg2mIYqJTrINPMUE7CY5sLiTed7esudWyEHMBm4hSvA3Cdu8eRgi8B8EMArnf3h3BCM5sDuKb+/G8o9sQrKMkVAN6/CaaAw6Q+80ObWuEStcpDKE0JB95HrJMoK5dIX+jUpdfJzWofkTsdLlvLabZKpkTg6L1Atp33OTEfDAiEVLUuLEhtbrW8bjKxB5RUQd8n25a2IUNIGWpjVMVpp4L0gdEJtyhH+rdj/UkQpe2jfuyEsXrqte0Laa84MIbnkvBtgTIzh2Hlh/kcvObKV5SfaBZ25EjvtNsQUix8rf+sL8M/BOB9ZvY59dJjALwOwN+Y2dfWdpqZPdjd9939IfXz/e7+QQA3mNnDq9f5SQD+66bad3iEIjDaPTToeLY8uzAGYKh7Q4hHUKRTW1/y5FI71cQE79Jm9h+ajO01oTV9Z9jPhO2U8AjzQOKwacKPBSB5sTubWK2L28C20hAa8UKwEID8Miubi4BlPqWPuv/c1sgj/dely1AjCZBSNqmwilQTmvR2h0DOEGSmIYidsROMJ05Mxx7qgipIWRclM8t3c03RYtG2xU4uSKdCG3C0VHoagBeY2ZsAPATAjwD4pwCeYmZvBHAtgK+cyPsvAPwCgHcAeCeA3wUAM/sqM7sOwCMA/Hczu/p2tu6Qqc/sMQxy7+xe2e6MDgEK6usCkTNbTdRttnzHh674ido3qDhhG2OnROQN9MJq8ZQ6GYKPy2ZhwAdjyCEZYdPM2snlZSEq7T7XEWnm1lT/jviaTmoJ7h4C01kNVIQe6fR/kC6CM0J+TLKNLkwHUWdnVsgCn6MvgPH0bSYWjFlZ2UKnFPzHOFyQM2aRRC5MUBd4fhAzy+2mzb1uwN3fAEDtjp/EMrxmVd7XAXhQcv03AfzmqfB1KJBiGkZCalQ34AMNTGH1SLcqljDKAaGfKVsY8zm18lN5LS+r3opAM+TB6TPzQeThhYMdUBEwnNnisjI5BlERtDoeFAmqcFgsOkSn8X9sA+T+ahOXaeo5c58tlp7qzNPsvLgFXxU1c1SC8qgItW0UyJw42h5ul4Y/TdnCZfHv/vO9eFYrbLuN7eiLZKHdGG0OKR5KOhRIMby+jWZLY7/F72ylVPTH+VVtFgTZ6jZ5F3OGAkADb8puVhIt0QmnDwePICVMCWtGl1H0/n7eT7xgsMBUhKJ9mtm/tE/ZPgf0yFQRNE/mZHEZnBpRPlNmqlAiwTw4woDlwR7twqhO+/4C8b4cJvXWegjVQM1Hdot9OFPRI111gHRRBixUp9rF92a9Q6kJ90CrjIDde6ch93OEdqkn+lTIAT+Y9/mspUMhFBvJRO7UTEVpqkqyIEkG3RB7SHVmnrnOQQE0tUgRkqZL64+6aTKu9QayUA0eecDLQQYpJUI+RWfaFm1HlKM2wQzxZsRlJ+aR9nVfXlg1IWCbmhkLHqWz+axvh9TB70NmARp26tY/EccItHHmx08s267E6rUKKu5Tbof0UWZbtOwe74rS56UIUZHrRujcFoqHQ30GetTBH6CbxN2uAXUeZM4T2XnReZ/V1rOQOMIMOQavibdwSCO/bXen9wZmbVUKfuPlRYyMp2xbfH9GwcHcT3pNkcxEv/D3SXNCrbPd537P+jTscYyoVVVnbUFsqdqOzuEhaVi1DwHaxaxm/aFOPDbLcFqtjwWWLijZ82N1nrd6Sp/ZfA47enQoe4hKYM1kk7RVn88MNaN7tqqyZzQesDo6gsROFmVmsXYDehRVuQtY1j3WLMSVB55YPDh5Zc9412sTE7NzUGQqp4StdG2fUtchtlvhhVFetKB7YTu3T/ozE7CdCs+2wileM/UzWTT1OTQVW1Bm8Btxoz4jIUfCfCh/3XNnVJ3Zyqkdk+hZgtE1wL6dlp61mX+z2r4p9Rk464XeOjo0QrGLvUqoC1yOQUMBrG3gsHE9yps6h06RzoQqPSC5usoP6k1853Qh1Ol7lw80OZS/cjM6KH9bILc1+kXV/JhcLgflstDWQOCpZ6EmCra7BQ9ZeVEfeVfT9q7qSyknfckZC3H639lFg++axxh5TtSp8Z7d/WyxqdeHfPrMOZ/2i9iVB7MD9cEwLjJVexPkwIaCtw8tHQr1uTlTyCM6qHxBPCBp4nUIg9VqFSCL5cudwlPXDiQgA3mrSycg0O2FXolosslNqDXi5QYVXBFYSTwKdt65wcIz/vNEZP6YLy2TQzk4baBPOX+wszXys5Dn1pk9qH8ifnJ4adOMDnZg/hiJZQIqxk48t8xEwQHYMQaOHh3HGtGwEPEixnXFswtUrUJZUXm2wC07bezPBH2nu2mojzcmEImtdZ+zmQ4FUnRgtd0jJmPLsEIt4cmtKkZMgKwOs1w48aBSYTUlAEVFSwftguIVa/7uBVZcZ1YHsDx8QhCp+/L06NTDXXlp6Du2njHaZJQypSpGXaqyCophBK6e0GGfOfVPYXWFbVd5CkEUaDUbJ7x4Rt3zObCzU17RUO8NWyOFP5vP4TwuWUNh5xELNqKmLbBZZG9vFL78TCcWys48EG0joeyr5tbJ0DnufT4cSDFBJm11rA+47bBQlYkfuMQGDtukMkGTqdg5k/3vEHwTntxW5qIPwlVBw+ioCcnMCRRsRL/o4QGZ+s2oGQAWvXrWeZqBvv/UtBDoi+PlgFyYxStSE0Ta+oIdXtoXTCIQeVdNyxPlcX2sDnPfcB3iMBpeScH9NCPn0f5+setF2dKP7NyLcdvSTdk/6/N35U8Rp+ZL+ohp0H42QObrP2czHQqkODggaFDqSssqa/ewk90Rk8gvQ3h6TVdfXv2zuC+a2M22OZ+3U0umBu3g3EhUzJpwsCO1+yLMG29iT+qCkev9QaQzMmGVj9Eoq9jBM/PKbZmixK6oITarYkgDaQ0OJ1Yvw+YsavnABwC/7bZOgA5qr6JZ1Qh43GVB1mJ6mErTakgW7gEtUnva86rjQWM4NyYYzwHv8jo6HEIRyJEaD0QWaDywMtTHjhi2PVo5568F4AKjtxvi7WOBIJMqM5Yrv6l3kXnl61OokxEOIzDOt8ZRwh525Ts10mfhNmqmUAEevJAqqVEFvr/I340t9XCYjO8vyvu5SQ3v+pWFU4acxUQByNiJ9vJYUPQl5QzPTxd2fh6E/NtzyNT76GOyOXexsjIOA7V2p8+zUI98wHT41O0mO+cdLYdDKGZ2lFXpCBVMevHq/UaxDc5lIvP3Wo5lk4jQk3E6FhI8qVbZAlvV3pfB9qXEFjW89jLKpMWii7PketiOJ+hmCIXKhIGqwiEQGfGwoFA+WnsSRL5kbBQACU+pB16fuQToD4tWVj/XE2aCQM6sHUQ/Bxpb5+EVdTpdJJXm87KPPUgXZraB8xhQO3Drog0Ksi1SPAPE6GfViqZqUAimUJOYEpTXXusZ9w9Ag8ClnSSdUGOEFOzKRCrXFu39J11oiqijnfCrZbeYSm5fUMRjJjbXFsaU2PjU5mZHj8BvuXV8F7WaHmpb7YLzgWPH4Lcdq5crf7K48VvnGq0wUQz86zOJNoc5QPnLFhitrz0TQo9cLqXJeGgqKYeF6b7nKRNDspi1xVgPn8jGqoZnadtkkduoXfFgU+espTXL1RkiXYGD+NACNlbHw5/PYUeOlO8xiNjGxGmrqtEZ6Gu+ZqxP7HONQhUmG1XqQCFBlwVDt2PE2LkxJSgS1X5AGjxxwj6mpoGojydYNtFOnIDfcusyXby4KdqROYCOnwAWPjhgJpEKX2d0zf8zey23N8j7XTOB6JSG3SosoGZJSNaUw0d40Oer20d5LHcOv+QZDM5GRevxHBIvNhXS9Q3n3+TeZ7it/5zFdNqEopk93sz+2sorCJ9+MG6SFwwBue1qsSiTOIRdXOPvCfKIeLlUnSw/xgkQ9jCXuDlVLRWRqQDgwc2CrKpi3WSJyVo97n5ir+6dFsHNNkcWxuyQUbMC88aLhvA7u+wusKNHBrXyx4+9Aj9+68uLUDp2rB3wars7y/AeiWfUGMQuyJtRJQulVQIc0eW9Ojy1CaDbAgcsz5esnxCcqfDIBPmMztTMEF1mPtEyMwSueVU1l62TsTUSs9nyXTJ9B411nyKd697n0yIUrRwf/mwAXwzgAQC+zswesDZjMgl8b69XUyLY1+nNbzroVPUgI/eQBujtYixIQriQfS+bNJPHj+lEUXspIyUW2tR2nDhRFwrr7XFTqE/Vx3VqWAjm3d2BL7/l1qVJYHe3F1jhLmFkt3DYBRcUwSjqYQtNYfWbBSHzA0JtUwIj2sRIW55x+x6OiMgfwjP44QVCx038ZjXV+/Mc40VpLT6TnR3RWzR27LyjfVvWPcdAwLpYaD28KMRY1mD7TZAf4HMW0+myKT4MwDvc/V0AYGa/hnKC7lsmc+hgB5odrbtfqThNZKXmAa2ru6DD9MU+ql7RdrnmNYwwG5pUXWByZuBfiENIBTmjOvQOmO6wAiXlV6/rtUw1j9/0drvKBPzGm5aCfH8fP35rOcT4wYsPA0D7/fS7fQ382PEixG+7bfJtg921zCHA6mQEVu/stJjArg+Dfz3wIzGPDH2kyFlVa3rmLdRLEFhvG3W0g26rYEz3NVe+7fzzSwREXMvGrLS31ZY9W96mmPWpquNbWkmnS32+HMD76PfwCkIzu8rq6w+P+225OgqMq3h8Z4EY19j4rWqJIo5KmQG6C/oOlTRslSEcS8JlHkaYuvJnNh0epGIL5IWAbU0uCKWhEonvG/jg+lRtD5roHztypKnwUxDAjx1fIpJAbxp2ktkKGY0zHyywoyxG6YzoWZhF3bPkRfdsrwtUHmUwmpLA9NRbTIK3oFlbLpLcb4G+ue/397G44cblwsh2cym78T1lRojxqVpBTTscdbchOtfV5zvM++zuz0F5xyvuNL/MAbTVsw3lzN4S18Xg3Dy0fF08cH1ya8ii8yKzEAwKpCADrCHDWMll8nRvjZtCN+xECfuh9/FpU5MyfcGTpmcUKhNsOIiC+QxVMNT2+RzffX45JT4Q4ndf9CUVHfX7yZvaCxQEGotYEB9ioItY9Enlob1rJMpUNKT8Z04SrSP6NuyZXJYKoxltvxSv9MpzCg8SRRECV9ViQrwdrUOErXgJWcqiM06WHDjXt/mdLqH4fgD3oN/rX0FYH9wwSVWtztQHRUkqcNRWkw04VsG4TLNy4vLNJ5bq7Jzeq5upgiFYMy8qIUpXIQ5CpMxD8EEOn5YWSE+g7vpCBQjXFWlZYJDQbgerGm2zjDYvFsCRIyUs55Zbi3A8fmJp64r9wRPU4iN5wVFzRSIwlgXYKPDiP6N5XkQpXVd+pI1ySYgMzht6Jp3GQsK845mJ86ktVVGymgq4jBXlxyLccm0scLvSWY4E19HpEop/AeC+ZnZvFGH4RABfvzJHHEwQgyxzaISNh7x/3ZmHQI6YVMWKdHE/1Ke4xurM/j785ltakc1YvkNdVydRrNCTp5ZoPRDb5oQtVM/XM5D3dDbr61xlQwtBx/Yu6tvh5V595zd17N9c8qXFjra/v0RyM4PfttcJoKaexvM6frw7IchquZ2JI5vo2bONvNxvkS4TIiyAwnGkpogQalzX1Lja3cXskoux+OT148uttG7+rmg8iGMuebGtHuXOlqqmA+ZLyh9iMDdAZ7t6vI5Oi1B09z0zeyqAqwHMATzX3a89UGa1EwJL1bgzbpPjhNUKHtxxndMBuToZpBMoE6gsqChdIMx010QMZlJl0kNr+XsgpuCX6mx1sCo4hcrYzhQIbNHvful4DltYeDzZbrlY9KeuVBXXdneWAk4QcaTjiTkEObfi8gVl2ItN+YaTejIVmsuZzQr6v1V41H5UrYXr3t8v3vmoj9Ho1G9eoFmgh4CWfigsSIwlLyIZehVK33N+qrQViidH7v4yAC87YOLyP5vYJCTZKTIMlHCyRMxcoi619CwsZ7PRhrkKMYFQDt9ndQfiiQ5BywJ3d7cXrizEIw+Qb3cL9TBUTd0OxsSLxO5uKU8Fa9THQl/slXakTNzYuaIIyG87thRomboZ/Z7xpypztLvykh2EwNcGdD71/ELonTjRC/ZszAW/7Cyhxdr39+G33NIfSMHqN7DaLBDEe5YVCQLdIjao7My/ahlqm16xUNxu2go0RZZmAAAgAElEQVTF009pH2fb9ipN7U7obG66MnN6ttNkglDq6/LGntRaxxB6QQhsGMiZl1gHq05oqVsH/tAfmfq2WNo4m42PbXkicLpjyQKN3Xpb/6oHQbWlz70c689ohoQeQiCr9z7MAFXAd/GMs8SJRWg8tk2myDN5hhxF0PiLvlQvcFVfu77QcrgOduyp2tySJAsGo31us7ZBtShup47d2XIL7HBO4ynQueBdXkeHQih2gz5Bhgb0k4P/J6pN51FGr5INtrOgKVWa7TSM+Jgn5qcKzSFQfH/5prrKSHFe7OwsA3OnbIJUxkrkwX3I10KAcznxPbNRsbBkoZ4JCLWxqt10NiseaBZa0VZF2cASaUY6fRaBkiuP3VRXFVe/68KifGaecO7LiuSHeFRNl/HEbeDvujgyD7wIybjIwreG9tbyNnoYBICt9/lMkD7YGASZWsjoRR0tyb5Qm/dnBg6OGSAf5KIOt3QqKDh9qbTnnQTfgCD394ugysJwuMxMHUwQRNcfjIQYiaiHdUavBOU+ZlV3lVeVUTn3ezwjdiSJ2tepvNyXbEbJFrDoY76naJeT83OTYO9UnWcnjApLENrTZ8L263q9U3/1mTEtehvv0Obo/1j4k7yduUa1B63vFGiLFM8EqReuUjeYAh3wYFQVm9WPKdVCkFYj/s6DixwGqc1KwzLqtRjYsztfWoTesWPj/mZtO5Oon+kpM4yemT9uL3szpc3dUWjaP9pvasfSftNJDSy/T6FfLiMzcxCPDakvC+7VcFVZqa615gX1PEf52SKgamvWF2yy4G2IQqrBYOHLA0OYDxbMXBePZfb0R5/Q2NzoKTlboXiGiJEWkE8gsYelcYDyfXjFpaaDeIGZBw0pYcHEKpeiFJrYfsON4yZ9Rrsq7FQ1pgndHRLLXsVMLWtNSbbGEakTg3kdzgokvtp9tpcF/0m5w1FXqvpm/cc8JuaVYayowOPyZCHs+lBV19uzXzhTbzNNhseECKjl7hXRQg7Svno/zA7tEOVoz4Lsi5ugrU3xDBHbkTIKlBC7I4hStZPsXcMBqhmyA/qwnszWtEjeSczfJ1QU39vr39MceVrTCSnE74qC+JDR1C7Ek0URHi0c2V7qVqZM6kH4JItUqt4JDTZcrW+KkpjMVv8KUwlV3DtzgoLfWv5UsPvSxOEou3kkbEqeA5/SPgh9HWMS7ZCeAs4aEe+kYj4ZNbOw3dvLx/YGVefC+GaLO2x0OIQi0D/MhqIWJf4NWB5YEGl50qgNSZ0pvK0sM5DHIAsPZDawpiYz228ytRyJQEv2Tg9b36aENPOdfM+ColnQtQD5mq8EgRdUEgdedKpohjq13ggvonwhLDrhyGYJJhY66hUO0kk+tcUx7iv/iuho4eBFpQmoOdW/YkywcG0L0A7FbCbIt6WPnVH8fDOtQu3XGfoNXnSsrFm8ToZswzL2sNHhEYrAiFqAMUQD6NGQ2pIgqqhs8O/q4joFGXWxiLWsLLC4IYXKhzM6C5TKCEPtQ2x7SgS5UuNhYh+slYbmQkORX0ykXUIq2aTkdnJ5jJylrxhVD/3GNrxssgcfbC6RNg0aAquciuLr851aLBtlC4IKGRU2IoAsEBtISC4SJ4qi9SiL+0QdTlzfKjQo2sKU6WRLOW1+GTlZYu9n0IzOh8ucAZn9kb+LkEsHtw5yRUI0AAehsFgUlEVp20ETkTZUIZ78mdDO0I3SbLmlLz0lJoRwIE1Ws7gtnE8RaIaOF8nhBxnajXaJsOwWqUBDGnGAcDZI21lIZKhnsei3E8Z4UTQ1m/XmDxEq7bBi6YvhCLKkflVn24ukqO3pzpQol9rQXp3K5SeL6KRwZHI6FHmT5Af4nMV0eJAib4FjYcFqxZSaxDtBMrUTghKnhM5i0VTBAbmpQGxocNZPvil1TeP0mIfa7qnXeSrPnZ1wYrvcIKxi8qnTgngd9skm6mIjRllVZeyO4tL0xGO3dZGoQ1EsdBOEz2aThtR5QVABlyy4XG56WhD3RWZH1X7VZ6dbQaUfOsRIdsXUTKELZvQRm124/ljIN7jnuVR07jtaDgdS5IFTV+dmNK7XxixlFZzd+VLMLr1TuRir6ip7lP6WVT4E8HA+oiJNXb0DpfEg5PR1krqqteRl1J0czFuHWBiJ8vUoS4UBhTI1pKkCUdvDE0xV3/jNKH6X1tds8SLqbIv80fAV9Qor4qR+6U6mnhBETVtg/pSXKJeFKgtERtxTjjBFg6yV1DHg+2MZA69MMcZ2d3szyIxsoVm/8/jdFG0QKZrZ3Mxeb2a/c3tYsInXnZjZU+s1N7PLbk+ZQYdDKAaFl9csXynpIUcaP3ZsuTG/DvCw6bS8cU/RGd9nYcdCZ2BxKVBMJzALJc3PNiSebIpAmQQdDrYhRaOal4QdC9U2gaI8UvnSvbIhIHRBCPVsf79sA4wFSVVjVfWCLy5T+1DbIuFRkyqj1iv92KXR50XCb/Cw1/cBrVVHKXpgqK8K/lI27W6ixWAQrvw5caKcSqRCOZ5jNmYzs9Cp0mbV528H8Nbbk8FWv+7kTwA8FsB7bhcXRIdGKHb7lTNiNENoym+9rQRGk+3JeOKo7UuN0yQYUmRQ89j555UQH57wq2x2HGjOjpgQJJEvi4nLBEjwJ5OyEy6Jiqi2rda++EyVF3xEnxG/trNTDoigcj3QtTq+1gkoDVXiCU/I2llY0sIyLJ5UfofKQxgJP937dWSMNLTI4ypIkCV14KTzqNktZ8uTwZt9WBYtM1v2afZsuU+5bp1DOk5PkQzF+7zuc6CyzK4A8KUAfoGuPdTM/sjMrjGzq83s7knW9roTdz8OIF53And/vbu/+1TaeGiEYqCX7pWVQG/fUQSyWPSGaR48oabEfRV0PHmnnCBEftuxJTKaUOebENZByY4FAM05o4JaHRUQ+9JEn9iRI5h92mX4X+19bax1W1XeM/Y+53JvuReI3hpbLinUAgZRUdGKpNYPYm8tkbbBBNPaWpsSLUVNMEYkbdIfJm1t/EBrDQGaNiWxVqXQ1i9M0X6koEC5Ily1VIlcgkFMqVy8933fs/fsjzXH3M945phrn/ecfc7Z53WP5OTsvfZac40515zPHOMZY85l99zT809utahlnLnaIBdewYW/n5xMCcLCyTXQkShxNqCNdaL2YkvM697+O3joUj1seFY91j1jeTbdTu0K4t5mC8rxlAkz9ImRVZp5H8r5ucVoNu1oxJY7tUf3fJk2ailFssZ74KbftlRO8RSvI7jf6utG6t/Lk9J+EMB3AdObpM3sGMAPA3hpKeWLALwRwPcm12193cl5ZD8CLf5gvZPwAJTAgJLQKbcGhNSYdObke6urWMvp1u36/dQCZD04743vybyhrePgG30GDWzWW8CxrFbAo58CyroDGMvSOoTD9d98xxm2grKBWG6dTPmjrCvXn463J7WmlThzwMFl8qYP63UPKiSqZ+gnuibYgXoUpKqiVqgmwrfjGpSbFMitYqIzwvJFpXbEOGgWZX1TYntWzH07EKrVm3gE55LTFfXxUsrzRz+a2YsBfKyU8m4z+4p6+NkAngvgbbXtlwA+ei5dzyD7AYqgToLYYVqun3N9QOQYgfiZZ2X/P7cJgFpf9J/XXnf5gZOiTV8Aw6husAQ9OloHjSdP23LqvO0+rjfLDJiUGzfCsQAKbOmMuFWVWtcAqq3NLJTR7R3p95T2YPc1XcutKnCep04K8gzD5GUW+43m6mUTKbcZu7wlRnV1c4csnzVI0l8CwIsBACBf/ljLcn34dbdNR68r6VJKmV6ilfWds8pu8PWFAL7OzL4WwN0AngTgHwN4fynlBXyimT0NwH+sX38MwEO43ded3IbsDSiyO2W+eiUDL/48Ajig7SPXylfACjcfdMPs3CQFIqzaYD1G1gMNcpuUa1Zds0YIPLtdfBQwM8tVUzx40tD6uFW6FPD2uvh3r2/SLvrK2G4Fi7SFT24deHP7sVXXAbNQDmwJaTtrXZVDdp3UNeVJ1Qx2333AwlAe/VTe/uv1ZnVPjTBP9s5y036su/YF18f1ZIubVyJxW7He/Lx9UsPuZRcpOaWUVwN4NQBUS/E7AXwDgA+Y2QtKKf+zutPPKtOu/c9r9zc7wu2+7uQ2ZG84xcZ7MLGuv5PwO0rCOf5Hg5QHnQc6QvCiduDwgnP/XfnMKlk+YbMyyJJAdmxzQYw61vvpRqghQODiujmnpcDu5/OKnpFVuJYglgZv1NLhSWlBycoU1TaqW7Di1+uOi+xSZRgc+L5JG4bULW/Tqlv3wnqvV8Yv+vc19QM/18/znY4Y3IXLtKMjejYJpUBt3srxrIksQKP5mmyBs95eZ5+M/DlyPXYlu40+b4qdgiYvBfBPzewhAO8F8GXJeScA/HUnDwP4iQqcMLNvM7NHMFmPv2Zmr9frt8neWIoAuWE+a/IMK67kiJtSbqnjf5Qod1l4FHXqQB0Q6YD1z1lZbPFoZ0yspo7rwoZDaoMM6N1Bphj0Hq6HWjSDumRBgGGd1KXN2mmu7kJnDK3YjXLD9t+sVSY31Sca5vpWyQvCpG9AJ1Gpd7l5E7h5swtiBAv5scfyOnB7cvk8QWzjJh1s1ZLnptKAGlEDO5Fy+ujyqYss5ZcA/FL9/F4AX36Ka9LXnZRSXgvgtefRZ38sRRbm7nS29t8zl4u5pfWGC5zlUxyEKX9sqJNLZtXQ52aNatS5ntN+Y0K8HuM1su66aSJxmjbj5ak1xW3Cevt3kS5pmXRu55Mr3QbxaHLgNslWumSAyMCbcG+ZFdyE25Wef8pduv5enlMW3hZls0wu6J60s6fPjPJBWzqO/86cZwZYcxNE5kbrf/cglKLYhVyQpbgvsj+gmFk8owfu1hLP2IOH3q0SaT/k7mYDIx7omVUHV6UH3S6pO5nN7egopoAkbjQDZCtLI8es44xbxXUM+kpKSLBMCSBDgGHgjrW2UO5ssYjtSvpNYFkH7j335Mv7GEhGz1uBgp4ntxMDXerSkt5tovRbrKj+Up7XTyflFOxdX3aJ+Rw9N7P8tY9qfTzYxG2wIzllSs61lf0BRRYa/KEDi6sWchDZnXAuhvO+XNjamRMv18tRfkYBSPQLeXsqPDBIH50UWmeW3Mc0YZkBLwMOBUb/jVOI2JJJ2iPwjnpPBoRqIet7nzPxSKpZXZ3kwTEHAO8HJydx+zjin1OLN3MvK8h1eZ8MLlx/Dv7ZtIFrahWvB/mr7GJna9/ZC0r41vRZMAVRyx7K3G/nkYOleEmi7pEMpI6EZvHVD7dONp2L3ahsRmarE4gDhK0VdVUFCDur0stobhOBB1tgJG3lAp/D7mMWaGAXU/VWCytz/dn6OqW0KHuWQ6eiIM/udmZBHR83PUsFal6t04CM25uBni24Ba3x5t8ZuJTX42csFnXX9+gZ2HIJOz6Kmxm7bjyRiBU3pDLU2q7XZqDbraP356GTmxgT55LTAOI1B8X9CbQMOCl327I8uJa+4rP5tvKALnWk00GAMwRuXAdPKj46agGR7h6uE+WTdW/Yq5xhGhXWiC82QaSW6qHiVp+S8yt6FYACNA3akFDsAMXglliLhfVUkHILVwcvtW+YKMwCmKW79XjdPMFZV+DQucMNhTmAxxOE6qhW4yq+UsGfc9jdqPKHafBL217rBuRW9aIPRnVATe8Rt5GFuQMxXH/3eJvsByiyO+r/fZApoe/iHR/YWFPqIrrQ712OnUa2/Z6cXOyDiAeSGexJ9wGPPYbyqVtB15ZIS3p2QlZNqJuCCbdLfa8ygLD6pO3ezOWKDJcKstUDxJ2gF/Ja2WTAKUA00dUbbEmpG75YBGBp739WgOBJStN1uDz/nDyHsAmwWtpzbSiA2KxX/0+rUDgBnJ9Lmpzt17l+tEFtqNeMtJVL04XxPTtVl51Fn3EAxUuR1saZy5gsm5vdd1BlRGK3wkq/izUPpCz6WUG7PPro5iVBbGG5XhosYE5MBqmX0RK3nRIoyc7Ny+iKhlSOLDo5B8zeDKsV7K674qShZXl9FNRcrwFYNWsqsziBuo1+merF5bmOnGuZBciyyVBBxT0OAt5Okvbj5xPAT+854jfpezdx8UQleoXVR9o/lsuwHLK9hsC/e6I3W819bc8udzgo7gWnGAh6l+USi0//tI7zGgEicy7MRTVezyUDOuUN1V1TqZZMuXlrU+ZUEbTIuJfpAYEMnIG4MQLioOdE3LC0KwM51pdc385V80FG6SLtHH8xmNadAz2ZqzmgPvyv7XRDwOGR3QY4y76M4Jr7vbJEZ7PNPoN+b7+O+UZd7qn6zvG3Wf1HNEwGuHrMJ03uJ+vN4oJg2Xuqj3sh6mVoe7huTnFkyePnkTucU9w5KJrZ95nZb5jZr5nZm83sKafTZLGJ9nlHePRTsWNDeKJsQC4W0ZpYrxH2rstcW+XK/DqXUrrfW+fU/DXwZZtrwioa1q2WHwagWl/q5jHoAa3TK/HeBTqoPjq4001K1d1lYbBx0bXaoSiLHCHpOXqOrY4OIK6LnsvPlnlfKasLRPFkxhkGfq6k5HSeSVZfL4v52yw4xLrwSiBfHSR17HJNE5Brq7KY6vFrt7jgp5ayPR3nurvXF2Epvg3Ac0spnwfgt1DXN55K3MX0DT1v3MhnuARYutUKwCad4uho8zu7tNTh1YXlsgNvBDSAmzrhOgAEz8rZNlyt3PUmsTit3zZrlQMhVbcOqAUIhzl64i62Nrv3iZutqPj8NUVDRddgsdSJLrTD3AAVEA9bZzGgriXwVuucArtfw1aZtieDW8Z78m8MbKwXnzvi77gdmWbgZ+KTIF/D9+RJUrlCfSUr0G9ltws5WIq3J6WUX6hrEwHgHZjWIG67SLTa8G1pOoQO0sy6447rf+ymk8sSzsksTV5L6vwXgOnl5bL2NQkObNxfyUMksGggw+4z56IpiK3iK0TBIMWcZtbcUi6vvHEdy8kJcPNWB4ZBGNy5/uFmApJ6HVMgJyfRmnad7rkb9ifu6QGC64IKyGy5q77MWesxtkQzy5HrkvW/rM1dzySFKVA8Wh7Vzde1+6a+jZ5hY4AtTDYO+DnfRurVNtnVJrP7KhfNKX4zgJ/NfjCzl1vdgPImbggHuBnwaR4g83Bsqahr4R3v1q1Nh814IHZz2DLJeCdgmpGZDxMXK3A4PgjZpRq4inqsc5mUQ+J2UU6Uj9XPukIjUAHeNmwx37gxjO62qPMcx8buv38nizVMINqW7gYuFrD77p2CB0DXtmzBz3G3ASS8X2QTEP/ueiVeQHuuSn1wGQPxvhMsd9S2HIAe1tJ3lRulQFcol7nIHcmd7j6fKfpsZr8I4DOTn15TSnlLPec1AE4AvCkro5TyOgCvA4AnLz5904xr2bRBZtAWnaXfWn6df/YBKzmJ/L1bD8oW4RY+itNSgqXmv/OM7b8pj8UJ4iCulM/N8vQcyHjtLE8EGUjy8cVivI1aKXGDU7XM9LqFRHFZFwYIzp9U3Ws9QxoT62uG9e//QV29VKp1LnXMaAj9LuAb9Fqvo9upliP80j7ABz1Glv9wVyJ/hqKjroBqfevWrckj9XbTnN1MT36OwLSZxS7kDnCPt8mZQLGU8qK5383smwC8GMBXl9MkSLG7oh1OjmczHhP4DKjqyvLACx2S3RzPR/NBkgUrdAC7MGgokHgdCNSa28q/sdXmOvH9V3HTiBSEgD51hHXMBpX/xjrT8S75m9uDr3UQ9fYj4An35PZbrTbjTO9Ndd3kaCYb1Gp7q078HBg0fDLVHEMgfy0E9ykvi9rHU4/MLdFR9x9M/GFSZ2uaPR3pS52FntA3u8xTPIDibYqZPYjpvQt/sZTyR6e+0NM2GLSY/5EZOliB2vmUB8v4IS+Tri2rNZisdqvUnvAEoKxRbt2KYKPBG+2M2FixYdCYbdb4av0S0GiWsKfnZFZQBcxO2Grx8ll/BS3n+iTJObjy3mbK//GgJWk5iuzaJ23F509FRvrA1mKBCQCGeuo92BL3lTDUfpwLmnGB5vcjj4Dbu+nkxzIrcTTp+D15Ul/L7tyD80N/Uaud7rsr99lw/d3jbXIRnOKPALgP03sW3mtmP3aaizyay3xL64SJVRce8mrVB2Xo/AZAmVvMbscyci8t6HDzJsrjNxpYdPlzNHuHdcHE5XT8owMEu6ouGSfo98vcxMT1DMnRDBQ02YRXwZKbzEGhjtccCT8vqzmVymdlQEiTCS/vC69hZR25zlJWuA+nvJhhcf+nxW35lZNdLAIgdXVmtxXIrXBQWzsnqvXlOnB/SWgOf44hICXUUTumludpn9sZxNZl6991lp1biqWUP3eW69pKDf7ukpHKc9wWz6I8S2Y8pVgaALnWyQytbqQCcZcSkblyXge2rJgmEF5K9QmWlLqNA26OLazunSrsquu6YJU5i4PODwnztcymt4JLvTZsLMztRRZ5q/MAlLhcrvP69/9go0stK9U/owfYctRJSa+jvto9Jwcrek6zr8nQewA9z+v9ja1cdf+9f+1CDpziJQlFIzsXiAesP3DKF5xNrvVjQBok0PNs5ILOcFgm4Njceu6Y7GZSnYe6rGP+WRep5eMDN10TpZkS8OMt4KGRV/7PLra3xWiA8bmc7lTrmrqD/kyY62QQ0TY6Pu7vK9aXHzMvSwNdM1Hqjr6hNuz4ysytpXP1mYU1yn7fzBJ2MYuAmtEGGhRTi5LbYEdycJ8vQQoQ+CsfvC3vTK0CjmJmFs9ots86Ol+jLlWmK4MFzdopZ+OAoMKDV+rk1zUXMrlvSFeakwR0Q5oTc61ebwfZFb3cnvlABxUPOvB9xLULK25Ej5Rq8Do7FVEo4ZyB6dateH9qt67uDPKrVYjCponsWfuxnnTftoJkxPl5GVzXet+uH7FkXhI9qy6/09NwmNrhDAieFHYh5RR/11j2wlL0xxX4vFFklT/rg1a3MesIyt2Jq92u9esFkE0HCpXZuc48u+sA5PsJeDVrsw6EAH5cRxooxoDFbcBgppY0JW2H9vC6sGSWrtwnRMXXtIsMWTi8u09zl+uEVE5OwmTjmhV+Bj5xJQDcuOgqRr/zOe13pWGAPlCTURTt9+heh/LFEjW6vn328hjwJI2s6z9MFaneU6XiBrx+bIe84sFSvEzRhOnM4lvT7KxEdjbj+6y6TqKLeq+5JVHqbrFFOcdVqrWl4JTcr8uf1M8+kPTYep0O9GARsd7s4rPVPSkR/9d2D0nOLAt65QC3kd+nTRzSbgw0zK9SncOA5/IYFGq904R6OSfUOetj3GbefwaW3ggEw58+a/rOG8JuNs6g5+RRfv+d+jEHplp7bKNqdiEHS/ESRM185lD02IJSM/wancGBtBM0Fz2ZaRvfwwPVuU6/H7mWwIw1xVYADd6mc2atchns+qn1qPWac/V1KeCibozLkVFvX7YaXSfhongVRoActojVRVY+lCL1ob4ja2yLZLxy4BJHbUQAFHYfunUrn0QSFzkDxK5tarkhLYnqzWlarbxlD+jdBAB09Q79jvWfFOvLPIsUXPtlfNtkP0CRHxiDHItYhgEgEt4lXEcgaxzx3RQGq5xMt8qgdmTu0EMuT/UmMIS6RRV0ysnJOMBTdcuiycP70v92N7aA/D4MVBwcYX3Vuq3ndmDAwSTRJd230O/h5XFe6Jzb7kEeWq3RWdUMAqcA27Bl2WrVg+y6z0pIV7awW8yTfD22AbxlB9Zdb8oCI9yW2UIH0ASxxTg4jxjufPd5P0BRLZRRx3bR9BXm/hIesJ2XEP7dJqJ8jd+Lzx/xRlSuRtCV60qFeLlm0VIbbAuqpK9ZIIuX9U1XTZCEPNGEGmi6KS/KXKzrPrBgmDN1gO5SoRKLPrQX11OFr+VUJI0I3zoBFtY42S7NiduDjzNNoML9T/k8DrAwrcK/FUpf4uc055mg57TbNbsMsrCed6jsB6eojZxZgS5uMTJHqNdnHWEwcEIkFhQxVU5JrUQHAAJN53wYZHTRf9ODgTjpZOm+kSOgIC6tDVi3qMQ67fgosXoCT8WDrpbXNpVgTk95QAZhndRWmw1T2+BPggIhqXkEKv6fXffsOVf9QpqSR5Gz93y73t4Gem1WZ2m/wPHJMw/fuZ/L8wq0h64IGgUReWL28jWR/Jxy2BDiMkRnXeV0/ByyPLq8QHXdsrQGP+5uWGJZNqvRt3Sv0pYA6syvVimVmVlh7ZyBJWY60DKLN+P8uC7V/VX+Ti1EtyQ6azkr2y06nax4oAofm1lRXbI7tydZtbPBpsw6m7Ms1XqqdeP3tbQ6s3ubTbrqpdTynXPO8knneD1bLqeoO+dOKhXifeOz/gzWdx3BfutD4b5ATunMejVnlTsgkLJN9gMUgdipyaIJpPlakpMFHAFKTeEZWrmmzPLQga3cowKiAjkT/PzdyyYerKMJGLT4mqyNzIDjY5SbNzduNvNVLARq6joy19cB4RxP53qOLHlvA+Fs+bzUzdeykhU+6YQnZXRupnKIBBJtIshWuWTgm93XXf/B/dI6ev9cLuMkoZau6FDMUI4X+S5NTEVcsBwCLZcl0hlCGoYPjuVysyJgHSObw+VSfowHCr1AimfaNtiTmT5EvBU85gBiZNno7+t1xwmlYCBWiyWEeyZdezpIs9Wpkq2gWC43XKdfw0vq2KoUS6ulnHDgKeEsN0qL5eYy4wp29QTiPbSPaBurLvw5W8njdc2WHW579hJcaxOX9RsSAwB+58NYqBHA/K1SIFlAaAdyp4PifnCKQOTNXJwP8d+B7durj1yr1WoCQ+mIvOV9I//93mpF6ufM5QFy939Fm3/ODR7h8gJHRZ+zlR/BClb+aU3v7/B2VRpAwcYHHm+koL/z4GTulfmxhBpIaQVt26pTt2kC/1cLjH/X50nSrWbxsrJJ1T0X10nB5hSA2FbnZOuns0Bf1ibEb6crmqqu3W86rs4jpeqx7e8ay/5YilnEkiXprP7Y26DJ8gfJdRFctHUAACAASURBVB1aKG4pqdtN9wLI5VQLhEXdIa5ixhdmXNhiAbvrro3lqm43D8LMSla9PNBAyb+pLn49W2ccOeffiSYIOXhO9PMOPCLprj3DkwdWE+ss3GnXRtQmbGGFdphzeyvQNM/Eyz4+niZabi9NgOfAlT4frTunNmUpWlxntaDZc6oSPJ8dutXXPZCyTfYHFDOrzEX5PvktRISzwZbltHG53HGydbqj9Istrq1vaWbL5WRtoYJFYrWo297W6A4AI7zUXfjVNDeQEtO7yDa3BesgAOKfq9LtusBr1ft4vTXpnNe2h/JclCPL0qu4zVyXkatI3zVR20YTsYLpgEu0u++OifDcjxi0MsuQrXqmGziqnhkI8rkBoVIp/H3H7vMh0HIJ4m5FIPxlWZWNgImPMbBocGEOWDc3it99kEtna+eqdSXAaPfdCzz6qenLPXfDjo9R/u8nNvVht5I6bwOMGgHXKGJY3ys8ZFdfrmsWGfVIrAIqWUfts+5LOQNWYT/GhBt1bth52s4KU1EecCU7pGecoICc8fX6PLk+al1lHshqhfL446ENw7Xr9WRJZhOnrE/vckw5wMR6cfnUB9pqLK3HBYjhYCleqpTVxJt1r6/0z9z5MuvQj2/jNITrSS0nYNN568Auq9VmRQzfkyN/5K7jscc3q0Q++ShK1a0LqPC91xRwqQDBaTSTOpvvbGFy8nOa7Mvt4//VYhvxmaLHcIdrvY/ULQzoxWYVyHBlilptg4mktXsGIHMTIrdFsCpXm12YHKj9swO99w+mU1wn/j2rG/GeacBpBGx6nv7G978IKdd/E9ltsjeBlinNZhEBsf4PS5fW6+iyeIdiMpzI6OGApfvqsY5051l4vY6DwAc5g4oHQKoLXEqZVk5ks/nxccyZZKmdnK013i6tq5+Xs1j0wKvtyq5V5t4tTkHY83cHZnWvdRDzoB0M7MIAVaOx4T5zlr7WKQMIBUmdRBfTGvFWDvOOfD4HkrhMpV44SFMn1/Ac2JJkHbPnonryM5qhDnYq5RR/W8TMnmZmbzezD5jZ+83s229HBTP7IjN7n5l90Mxea7URzOzra3lrM3v+7VVskr0AxY6x0QHkIsRys3QYBNupJZaRcX/aGcV9bcd0wAP9agbv9BW0W7SxuoV2fBTPr/f3yaBblaEgts0ymLGQ26oR0jFYgj6QnQIYrcOueoSJhMErc81J13A/P99Bnq4Nq2kYgBKA7Jbfcdvwf9aRy8+saS7DAU3bN5sc1W33fqG0CHsUDOIEll1wi9tiVD+pW5cDuSPZ0YqWEwCvKqU8B8CXAniFmT3nNtT4lwD+HoBn1r8H6/FfB/DXAfzX2ygryF6AYmoh8Wd5oF3EMLEuAWwGDB9XIPSOzzyd6JNumeXHFKzIouN0H7aO2BIqJycbK5K5usxNZ2uSgUX14XbRyYKW2TV3eBTp5Pv6scVism4reBbXk9so00ndRLaup0bv76uyWm14zbmJ0NsmS2ZX/o8msnCtTjAKiPp5ZP0qeAnITXw6bajLQcGsv44mR13S6TTM3KR6FikA1mX737ZiSvloKeU99fMnATwM4Klm9llm9nNm9m4z+29m9tl6rZn9KQBPKqW8o0wP/d8A+Ku1rIdLKb95niruD6fInZEJd+9sDCgr2amFXEo/NnznLrvEyewZorr1HmEH7NWq3xQ3s0L5XnyORAU1ehzSa9wVVr29PAZX11fdcB8wNWE94wLb6zinguK1av2wDl6eRvdVJ7kmiCfSq7XVbhD5vs6dp77S2s7rOIr+MneYRcG5jdnLGJXDz86vY6nldLrXYwWkQ2YBTopu9PX+7/1yFd9H0/Tx9ssS4M8jpyvqfjN7F31/XZne9d6JmT0dwBcAeCeANwP4llLK/zazPw/gRwF8lVzyVACP0PdH6rGdyP6AIgNKKZHgdlks+g67rUwgbInlHU87se5pFyw8HxhmMdrJrqPfRwEc6FNkGCD9mOcRZoNCgwpJ3e3oKBL/3J51j0DTyYP1SXi1VA8vPztnvc4nk9NwXTL5DfVwAFdOTrlDP6Z8sLvc69Le3li4DH62Du5ZShfrWful173LpPBzEmpiei5JvXVyrMeCvvwc3CLniTTTdQdySvf446WUrZyemd0L4KcAfAeANYAvA/Dvqa2fcDYtzy77A4oahGCuiCR0Np8BNdCR8DCa1qK5jR1fqGkRVJaX0a7PZBHfJdwAUlNIGMC0LOYWGcxcj4EFEFI8fAJwV1zrNADsdCuyzPKS9k5XZWidEjohuviUTziaUDIQ9OPcbomeZoaihph6H3J+Wg/6zOd0Vqe65NR/eR1zWjf1llAnWZqo256UGR/s1+/QUtxV9NnMjjEB4ptKKT9tZk8C8IlSyvPkvCWAd9evb8XEJz5ApzwA4CM7UQrYE04RiPxJZlmoqyo8YFeWn1O/D60kOj+Q9hknlQGhd3YHP+fRxFpqkWMWB14HRu/Y5Fq35Yl+D47MModIAYkugn8aEWsvtMdcMIJdVrpfiCCrLvUeXXK4gxDvhE3PpnB71bo2Hfy7W4L1r9WBg0O1Puw6d2lZBExdwILrmgFz1kbaHlPBm/Zgfnpk2am3oIaAWvAXIaNos/5tkRotfgOAh0sp3w8ApZQ/BPA7Zvb1fo6ZfX4pZVVKeV79+0ellI8C+EMz+9Jazt8C8JZdVXF/QPF2RDtRMhOmwQMKMnRlLRZxYDggJQBtdx3XUwik1Kpgl0yuD8JWhLrESrazu18G66hjAdM/H8T+P+MJqWwOFKXiZS2XPVB5GcJvZXXTlKLsPAenAFwOWEqlMMjWZ9csZGmrFgRTF5+pj8xSzADO28MDV0X20JwKid+1vRjg+f6ZG8768udBW+9SDICVsvXvFPJCAN8I4KvM7L3172sB/A0Af9fMHgLwfgAvGVz/9wG8HsAHAfwfAD8LAGb218zsEQAvAPCfzeznb7eO++E+zzWi2RTtFPeaAxKNH6QgwmhAG6d1VDCyzY9jC5RdKs4TJJ38vM5FZndc68b3W602QSOyPDP3rctDzAIaWg/WgTlPo2V8bCHdOplSibIJhspu7i9bcMrjjoSBRyeIzI30z0C35njxlCdj/Yn/13NqdJ/OHdWAHrvztfxsbXdXrpfNtElC43QuMutC927twTx7aH56p9CCOPHlMu8z22iN25EdYG4p5b8D/WOq8uDgOF//LgDPTY6/GVOw5syyH6CowrxWFmQAQqCkRYuBzSJ9GvTAACSVNxT+h4GC71cee3x+cOgxT9otFCHkCDufy4PJLYiRdeeTBXOwc1wn13UdedWsd9pykYMJ61cHYQDOJLDFOnQrYrT8+p3LNbLC0rXdi8WUrkMWpO4hGSK3yjsy+HJdHKDW9fWscttQF9Kre25eNvN+8jzaZE9tFYCZn5c+i0R0xdOu5JSW4LWV/XGfdfCxxaUBjmx2Xi5h9z5xisL6MTNgnWzCqgED4Wns+GjjetbOnS499MHFA4r5qoR30h19dAYPbj/VU12yYIkk7dF0YyHdO/cus2L19wW9JoAAZWQJTm4sUR3kCodVHa6vuLIhii0DO7SHg94nH20TYgB9f4aUXcATieabdkBeP3dcp0xoZSUTkk5QHEzi9uf7KTUycoelPbJAz2we6lllR5ziPsuFgaKZvcrMipndv+3c0IY+OzKP4p2CBodyMFivp/XFN26IZZlEAqcCog4+MEpBuTkFNuzoaDNQVRcgB1cPmrhkndcHQcZntaJjz3IOjBOvw447zj/55ywPk84dRlTZYuX2KaWzvAAacA4StNRw0nsxBM2u3dhy18mIgC0AXGYtCcDM7nyEhDN0fpB06kCM21FBM9Nl05D9bzrp8318IuF2kMmyu8/xcTg25IbPJNPa521/11kuxH02s6cB+BoAv3vK8+ddAe08ystUmYCtdJ0zuBHEvbnrk1pnalUSdxWAT6wv004s9WxcENWjuZkeTDDrd9TmFA6xpFoZrONqs21Zu96F6+aDX91XsdxDfcKEVIBjAhFZqpYCorqJ/FxcBsn17b7aX0Rfp1OCa8uWqdAU4f5JIKVx1swB+72y5PWsP9NzNz5vungqd0k7m2f9Uo5192TvJWufXcjBfT6T/ACA78JpDekBKW3HR7HT+J+LdOYpYkr8UTuWJyw3N65TJ3bIZkXyYEo6ht6ni5Bj7GpyvRtYsYvGlpTXne/tg4F0GwFSGNhq7epnvpeUP9U54c3quV0gSjlPdouZl2zuKH1nIGJw8+eSuaj83Pw3niQ4Gu9tojpyH6O24mebRdKn06mdjo5iO7L+XLb3vaS8zm3X4yPZJa9YAFtv/7vOsnNL0cxeAuAjpZSHtibyZuKdcLWa3FgFIxV2dZgXkmV/wWqhDtJmbDo+eoF7reCmMxPP2Jb/cdkKklngwc9jayapT7s+Cw6wpZpYsMEa1PqqpSh1UF1aQIujttMP/eDLEs4zK9/ie7Gbi6z11Ou0Dn6fRRI4GtWJ81d9MsosaVDgLytfgW6bZcY8ZKF3fZMn0+k/4oi5TXiSW9BejQdL8dRyJlA0s18E8JnJT68B8D2YXOdtZbwcwMsB4G574uYHdeFGA8cDGicn3dK7NnC9DNSB5p2cLQ/tzKNZNzsmA1kjhJ1LyMfZQuL713Lt6Chyhuxycrkc6eUJguufgQqDmYKmtgs2gNvuoe72SLLBzf/hVbI4aWTW50jYss7qy1Ij9g3gtK0y2kSt/JElLRLWzDPg0gSSmg0j/bnvZNakL42l++yWT/TK7L7IfZIzgWIp5UXZcTP7XADPAOBW4gMA3mNmX1JK+T0p43UAXgcAT17ePzWzzrLcyQf5evXGYwttDswyqwqIA2ZNUXAeRIlbFaxAjZzzAJiLKKrrSDoBSFMsgoWRucP+XV1PbgMGy8Q9Td+SJ/fxyahNDhlXyVaaWqqcEqVWLVv7I+5T68p6uf41F1S9gda+OsmwSLuEiUj1GX3W3bi5nf2+65r+49vNqSVIK6aCMF2jHPAOrbtZT+oOkJ26z6WU9wH4DP9uZh8C8PxSyse3XDj9HxHUQOh0bAl2s/3IQuN7KWckIMypDMH1UPJfy3JAPz7eRIoZTCstkFqJXn8X2VyCE6wDoIDAeAQ4Sh/wfTMyf5t1yTon0k1e/IwkSMHldoDH5+l3BnnWm0UtPD7uVpXq7BTK5of0JVFduQrOCXC34M/dT0B5/EbkKAW0mBvvglya35pxjOwB7Dol587GxD1J3nYQcrBzLmRg7YTBzSDIicx8TS0zfNfP9D0k+QKbAeTEvAQEmviAWi6mPRL9mPJtIxdv5JJ6uQqwrpsCjotaqeTClZOTDQcq1AOAfpcV/q8cpgJJBtBJPfl5B2Frh9z6bqqToEqw1AnMPAm8/U6TQhc5ZqvKDHbvE4F1QfnkJze663NSvlSFJ9vVCuWxnlIYJVmnSe7eRgva7UknBx0HOxLDqZfxXVu5UFAspTz9Ns6dPugsrG6Quncs3jHECuBy9X7ZfcIAWq83Za2mTU7D6grlmRaLKUDkx0aWTtYGPkD5ILvF6vJq56S26ZamuYUBTO6j5kiSTi2XcrXqJ4gRTeF6EUAHS1tdcq0f0D8ztnYy654sps4aHP3GVtmW5YtYraaEcAZptawV4L1OPFkJNdN04iAPc5naX0bWMFM7pF86Tnbp8h5A8XIk7AfIAEUAsJXodiArJVoJrajN4Ojee0yfg8XBnV5SL7pcOymvcWzMhWV8Zx345tcqmOpGCKq3ThKLAYGvQQQWrpO34SLJJ2RhcM4sQq+7BMLa5Q4coyV/CoRqffI5mX6sCz9Hf0viAOwhIIrFYqOjuv3uvbAe2r4MVlw2n8cTiuqi9BDXaxvYsfezKzmA4sWL5gUG0dQPIG4QMdM5sg0POrLeOw2lWzQh6yoA5dwABl9uwlGuI7iPgiNs8fLveo62C39Wqya7h5/HvGMGwAP3t9t3kq1RkOunkwzXW65JdczqTZZ85/IrADOQrVYxMKbXMAixtcc5jyOaYzAx+DU2OidzdVkPmWRDuSw+kYVDlDWwC/ljwCkOSJDLF1++5svJAATAap0RCNG3QgO6zeiQzrCIy8PYVS8KONyZ+b48mPmYf5f/Ye1zFX5bYdvsYM4KIG4r3DfT1UGK20jLzAA0c8H9Oh3AXB8fyGvaEWZUlgOKu+U8CfK66DkZ8XUjGQBrlrPX1hp7HVxX7RseaR7pl1m42VJO9gb4N5eRNennuustxz2hnN8plOXFnldsvd76d51lLyxFoLeqAPRuqtlmF5zaKdwVCukucn2wJhebrZ04765tw5QNhpE7kz18zw/LIql0TcvLyyxCHzgj7jADqrkB5gCdRS792uUyvlBd9A33ZsulUgrd7jUjfVx4g41MGDjUxR2lmHB7sBXLk4VylBgEMzLr2S1bBZrMxXZJ9PR2C7yv6zWTPtO53mo9a9Aoa5dzSzm4z5chXQRyTetJuUO7heGddpRyopxaAgIBaAccUFQyRiVD4ENAN9xTryULrOmgmyp4fTzirW6xt4W3FdW94/CUfshcPwUudRP1PAeqm7eApeX3YTBiqzQDetEnXbM9AkEg5yyVOmBg3VZ/F35GI6BTWkU5R7+Wrm9eC983C7Cwbtn95fmEF6Bty9U9qxQcQPHSRDsSaMbTAc3ifE1iAbTPktrRrgPy4E7GcXFn18ACW2HZgGML1F2fjHivEvhP7YAZN5gkHGereoI+fp+6gUaniw5GbffVClhs9rF0yzik6YyeB1t/+lzUW9CJyIWzAwiIy2oN3hnJkkmkK48nIi+H68GilIbyoou4M0/nrmcASr/x1nIpj81tp3mzrWloTf+6RB57F3K9veOtsh+gyB2cXOaW3Z881HAsc3mBCHyjmZddVQejbPUMW0B8HZBHlEfXOGdKu/lo8CfwnirsulbLIFgeGQDqAFwTl7qwjWU952apFYgkJUeOdbyWW49zPCbfR6PNy+Xmlah+PiSyuiB3fg4IFLBBk/DxcaQSZlzaVtZiEfMhWUYTS2Khdssd3WvJKBMNkgWevUxtcXw0toTPKIc8xcsUseR4Z209z6Yfpk6RgRml5MwuiDfatbte26VpSIS04894gPHAV1eNXJuCPqrOyczplvbsevM1CSfWvT9aJp7gFiodoNwjTxh8vwwseGIbLZPzgZ5NeHPt6lZ2AmipjDwL7iN+H63DSEbljdZ4Z1Yi66Y0C09QI1pDPytdVK34TsddyQEUL1HYtUtm3C4ZWN3YdZk4LiBPu3ARvqflyzE46OCkzhiW1WUL892lHVgCAAEGd2i5pnOfGCAceBTA48364xlQZG5Ydp4Lc4Sui1qSDLYsGuTI3ODErWx6JsGMzjLMuEC+RvsNTwJzOiQWvGXPot5vFvCTZ9kWBpAFHOqhz0gtb5rUg/exSykFWO24zD2THTKw55CEOPYAS6Hf7Oionp5HmJ0bC6IdKUllaB1HgcytE//z8nQ3Gh9k/tdSI+hVngoeXibpw3xgpzt9nnWt/c8nDRWuhwOXpo1IYCI8A40aK6Dxb9nE4O2jbcnljXQm0TbN+Ldur8xMR43OZ+DjUifllsTvO0wLSHuKT1j/zvcclK9vXQwTHlMK/Mz4GVcus1nzmne7K3HvbO7vGst+gCKQm/fEn22d8bLZ2oWBkB+Ydxp267KO6w+aQUOBPFshszDY0dGms6v1Cdz+Yn12xZX3yhJ0R9YCTwRsOXEd6/eQASC8IoDNcrXs/pmlqzroZxUKJDUPINNJpFlK2e8cKPFzOBgzcF05/9SWAqBivXY71YxEf0tevNaeDee36qRWz+ssxF2D4wEUL0myhlzQ+37d8mFrBWINOMCxNcZla8fQiGRmZbJrkrkvan1QB/QcynAPP280SJRj0t/WsmOO1oVBnq29k5O+jjpw9HgW/WQdGUx0MGTt1VQlENYJhu8h1wcrlYNspHeXxL/ioNKiB8A5SfQP9E0tT61B3qasPatRu9fzwm9cd36WmXVercPu+XjbZu10HimYLORtf9dY9oNTdJdAO6pH0zhvjWUt6Rv82TvFiFvyfe1GgQKyUp1zZJ5T+SLeKGJzUNyrzOWuutnREcLmC85LJvdqnBkPWv7O7nDlyobJ5Hqtt1UGNllkVaOj3nZ6H6Y5jJK91YJmnTIddfDT82vtxIGPhXWvIJ3l6shCMy2X60E6tr452NEmtJhOMskElraDlB1yEWcmoNZmO7PeClDubE5xP0ARiNFN6TAjwjisDFCrhgeQdkQvT93ygRvmoKybQfDa6m4T1tMQ3DrINc+R3VeuC38nTqysVpN1mrjAWXuElBmnAAag1uXs+f10Dbo+AxcZ0Hy8bdLgOst5XfDBAUT7Cbcj601t2fG1Xm+6bwB+X+njE0WSaN/KGLV7NuEj9vUOeBsvOnHloUwQEK/7hQvTdbTj0k4jz7jjAy37A4pVQmQ341QG1wTLQ10P77Ru1RBopVaAWokDLqwDCgUf7rBcn6S8cvPmuBzXjQGXV93w5MH3GAVFkACOTTlt5fEb6BK/VRePcjo9wOBUywv7YY7cUAb9BChVP/DzmNssg+vq/31ZqFiNrQwuc6DzbFI/0yz+38sY8Mbdrk3UbzbHyZIceFJTFeLrKExB+jST9GnlmnOG22R/QNEfepLb1m04m7gl4RWmzFm5cMdgtzU51n5zvRRkM70TCy+UPwJ4smDSXVv4Hux6Y2A1KWDMWNjhXoiWY2dFsk4KAMmEFSYMmTxm3Wa+PmuLTSXixEDtmb1ErKMyEiBXayzUzwMaymGr8DJU0itdITOwINvPPBY8HUl1XyTLBV2U2tiV3OGguOPWOocwWcwAsF5HzoZdFhYHSk5FULfh+LgPKIw4nOweGQ/lZLaU46k4vCMPFMBIbLls54Zr3ArhSOpaIs9i/bXBzboyKI2ikXQub/TA17dUGLaOuM1K6VNhFJxVt5EemSjdUI/FqHDO7aVlZMDMzxbIA14a2abPxXNXHUiXyWt2+bmO+EQV7Wt8f9dJqKaSPY9zySkiz9ccNPcDFHkQA5tBrx024QAVHNqWXECbXVt0kLbqKjqgGTwSC0O3M9PByds2Yb0JAIV8Sy9DU4HEagoWRbV8UEoAzRDkKZu8uGCd+QTBUegs+MTCA3SxCdD49m0BKLMBrQBN13btxr8rwLjwdXOW6SAQ1/6rXiw6Cc2JAqYf43xTznpQkPB6UIRawbdon/N7eB/menD7eB8jHcPktgspCJPz8O8ay36Aog8ucmFC+gW7xgwkHvxQ64EfzK1bgb9qrmo2c6uQRRB4x2RwZu945nK7RF6+J53TyqX/JjqE/3yOgoUPIrJcOheOAZqsv7Jaodw6CRZYV7ZY9K08PmexCMGf1gb0jNhybdZ19sx18sxkZCEyyOsEqLprGQpumWWnG4t4HSRrAItFv48mW96LzSKFkGfJ1p8HFkNq2jr3fLjeu5I73FLcC06xZB1uhoNLo3bCFZVSYHfd1eVwKVh1FpZ2LLbC9LqkLmFzBhcv0wcOz/IMZmItzwqBvFqWoaw5isAHfD2Hg1UTaU/WK1tCqp+sR+5e9UCJ3cNXpXrZGmhjgFSwGrnBPNko8CXtEp5l9gwUNDOPRXlal1F/83bxY9k6cdaZ2qDRA9SPTCcuuWZ3Uu746PNeWIpht5vR0ix66GGNq39mN9HLvHWr41Naao3eO3PfBhZc6MSLuAwunalVePYnYAguurtkWqa63VouyBpTK0v4waYDW78qyT2CNedl++cskdjvzc+CLUY+V55jJ8KHtrbfYqm342zRieWdudHhfK/Deh3vy/2QLbgMyKOSwTJs99bULC+Tj3HZHhib8UK2WtinlQKUst76d51lLyzFIN75+DuQWnthdufOS+dy2khLQGar5+go8pDA1ujt5mJrvM3QdQmni9vN/xeyiYEO5rJJs2B+KEsZ6vYQdGHQRG+Vpfl/DJpa3jZrVN30zHqeFIkWYVYOgy6dk1nJac6fW1lKt2TUAP3Wba2/WMSJVifxLNpLbRhWOWXnyLFSZE9NfjZuoVN9ec317GYh55FrvmJlm+yFpfjML3h6nCFdNG/MI3m+nphm0GAZ8CBhAHFObfqymWG9I/vflsEfLAsCtTBbk3vb6TVyA71O2QRA1kSox6gsPu7BEgLDzjpPBn/TwduI6sDgMrs2na0cLV8tNGqzohZZZoXzvdl914CTtkcmaoFnbifRFX4PzRgIG4DotRzoW60iuGbuMyIV4QG1JupZKHeoXsGuZEecopk9aGa/aWYfNLPvPq9aZvbqWtZvmtlfOut99gIUm2S8iH/3jq+dQjsKd+zRgKhAUE5O+nWn7nbyQB/xNCM3O6sX84cO8FrPaskWAfK0XSoNEAIhajX5NWwZwm9VNvfUtuN2WshuLSQOSEMdHQi8rRVo1IUWSzbcm1199gSyOvNz4vw+oTzK6Bnzd9GnAybXwzlTs5jypDQI9ZXADaqlzPenjIYuoKbniT56zrlFA2Wjvy1iZksA/wLAXwbwHADfYGbPOY0KZvah5NhzALwMwOcAeBDAj5rZ8iz32R/3WTu2pxIA/S430pmzROQOKBjonMtR11A6W3BB6mtVO1J+m+uolmDWYeReXYTYy2GusB5ru22PXKRMD/5NrSoG/BE35vf2umodVB99RtxGXLa7mBx0q+WMEs6DjNot+R7czeVyk/DNdaG0GY/4hvSWNa3j5ntxXxuBU7KKJ/wXyiIN5Ejdwutps3vuSk5pCW6RLwHwwVLKbwOAmf04gJeY2T0Avh/AvQA+DuCbSikfPUV5LwHw46WUGwB+x8w+WO+B7D4APjAqaH8sRaOcN2CyAE9O4kxKs9Roxg7lsWU54mXY9fKB6qDJM35zVQbcF0VYm/g5zPERDcCAFZLO1fXz+yeDvgsojGZrubaLYKqbzyDMdfFr9XmhuvgDNze14tytX1Gu27ArRwAACJhJREFUoj5Tmijcsgsua9ZWIqEfSR3DKx1mrMVuSR4/uyOxLdjl5+9c/8zF5nMkGKgTQuem6/Mg0B8+kzNJaa7/3B+A+83sXfT3cinoqQA+TN8fAfAAgB8G8NJSyhcBeCOA7z2lYll5T505PpT9sRSBbglV9iqBoUWlQFePddYec2h0bWrlJMGFthaVXD1z3dW1qWVMrx9ABOjMJed7qlXku/r49dIGbXOGQbsFgNZBm0VffWDWpWJhHbi7pTrQ6kAN53GduM24/is5rmDidSVLe25jCBXtJ+H+TfUYvOMgVAhIcd/J7pdYmf6Mu/tyKlMCxE2SPs4bkqR9n6iN4S5TZ5ECnDLQ8vFSyvNvs/RnA3gugLfVZ7sE8FEAMLPXAPj6et6fNrP31s//o5Tyitu8z6zsDygm7mv8Wd4xnJHp7A4qj+b8FK9woc5oOmDrwO/cFl5Pym4hA4lLdZ/aXobKVWbWk4uS8A66bM3Sb80dpOtm11O7a1c2qTLpC7ta+k2B5y22ddHJYDa3WLK6LDwYVTavRgX6HdO3uHydjjphjK7Pnis9x+FqFwbpbCIRYNMgVADa2g5tssyW4IklHUCVymadM963/bYtzel2ZTcpNx8B8DT6/gCAXwbwxFLKC7pblvK9qFajmX2olPK8U5T3kfp5dDyVC3GfzeyVZvYbZvZ+M/tnp7nm5x/918E1quWEh54uZwNiZ6Uld028U928Gcl1F56Bhbgf8pXaCdkV906oA4gijoWtwcxS4HLdVeJ6ct3reTpA7OhoAozj4+mcJLrNbdS1q4Pfcjm9FY4GfSD8ZeC2eqzXMbhUwadZ21Rf8/XBi7i7dWgL/8xt49abcs2kQ1ghs56Sn5urT+AF0KSi7azuKa9KEc+ka1MvjykWp4co4DO5nqSvGgokXTBO25yf0w6lACjrsvXvFPKrAJ5pZs8ws7swBUneCuBPmtkLqv7HZvY5p1TtrQBeZmZPMLNnAHgmgF+Zuc9Qdm4pmtlXYiIyP7+UcsPMPuMs5bgLYmpduZVBnRxABDK1ujKg4XJ5cPH5zB25ZcWf2RVlQKEB1vLSxN1lUj+UOXCjgmufuJCdxeEubmxUwC1tdaO57AyQavuohZVairzG2+uube7XJC5oN5BHA9vbTN1/aQeT7x0dgSRIoXVTKoDrQ+d10WEuk60/AsdABy1JWwm0tDKwsdLTF3Rxm6v+55VSdmIpllJOzOwfAPh5TG7yG0spD5nZSwG81syejAmffhDA+09R3vvN7CcwBVBOALyilLICgOQ+s+VdhPv8rQD+SY0CoZTysdNe2AbccjlFIBW4bp1srAx2lVW0M4w6RMZfcYdm3kyDJeRCp53PraLjo2ldKt9jU+FNGQ6II0vF2wfRSmBaIW0HdUu53IQz7Taa5UG3WPScmN6P78N14fc1a/0cSJ0mkLp1aVP8G1trrLuXsVp1IF5Kgen7o7UsrYvWSUVW2oQJmMWf+RwfKBHrlpC9Wk/DekAxhfEwx2WeU3a1604p5WcA/Iwcey+AL99y3dMHx5uLve0+c2KzEdwzSCVA34IpV+hxAN9ZSvnV5LyXA/CI1HMB/PpOFTm/3I8pJWBf5KDPvOybPsD+6fTsUsp95ynAzH4OU722ycdLKQ+e515XJWcCRTP7RQCfmfz0GkxI/XYA3wbgiwH8OwB/tszcyMzedYZI1YXKvul00Gde9k0fYP902jd99lXO5D6XUl40+s3MvhXAT1cQ/BUzW2OaWX7/bCoe5CAHOcjlyUVEn/8DgK8EADN7FoC7sF8uxEEOcpCDDOUiAi1vBPBGM/t1ADcB/O0517nK6y5Aj/PKvul00Gde9k0fYP902jd99lJ2Hmg5yEEOcpDrLBeSvH2QgxzkINdVDqB4kIMc5CAkewWKZ1keeNFiZq8ys2Jmp8nNumhdvq+2z6+Z2ZvN7ClXpMdONwc9py5PM7O3m9kHar/59qvUx8Wmvfz+l5n9pz3Q5Slm9pO17zzsy+gOksvegKIsD/wcAP/8ilWCmT0NwNcA+N2r1qXK2wA8t5TyeQB+C8CrL1sBO8fmoBckJwBeVUp5DoAvBfCKK9bH5dsBPHzVSlT5IQA/V0r5bACfj/3Ray9lb0AR51geeIHyAwC+C9M6+CuXUsovlFLqljt4B6YdPy5b2uagpZSbAHzTziuRUspHSynvqZ8/iWnAz+6Xd9FiZg8A+CsAXn+VelRdnoxp2dwbAKCUcrOU8omr1Wq/ZZ9A8VkA/oKZvdPMftnMvvgqlTGzlwD4SCnloavUY0a+GcDPXsF9b3vTzssSM3s6gC8A8M6r1QQ/iGky3d2C47PLMzAtnPhX1Z1/vZk98aqV2me51P0UtywPPALwaZhcoC8G8BNmNrs88IL1+R5MrvOlypxOpZS31HNeg8ltfNNl6rbPYmb3AvgpAN9RSvnDK9TjxQA+Vkp5t5l9xVXpQXIE4AsBvLKU8k4z+yEA3w3gH16tWvsrlwqK+7Y8cKSPmX0uphn2obr7zAMA3mNmX1JK+b2L0mdOJ9LtmwC8GMBXX+SEMSNzm3leiZjZMSZAfFMp5aevUhcALwTwdWb2tQDuBvAkM/u3pZS/eUX6PALgkVKKW88/iQkUDzKQfXKf92Z5YCnlfaWUzyilPL1uU/QIgC+8aEDcJmb2ICa37OtKKX90RWrc9qadFyk2zVpvAPBwKeX7r0oPl1LKq0spD9R+8zIA/+UKARG1z37YzJ5dD301Zl7adJB9eh3B2ZYH/nGTHwHwBGzeYfGOUsq3XKYCg81Bt24CeoHyQgDfCOB9tnlvx/fUPfQOMskrAbypTmK/DeDvXLE+ey2HZX4HOchBDkKyT+7zQQ5ykINcuRxA8SAHOchBSA6geJCDHOQgJAdQPMhBDnIQkgMoHuQgBzkIyQEUD3KQgxyE5ACKBznIQQ5C8v8BuoIhBCGh4vUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# stacker\n",
    "s18d_stack, stack_num = s18dStack(mdcwRA, mdcwDec)\n",
    "\n",
    "out, cen, cen_d, out_d = 0,0,0,0\n",
    "\n",
    "#plot stack\n",
    "plot = plt.imshow(s18d_stack, extent = [w1.value,x1.value,y1.value,z1.value])\n",
    "plt.scatter(0,0, marker = '+', color = 'r')\n",
    "plt.colorbar(plot, format='%.0e')\n",
    "plt.title(\"VLASS Stack on 720x720 MaDCoWS Clusters\\n Total in Stack: {}\".format(stack_num))\n",
    "#plt.savefig(\"vlass_stacks_on_mdcws.pdf\")\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# give units back to vmap\n",
    "vmap = s18d_stack*(u.Jy/u.beam)\n",
    "\n",
    "#Convert from Jy/beam to Jy/pixel\n",
    "conv_factor = (np.pi*(2.5/2)**2/np.log(2))*(u.pix/u.beam)\n",
    "vmap = (vmap/conv_factor)\n",
    "\n",
    "#Smoothing map\n",
    "filt_scale = 1*u.arcmin\n",
    "beam_fwhm = 2.5*u.arcsec\n",
    "smooth_scale = np.sqrt(filt_scale**2-beam_fwhm.to(u.arcmin)**2)\n",
    "\n",
    "filt_pix = np.round(smooth_scale/pixel_size.to(u.arcmin))\n",
    "kernel = Gaussian2DKernel(x_stddev =  filt_pix)\n",
    "\n",
    "hdata = convolve(vmap, kernel) * (u.Jy/u.pix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# saving stacked array as text file\n",
    "# np.savetxt('/home/alex/Downloads/stackedSmoothedVLASSdata', hdata.value)\n",
    "# b = np.loadtxt('/home/alex/Downloads/stackedSmoothedVLASSdata')\n",
    "# plt.imshow(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'Quantity' object has no 'get_array' member",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-13-02f59ced6809>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtitle\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"VLASS Stack on MaDCoWS Clusters\\n Total in Stack: {}\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstack_num\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmarker\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'+'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'r'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolorbar\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mimshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextent\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mw1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mx1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0my1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mz1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;31m#plt.savefig(\"/home/alex/Downloads/vlassStackOnMadcows.pdf\",bbox_inches='tight')\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/my_project_dir/my_project_env/lib/python3.8/site-packages/matplotlib/pyplot.py\u001b[0m in \u001b[0;36mcolorbar\u001b[0;34m(mappable, cax, ax, **kw)\u001b[0m\n\u001b[1;32m   2007\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0max\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2008\u001b[0m         \u001b[0max\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgca\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2009\u001b[0;31m     \u001b[0mret\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgcf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolorbar\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmappable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcax\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcax\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkw\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2010\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mret\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2011\u001b[0m \u001b[0mcolorbar\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__doc__\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolorbar\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolorbar_doc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/my_project_dir/my_project_env/lib/python3.8/site-packages/matplotlib/figure.py\u001b[0m in \u001b[0;36mcolorbar\u001b[0;34m(self, mappable, cax, ax, use_gridspec, **kw)\u001b[0m\n\u001b[1;32m   2236\u001b[0m                              'panchor']\n\u001b[1;32m   2237\u001b[0m         \u001b[0mcb_kw\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mkw\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mk\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mNON_COLORBAR_KEYS\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2238\u001b[0;31m         \u001b[0mcb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcbar\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolorbar_factory\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcax\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmappable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mcb_kw\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2239\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2240\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msca\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcurrent_ax\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/my_project_dir/my_project_env/lib/python3.8/site-packages/matplotlib/colorbar.py\u001b[0m in \u001b[0;36mcolorbar_factory\u001b[0;34m(cax, mappable, **kwargs)\u001b[0m\n\u001b[1;32m   1683\u001b[0m         \u001b[0mcb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mColorbarPatch\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcax\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmappable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1684\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1685\u001b[0;31m         \u001b[0mcb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mColorbar\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcax\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmappable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1686\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1687\u001b[0m     \u001b[0mcid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmappable\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcallbacksSM\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconnect\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'changed'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcb\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mon_mappable_changed\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/my_project_dir/my_project_env/lib/python3.8/site-packages/matplotlib/colorbar.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, ax, mappable, **kw)\u001b[0m\n\u001b[1;32m   1203\u001b[0m         \u001b[0;31m# Ensure the given mappable's norm has appropriate vmin and vmax set\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1204\u001b[0m         \u001b[0;31m# even if mappable.draw has not yet been called.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1205\u001b[0;31m         \u001b[0;32mif\u001b[0m \u001b[0mmappable\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_array\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1206\u001b[0m             \u001b[0mmappable\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mautoscale_None\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1207\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/my_project_dir/my_project_env/lib/python3.8/site-packages/astropy/units/quantity.py\u001b[0m in \u001b[0;36m__getattr__\u001b[0;34m(self, attr)\u001b[0m\n\u001b[1;32m    826\u001b[0m         \"\"\"\n\u001b[1;32m    827\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_include_easy_conversion_members\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 828\u001b[0;31m             raise AttributeError(\n\u001b[0m\u001b[1;32m    829\u001b[0m                 \"'{}' object has no '{}' member\".format(\n\u001b[1;32m    830\u001b[0m                     \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__class__\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'Quantity' object has no 'get_array' member"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV4AAAEXCAYAAADlfCy+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df7hVZZ338fdHT2ilCQiZngMhHSPBmZzmoPlUxqSJngp8zBwsE3/lzITzzGR1RekYmTVYM1kNNj2UXGolaFYDTykOqeTUaIjmlOCPg2JxjpgISKEpevw+f9z3kcVmnx9w9lnnB5/Xde3r7HWve631Xfvs/d33vtda91JEYGZm5dmrvwMwM9vTOPGamZXMidfMrGROvGZmJXPiNTMrmROvmVnJnHhtl0k6S9LP+zsOq07SHEnf7e84rHNOvN2QtFTSpVXKp0t6QlKdpKslXdbNeq6W9KKkgyvKh0takNf1R0kPS5pdsZ37JP1B0lOSbpN0aCfbaJD0g1xvi6T7JZ2V542TFJLqduuFGAAkLc/78OaK8h/l8ik9XE9IekbSVkkbJd0q6a+r1Jsq6Y78f9kg6WeSpvVg/bdI+lRhuj5vs1rZ6/L0ZyStzTG1Srq+m218UNLKXH+9pJslvb0n+98TQ+H9MpA58XbvGuAMSaoo/zDwvYh4sbsVSHo18H5gC3BGxewrgP2Aw4EDgGnAmrxcI3At8PE871DgSqC9k019B1gHvB44MMf4++7iG2QeBs7smJB0IHAMsGEX1/PmiNgPmABcDcyT9NnCek8Fvk96/RuAg4BLgPf1YN13AMcWpo8FHqxS1hIRT0iaSfpfHZ9jagJu7Wzlki4Evgp8Mcc1FvgGML0HsZXCCbsbEeFHFw/glaSEeWyhbATwHOnDC+mDe1kX6ziTlBD/Abi/Yt79wMmdLHcqcN8uxLoVOLKTeb8DItfZSkpWbwBuAzYCTwHfA4YXlhkD/JCU1DYC83L5WcDPC/W+DPwcOKDKdvchJYnH8+OrwD553hSglfTF8iSwHji7i/1bTkp+rcDeuewC4N9z2ZRcdhRwJ/B0Xuc8YFhhPQE0VnmtnyN9YSm/Xp/sIpa9gIuB3+bYr+3Yf+Adedt75elvAH9D+hIsln07P58HfLWH/+MD8v/vA13UmQN8t/gaV8x/jJTkO16rlcAfcnxf6ez9ksvPAR4ANgO3AK+veF1nAS3A2vw6XpFfnz8AvwGO6O/P9EB4uMXbjYj4E3ADhVYWcBrwYET8Tw9XMxNYCCwC3iTpLwvz7gK+IOlsSYdVLHdvrn+FpL+StF8327kLuFLSDEljK+Z1tLaGR8R+EXEn6YPxz8AhpBb3GNKHFkl7Az8mJZZxQH2O/2WS9pL0LeDPgRMiYkuVmC4C3gocCbyZ9EG/uDD/daRkUg+cm+Mf0cU+Pg6sBk7I02eSkl5RO/AxYBTpC+Y44KNdrBNgMVCX45tAei1u7KL+WfnxV8B40q+WeXneCtIXTkeXyLHAMtIvmWLZHfn5XcCZkj4pqSm/9p05BtgX+FE3+9NTXwO+FhGvIX0R31CIDwrvF0nTgc8ApwCjgf8iva+LTgaOBiaS/kfHAm8k/Y9PI32BW39n/sHwAN5OasHsm6d/AXysMP9qOmnxkn4GvkRuiZJaCV8rzH8l6c18D/AC6cN5UmH+W0kfhg2kFtnVwH6dbGsEMBdYRUo+9wGT87xxpBZJXRf7eTLwq/y84+f7TvVJCeeXwPXADyi0JqvUfQRoLkxPBR7Lz6cAfypug9Q6emsn61oOnEfqrlkIvAl4OM97ucVbZbl/BH5UmN6pxZvLnwA+BLwt19m3i/26FfhoYXpC/v/VFWL9B2AksC6XzS2UvcSOrcUPAT8FniElp091st0PAU90836dQ89bvHcAnwNGVdTZ6f0C3AycW5jeC3i2Yz9y/XcV5r+L1DX0VnJL34/0cIu3ByLi56Sf4idLegOpVXRdDxf/MPBARNyXp78HfFDSK/K6/xQRX4yIvyT9zL0B+L6kkXn+XRFxWkSMJv2EPZbUiqwW5+aImB0Rk0h9f/cB/1GlfxoASQdJWiSpTdIfgO+SWomQWny/jc77sBtJfYqfi4htXez/IaRWc4ff5rIOGyu28Syp9diVH5I+1BeQ+rV3IOmNkn6cD1j+gdQXOqqyXsUyryC14jaxvVV2cOdLVN2vOtLrDtv7ed9B+qKG1B3TUbYuIl5ePiK+FxHHA8OBvwU+L2lqle1uBEbVsA/1XFKL9EFJd0t6bxd1Xw98TdLTkp4mvVYi/VrpsK7jSUTcRvoVcCXwpKT5kl5To7gHNSfenruW9LP2DOCWiOjpQaszgfE5CTwBfIWUBJorK0ZER5J4NelAWuX8u0lJ54juNhoRTwH/QkoQI0mtkUpfzOV/Fumn5hmkDxKkD9DYLj7gDwBnAzdLmtBFKI+TPrAdxuay3RYRz5JaX39HlcRL6vN9EDgs79dn2L5fnZkOvEjqJniItP/v76J+tf16ke0HM+9g+xflf+WyX5Ba08Vuhh1ExAsR8X3g11T/P98JPE/6ddITzwCv6pjI3RijC9triYjTgdcClwM35oPB1d4v64C/iYjhhccrI+K/i7tQsT9fz42KiaQE/8kexj2kOfH23LXA8cBHSGc6VNpb0r6FxzBJHQewjiL1cR5J+jBdR+4zlvRPkibn+vuSfoo+DTwk6e2SPiLptbnum0hnPdxVLUBJl0s6QukUt/1JiWlNRGwkdRu8ROqP7LA/6cDJFkn17PihWEE6MDVX0qvzPr2tuL2IWEhKaj/NvwSqWQhcLGm0pFGkg2O1OMf0M8A7I+KxKvP2Jx3M2Zpfs7/rbCWSRkr6EKlVdnlEbIz0O/lC4J9y3/trcn/22yXNL+zXxyQdmvvevwhcX2i930lqvZ5BTrwRsZn0fziDQuJVOi/6PZL2z9s5CZhE6s7ZQaR+9EtIfeEnS3qVpFdIOknSl6rs4sPAvnn9ryD1r+9T2PYZkkZHxEuk9x2k90m198s3gU9LmpSXPUDSB7p4bSdLOjpv9xlSV9lLndXfo/R3X8dgepD67TaTj8oXyq8mfdMXHz8nvVF/UGU9R5FaLSNJH4T7SYliU97G/8r1jgD+H6kVtZXUN3c58IpO4vs30hHlraQPzo+BwwvzL83lT5P63SaR+pa3krolPk6hP5DUivsPtp/18PVcfhY7ntXwEfJBuCox7Qt8nZTE1+fnHX3lU+ii/7GT1/+8TuYVz2roOH1rKynpXVoRb5ASwdb8mt8OfLDKOk/My3e8nsuB9+R5e5ES4Lo877vAiIrl78zvl70KZd/I259QKDuF1BrezPaj/2d18178EOlshGdIfdM/Kbxv5pD7eAv/r/Wk/vNPsGMf73dz+VbSsYGTC8vt8H7JZR/O8f0h7/uCite1sTB9HKnlvpXtZ81UPT6xpz2UXyAzMyuJuxrMzErmxGtmVjInXjOzkjnxmpmVzInXaqK70awkrVIPRw/rL0qjn53X33HY0OfEO8RJGpuHDux4FIdE3CrpHZ0sV9MxdyNiUkQs351l1cXQmBogY89KukBpmMbnJV1dZf6rJH1D24fsLJ7HO0fSCxX/p/F53hslLVYalnKT0pCTXV2wYoOAh24b4iLidxQuwZUUpFHV1vRfVD2n7UNjnkIaSW0/0uArnQ2N2V8eBy4jjUXxyirz55M+b4eTzh0+smL+9RFROWQopIswlpCuEvwj6dzhxaRxKmyQcot3D5avPLo2t6Z+K+nifOXU4aSLP47Jra+nc/33SPpVbnmukzRnF7b1mKTj8/M5km7I2/5j7oZo6mTRI4G1EXFrJH+MiB9ExO8knUi6gu2vc5z/k9d/tqQH8roflfQ3FbEUW9CP5PVUxnuwpF9L6tElrhHxw4jouNikcl0dVxyeHxEbIqI9Iu7p4XpXRMRVEbEpIl4gDbM4QWkcYhuknHj3bP9GGq5vPPBO0mXMZ0fEA6SBWu6MNCTg8Fz/mVxnOPAe4O8k9XTMgErTSMNMdrTo5nVSr9OhMSNiKdsv1d0vIjqGXHwSeC/wGlJL8QpJbwGQdBSpBf3JvO1jSVdyvSx3Y/yMNP7wl3PZbEk/3s19PYp0Zd/nclfDbyRVjgPxvtyVsEpSp5c453ifiHQZuA1STrx7KKXBUmYAn86tyMeAfyVdElpVRCyPiN9ExEsR8WvSeAXv3M0Qfh4RN0VEO2mgmzdXqxQRj5IuLa4njdz2lNJtlDodwSwifhIRj+QW8s+A/yQNWANpNK4FEbEs70dbRDxYWHwi6RLiz0bE/MI650ZEVyN3daWBdPn3FtKgRRcA1+RfFuT9Opw0eM1HgEsknV65EkkNpDElLtzNOGyAcOLdc40CXsHOQxvWV68OecCT23PXxBZSq7jL4Ra78ETh+bOkgVyqHnOIXRgaM8d5kqS7cgvyadJIcMXhLh/pIq4PAW10PQj6rvoTaazeyyJiW/4yuJ08mHtErI6Ix3MXxH+TBic/tWKfRpO+QL4RaXAiG8ScePdcT5GSQeXQhm35ebVBPK4jdQuMiYgDSP3A3Q23WFOx89CYO8QpaR/S4Oz/AhyUu0luKsS5jjRiXGfmkF6b69T1nSB2xa+rlHU1SEpQeF2V7sjxn8CSiPhCjWKyfuTEu4fKP/FvIN12aH9Jryf9hO04Nev3QIOkYYXF9gc2RcRzua/0g30dp7ofGvP3wDhJHe/lYaRhDzcALyoNsXhCYZVXAWdLOi4fSKzP6+zwAvAB0pjI1xbW212cdUrDeu7N9iFCO1rwd5DuYfbpXO9tpFsG3ZKXnS5phJKjgP9DOnMBpYHDbwF+ERGzsSHBiXfP9vekA2aPkoaxvA5YkOfdRhom8AlJT+WyjwKXSuo4rekG+t7TpET7G0lbgaWk+411jD37/fx3o6R7I+KPpMR1A2mYxQ+SWulAOkuAfMCN1Of6M3Zs9RPpjhqnkO4msSAn6M9IurmLOC8mdSnMJo23+6dcRj4bYTqpy2ML8C3gzELf8gzSLZ/+SDrwd3lEdIz5/L+ByaQvi+J5vpX31LNBxMNCmpmVzC1eM7OSOfGa2ZAhaYGkJyXd38l8Sfq6pDX5Apm3lB0jOPGa2dByNemWTZ05CTgsP84n3Ri1dE68ZjZkRMQdpLEwOjMduDZfXHMXMFzSweVEt92gHCRn1KhRMW7cuP4Ow8wGoCOOOII1a9bQ1NQUAPfcc89T+eIbSBcIrStUb81l64vrkHQ+qUXMq1/96r9805t6NiZRxbY6NSgT77hx41i5cmV/h2FmA9Bjjz3Ge9/73pdzhKTfdrPITvLl4vMBmpqaoqf5pqfbcleDme1J2kiXjXdoYPvVmqVx4jWzPckS4Mx8dsNbgS0Rsb67hWrNidfMhozTTz+dY445hoceeoiGhgauuuoqgNGS/jZXuYl0peYa0hWEH+2POAdlH6+ZWTULF+48cNt55523ISK+CRDpUt1ZZcdVyS1eM7OSOfGamZXMidfMrGROvGZmJXPiNTMrmROvmVnJnHjNzErmxGtmVjInXjOzkjnxmpmVzInXzKxkTrxmZiVz4jUzK5kTr5lZyZx4zcxK5sRrZlYyJ14zs5I58ZqZlcyJ18ysZE68ZmYlc+I1MytZTRKvpBMlPSRpjaTZVebvI+n6PP+XksZVzB8raaukT9QiHjOzgazXiVfS3sCVwEnAROB0SRMrqp0LbI6IRuAK4PKK+V8Bbu5tLGZmg0EtWrxHAWsi4tGI2AYsAqZX1JkOXJOf3wgcJ0kAkk4G1gKrahCLmdmAV4vEWw+sK0y35rKqdSLiRWALcKCk/YBPAZ/rbiOSzpe0UtLKDRs21CBsM7P+0d8H1+YAV0TE1u4qRsT8iGiKiKbRo0f3fWRmZn2krgbraAPGFKYbclm1Oq2S6oADgI3A0cCpkr4EDAdekvRcRMyrQVxmZgNSLRLv3cBhkg4lJdgZwAcr6iwBZgJ3AqcCt0VEAO/oqCBpDrDVSdfMhrpeJ96IeFHSBcAtwN7AgohYJelSYGVELAGuAr4jaQ2wiZSczcz2SLVo8RIRNwE3VZRdUnj+HPCBbtYxpxaxmJkNdP19cM3MbI/jxGtmVjInXjOzkjnxmpmVzInXzKxkTrxmNmQsXbqUCRMm0NjYyNy5c3ean0dCvF3SryT9WlJzP4RZm9PJzMz6W3t7O7NmzWLZsmU0NDQwefJkpk2bVlntYuCGiPj3PIriTcC4smN14jWzIWHFihU0NjYyfvx4AGbMmMHixYsrqwXwmvz8AODx8iLczl0NZjYktLW1MWbM9mFjGhoaaGurHDaGOcAZklpJrd2/r7auvh4N0YnXzPYkpwNXR0QD0EwaymCnPNjXoyE68ZrZkFBfX8+6dduHBm9tbaW+vnJocM4FbgCIiDuBfYFRZcXYwYnXzIaEyZMn09LSwtq1a9m2bRuLFi2qdnDtd8BxAJIOJyXe0u+s4INrZjYk1NXVMW/ePKZOnUp7ezvnnHMOkyZNAjhE0rQ8UuLHgW9J+hjpQNtZeYjacmMte4NmZn2lubmZ5uadTs19PCddImI18LbSA6vgrgYzs5I58ZqZlcyJ18ysZE68ZmYlc+I1MyuZE6+ZWcmceM3MSubEa2ZWMideM7OSOfGamZXMidfMrGROvGZmJXPiNTMrmROvmVnJnHjNzErmxGtmVrKaJF5JJ0p6SNIaSbOrzN9H0vV5/i8ljcvl75Z0j6Tf5L/vqkU8ZmYDWa8Tr6S9gSuBk4CJwOmSJlZUOxfYHBGNwBXA5bn8KeB9EfFnwEzgO72Nx8xsoKtFi/coYE1EPBoR24BFwPSKOtOBa/LzG4HjJCkifhURj+fyVcArJe1Tg5jMzAasWiTeemBdYbo1l1WtExEvAluAAyvqvB+4NyKer0FMZmYD1oC42aWkSaTuhxO6qHM+cD7A2LFjS4rMzKz2atHibQPGFKYbclnVOpLqgAOAjXm6AfgRcGZEPNLZRiJifkQ0RUTT6NGjaxC2mVn/qEXivRs4TNKhkoYBM4AlFXWWkA6eAZwK3BYRIWk48BNgdkT8ogaxmJkNeL1OvLnP9gLgFuAB4IaIWCXpUknTcrWrgAMlrQEuBDpOObsAaAQukXRffry2tzGZmQ1kNenjjYibgJsqyi4pPH8O+ECV5S4DLqtFDGZmg4WvXDMzK5kTr5lZyZx4zcxK5sRrZlYyJ14zs5I58ZqZlcyJ18ysZE68ZmYlc+I1MyuZE6+ZDRlLly5lwoQJNDY2Mnfu3Kp1JJ0mabWkVZKuKzlEYIAMC2lm1lvt7e3MmjWLZcuW0dDQwOTJk5k2bdoOdSQdBnwaeFtEbO6vsWGceM1sSFixYgWNjY2MHz8egBkzZrB48eLKah8BroyIzQAR8WS5USbuajCzIaGtrY0xY7YPDd7Q0EBbW+XQ4LwReKOkX0i6S9KJ1dYl6XxJKyWt3LBhQ81jdeI1sz1JHXAYMAU4HfhWHhd8B3194wUnXjMbEurr61m3bvvtH1tbW6mvr7z9I63Akoh4ISLWAg+TEnGpnHjNbEiYPHkyLS0trF27lm3btrFo0aKdDq4B/0Fq7SJpFKnr4dFyI/XBNTMbIurq6pg3bx5Tp06lvb2dc845h0mTJgEcImlaRCwh3SnnBEmrgXbgkxGxsfRYy96gmVlfaW5uprm5ubL48Zx0iYgg3X7swrJjK3JXg5lZyZx4zcxK5sRrZlYyJ14zs5I58ZqZlcyJ18ysZE68ZmYlc+I1MyuZE6+ZWcmceM3MSubEa2ZWMideM7OSOfGamZWsJolX0omSHpK0RtLsKvP3kXR9nv9LSeMK8z6dyx+SNLUW8dgeaMqU9DAbBHqdeCXtDVwJnARMBE6XNLGi2rnA5ohoBK4ALs/LTgRmAJOAE4Fv5PWZmQ1ZtRiP9yhgTUQ8CiBpETAdWF2oMx2Yk5/fCMyTpFy+KCKeB9ZKWpPXd2cN4rI9QUcr92c/23F6+fJ+CMasZ2rR1VAPrCtMt+ayqnUi4kVgC3BgD5cF+v6un2ZmZRk0d6CIiPnAfICmpqbo53BsoOho2bqla4NILVq8bcCYwnRDLqtaR1IdcACwsYfLmpkNKbVIvHcDh0k6VNIw0sGyJRV1lgAz8/NTgdvyvY+WADPyWQ+Hkm6zvKIGMdmeZvlyt3Zt0Oh1V0NEvCjpAtLdO/cGFkTEKkmXAivzTeauAr6TD55tIiVncr0bSAfiXgRmRUR7b2MyMxvIatLHGxE3ATdVlF1SeP4c8IFOlv0C8IVaxGFmNhj4yjUzs5I58ZqZlcyJ18ysZE68ZmYlc+I1MyuZE6+ZWcmceM3MSubEa2ZWMideM7OSOfGamZXMidfMrGROvGZmJXPiNTMrmROvmQ0ZS5cuZcKECTQ2NjJ37txO60l6v6SQ1FRieC8bNLf+MTPrSnt7O7NmzWLZsmU0NDQwefJkpk2btlM9SfsD/wD8svQgM7d4zWxIWLFiBY2NjYwfP55hw4YxY8YMFi9eXK3q54HLgefKjXA7J14zGxLa2toYM2b7LRwbGhpoa9vxFo6S3gKMiYifdLWuvr6ruROvme0RJO0FfAX4eHd1I2J+RDRFRNPo0aNrHosTr5kNCfX19axbt+7l6dbWVurr64tV9geOAJZLegx4K7CkPw6wOfGa2ZAwefJkWlpaWLt2Ldu2bWPRokU7HFyLiC0RMSoixkXEOOAuYFpErCw7VideMxsS6urqmDdvHlOnTuXwww/ntNNOY9KkSQCHSNr59IZ+5NPJzGzIaG5uprm5ubL48YhYUlkYEVNKCaoKt3jNzErmxGtmVjInXjOzkjnxmpmVzInXzKxkTrxmZiVz4jUzK5kTr5lZyXqVeCWNlLRMUkv+O6KTejNznRZJM3PZqyT9RNKDklZJ6nzUYjOzIaS3Ld7ZwK0RcRhwa57egaSRwGeBo4GjgM8WEvS/RMSbgL8A3ibppF7GY2Y24PU28U4HrsnPrwFOrlJnKrAsIjZFxGZgGXBiRDwbEbcDRMQ24F6goZfxmJkNeL1NvAdFxPr8/AngoCp16oF1henWXPYyScOB95FazWZmQ1q3g+RI+inwuiqzLipORERIil0NQFIdsBD4ekQ82kW984HzAcaOHburmzEzGzC6TbwRcXxn8yT9XtLBEbFe0sHAk1WqtQFTCtMNwPLC9HygJSK+2k0c83NdmpqadjnBm5kNFL3talgCzMzPZwLV7ix3C3CCpBH5oNoJuQxJlwEHAP/YyzjMzAaN3ibeucC7JbUAx+dpJDVJ+jZARGwi3dXz7vy4NCI2SWogdVdMBO6VdJ+k83oZj5nZgNergdAjYiNwXJXylcB5hekFwIKKOq2AerN9M7PByFeumZmVzInXzKxkTrxmZiVz4jUzK5kTr5lZyZx4zcxK5sRrZlYyJ14zs5I58ZqZlcyJ18ysZE68ZmYlc+I1MyuZE6+ZWcmceM3MSubEa2ZWMideM7OSOfGamZXMidfMhoylS5cyYcIEGhsbmTt37k7zJV0oabWkX0u6VdLr+yFMJ14zGxra29uZNWsWN998M6tXr2bhwoWsXr26stqvgKaI+HPgRuBLpQeKE6+ZDRErVqygsbGR8ePHM2zYMGbMmMHixTve+Dwibo+IZ/PkXUBD6YHixGtmQ0RbWxtjxox5ebqhoYG2trauFjkXuLnaDEnnS1opaeWGDRtqGyhOvGa2B5J0BtAEfLna/IiYHxFNEdE0evTomm+/V7d3NzMbKOrr61m3bt3L062trdTX1+9UT9LxwEXAOyPi+fIi3M4tXjMbEiZPnkxLSwtr165l27ZtLFq0iGnTpu1QR9JfAP8XmBYRT/ZLoLjFa2ZDRF1dHfPmzWPq1Km0t7dzzjnnMGnSJIBDJE2LiCWkroX9gO9LAvhdREzrYrV9E2vZGzQz6yvNzc00NzdXFj+eky4RcXz5Ue3MXQ1mZiVz4jUzK5kTr5lZyZx4zcxK5sRrZlayXiVeSSMlLZPUkv+O6KTezFynRdLMKvOXSLq/N7GYmQ0WvW3xzgZujYjDgFvz9A4kjQQ+CxwNHAV8tpigJZ0CbO1lHGZmg0ZvE+904Jr8/Brg5Cp1pgLLImJTRGwGlgEnAkjaD7gQuKyXcZiZDRq9TbwHRcT6/PwJ4KAqdeqBdYXp1lwG8HngX4FnKxeq1NejBZmZlaXbK9ck/RR4XZVZFxUnIiIkRU83LOlI4A0R8TFJ47qrHxHzgfkATU1NPd6OmdlA023i7eoSO0m/l3RwRKyXdDBQbdCJNmBKYboBWA4cAzRJeizH8VpJyyNiCmZmQ1hvuxqWAB1nKcwEFlepcwtwgqQR+aDaCcAtEfHvEXFIRIwD3g487KRrZnuC3ibeucC7JbUAx+dpJDVJ+jZARGwi9eXenR+X5jIzsz1Sr0Yni4iNwHFVylcC5xWmFwALuljPY8ARvYnFzGyw8JVrZmYlc+I1MyuZE6+ZWcmceM3MSubEa2ZWMideM7OSOfGamZXMidfMrGROvGZmJXPiNTMrmROvmVnJnHjNzErmxGtmVjInXjOzkjnxmpmVzInXzKxkTrxmZiVz4jUzK5kTr5lZyZx4zcxK5sRrZlYyJ14zGzKWLl3KhAkTaGxsZO7cuTvNl7SPpOslrZH0S0njSg8SJ14zGyLa29uZNWsWN998M6tXr2bhwoWsXr26stq5wOaIaASuAC4vPVCceM1siFixYgWNjY2MHz+eYcOGMWPGDBYvXlxZbTpwTX5+I3CcJJUaKFBX9gZr4Z577nlK0m/7O46CUcBT/R1EL3kf+t9gjx/6dx9GAK8p5IaRwH7A6wt16oF1ABHxoqQtwIFUxCzpfOD8PPm8pPt7GMOEnlQalIk3Ikb3dwxFklZGRFN/x9Eb3of+N9jjh/7dB0mnAidGxHl5+sPA0RFxwa6uKyLmA/Pzenq8T5JW9qSeuxrMbKhoA8YUphtyWdU6kuqAA4CNpURX4MRrZkPF3cBhkg6VNAyYASypqLMEmJmfnwrcFhFRYozAIO1qGIDm93cANeB96H+DPX7ox33IfbYXALcAewMLImKVpEuBlRGxBLgK+I6kNcAmUnLuzq7sU4/qqh+SvZnZHs1dDWZmJXPiNTMrmRNvD//uX9cAAAOUSURBVEkaKWmZpJb8d0Qn9WbmOi2SZlaZv2QXzgmsqd7sg6RXSfqJpAclrZK08/WYfRf3iZIeypd5zq4yv9PLQCV9Opc/JGlqWTFXiXG39kHSuyXdI+k3+e+7yo69EONu/x/y/LGStkr6RFkx91Z3+1yot0DSkz3+bEeEHz14AF8CZufns4HLq9QZCTya/47Iz0cU5p8CXAfcP9j2AXgV8Fe5zjDgv4CTSoh5b+ARYHze7v8AEyvqfBT4Zn4+A7g+P5+Y6+8DHJrXs3c/vO692Ye/AA7Jz48A2vrpvbPb+1CYfyPwfeAT/bEPfbHPhbrHAm/p6WfbLd6eK15qeA1wcpU6U4FlEbEpIjYDy4ATASTtB1wIXFZCrJ3Z7X2IiGcj4naAiNgG3Es6T7KvHQWsiYhH83YXkfajqLPLQKcDiyLi+YhYC6zJ6yvbbu9DRPwqIh7P5auAV0rap5Sod9Sb/wOSTgbWkvZhsOjJPgMQEXeQzpLoESfenjsoItbn508AB1Wp8/LliFlrLgP4PPCvwLN9FmH3ersPAEgaDrwPuLUvgtzVeKi4DBTouAy0J8uWoTf7UPR+4N6IeL6P4uzKbu9DbnR8CvhcCXHWUp+9f3web4GknwKvqzLrouJERISkHp+HJ+lI4A0R8bG+Hoaur/ahsP46YCHw9Yh4dPeitF0laRJpJK0T+juW3TAHuCIitvbDeDQDkhNvQUQc39k8Sb+XdHBErJd0MPBklWptwJTCdAOwHDgGaJL0GOk1f62k5RExhRrrw33oMB9oiYiv1iDcntiVy0BbKy4D7cmyZejNPiCpAfgRcGZEPNL34VbVm304GjhV0peA4cBLkp6LiHl9H3av9N37p787sAfLA/gyOx6Y+lKVOiNJ/Vgj8mMtMLKizjj67+Bar/aB1D/9A2CvEmOuIx3gO5TtBzgmVdSZxY4HdW7Izyex48G1R+mfg2u92Yfhuf4p/fGeqcU+VNSZw+A5uNbtPlfU7/Fnu993brA8SP1ttwItwE8LyagJ+Hah3jmkgzhrgLN7888ZSPtA+rYP4AHgvvw4r6S4m4GHSUeYL8pllwLT8vN9SUfL1wArgPGFZS/Kyz1ECWdh1HofgIuBZwqv+X3AawfTPlSsY9Ak3s72uZN6C4H1wAukvuBzu1qvLxk2MyuZz2owMyuZE6+ZWcmceM3MSubEa2ZWMideM7OSOfGamZXMidfMrGT/Hzv4fz0xDp0NAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot smoothed stack\n",
    "plt.title(\"VLASS Stack on MaDCoWS Clusters\\n Total in Stack: {}\".format(stack_num))\n",
    "plt.scatter(0,0, marker = '+', color = 'r')\n",
    "plt.colorbar(hdata)\n",
    "plt.imshow(hdata, extent = [w1.value,x1.value,y1.value,z1.value])\n",
    "#plt.savefig(\"/home/alex/Downloads/vlassStackOnMadcows.pdf\",bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "emissionStack = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fill emission stack with flux values for histogram\n",
    "# emissionStack = []\n",
    "i = 0\n",
    "while i < len(stackArr):\n",
    "    try:\n",
    "        cent_em = central_emission(stackArr[i], int(100), pixel_size.to(u.arcmin).value, 1.0)\n",
    "        emissionStack.append(np.sum(cent_em))\n",
    "    except:\n",
    "        pass\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1378"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(emissionStack)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "emissionStackArr = np.array(emissionStack)\n",
    "a_file = open(\"/home/alex/Downloads/test.txt\",\"w\")\n",
    "np.savetxt(a_file, emissionStackArr)\n",
    "\n",
    "a_file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fill values for histStack ... cut out values below VLASS noise floor\n",
    "histStack = []\n",
    "i = 0\n",
    "while i < len(emissionStack):\n",
    "    if(emissionStack[i] > 6.9E-5):#combined (3 epoch) sensitivity 69microJy RMS\n",
    "        histStack.append(emissionStack[i])\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'histStack' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-18-40e2635d070d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhistStack\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'histStack' is not defined"
     ]
    }
   ],
   "source": [
    "len(histStack)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAElCAYAAAALP/6mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debwcVZ338c+XAAlrEBIJIcBFNgUGl4ngOOr4KCqLERccFhFRFHnmQWV8xA0c0QdHZMYRBH0pI8o6sipDFAbGBRAGlQRBEhBlSUgAQwgSSFgTfs8f53SlbtO3b93crtu37/2+X69+3e5T26+r69avqk7VOYoIzMzMANbpdgBmZjZ6OCmYmVnBScHMzApOCmZmVnBSMDOzgpOCmZkVnBTGCUlvlLR4pKaz0UPSAkl7dzuO0UbStpJWSJrQ7VhGEyeFDpJ0qKQ5eUN7SNJVkl7XgfmeKOn8TsTYZhkhaWWOfYWkx+pcXoV4rpX04TbD+3LMv2sqnyLpWUkLKi7nCEmrS9/7Pkk/kLRz03jr59/hT3k9LZD0fUl9g8x/qxznlqWy4wco+6/8foakyyQ9Imm5pHmSjmizjE0lnSrp/vwd7smfp1RZB1XUsQ2OxHbdtLx+yTEi7o+IjSNi9UjF0AucFDpE0ieBU4F/BrYEtgW+DRwwAsuWpE78li/P/yQbR8RmHZjfSNhQ0u6lz4cC9w1xHjdFxMbAZGBv4ClgbtN8LwXekec/GXg5MBd4c7sZR8RDwN3AG0rFbwD+0KLs+vz+PGARsB2wBfB+YEmr+UtaH/g5sBuwD7Ap8DfAMmDPdrGNJEnrdjsGqygi/Brmi7STWAG8t8046wCfBe4h/cNeDGyeh/UBAXwAuB94BDg+D9sHeBZ4Li/jtlx+LfAV4EbSTmxH4IPAncATwL3AR0vLfyOwuE18AezYorzfdM3jAWcDJ+X3nwF+A6ybP/9vYD4wqcV8XwT8BFgK/CW/n5GHfQVYDTydv/MZLaZvrLMTgH8plc8BjgcWlMoa6/0J4A7gXaVhRwA3tJj/T4BL8/tGotimzfqbDlwBPEpKAh8pDTsLOD2/n5C/89FNZY8Dr8ufVwCvqLjtfZiUMDZuM84CYO/m32uA3/czwAN5Xd1FSnoDbYOT83d7KE9zEjChtF5vBL5B2t5PahHXicD5A8T8jrztPEba1l9WGrYN8KO8Hpc1tg9gB+AXuewR4AJgszzsPOD5/DuuAD5d2oYa22u73/BE0v/suXndzAdmtltv3dofDffV9QDGwiv/06xqbFwDjPMJ4NfADGAi8F3gh3lYY+P8d2AD0lHoM41/hFb/PPkf5X7SEeK6wHrA/vkfQ8DfAU8Cr8rj9/vnbxFfJ5LCOqSj3ROBnUg7+1cOsLwtgPcAGwKbAJcAlzd9vw+3ibexzvpIR9UTgF1JR+B70z8pvDf/w68DHASsBLbKw46gdVL4ELAkvz8ZuG6QbeB60pnhJOAVpB3Wm/KwD7BmRzozj7tTU9lTwPr5889IO9SDgW0HWe6FwDmDjLOACkkB2CWvy+mldbxDm23wx6TteCPgxcBvyQcieb2uAj5G2j43aBHXC+aZy3fOv9FbSNv1p0k76fXz73wbKdlslNd3I5numKeZCEzN6/nUVuuhaRtqJIV2v+GJpIOU/XIMXwV+Pdh668WXLx91xhbAIxGxqs04R5OO/hdHxDOkjezAptPqL0XEUxFxG2nDf/kgyz07IuZHxKqIeC4ifhoR90RyHXAN8PohfI9bJD2WX98cwnQARMTzwOHAx0lHXKdExO8GGHdZRFwWEU9GxBOks4O/G+oygcWkI7O987LPa7GsSyLiwYh4PiIuAv7E4JdWHgQ2z++3IB0NtyRpG+Bvgc9ExNMRcSvwvRwPwHXA7pI2I/0ev4qIPwFTS2W/john8/jvBX4FfAG4T9Ktkl49wOLbxjZEq0k71F0lrRcRCyLinlYj5vqQ/YBjI2JlRDxM2lEfXBrtwYg4PW+fTw0hjoOAn0bEf0fEc8C/kg6WXkv63aYDx+XlPh0RNwBExN15mmciYinwb1Tcpir8hpAOHq6MVAdxHmv+Pyuvt17gpNAZy4Apg1w33Q74cWOnS7rMs5pU/9Dw59L7J4GNB1nuovIHSftK+rWkR/My9gOGUtn4qojYLL8+PoTpChGxAPgl6WjpWwONJ2lDSd+VtFDS46SjtM3W8k6Qc0lHpofQIilIOjzvWBvrfncGXy9bky4jQPp9t2oz7nTg0ZzcGhbmeTTWyQOknf8bSDt8gP8plTXqE4iIv0TEZyNiN9L2cStwuSS1WPZgsVUWEXcDx5IOWB6WdKGk6QOMvh3pKP6h0nr9LumMoWFRyykHN520/hpxPZ/ntTXp0tHCVgdgkrbMMT+Qt6nzqb79t/0Ns+b/z0mS1h3iehv1nBQ64ybS5Z53thlnEbBvaae7WURMiogHKsx/oKZsi3JJE4HLSEdVW0aqKL6SdCmpk54kXfJpmFYeKGl/UkXnz4F/aTOf/0s67d4rIjZlTaVrI96hNN97GenS2b0RcX9TPNuRLssdA2yR18s8Bl8v72LNzvtnwJ6SZgww7oPA5pI2KZVtS0oEDdeTvuPfkJIBef5vAF5HKSmURcQjpN90OmvOXMp+BrxN0kaDfJ+GlbT5/SLiPyLidaSdfgBfawxqms8i0jY/pbQ9b5oTGQNMU9WDeflAupGClAweyMvddoADsH/Oy/yrvE0dRv/fuV08VX7DAbVZbz3HSaEDImI58E/AtyS9Mx8Fr5eP3E/Jo30H+EreSSFpqqSqdyYtAfoGucNofdIp7FJglaR9gbeu1Rdq71bgUEkTJO1D6fQ83wL5PVLl5weAWZL2G2A+m5Cuoz8maXPgi03DlwAvqRJQRKwE3pSX22wj0j/p0hzjB0lnCi+Qv9P2kk4nXWv/Up7/z4D/Jp3p/bWkdSVtIuloSR+KiEWkHf1XJU2StAdwJOlIteF60qWIByPi8Vx2Qy6bTDqwaMTxNUm7N5ZDqrC/OyKWtQi7cafSZZJeKmkdSVtI+vwA6/5WYD9Jm0uaRjrCbSx3F0lvygcYT5N+n+fz4H7bYKS7qq4Bvp5viV1H0g6ShnoJcJ28zhqviaQK3f0lvVnSeqQDiGdI6/i3pMtlJ0vaKE/zt3lem5AqkZdL2ho4rmlZA25TFX/DlgZZbz3HSaFDIuLrwCdJd8MsJf2jHgNcnkc5jXSd/RpJT5AqnfeqOPtL8t9lkm4ZYPlPkK7lX0yq4D00L6/TPgHMIt0V8j7WfD+AM4H/zNddl5H+qb4naYsW8zmVdJ34EdK6+K+m4aeR6lz+UqV+IyLmtLqOGxF3AF8n7XSXAH9FqsQt+xtJK0h3AF1Luq3z1RFxe2mcA0lnXhcBy0lnGzNJR+qQLl31kY44fwx8MSeThutIl1ZuKJXdSloHcyPiyVL5hnkej5HuItuOdDdOq+/9DKk+5Q+kxPU4acc5hXQnWLPzSPVVC0g79YtKwyaSKtUfIV0qeTHwuTys1TZ4OOlg5A7SNncpQ7+UdQhpJ9p43RMRd5GO8k/PscwCZkXEs/l6/ixSpfL9pDqlg/K8vgS8ivT7/JR0h1LZV4ET8uWuTw0QSx8D/4YDabfeeo4i3MmOmZklPlMwM7OCk4KZmRWcFMzMrOCkYGZmBScFG/U0SIupNS97A0mzlVorvWTwKfpNG5J2rCs2szo4KdiQ5Z3d7eXnJiSdJOnsLoZVlwNJTxVvERHvHemFSzpb0kkjvdy87PUlXarU5HRIemPT8OOUmvV+QqnJ8eOahv9S0lJJj0u6rfm5HKWm5hcqNUV+eX5exbrMScHW1nT6t3Mz6ikZ6ja/HfDHQdq1GrUGePJ3KG4gPTPw5xbDRHpW4UWkRiGPkVTeJj5BanhwU+Ao4HxJW+W4diM1i/F+UtJ9ktQYnXWZk4KtrVOAL7Xa6ahFb20qdXCi1LnKJZLOz0eZt0vaWdLnJD0saZGk5qexd5D023zU+Z/lo0pJr5H0P/mhpNvKR7T50tNXJN1I2vG84IlWSS/L4z0mab6kd+TyL5GeVD9IqfOaI1tMOyE/PXxP/i5zlRpXax6v3yUwpc59bsjvJekb+bs/ntfH7pKOIj0g+Om8/Nl5/OlKnfAszUfoHy/N98R8dH++Uvs/R0jaU6nzp8clLZH0b83xtZIfFjs1Nzj3go5oIuKUiLglN3h3F/CfpEblGsN/X0qmQWorqbFu3gfMjojrI2IFqfG/d6t/MxPWBU4KtrZ+RHp69oi1nH4W6enaFwG/A64mbY9bA18mHUWWHU5qznorUpPM3wRQas7gp6S2/DcHPkVq8mFqadr3k45UN6HU0Fqefj1gNunp3heTmnq+QNIuEfFFUns6F0XqeOisFt/jk6QnYfcjPQn9IVLyGYq3ktpA2pnU5MXfA8si4kxSnwCn5OXPymc6s0lPJW9N6u/gWElvK83vANLTxZvl6U8DTstH7DuQnnpvfP/fSzp0iPG+gCSRGveb31T+E0lPk56uvpbU3wWkJt9va4yXn0Z/Nq8D6yInBVtbQTq6+4JS719D9auIuDofSV5Cav/+5NxU8oWkdnbKvb+dFxHzcjtHXwD+XqlF1cOAK3PTGs9HxH+Tdjzldn/6NTHeFMdrSK3RnpyPjH9B6mDnkIrf48PACRFxVyS3DdBGUTvPkRLWS0mtDNyZ2xZq5dXA1Ij4co73XlKDf+XLNjdFxOV5fTyV57+jpCkRsSIift0YMSL2iIj/GGK8rZxI2p/8oFwYEW/P320/4Jrc4imkdb68aR7L87jWRU4KttYi4kpS2zMfXYvJy91LPkXqj2J16TP0bzq83AzzQtKliCmka/7v1Zp+IB4jtTq61QDTNpsOLCrtrBrz33qA8ZttQ+rVba3lRHQGqanxhyWdKWnTAUbfDpje9H0/T/8m2Ju/75GkI/A/SLpZ0tuHE28zSceQzuT2z20x9ROpr4+rgLc2Ls2RGq5r/o6bknousy5yUrDhOp60Uyo3x9yveeZ8RD+V4Slfp9+WdPT7CGkHeF5Tk+QbRcTJpfEHazJ5m6YK6MpNJufl71BhvMGarP5mRPw1qfe4nVnTwmerJqvva/q+m0RE+cyo3zQR8aeIOIR0eexrwKWq3tR2W5I+ROru9M0RsXiQ0ddlzbqaT6kTKUkvITUs98dOxGVrz0nBhiUiriW1GPqBUvEfSR2Q7J+v2Z9A+ocfjsMk7SppQ1Kdw6X5zOJ8UhPdb8uVvpNyRfdAfR80+w2pDuDTSs2dv5FU33Fhxem/B/w/STvlCuM91LpV2FtJFakbKj27UFRaS3q1pL3yulpJan653GR1uXL8t8ATkj6j9AzFhFwpPVDPbEg6TNLUfDb0WC6u1LSzpImSJuWP6+f1qzzsfaQ6l7fky1jl6V6q1HT8Bnm9HkaqN7kuj3IB6Xd7fU5QXwZ+FP07ubEucFKwTjiBUgcwuX+JfyDtMB8g7egGO4oczHmk/oX/TOpD9+N5WYtIFaufZ02T5cdRcduO1AXmLGBf0pnHt4HDI+IPFeP6N1LF7TWkivezSM1hN/sGqSJ1CXAOaafYsCmpXuAvpEtXy1jTQdFZpG4eH5N0eU6Ebyf1IXxfjvl7pArqgewDzFdqHvw04OBc10C+2+p9baa9i3Q5b2vSzQBPsaYDnJNI3YHenO+OWiHpO3mYyD2RkX6XTwAHRcQtABExn9RF7QV5nE1I24x1mZvONjOzgs8UzMys4KRgZmYFJwUzMys4KZiZWWG4jWV11ZQpU6Kvr6/bYZiZ9ZS5c+c+EhEtnx3q6aTQ19fHnDlzBh/RzMwKkhYONMyXj8zMrNCTSUHSLElnLl/e3J6WmZkNR08mhYiYHRFHTZ7c7iFOMzMbqp5MCmZmVg8nBTMzKzgpmJlZwUnBzMwKPZkUfPeRmVk9ejIpdOLuo75p05DU79U3bdrgE5qZjWE9/UTzcCxcsuQF/RxqyZKW45qZjRc9eaZgZmb1cFIwM7OCk4KZmRWcFMzMrNCTScG3pJqZ1aMnk4IbxDMzq0dPJgUzM6uHk4KZmRWcFEomgp9yNrNxbdw+0dzKM+CnnM1sXPOZgpmZFZwUzMys4KRgZmaFnkwKfnjNzKwePZkU/PCamVk9ejIpmJlZPZwUzMys4KRgZmYFJwUzMys4KZiZWcFJwczMCk4KZmZWcFIwM7OCk4KZmRWcFMzMrOCkYGZmhZ5MCm4Qz8ysHj2ZFNwgnplZPXoyKYwk99tsZuOJ+2gehPttNrPxxGcKZmZWcFIwM7OCk4KZmRWcFMzMrOCkYGZmBScFMzMrOCmYmVnBScHMzApOCmZmVnBSMDOzgpOCmZkVnBTMzKzgpGBmZgUnBTMzKzgpmJlZYVT1pyDpncD+wKbAWRFxTZdDMjMbV2o/U5D0fUkPS5rXVL6PpLsk3S3pswARcXlEfAQ4Gjio7tjMzKy/kbh8dDawT7lA0gTgW8C+wK7AIZJ2LY1yQh5uZmYjqPakEBHXA482Fe8J3B0R90bEs8CFwAFKvgZcFRG3tJqfpKMkzZE0Z+nSpfUGb2Y2znSronlrYFHp8+Jc9jFgb+BASUe3mjAizoyImRExc+rUqfVHamY2joyqiuaI+CbwzW7HYWY2XnXrTOEBYJvS5xm5zMzMuqhbSeFmYCdJ20taHzgYuKLqxJJmSTpz+fLltQVoZjYejcQtqT8EbgJ2kbRY0pERsQo4BrgauBO4OCLmV51nRMyOiKMmT55cT9BmZuNU7XUKEXHIAOVXAlfWvXwzM6uuJ5u56Pblo4kphhe8+qZN60o8Zmad0pNJoduXj54BosVr4ZIlXYnHzKxTejIpmJlZPZwUzMys4KRgZmaFnkwK3a5oNjMbq3oyKXS7otnMbKzqyaRgZmb1cFIwM7OCk4KZmRUGTQqSfl6lbCS5otnMrB4DJgVJkyRtDkyR9CJJm+dXH6lDnK5xRbOZWT3aNYj3UeBYYDowF1Aufxw4o+a4zMysCwZMChFxGnCapI9FxOkjGJOZmXXJoE1nR8Tpkl4L9JXHj4hza4zLzMy6YNCkIOk8YAfgVmB1Lg6ga0lB0ixg1o477titEMzMxqQqnezMBHaNiKg7mKoiYjYwe+bMmR/pdixmZmNJlecU5gHuPaaCVp3vuOMdM+slVc4UpgB3SPotqX8ZACLiHbVF1aMane+UyR3vmFkPqZIUTqw7CDMzGx2q3H103UgEYmZm3Vfl7qMnWHNVZH1gPWBlRGxaZ2BmZjbyqpwpbNJ4L0nAAcBr6gzKzMy6Y0itpEZyOfC2muKpxA3imZnVo8rlo3eXPq5Dem7h6doiqsDPKZiZ1aPK3UezSu9XAQtIl5DMzGyMqVKn8MGRCMTMzLqvSic7MyT9WNLD+XWZpBkjEZyZmY2sKhXNPwCuIPWrMB2YncvMzGyMqZIUpkbEDyJiVX6dDUytOS4zM+uCKklhmaTDJE3Ir8OAZXUHZmZmI69KUvgQ8PfAn4GHgAMBVz6bmY1BVe4+WgiMqhZR3cmOmVk9BjxTkPQvkj7aovyjkk6uN6z2ImJ2RBw1efLkboZhZjbmtLt89CbgzBbl/w68vZ5wzMysm9olhYmtuuCMiOcB1ReSmZl1S7uk8JSknZoLc9lT9YVkZmbd0q6i+Z+AqySdBMzNZTOBzwHH1h2YmZmNvAGTQkRcJemdwHHAx3LxPOA9EXH7SARnZmYjq+0tqRExD/jACMViZmZdNqROdszMbGxzUjAzs0LbpJDbOvrHkQrGzMy6q21SiIjVwCEjFIuZmXVZle44b5R0BnARsLJRGBG31BaVmZl1RZWk8Ir898ulsiA1g9EVbhDPzKweatGSRc+YOXNmzJkzZ62mlUTzNxesddlA5ZOAZ5rKtttySxb8+c9DiNbMrHMkzY2Ima2GVemjeUtJZ0m6Kn/eVdKRnQ5yrHqGlCjKr4VLlnQ1JjOzgVS5JfVs4GpS/8wAf8TNXJiZjUlVksKUiLgYeB4gIlYBq2uNyszMuqJKUlgpaQvy5XJJrwGW1xqVmZl1RZW7jz4JXAHsIOlGYCrw3lqjMjOzrqiSFOYDfwfsQrrB5i7cPIaZ2ZhUZed+U0Ssioj5ETEvIp4Dbqo7MDMzG3kDnilImgZsDWwg6ZWs6YJzU2DDEYhtzJpIek6ibMN11uHJ55/vV+bnGcxspLW7fPQ24AhgBvB11iSFJ4DP1xvW2NZ4dqFMzz//wjI/z2BmI6xdz2vnAOdIek9EXDaCMZmZWZdUqVOYIWlTJd+TdIukt9YemZmZjbgqSeFDEfE48FZgC+D9wMm1RmVmZl1RJSk06hL2A86NiPmlMjMzG0OqJIW5kq4hJYWrJW1CbvLCzMzGlioPrx1J6lPh3oh4Mjd58cF6wzIzs26okhRel//u0XxvvZmZjS1VksJxpfeTgD2BuXSx5zUzM6vHoEkhImaVP0vaBji104FIeglwPDA5Ig7s9PzNzGxwa9Ow3WLgZVVGlPR9SQ9LmtdUvo+kuyTdLemzABFxb0S4Rzczsy4a9ExB0umsaZVhHVKl8y0V5382cAZwbml+E4BvAW8hJZibJV0REXdUD9vMzOpQpU5hTun9KuCHEXFjlZlHxPWS+pqK9wTujoh7ASRdCBwAVEoKko4CjgLYdtttq0xiZmYVValTOKfDy9waWFT6vBjYK9/q+hXglZI+FxFfHSCeM4EzAWbOnNnchpyZmQ1Du6azb+eFjXlCepo5ImKPTgYSEcuAozs5TzMzG5p2Zwpvr2mZDwDblD7PyGWVSZoFzNpxxx07Gdeo06rfBfexYGZ1anf30XrAjIhYWH6RduJV6iIGcjOwk6TtJa0PHEzqA7qyiJgdEUdNnjx5GGGMfo1+F8qvhe5jwcxq1C4pnAo83qL8cSo+pyDph6SuO3eRtFjSkRGxCjgGuBq4E7g4N7JnZmZd1u6If8uIuL25MCJub3FHUUsRccgA5VcCV1aZh5mZjZx2ZwqbtRm2QacDGQpJsySduXz58m6GYWY25rRLCnMkfaS5UNKHSW0fdc14qVMwMxtp7S4fHQv8WNL7WJMEZgLrA++qOzAzMxt5AyaFiFgCvFbS/wJ2z8U/jYhfjEhkZmY24qo80fxL4JcjEEtl4+U5BTOzkbY2raR2nesUzMzq0ZNJwczM6uGkYGZmBScFMzMr9GRS8MNrZmb16Mmk4IpmM7N69GRSMDOzejgpmJlZwUnBzMwKTgpmZlboyaQwnu8+anTRWX71TZvW7bDMbIzoyaQwnu8+chedZlannkwKZmZWDycFMzMrOCmYmVnBScHMzApOCmZmVujJpDCeb0mtqm/aNN+6amZD1pNJYTzfklrVwiVLfOuqmQ1ZTyYFMzOrh5OCmZkVnBTMzKzgpGBmZgUnBTMzKzgpmJlZwUnBzMwKTgpmZlboyaTgJ5pHhzqemvaT2GbdpYjodgxrbebMmTFnzpy1mlYSzd9csNZlw51+2GVNv+OA36+Dv3cdyxiJuM3GO0lzI2Jmq2E9eaZgZmb1cFIwM7OCk4KZmRWcFMzMrOCkYGZmBScFMzMrOCmYmVnBScHMzApOCmZmVnBSMDOzgpOCmZkVejIpuEG8/ibCCxqR67RWDdWZ2djTk0khImZHxFGTJ0/udiijwjOkRvLKr05buGRJ7csws+7ryaRgZmb1cFIwM7OCk4KZmRWcFMzMrOCkYGZmBScFMzMrOCmYmVnBScHMzApOCmZmVnBSMDOzgpOCmZkVnBTMzKzgpGBmZgUnBTMzKzgpmJlZwUnBzMwKTgpmZlZYt9sBNEjaCPg28CxwbURc0OWQzMzGnVrPFCR9X9LDkuY1le8j6S5Jd0v6bC5+N3BpRHwEeEedcZmZWWt1Xz46G9inXCBpAvAtYF9gV+AQSbsCM4BFebTVNcdlZmYt1JoUIuJ64NGm4j2BuyPi3oh4FrgQOABYTEoMbeOSdJSkOZLmLF26tI6wx5W+adOQ1O/VN21at8PqZyLUHmPV9dBqvI0mTBj163A4euE7d3o7Hs786lhfI/l/qoioZcbFAqQ+4CcRsXv+fCCwT0R8OH9+P7AX8BngDOBp4IYqdQozZ86MOXPmrG1cNH9zwVqXDXf6EStr+r0HXA9Vx6sw7VCM1HIqL7dL62E06YXvXPX3G4n51bG+avh+cyNiZqtho6aiOSJWAh/sdhxmZuNZN25JfQDYpvR5Ri6rTNIsSWcuX768o4GZmY133UgKNwM7Sdpe0vrAwcAVQ5lBRMyOiKMmT55cS4BmZuNV3bek/hC4CdhF0mJJR0bEKuAY4GrgTuDiiJhfZxxmZlZNrXUKEXHIAOVXAlfWuWwzMxu6nmzmwnUKZmb16Mmk4DoFM7N69GRSMDOzetT+8FqdJC0FFg5hkinAIzWFMxyOa2gc19A4rqEZD3FtFxFTWw3o6aQwVJLmDPQUXzc5rqFxXEPjuIZmvMfly0dmZlZwUjAzs8J4SwpndjuAATiuoXFcQ+O4hmZcxzWu6hTMzKy98XamYGZmbTgpmJlZYcwkhQH6fS4Pnyjpojz8N7nzn8awz+XyuyS9bTTEJWkLSb+UtELSGZ2MaZhxvUXSXEm3579vGiVx7Snp1vy6TdK7RkNcpeHb5t/yU92OSVKfpKdK6+s7nYppOHHlYXtIuknS/LyNTep2XJLeV1pXt0p6XtIrRkFc60k6J6+nOyV9riMBRUTPv4AJwD3AS4D1gduAXZvG+QfgO/n9wcBF+f2uefyJwPZ5PhNGQVwbAa8DjgbOGEXr65XA9Px+d+CBURLXhsC6+f1WwMONz92MqzT8UuAS4FPdjgnoA+Z1cpvqUFzrAr8HXp4/bzEa/hebxvkr4J5Rsr4OBS4sbf8LgL7hxjRWzhQG6ve57ADgnPz+UuDNkpTLL4yIZyLiPuDuPL+uxhURKyPiBlL3pJ02nLh+FxEP5vL5wAaSJo6CuJ6M1Cw7wCRa95g64nEBSHoncB9pfY2KmGo0nLjeCvw+Im4DiIhlEbF6FMRVdkietlOGE1cAG0laF9gAeBZ4fLgBjZWksDWwqPR5cS5rOU7eeSwnHYlUmbYbcdWpU3G9B7glIp4ZDXFJ2kvSfOB24OhSkuhaXJI2JvU//qUOxTWPBO8AAAP+SURBVDLsmPKw7SX9TtJ1kl4/SuLaGQhJV0u6RdKnR0lcZQcBPxwlcV0KrAQeAu4H/jUiHh1uQKOmj2brLZJ2A75GOrobFSLiN8Bukl4GnCPpqoio40xrKE4EvhERK+o/SK/sIWDbiFgm6a+ByyXtFhHDPsocpnVJl0xfDTwJ/Fypg/mfdzesRNJewJMRMa/bsWR7AquB6cCLgF9J+llE3DucmY6VM4Uq/T4X4+TTrcnAsorTdiOuOg0rLkkzgB8Dh0fEPaMlroaIuBNYQarz6HZcewGnSFoAHAt8XtIx3YwpXypdBhARc0nXtHfuQEzDiot0lHx9RDwSEU+SOuJ61SiIq+FgOnuWMNy4DgX+KyKei4iHgRuB4beN1KkKk26+SEcY95IqihuVNbs1jfN/6F9Zc3F+vxv9K5rvpXOVW2sdV2n4EXS+onk462uzPP67R9nvuD1rKpq3Ax4EpnQ7rqZxTqRzFc3DWVdTG9s4qYLzAWDzURDXi4BbyDcNAD8D9u92XPnzOnk9vWQUbfOfAX6Q328E3AHsMeyYOvkFu/kC9gP+SDrqOT6XfRl4R34/iXT3x93Ab8s/LnB8nu4uYN9RFNcC4FHSUe9imu5K6EZcwAmk65i3ll4vHgVxvZ9UkXtr3rG8c7T8jqV5nEiHksIw19V7mtbVrNGyroDDcmzzgFNGUVxvBH7dyXg68DtunMvnkxLCcZ2Ix81cmJlZYazUKZiZWQc4KZiZWcFJwczMCk4KZmZWcFIwM7OCk4JZJml1U2uYfZLeKOknHZh30TJpqWzFINO8XtIdkkbLE7Q2DriZC7M1noqIfk0il5t17oB7muffTkT8StJ+wLCTkllVPlMwq0jSiSr1hyBpXj4DeLWk30uaJGmj3BdA5SY2JJ2bW1JtfL5AUnNLmWYjwmcKZmtsULq8c19EVOqoJyJulnQFcBKpCePzY2iNpp0F/COpYbrJwGuBDwxherOOcVIwW+MFl4+G4MvAzaT+Lz4+lAkj4jpJ35Y0ldQExWXRuWa/zYbEScGsulX0v+Ra7ipyC1JbNOvl8pVDnPe5pHZ/DgY+OIwYzYbFdQpm1S0gN+Us6VWkli0bvgt8AbiA1M/EUJ1NalqbiLhjOEGaDYfPFMyquww4PPfu9htSy5ZIOhx4LiL+Q9IE4H8kvSkifjHQjHK7+EWPdRGxRNKdwOW1fgOzQTgpmGURsXGLsmuBa/P7p2jd09wC0uUfIvUpvFeFxe1GaioZAEkbAjvR+U5czIbEl4/MRsZqYHJ+KO5o0s7/BABJewN3AqdHxPLGBLnv5NnAI12I18Yp96dgZmYFnymYmVnBScHMzApOCmZmVnBSMDOzgpOCmZkV/j+IDmifbF2yQwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(histStack,bins=60,align='right',edgecolor='black',color='red')\n",
    "plt.xlabel('Flux [Jy]')\n",
    "plt.ylabel('Cluster Count')\n",
    "plt.yscale('log')\n",
    "plt.title('Central Flux at MaDCoWS Cluster Locations \\n Number of clusters: ' + str(len(histStack)))\n",
    "plt.savefig(\"/home/alex/Downloads/madcowsHistogram.pdf\",bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1230"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(histStack)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pick out central emission >300mJy\n",
    "i = 0\n",
    "while i < len(histStack):\n",
    "    if(histStack[i] > .3):\n",
    "        print(str(i))\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# central emission #hdata\n",
    "# cent_em = central_emission(stackArr[0], int(100), pixel_size.to(u.arcmin).value, 1.0)\n",
    "# print(np.mean(cent_em), np.std(cent_em))\n",
    "# print(np.mean(hdata), np.std(hdata))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# emissionStack = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# i = 0\n",
    "# while i < len(stackArr):\n",
    "#     if(stackArr[i].shape == (720,720)):\n",
    "#         cent_em = central_emission(stackArr[i], int(100), pixel_size.to(u.arcmin).value, 1.0)\n",
    "#         emissionStack.append(np.sum(cent_em))\n",
    "#     i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# len(emissionStack)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# i = 0\n",
    "# histStack = []\n",
    "# while i < len(emissionStack):\n",
    "#     if(emissionStack[i] > .0025):\n",
    "#         histStack.append(emissionStack[i])\n",
    "#     i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# len(histStack)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.hist(histStack,bins=60,align='right',edgecolor='black',color='red')\n",
    "# plt.xlabel('Flux [Jy]')\n",
    "# plt.ylabel('Cluster Count')\n",
    "# plt.yscale('log')\n",
    "# plt.title('Central Flux at MaDCoWS Cluster Locations \\n Number of clusters: ' + str(len(histStack)))\n",
    "# plt.savefig(\"/home/alex/Downloads/madcowsHistogram.pdf\",bbox_inches='tight')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the rest are just notes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pixel_area = (pixel_size**2).to(u.sr)\n",
    "\n",
    "# print(pixel_area)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# vmap = vmap*pixel_area\n",
    "# print(vmap[0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(vmap)\n",
    "# stamp = Cutout2D(vmap, (int(723/2),int(724/2)), 200. , mode = 'strict')\n",
    "\n",
    "# fig, axes = plt.subplots()\n",
    "\n",
    "# plt.imshow(stamp.data.value)\n",
    "# draw_circle = plt.Circle((100, 100), 50, color = 'r',fill=False, alpha = 0.5)\n",
    "\n",
    "# plt.colorbar()\n",
    "\n",
    "# axes.set_aspect(1)\n",
    "# axes.add_artist(draw_circle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# fwhm = 6*u.arcsec\n",
    "# omega_B = np.pi * (fwhm/2)**2/np.log(2)\n",
    "# print((vmap*u.Jy/u.beam).to(u.MJy/u.sr,equivalencies=u.beam_angular_area(omega_B)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# vmap = vmap*pixel_area\n",
    "# print(vmap[0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cent_em = central_emission(vmap, int(100),pixel_size.to(u.arcmin).value,1.0)\n",
    "\n",
    "# aperture_area = np.pi*(1.0*u.arcmin)**2\n",
    "# average_flux = cent_em[0]/aperture_area\n",
    "\n",
    "# print(average_flux.to(u.MJy/u.sr))\n",
    "\n",
    "# map_var_Jy = np.std(vmap[150:199, 0:199])/pixel_area.unit\n",
    "# map_var = (map_car_Jy).to(u.MJy/u.sr)\n",
    "# print(map_var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# how to read text files I sent to Jack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "# path = \"/home/alex/Downloads/\"\n",
    "\n",
    "# mooFile = open(path + \"madcowsCentralFlux.txt\",'r')\n",
    "# mooLines = mooFile.readlines()\n",
    "\n",
    "# actFile = open(path + \"actCentralFlux.txt\",'r')\n",
    "# actLines = actFile.readlines()\n",
    "\n",
    "# mooNameArr = [] # moo names stored here\n",
    "# mooFluxArr = [] # moo fluxes stored here\n",
    "# for line in mooLines:\n",
    "#     mooNameArr.append(line.strip()[2:15])\n",
    "#     mooFluxArr.append(float(line.strip()[18:len(line.strip())-2]))\n",
    "\n",
    "# actNameArr = [] # act names stored here\n",
    "# actFluxArr = [] # act fluxes stored here\n",
    "# for line in actLines:\n",
    "#     actNameArr.append(line.strip()[2:20])\n",
    "#     actFluxArr.append(float(line.strip()[23:len(line.strip())-2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# old central emission\n",
    "# def central_emission(hmap, map_w, pixel_w, max_sep):\n",
    "#     #Returns the average surface brightness in Jy/pixel of pixels within max_sep \n",
    "#     # of the center of hmap (assuming hmap has units MJy/pixel)\n",
    "#     center = map_w + 0.5\n",
    "#     emission = 0\n",
    "#     div = 0\n",
    "#     hist = []\n",
    "#     for i in range(2*map_w):\n",
    "#         for j in range(2*map_w):\n",
    "#             if pixel_w*np.sqrt((center-i)**2+(center-j)**2) <= max_sep:\n",
    "#                 hist.append(hmap[i,j])\n",
    "#                 div +=1\n",
    "#     return  np.array(hist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#old attempt at naming central emission clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# nameArr = []\n",
    "# j = 0\n",
    "# while j < 2603:\n",
    "    \n",
    "#     if(j == 1059 or j == 2544):\n",
    "#         hdu = fits.open('/home/alex/Downloads/moo/'+mdcwName[j])\n",
    "#         header = hdu[0].header\n",
    "#         obsRA = header.get('OBSRA')\n",
    "#         obsDec = header.get('OBSDEC')\n",
    "\n",
    "#         done = False\n",
    "#         sep = .8\n",
    "#         while(done == False):\n",
    "#             i = 0\n",
    "#             num = []\n",
    "#             while i < len(mdcwName):\n",
    "#                 if(mdcwRA[i] < obsRA+sep and mdcwRA[i] > obsRA-sep and mdcwDec[i] < obsDec+sep and mdcwDec[i] > obsDec-sep):\n",
    "#                     #print(\"real one is \" + str(obsRA) + \" \" + str(obsDec))\n",
    "#                     #print(mdcwRA[i], mdcwDec[i])\n",
    "#                     num.append(i)\n",
    "#                 i+=1\n",
    "#             if(len(num) == 1):\n",
    "#                 done = True\n",
    "#             elif(len(num) == 0):\n",
    "#                 sep +=.0001\n",
    "#             else:\n",
    "#                 sep -=.0001\n",
    "#         name = mdcwName[num[0]]\n",
    "#         nameArr.append(name)\n",
    "#         print(str(j) + \" complete\")\n",
    "#         j+=1\n",
    "                \n",
    "#     hdu = fits.open('/home/alex/Downloads/moo/'+mdcwName[j])\n",
    "#     header = hdu[0].header\n",
    "#     obsRA = header.get('OBSRA')\n",
    "#     obsDec = header.get('OBSDEC')\n",
    "    \n",
    "#     done = False\n",
    "#     sep = .8\n",
    "#     while(done == False):\n",
    "#         i = 0\n",
    "#         num = []\n",
    "#         while i < len(mdcwName):\n",
    "#             if(mdcwRA[i] < obsRA+sep and mdcwRA[i] > obsRA-sep and mdcwDec[i] < obsDec+sep and mdcwDec[i] > obsDec-sep):\n",
    "#                 num.append(i)\n",
    "#             i+=1\n",
    "#         if(len(num) == 1):\n",
    "#             done = Truepickle.dump(favorite_color, open(path+\"save.pickle\",\"wb\"))\n",
    "#         elif(len(num) == 0):\n",
    "#             sep +=.005\n",
    "#         else:\n",
    "#             sep -=.005\n",
    "    \n",
    "#     name = mdcwName[num[0]]\n",
    "#     nameArr.append(name)\n",
    "#     print(str(j) + \" complete\")\n",
    "#     j+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#old attempt at pickle\n",
    "# cutout = fits.open('/home/alex/Downloads/moo/'+mdcwName[i])\n",
    "# vmap = cutout[0].data[0][0]\n",
    "# #Convert from Jy/beam to Jy/pixel\n",
    "# conv_factor = (np.pi*(2.5/2)**2/np.log(2))\n",
    "# vmap = (vmap/conv_factor)\n",
    "# cent_em = central_emission(vmap.data, int(100), pixel_size.to(u.arcmin).value, 1.0)\n",
    "# toBePickled = {nameArr[i]: np.sum(cent_em)}\n",
    "# pickle.dump(toBePickled,open(path+\"save.p\",\"wb\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}