specreduce_extract_jwst_OLD.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from astropy.io import fits\n",
    "from astropy.table import Table\n",
    "from astropy.visualization import simple_norm\n",
    "\n",
    "from jwst import datamodels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# useful function that work for boxcar, boxcar scaled with wavelength,\n",
    "# and psf-weighted extractions\n",
    "import numpy as np\n",
    "from gwcs.wcstools import grid_from_bounding_box\n",
    "\n",
    "\n",
    "def get_boxcar_weights(center, hwidth, npix):\n",
    "    \"\"\"\n",
    "    Compute the weights given an aperture center, half widths, and number of pixels\n",
    "    \"\"\"\n",
    "    weights = np.zeros((npix))\n",
    "    # pixels with full weight\n",
    "    fullpixels = [max(0, int(center - hwidth + 1)), min(int(center + hwidth), npix)]\n",
    "    weights[fullpixels[0] : fullpixels[1]] = 1.0\n",
    "\n",
    "    # pixels at the edges of the boxcar with partial weight\n",
    "    if fullpixels[0] > 0:\n",
    "        weights[fullpixels[0] - 1] = hwidth - (center - fullpixels[0])\n",
    "    if fullpixels[1] < npix:\n",
    "        weights[fullpixels[1]] = hwidth - (fullpixels[1] - center)\n",
    "\n",
    "    return weights\n",
    "\n",
    "\n",
    "def ap_weight_images(\n",
    "    center, width, bkg_offset, bkg_width, image_size, waves, wavescale=None\n",
    "):\n",
    "    \"\"\"\n",
    "    Create a weight image that defines the desired extraction aperture\n",
    "    and the weight image for the requested background regions\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    center : float\n",
    "        center of aperture in pixels\n",
    "    width : float\n",
    "        width of apeture in pixels\n",
    "    bkg_offset : float\n",
    "        offset from the extaction edge for the background\n",
    "        never scaled for wavelength\n",
    "    bkg_width : float\n",
    "        width of background region\n",
    "        never scaled with wavelength\n",
    "    image_size : tuple with 2 elements\n",
    "        size of image\n",
    "    waves : array\n",
    "        wavelegth values\n",
    "    wavescale : float\n",
    "        scale the width with wavelength (default=None)\n",
    "        wavescale gives the reference wavelenth for the width value\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    wimage, bkg_wimage : (2D image, 2D image)\n",
    "        wimage is the weight image defining the aperature\n",
    "        bkg_image is the weight image defining the background regions\n",
    "    \"\"\"\n",
    "    wimage = np.zeros(image_size)\n",
    "    bkg_wimage = np.zeros(image_size)\n",
    "    hwidth = 0.5 * width\n",
    "    # loop in dispersion direction and compute weights\n",
    "    for i in range(image_size[1]):\n",
    "        if wavescale is not None:\n",
    "            hwidth = 0.5 * width * (waves[i] / wavescale)\n",
    "\n",
    "        wimage[:, i] = get_boxcar_weights(center, hwidth, image_size[0])\n",
    "\n",
    "        # bkg regions\n",
    "        if (bkg_width is not None) & (bkg_offset is not None):\n",
    "            bkg_wimage[:, i] = get_boxcar_weights(\n",
    "                center - hwidth - bkg_offset, bkg_width, image_size[0]\n",
    "            )\n",
    "            bkg_wimage[:, i] += get_boxcar_weights(\n",
    "                center + hwidth + bkg_offset, bkg_width, image_size[0]\n",
    "            )\n",
    "        else:\n",
    "            bkg_wimage = None\n",
    "\n",
    "    return (wimage, bkg_wimage)\n",
    "\n",
    "\n",
    "def extract_1dspec(image, jdatamodel, center, width, bkg_offset, bkg_width, wavescale=None):\n",
    "    \"\"\"\n",
    "    Extract the 1D spectrum using the boxcar method.\n",
    "    Does a background subtraction as part of the extraction.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    jdatamodel : jwst.DataModel\n",
    "        jwst datamodel with the 2d spectral image\n",
    "    center : float\n",
    "        center of aperture in pixels\n",
    "    width : float\n",
    "        width of apeture in pixels\n",
    "    bkg_offset : float\n",
    "        offset from the extaction edge for the background\n",
    "        never scaled for wavelength\n",
    "    bkg_width : float\n",
    "        width of background region\n",
    "        never scaled with wavelength\n",
    "    wavescale : float\n",
    "        scale the width with wavelength (default=None)\n",
    "        wavescale gives the reference wavelenth for the width value\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    waves, ext1d : (ndarray, ndarray)\n",
    "        2D `float` array with wavelengths\n",
    "        1D `float` array with extracted 1d spectrum in Jy\n",
    "    \"\"\"\n",
    "    # should be determined from the gWCS in cal.fits\n",
    "#     image = np.transpose(jdatamodel.data)\n",
    "\n",
    "    print(dir(jdatamodel.slits[0].meta.wcsinfo))\n",
    "    print()\n",
    "    print((jdatamodel.slits[0].meta.wcsinfo.waverange_start))\n",
    "    print()\n",
    "    print((jdatamodel.slits[0].meta.wcsinfo.instance))\n",
    "    print()\n",
    "    print(jdatamodel.slits[0].meta.wcsinfo.crpix3)\n",
    "    print(jdatamodel.slits[0].meta.wcsinfo.crval3)\n",
    "    print(jdatamodel.slits[0].meta.wcsinfo.cdelt3)\n",
    "    print(jdatamodel.slits[0].meta.wcsinfo.ctype3)\n",
    "    print(jdatamodel.slits[0].meta.wcsinfo.cunit3)\n",
    "    \n",
    "\n",
    "#     print(dir(jdatamodel.slits[0].meta.photometry))\n",
    "\n",
    "\n",
    "    grid = grid_from_bounding_box(jdatamodel.meta.wcs.bounding_box)\n",
    "    ra, dec, lam = jdatamodel.meta.wcs(*grid)\n",
    "    lam_image = np.transpose(lam)\n",
    "\n",
    "    # compute a \"rough\" wavelength scale to allow for aperture to scale with wavelength\n",
    "    rough_waves = np.average(lam_image, axis=0)\n",
    "\n",
    "    # images to use for extraction\n",
    "    wimage, bkg_wimage = ap_weight_images(\n",
    "        center,\n",
    "        width,\n",
    "        bkg_width,\n",
    "        bkg_offset,\n",
    "        image.shape,\n",
    "        rough_waves,\n",
    "        wavescale=wavescale,\n",
    "    )\n",
    "\n",
    "    # extract the spectrum using the weight image\n",
    "    if bkg_wimage is not None:\n",
    "        ext1d_boxcar_bkg = np.average(image, weights=bkg_wimage, axis=0)\n",
    "        data_bkgsub = image - np.tile(ext1d_boxcar_bkg, (image.shape[0], 1))\n",
    "    else:\n",
    "        data_bkgsub = image\n",
    "\n",
    "    ext1d = np.sum(data_bkgsub * wimage, axis=0)\n",
    "    # convert from MJy/sr to Jy\n",
    "    ext1d *= 1e6 * jdatamodel.meta.photometry.pixelarea_steradians\n",
    "\n",
    "    # compute the average wavelength for each column using the weight image\n",
    "    # this should correspond directly with the extracted spectrum\n",
    "    #   wavelengths account for any tiled spectra this way\n",
    "    waves = np.average(lam_image, weights=wimage, axis=0)\n",
    "\n",
    "    return (waves, ext1d, data_bkgsub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data is taken from s2d file. x1d is used for comparison with pipeline extraction.\n",
    "from astropy.utils.data import download_file\n",
    "\n",
    "s2dfilename = \"nirspec_fssim_d1_s2d.fits\"\n",
    "x1dfilename = \"nirspec_fssim_d1_x1d.fits\"\n",
    "mainurl = \"file:/Users/busko/Downloads/NIRSpec_s2d/\"\n",
    "\n",
    "s2dfile_dld = download_file(mainurl + s2dfilename)\n",
    "x1dfile_dld = download_file(mainurl + x1dfilename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rename files so that they have the right extensions required for the jwst datamodels to work\n",
    "import os\n",
    "\n",
    "s2dfile = s2dfile_dld + '_s2d.fits'\n",
    "os.rename(s2dfile_dld, s2dfile)\n",
    "x1dfile = x1dfile_dld + '_x1d.fits'\n",
    "os.rename(x1dfile_dld, x1dfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use a jwst datamodel to provide a good interface to the data and wcs info\n",
    "s2d = datamodels.open(s2dfile)\n",
    "image = s2d.slits[0].data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'slit[0]')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAABwCAYAAACacAjbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAdpUlEQVR4nO3df6w0V33f8c939t7neQx2cAiuRW0KJLESOZFwoialIq0oqI1Jq7ptaGRSpbSicivZEkhRCgS1pVKRyB8NTZQ0klMotE3q0AQUFKGmCBy1/aMEp3FCDFhxYtPgGgwpP8wP+7l359s/5pzdM2dn9vfOzJ19v6Tn2d2ZM3POnP11PvfM7pq7CwAAAADQr6LvBgAAAAAACGcAAAAAMAiEMwAAAAAYAMIZAAAAAAwA4QwAAAAABoBwBgAAAAADQDgDAAAAgAEgnAEALjwze6uZ/adw/c+Z2VfNbNJS9kVm5qHMXWvu/91m9g0z+8w+2w0AQIpwBgAYFXf/P+5+rbtPJcnMfsvM/lFD0evd/d54w8xeaWafMrOvm9n9ZvbCZJ//QNKrDt54AMBRI5wBAI6emT1P0vsk/TNJz5X0gKRf6bVRAICjQzgDAFwoZvZGM3vczJ4ys4fN7JXZ+nja4omZvU3SX5L0c+E0xp9r2e3fkfSQu/8Xd39a0lslvcTMvvOgBwMAQIJwBgC4MMzsOyTdI+n73P06ST8o6bG28u7+Fkn/Q9I94VTHe1qKfpek30u2+5qkPwrLAQDoxEnfDQAAYANTSZcl3Wpmn3f3xyTJzHbd77WSPp8t+7Kk63bdMQAA62LmDABwYbj7I5LeoOq0wyfN7D4z+7N72PVXJX1TtuybJD21h30DALAWwhkA4EJx91929x+Q9EJJLumnVm2yxm4fkvSSeMPMni3p28JyAAA6QTgDAFwYZvYdZvYKM7ss6WlJ35BUrtjsc5K+dUWZ90v6bjP7YTO7IumfS/p9d//Uzo0GAGBNhDMAwEVyWdLbJX1B0mcl/RlJb16xzc9IerWZfdHMfrapgLt/XtIPS3qbpC9K+guS7txXowEAWIe5r3O2BwAA4xB+XPphVTNvP+Huv7jGNu+U9HclPenu337gJgIAjhThDAAAAAAGgNMaAQAAAGAACGcAAAAAMACd/gj15Lpn+8kN13dZJQAAAAAMxtVH/+8X3P2GpnWdhrOTG67XTW+7u8sqAQAAAGAwHv3Rt3y6bR2nNQIAAADAABDOAAAAAGAACGcAAAAAMACdfuYMR846+k09t27qAQAAAPaIcLaLrsIGNrOP+4WABwAAgI6NN5wNIDgxvO/HXu75ATx+ZgiKAAAAR2G04ayL4aztOoDvIgAc4cB+1yP2nvosPhosuz2ooIiL7whfEwAAuCjGGc7MNTmZypaMQdqCVdPytv3EsrPLZF2R7ScvOy+3uv58X7FcGiLKhgHXNiEjrd/d1g6gsY15O4qW487bGNc3HWubWFfcR1p307GXya7z9Xlfnk8LlWUhM5eZazIpW8svW1Ytb26/r2jvqu1XbberZfVutqNhh4G+wviutr57CPvjd0Ef0wCAkYazwlzPunI1GfDX1+dBIA0OaTiwZFlT2IrLJ1YN3E+KsraPQrGecla2MFeh+v4KKzVZ2Gb5+qiUaeqm0udfvNkU1PJt0n2Vstp+C/PaPmbHY2WtnsLmgWUSjiu2J7b3xKa1fUxUbTNVoTK0O/bPRGWtr2PZWL52DG6zNp75RKUXC30R18fjjetKt9myMgl28fq5T/Sn33iWrp6faFKUmhSlrj29qklRqnSrBzlZbR9RU3Cehcls+bQsFrap1s+vNwW5VfdzWldt2ZoDt1XlmsLbvva9bn3bSuu3DWLOPtuwvKLV/dPn8HvXQJt247HEiE4jMQEcFx1/YMARG2U4k6RJ4ZoUZeOMTBq64mVTQJuFBKvPmsR1J0WpQq6TYroQvPJAdVJMNTHXiU0XwkdhrtMsxMR6Y2iKbWgKZ7mpt/9CQlo+31dbHc3HXy07takmCserchai4rHF4zq187Cv+b6nbipVqFDVV6d2rkuhfK2usK8yCWhTLzRVFbZinXF/V/1kIbhUYTDdJlzOQl11eVZONFWhL11zja6WJ7P+uGZyNgunMQRKqoW9pvBXlWkKgiFAZuGuKeCtCoD5sdbKZH2Zzza2WTULWZVp3nb5LGD7urawuWyYuW1I2Hamctt9bl52eSMOHRJXtW/XQDu4YVcHA8FD1HBRZ32Hhig7QPyBoX+8vvRmnOHMXJdOznXaEM7yWa10BixdF2dz4jbpbFIMYSfFVKdW6vKkCh4xeMUZoNN420qd2nT2rwoy80AjSVeKs1kISQPMLNQsefuoQtxigMqVKmrBbR7GSk1Cm+LtuD6tdyqblwltf7ad6TSEyIn5LLRcsipKTaz6Mb1LZtVtmQqLocJ1Jlch6dQKnWqiZxWXNPV0Rm4xaE69VCnXmS8GuVKlvu5TTT0NuXHmzHUm6cyrYzkLIe3MC535JFw/0VSmp/1UUy9mtyfyWhg884kkzYJgLBdn8eL6tL9nIdBtFijTQFZfVswC3lk5qQW+NBjOyq8Ig/N18/2m69M/CtTLL4a5vK58Vrlp+9TiDGHTKaj1Wca29Zvss8mq01+32X6bALqwjy3q3UfZ+nZbbbZTnfvafv16Ztf2sK/Dtzm/T3b+3PMuLvDALT1l392G98eCC4goNUJjCsgX7PVqnOEsKN10Np3Mbi87ZbEpxKUBTqp/PmrdmbN4+7SY1kJfDHCxrlObNs5QzWbUkpe+fLasCkarw9k0GZin+5yFtHCZBrB4e5rVGcPlFTubhcN6OJvO9jEJ7YshMN1nDC8TK3Wqqa4UZ1KY5Zq1Mzu2MgSh6njq4W0q09PlpVp74/HEcHXVJ7UZtjirFgNaGq7OfFILWHF5DFXx9MzzWD4JUtNaSNksSOUhauHUyOyUzbbZs7bQsir8rNpPU51NNpl12yTorKp7k7eUZfXmn8Fcb3+71X2Iberb77T5GvvvOlgNwAV744+OeeZtoz9wHLAdwKAd8WtEn0YdzqZuevrqafU5IS2eVpJ/5qx+auO8XFM4i9ubeW2mK52Bs+R6DG5Np1k2faZttq7hbaEpnOUBplY+OYVumbyuxVC0GJhOinKttqefW2vaZxpmm0/VzENG0TzjItN5uTjbloee/FTCplMIpfnnwVZ9biw/xbD+mbN6W5Z9GUm6j9n6lnJty6rljYvX+BzZqvVLV2+1z33UMd+4i9mL/t+sGCzioAbwGAeAYzXKcOZloaeevixJOj+fLC2bBrKmb2VcdtpI/i2NTV8wkn8pSeuXj+w4fbzOtxzmwWFVnctOVytajj2Wy7dt+kKVfJC7bh+s+kKMGMbb1jXtY357sa5l4cmzsk3X58taGrVkm+XLN9/Xpm1avvHhBnB9BqCjDT6xz/d5KguDfAAANjLOcOamr3/tiqTtQs/WQWnN7ZZ9xT8Ob61AssOg8lDBYtShgUH8cHBfAADQm1GGs9Q2A9qFGZ21N1xecn565BaNwsHsO0xd+LuXwTkAAEAvRh/O9jHQXDrY3mCWbQifVUFPuO8BAACwwvjD2aF1NeiOIXAfnwvJ2zymr0vdFuEJAAAAPSOcXRR5eNhnmCCYAAAAAL1b/M5xAAAAAEDnCGcAAAAAMAArw5mZvcDM7jezT5jZQ2b2+rD8uWb2ITP7w3D5zYdvLgAAAACM0zozZ+eSftzdb5X0Ukl3m9mtkt4k6cPufoukD4fbAAAAAIAtrPxCEHd/QtIT4fpTZvZJSTdJukPSy0Ox90j6LUlvPEgrsRK/nTbHj3wDAADgItro2xrN7EWSvkfSRyXdGIKbJH1W0o0t29wl6S5JmjzvOdu280IjOHVr0/4mzAEAAGAI1g5nZnatpF+T9AZ3/4olI1p3d7PmH8ty93sl3StJl7/1ps5jCsEIq2zzGIkP/74eXwRKAACA8VkrnJnZqapg9kvu/r6w+HNm9nx3f8LMni/pyUM1clt7GzjH/VjDsgPz5DfIWvLvheArfkstHltbuUEce9K0vkN/3/WPEYEXAAD0bWU4s2qK7J2SPunuP52s+oCk10p6e7j89YO0cEc+XXPEtSI8pINhs2xwvM6POOeD6WWD65b9NW4SQ0ufPyS99Fji5Yr+3aEKtQW3WKWrHqzzMmkfrtrXquUt26czbWuFgE3C6DrHuaVBhOJl9vlb7AM/VAxH+h5AqAcA7NM6M2cvk/Rjkj5uZg+GZT+pKpS918xeJ+nTkn7kME3cnLskr2ZhvFzxztkUGlYEKbewbJNt1ymbLbOGbdy8ttzD4Lmp7Ma2GZx6drmwvmrXwhg/3S5vutfLtx7aQjDKyreFrnS7YrFPq02yBlt2udCWbHo1LWf17vFa+YY6sv1YEWYV2x7LGwaoTQaTnh/wvsLavga0G3++cEQJjFDQmzTIE+qxbwR+4Lit822N/1Ptw4BX7rc5++Glyb96Ug24J+m7aF4wCTX56DkJVU1BwVyy82y9z9el+7EyXoaytWX1uq2sX8agMssXTXW0dkRy3Ravz/YZ60z2GZebu6ysynphstJVTJN9lVXZ4lyy0uVmkknlpFo9Oau2t2nos9JrdRRnrvLEVJ5Ux1xMXXYuFee+GMqK2O4Y+OYH6GaaXjF5YVXdJpWTIvRlNV1VTkxeSF6EPjDp/IqpvBTaN533Q9yHLBz7SeyDej96LDNp6GfF9V7bblZ/rUyyLARWL1zlqaTnnMnPTcVXsqdrvB+LxcwU621qU6OGx0euMbDGP1S07au1vjVHtKv2tckgxnzFDOwG+9r3PrYMjekgbt2ZHDupnvA+jQ/gdSvrMIU0Pa7y5T0OYFedhr1sm16kry0X0ZGFFQI/sJv8vfGi/cFjo29rvDBCEGqalVmYYaoFsXiZzPak65XcwWkdnlSTBS4rs3Jl0/r6tmm9tffz9PZs/eKreC28ZH0Qw1NVINbpCwFxHs5CaAlhIQatuP94LDGEeVGFqvKkqmRytVpenLms9Pmxqqq3uFrKJ6bp5UI2dRVnruKsVHE2rbdVVTic3QGFpGmcNazWnU1PpEIqT2NIq9YX0yrAlJPq0idVeTepmLi8sCpcTn3ev6UlYbDKr5aGsziuzX8psGk2z5Kps6bAFPeR3F3ukslk7vVNs8DtVn/81Cbu3OazvI3pLVmePk6aZjNVtae2CyXbrrIwO2rzvzq0buPL9902gG+qL9ZZ23fDlO0y67y4t/1BZOV2WdvW3ayh6LKBXe0NatPZ9nXK7yuAtO0m/yNare6G5Qd6Q94klK2zjZkvrB/VLO+utu2K2uO9YRmAUcrfB/M/Xg79tPRRhjMrXOWVaZjqqK/zhXuscRQ3GyjHm/UKQpEkdKX7W5jVck8Cly0ErLjL+jZhedOb0so3qoYCSQX5IS895VAhXJons3fzdLdwLMlsk8xl06K+rziet3q9XmgePtN9tR7HojgjNZutagxKcf38uhfekF3qDa7Ndi3Menl7uxq3ycq3bluVM0kqpOl10+Z+Wbb9utZ9gdr0hWwPM0I76XqAO9AXeikZ7MeH9qRsL3yR7fs+2GR2uO11a51qunisDvjxeXDHfOzAEkMNKIc09GMeZTiTqTqdsWHmq8azd922kNK2eZGV8fBfmgpC/d6wvvXUxHQWa6twpvqAIp8sXDFZsVBPPvOS1d/URg+fk0pnyhb+ot30F+5anW2NbF6cn8rXGM6Sy6XBaqFsw4xXW1tWfp6s6XbLNnGCq/D5qZNt9a7a37Li275QpTNQ+x5cDvDFcxSzGU0zTIesZwwO9ccLAAex6uduhj44x3EbZziTVMTPVWx86s4WlS2pY9kpRr5i2+V1Llm37xedPCh1MT7d4S/QM3sIMPEFvHVX+/rs1NI2hKA7WVGwS7XjGUFgOUYDHJwwYAKwT7ym4CIaZTirfcRn04Hjpllu2devt+3OW2+MU1df87/trMYegtNB8eYyGE1fwLHrfgAAAKJRhrOaAw+CRnGa067aZrlaP7Q3AAyOOze2QDK24wEAAP0bbThrGzjxFbUHsOFnw9ANwgMAAMDFMtpw1qbvASvh8Dj0/TgDAADAxXN04axvqwbthLfhI3gBAADgEAhnA5MP/Alr/SCAAQAAoGuEsy0QmMaP30YBAABA10YfzghS2KdNH0+EOQAAAKxrtOFsr6Fs/qNpzcuxtqYfBR/zzxFs9DiMXdP2eMPRINQDAHCcRhnO3LWX4LQQJMI+Y5hoChpN5Y/aGj9A7cp+yHvoP1rdur8dt8+b0+HjZxQBeUSB5thn/M3mfUBQBQAck1GGM7lUnk2afxy5bdCzMDC25nUWbqYBsCXEySVr2o83lW2YUWoq1ySfaVl3YJdWuWwbUziWhrJu1fK8rWHfC8cQ95NmsaK5LevmhYWuM83vp3g7W1/bJq3Ikn1adjuWtaRsevfG/TS1Zx1NB2wt11dt11pH8+KN8vCWQe4gg+x9hcqBBYBRhOUdzF5jFULawO4fAKvF13x3/sgCbGKU4czdpKuFVEoq0hX1y9bglAWmWZjIXlzsXPVwEsftZdgmXqb7KVUPJ0lYqWWEMitTar6PhQNO2mfZMeTbhHKeBYumOmIZn0g2Df9Kn5WLfVBMk22n9YAyuarZK7OVrskzLitdxTNhg8JUnprOrylkZdh/qCf2Vbyb5uGw3gnlaVG1M4Q8L0xeSNPT6tInVj+eUNYn8SBCOTOpmK9XIZUTzZal/eaTUG62D6/WFYt96+n9EvsmXx8PsGl5kd1WVj4N5WlwXKEWJvcRzlZNJG/z5rwypFjtYlWx3erazC6DEW9rcPqEyJd36dADrQPOIh978D0YBt/IpG/VXZ0NkAZCjM+xhPxRhrOZIrsdB8jh0pM3aWt8ZwmnL8qyGYwwQCokcw+nUVYhzT2pN8zemHtVn0tWzE+HjKHMVYWb2YtJw1+ZGmetGppUCwH5tnEXTYPxhhmt2THEsw7LKtCm4UylQqhKdz5XTOeDSXOTTas6ium8eAxKs/3ks3TLZj8tBKikrTEglSc2n0WzhjJp6MnDVxbC6ut9vm1yX89m6xrDWZJ4i3on5yFsYXkW2mb1zY5nxYxbm33N6G27/20GyWu3edP97u+d/OBvHoSLnbgbAe0Q8jM4gB4QyvZraGH3WGZhRxnOzFw+qWZeJNUDQ9tpXeGdxWIQS0NB06MypKH4+bZZSIv1pSEw7ipcl3u1eVgWGlCNwUubV9k2+7XE0tmJtlmVOIFSarWmdnlyPWuHlVZv+5IZQCvn6xvb3SQPVFJ99ivbvha0imwfkjwGp4bwNWtyXt9s+TxE1U+bTNvgzdvlx5ldXwhgG8yQNdr1VMh9zURtNZu2xTZL97efd529vmF0ER7SGdcjQ0A7oGXdegSDKmCM0s8B992OKG3P2ALbKMOZTNJJOQ9JtQ84NZRPRss+S1Dp+uZqZl8Iks1qrf05s6TuhVP2lsyQtVnrtLG2MiuedK31+vKKF7bL+yMbIC6UX3FMtcCUXK/NVmXr01m01tMLG9av9fmyNdq89kxUW4hbu54V6zexxUD2woUVaTCDx06Dw0COeW+a/vg0tmMEsLOxDegPbWj9NbT27NM4w5kkmyTnDK5rww/F1GbK0uULV5o2brm5bht2PYVkn9+O2HSc2/xVftPxaFuzG4POBqfmrbt9U7FdThVsLbt+0dZdXPRZghG/CI/ZIN48h9AGAADWNMpwZumsyUZvzL766/HTepLttjo1qLGuPgbR29fZNs3d13nKGw0GOw5Is10d2SwQNjeIUAMAADo3ynAmaeuBabczDG3TZxdHl2fW7cXgGnR8CB4AAADNRhvO8gHgED7IuBKDVuwJAQgAAODiGW04yw1hsHohAiL2YgiPNwAAAFwsRxPOhoABOwAAAIA2+c80AwAAAAB6QDgDAAAAgAEgnAEAAADAABDOAAAAAGAACGcAAAAAMAA7hTMzu93MHjazR8zsTftqFAAAAAAcm63DmZlNJP28pFdJulXSa8zs1n01DAAAAACOyS4zZ98v6RF3/2N3vyrpPkl37KdZAAAAAHBcdglnN0n6k+T2Z8KyGjO7y8weMLMHpk99bYfqAAAAAGC8Tg5dgbvfK+leSTKzzz/6o2/5mqQvHLpetHqe6P8+0f/9ov/7x33QL/q/X/R/v+j/ftH/cy9sW7FLOHtc0guS2zeHZa3c/QYze8Dd//wO9WIH9H+/6P9+0f/94z7oF/3fL/q/X/R/v+j/9exyWuPHJN1iZi82s0uS7pT0gf00CwAAAACOy9YzZ+5+bmb3SPpNSRNJ73L3h/bWMgAAAAA4Ijt95szdPyjpgxtudu8udWJn9H+/6P9+0f/94z7oF/3fL/q/X/R/v+j/NZi7990GAAAAADh6u3zmDAAAAACwJ4QzAAAAABiAzsKZmd1uZg+b2SNm9qau6j1mZvaYmX3czB40swfCsuea2YfM7A/D5Tf33c4xMbN3mdmTZvYHybLGPrfKz4bnxO+b2ff21/JxaOn/t5rZ4+F58KCZ/VCy7s2h/x82sx/sp9XjYWYvMLP7zewTZvaQmb0+LOc50IEl/c9zoANmdsXMftvMfi/0/78My19sZh8N/fwr4RuuZWaXw+1HwvoX9dn+MVhyH7zbzB5NngO3heW8Bh2AmU3M7HfN7DfCbZ4DG+gknJnZRNLPS3qVpFslvcbMbu2ibuivuPttye9KvEnSh939FkkfDrexP++WdHu2rK3PXyXplvDvLkm/0FEbx+zdWux/SXpHeB7cFr7ISOE16E5J3xW2+bfhtQrbO5f04+5+q6SXSro79DPPgW609b/Ec6ALz0h6hbu/RNJtkm43s5dK+ilV/f/tkr4o6XWh/OskfTEsf0coh9203QeS9BPJc+DBsIzXoMN4vaRPJrd5Dmygq5mz75f0iLv/sbtflXSfpDs6qht1d0h6T7j+Hkl/q8e2jI67/3dJ/y9b3Nbnd0j6D175X5KuN7Pnd9PScWrp/zZ3SLrP3Z9x90clPaLqtQpbcvcn3P1/h+tPqXpzvkk8BzqxpP/b8BzYo/A4/mq4eRr+uaRXSPrVsDx//Mfnxa9KeqWZWUfNHaUl90EbXoP2zMxulvTXJf27cNvEc2AjXYWzmyT9SXL7M1r+hoH9cEn/zcx+x8zuCstudPcnwvXPSrqxn6YdlbY+53nRnXvCKSvvsvmpvPT/AYXTU75H0kfFc6BzWf9LPAc6EU7nelDSk5I+JOmPJH3J3c9DkbSPZ/0f1n9Z0rd02+Lxye8Dd4/PgbeF58A7zOxyWMZzYP/+jaR/KqkMt79FPAc2wheCjNsPuPv3qpq2v9vM/nK60qvfUeC3FDpEn/fiFyR9m6pTXJ6Q9K/7bc74mdm1kn5N0hvc/SvpOp4Dh9fQ/zwHOuLuU3e/TdLNqmYhv7PnJh2d/D4ws++W9GZV98X3SXqupDf22MTRMrO/IelJd/+dvttykXUVzh6X9ILk9s1hGQ7I3R8Pl09Ker+qN4rPxSn7cPlkfy08Gm19zvOiA+7+ufBmXUr6Rc1P26L/D8DMTlUFg19y9/eFxTwHOtLU/zwHuufuX5J0v6S/qOpUuZOwKu3jWf+H9c+R9KcdN3W0kvvg9nDKr7v7M5L+vXgOHMrLJP1NM3tM1UeYXiHpZ8RzYCNdhbOPSbolfFvLJVUfQP5AR3UfJTN7tpldF69L+muS/kBVv782FHutpF/vp4VHpa3PPyDp74dvi3qppC8np35hT7LPD/xtVc8Dqer/O8O3Rb1Y1QfCf7vr9o1J+KzAOyV90t1/OlnFc6ADbf3Pc6AbZnaDmV0frl8j6a+q+tzf/ZJeHYrlj//4vHi1pI+EmWVsqeU++FTyxyFT9Xmn9DnAa9CeuPub3f1md3+RqrH+R9z974nnwEZOVhfZnbufm9k9kn5T0kTSu9z9oS7qPmI3Snp/+FzliaRfdvf/amYfk/ReM3udpE9L+pEe2zg6ZvafJb1c0vPM7DOS/oWkt6u5zz8o6YdUfQj/65L+YecNHpmW/n95+Npkl/SYpH8sSe7+kJm9V9InVH3L3d3uPu2j3SPyMkk/Junj4TMfkvST4jnQlbb+fw3PgU48X9J7wjdeFpLe6+6/YWafkHSfmf0rSb+rKkArXP5HM3tE1RcZ3dlHo0em7T74iJndIMkkPSjpn4TyvAZ1443iObA2I6ACAAAAQP/4QhAAAAAAGADCGQAAAAAMAOEMAAAAAAaAcAYAAAAAA0A4AwAAAIABIJwBAAAAwAAQzgAAAABgAP4/pR4E6eQ9Y8YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x1080 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# display\n",
    "norm_data = simple_norm(image, \"sqrt\")\n",
    "plt.figure(figsize=(15, 15))\n",
    "plt.imshow(image, norm=norm_data, origin=\"lower\")\n",
    "plt.title(\"slit[0]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<TableColumns names=('WAVELENGTH','FLUX','ERROR','SURF_BRIGHT','SB_ERROR','DQ','BACKGROUND','BERROR','NPIXELS')>\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAGDCAYAAADeRuzbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3debyMdf/H8dfnHLLvRGU5CskScsSdpWhBOZUsKUISFbcWFKnQopKEkkiqm1LRL5WU5M5WZAttlkJSiMp2V9bv74+5TvfkPnOcw5m5Znk/H495nDnXzFzXe2bOmfdcuznnEBERyUiS3wFERCR6qSRERCQklYSIiISkkhARkZBUEiIiEpJKQkREQlJJyEkzs/fNrEsOjKermS0K+n2/mZ15suMNJzObZ2bd/c4RSWaWYmbOzHL5nUXCTyURA9I/iMxsnZldGzS8offPeuywfWaWy8xOMbMnzWyr94G72cxGeffbH3Q5amZ/BP3eMYMML5nZQe/2X81sjplVBXDOtXTOvZzTz9s5V9A5tzGnx5sZM2tqZh+b2R4z2xzJaQdluMjMtoZx/EPMbEq4xh9p3v9AJb9zxCuVRGxZADQJ+r0JsDaDYYudc4eBgUAqcD5QCLgIWAl/fQAXdM4VBLYAaUHDXgkx/eHe/csCPwMv5dQTiyL/ASYB/f0Okhl9i886vVYnRyURW44ticbA4xkMW+Bdrwe85Zz7yQVsds7962RDOOd+B14FasDfF7l4i4w+MbNnvG/ja83s4vTHmlkRM3vBzLaZ2Y9m9rCZJWc0neBviN6czFgze8+bU/rMzM4Kum9Vb+7mV2+Oq32IcRb35qzSvN8Lmtm3ZtbZe25LnXOTgQznYMzsUu857TGzZwAL9TqZWZKZDTCz78zsFzN7w8yKe7eNM7M3g+77uJnNNbMCwPvA6UFzdqd73/6nm9kUM9sLdDWz881ssZnt9l7PZ8zslKBxVg96TXaY2b1m1gK4F7jWG/fq470vZpZsZiPMbJeZbQSuCPWcvfvf441jn/deXOwNT38Or3u3rTSzWkGPO93M3jSznWa2ycz6BN2W7OX/znvsCjMrZ2bpf+urvedzbfqcmJdjO/CiHbMo0xvnsX9fz1pg0el+72+4jJmNMrPfvPe8TmbPO24553SJ8gswD+gOVACOAsUJFPzPQD7gh6Bhe4Am3uPuIzCXcBtQE7AQ498MXHKcDC8BD3vXCxIoiYXB+bzrXYHDwJ1AbuBaL1Nx7/a3gPFAAeBUYCnQM+ixi4Km6YBKQdP/hcBcUS7gFeA177YC3mtwo3dbHWAXUC3Ec7kM2O5N/3lgegb3uQTYfMywksA+oK333O70nmv3ENO5HVhCYM4rj/e8p3q35QfWe8+5sZe3rHfbRcDWY8Y1BDgEXO29z/mAukAD7zmnAN8Ad3j3LwRsA/oCeb3f6weNa8ox48/sfbmFwBxrOQJ/Zx97702uDJ7z2d57cbr3ewpw1jHPIf316wds8q4nASuAB4BTgDMJFHVz77H9gS+88RtQCyhx7N9J0Ot3mMAXqDzea9WVoL+tEH9fu7zXNC/wby9bZyAZeBj42O/PAl8+f/wOoEsW3qS/fwhvBq4i8EH4iTfstaBhfwB5vOHJQC/gE+AA8BPQJYPxbyZrJfEnsJvAB+w7Qf/8wfm6etOxoMcuBW4ASns58gXddl36P9+x/8gZ/BNPDLrtcmCtd/1avMIKun08MDiT5/O096HzY/qHzTG3Z1QSnYElQb8bsJXQJfENcHHQ76cR+JDM5f1eH/gV+B64Luh+F5FxSSw4znt0B4E5x/TX9fMQ9xtCUElk4X35N3BL0G2XEbokKhH48nIJkDuD6Qa/fkkEiqyx91psOeb+A4EXvevrgKtCPJ+MSuIgkDdo2N/+tkL8fT0fdNs/gW+Cfq8J7D6Z/+NYvWhZXexJX+S0BVjoDVsUNGypc+4AgHPuCDAWGGtm+YBuwCQzW+qc++YEpj3COXdfFu73o/P+szzfA6cTmBPKDWwz+2spTRKBb55ZsT3o+u8E5mjwxlvfzHYH3Z4LmJzJuCYAvYFhzrlfsjj904OzOuecmWWWvQLwlpkdDRp2hMCH8o/Ouc+8xTenAm9kYfp/m5aZVQFGEljvlJ/Ac17h3VwO+C4L40zPmdn78rfnTeD9zJBz7lszu4NAIVQ3s9nAXc65n459Ds65oxZYQX86gQ/s0495D5P57994dp4PwE7n3J/ZuD/AjqDrf2Twe0ESkNZJxJ70kmjMf/+BFgYNW5DRg5xzfzjnxgK/AdXCnPEMC/q0AcoTmLv4gcA31pLOuaLepbBzrvpJTu8HYH7QOIu6wAr4WzO6s7esfQLwL+A2y/qWMdsIfFilj8eCfw+Rq+UxufI65370Ht+LwOKQn4C7gx4X6tDMxw4fR2AxUGXnXGEC6xrSX/cfCCyyycp4jve+/O15E3g/Q3LOveqca0SgfByBxT7pgl+/JAKL4tL/NjYd81oVcs5dHpTxLLLu2Of4HwJFmj7tMtkYV0JTScSeBQQWKzUhsBgJAotNKgJNCSoJM7vDW4mXzwKbxHYhsGz68zBnPBXoY2a5zawdcA4wyzm3DfgQeNLMCnsrds8yswtPcnozgSpmdoM3zdxmVs/Mzglx/3sJfIh0A54A/hW0kjbJzPIS+GZtZpY3aGXwewS+HV9jgS1m+gCZfdg8BzxiZhW8cZcys6u861UILOfuRGBR3N1mVtt73A6ghJkVOc7zLgTsBfZbYHPk4FKcCZzm/Q3kMbNCZlY/aPwp3oc0WXhf3iDwfpY1s2LAgFCBzOxsM2tmZnkILJ78g8B6tHR1g16/OwiU0xICiyT3eSub83krqmuYWT3vcROBh8yssgWca2Ylgp7P8fanWU3gvavtvb9DjnN/8agkYkdgIapz64GdwHbn3G5v2FEC/2SFgU+DHvM78CSBxTS7CKyfaOPCv+/BZ0Blb5qPAG2DFul0JrBi8msCczXTCSyrP2HOuX0ElpN3IPCtdDv/XWn5N2ZWF7gL6OwtjnucwGub/sHXhMAH2ywC35j/IPABinNuF9AOeIzASvTK/LeoMzKawLqbD81sH4EPw/reB+QU4HHn3Grn3AYCxTXZzPI459YCU4GNFthy6fQQ4+8HXE9gZfrzwOvHvCaXAmne67GBwJcIgGnez1/MbKV3PbP35XlgNoEP2pXA/2XynPMQeH128d+NAwYG3f42gXVIvxEox2ucc4e896IVUJvACuNdBIohvShHEiirDwkU4wsEVkhD4AP/Ze+1ynCrNu//5kHgI++1WJTR/eR/2d8XHUs08v6RH3TOzfA7y/GYWVcCK3Ib+Z1FoouZDSGworiT31kk6zQnEeXMrDqBxTXhXkQkIvI/VBJRzMweJzB7fY9zLuQWJSIi4aLFTSIiEpLmJEREJCSVhIiIhBRXe1yXLFnSpaSk+B1DRCSmrFixYpdzrlRGt8VVSaSkpLB8+XK/Y4iIxBQzC7lhjBY3iYhISCoJEREJSSUhIiIhqSRERCQklYSIiISkkhARkZBUEiIiEpJKQkREQlJJiIhISCoJEREJSSUhIiIhxdWxm+LFn3/+yVdffcVpp53GV199BUDBggU5//zzSU5O9jmdiCQSlUQUevzxxxkyZMj/DL/yyiuZOnUq+fPnj3woEUlIKokoNGPGDGrUqEG3bt2oXbs2uXPn5tNPP2XAgAH079+fsWPH+h1RRBKESiJKHD16lMOHD7Nz505WrVrFY489xp133vnX7Y0aNWLLli0899xz3HbbbVSvXt3HtCKSKLTiOgpMnjyZokWLUr58eV555RUAWrVq9T/3GzJkCIUKFaJhw4Y89dRTHD16NNJRRSTBxEVJmFmamU3Ys2eP31Gybe3atfTs2ZOUlBR27NjBwIEDqVOnDtWqVfuf+5YsWZJPPvmECy64gLvuuos2bdr4kFhEEklclIRz7l3nXI8iRYr4HSXb+vfvT548eZg9ezbNmjXj6NGjjBgxAjPL8P7VqlXjvffeY8CAAcyYMYPvvvsuwolFJJHERUnEqk8//ZSZM2dy9913c9ppp/Hss8/y/PPP06xZs0wfZ2b07NkTCKzkFhEJF3PO+Z0hx6SmprpYOce1c46mTZuydu1avvvuOwoUKJDtcdSpU4cCBQqwaNGiMCQUkURhZiucc6kZ3aY5iQg7dOgQ//znPxk0aBDz58/nvvvuO6GCAGjdujWffvopvXr1okmTJrz66qs5nFZEEp3mJCJs/Pjx3HLLLQCce+65LFu2jFNOOeWExrV3715uueUWpk6dSqFChThw4AALFiygfv36ORlZROJcZnMSKokIO/PMMylevDjnnHMOAwcOzHArpuzavXs3R44coW7duhQpUoRVq1aFXPEtInIsLW6KEj/88AObNm2ic+fOTJ48OUcKAqBo0aKUKFGCwYMHs2bNGubMmZMj4xURUUlE0CeffAJAw4YNwzL+66+/ntNOO40nnngiLOMXkcSjkoighQsXUqBAAWrVqhWW8efJk4e77rqLjz76iPnz54dlGiKSWFQSYfbbb7/Rrl077rjjDiZNmkTLli3JlSt8h8zq1asXZcuWpX///sTT+iYR8YdKIoz27dtHmzZtmD59OqNHj6ZkyZI8/fTTYZ1mvnz5eOihh1i2bBmTJ08O67REJP5p66aTsH79eubPn0+uXLmYO3cuI0aMoHTp0pgZe/bs4eKLL2bVqlVMmDCBM844g0aNGp3wPhHZcfToURo2bMjGjRtZu3YtxYoVC/s0RSR2ZbZ1kw4VfoL27dtHo0aN2Llz51/DXnnlFSpUqMA333zDiBEjWLFiBe+88w5paWkRzZaUlMS4ceM477zzGD16dIYnMBIRyQotbjpBkydPZufOnTzwwANcffXVvPbaazRo0IDvv/+eSy65hEcffZT27dtHvCDS1a5dm5YtWzJ+/HgOHjzoSwYRiX1a3HSCLrvsMjZt2sSGDRv+Guac49xzz+XLL7/k0ksvZdq0afh5ZNr333+fyy+/nKlTp9KhQwffcohIdNPOdDls8uTJzJkz53/O52BmzJw5kz59+jBlyhRfCwKgefPmVK5cmSeffFJbOonICdGcxAmoWbMmefPmZfHixWHdnDUnpB8rau7cucc9BLmIJCbNSeSgI0eOsG7dOi6++OKoLwiALl26ULp0aYYPH+53FBGJQSqJbPr+++85dOgQlStX9jtKluTNm5cOHTqwYMECnRNbRLJNJZFN6YuzcurgfJFw9tln88cff7Bt2za/o4hIjFFJZNNbb73Fqaeeyvnnn+93lCyrVKkSwN+2xBIRyQqVRDY455g/fz6XXHIJycnJfsfJsvRFYyoJEckulUQ2bN26lW3btnHBBRf4HSVbypUrxymnnMK3337rdxQRiTEqiWxI/5A9++yzfU6SPcnJyZx55pl/zUlcc8019OvXz+dUIhILon8bziiyadMmIHAK0lhTuXJlNmzYwI4dO3jrrbdISkqie/fuVK1a1e9oIhLFNCeRDZs2bSIpKYly5cr5HSXbKlWqxHfffcesWbOAwEEA77vvPp9TiUi0U0lkw8aNGylfvjy5c+f2O0q2ValShT/++INRo0Zx2mmnMWjQIN58800WLlzodzQRiWIqiWzYtGkTFStW9DvGCWnTpg0pKSmsWbOGli1b0r9/f8qWLUufPn04cuSI3/FEJEqpJLJh48aNMbk+AqBUqVJ8+OGHNG7cmB49elCgQAGeeOIJVq1axQsvvOB3PBGJUiqJLPr999/ZsWNHzM5JQGDl9YIFC6hfvz4A1157LY0bN2bQoEH89ttvPqcTkWikksiiWN6yKRQzY/To0fzyyy+MHDnS7zgiEoVUElmUXhKxPCeRkTp16tCqVSvGjx/Pn3/+6XccEYkyKoks2rhxIxBfcxLpbr/9dnbu3MnUqVP9jiIiUUYlkUWbNm0if/78lCpVyu8oOa5Zs2bUqFGD0aNH6wx2IvI3UVMSZnammb1gZtODhhUws5fN7Hkz6+hnvvQtm8zMzxhhYWb06dOH1atXs2DBAr/jiEgUCWtJmNkkM/vZzL48ZngLM1tnZt+a2QAA59xG59xNx4ziGmC6c+5m4MpwZj2eWN5HIis6depEiRIlGD16tN9RRCSKhHtO4iWgRfAAM0sGxgItgWrAdWYW6gw+ZYEfvOu+7fHlnIvpfSSyIl++fPTo0YO3336bzZs3+x1HRKJEWEvCObcA+PWYwecD33pzDgeB14CrQoxiK4GigBBZzayHmS03s+U7d+7Midj/Y9euXfznP/+J6zkJgNtuuw0z45lnnvE7iohECT/WSZzBf+cOIFAEZ5hZCTN7DqhjZgO92/4PaGNm44B3MxqZc26Ccy7VOZcarpXK8bxlU7CyZcvStm1bJk6cyJ49e/yOIyJRIGpWXDvnfnHO3eKcO8s596g37D/OuRudc7c6517xK1u87iORkX79+rF//35atGjB7t27/Y4jIj7zoyR+BIKPtV3WGxa1tm7dChCThwjPrtTUVKZNm8aKFSu4+OKL+eWXX/yOJCI+8qMklgGVzayimZ0CdADe8SFHlu3evZukpCQKFy7sd5SIaN26NTNmzOCrr74iLS1NR4kVSWDh3gR2KrAYONvMtprZTc65w0BvYDbwDfCGc+6rcOY4Wbt376ZIkSJxuY9EKJdffjkTJ05k8eLF2ixWJIFZPOxha2ZpQFqlSpVuTj+Pc0664YYb+OSTT/5agZ0onHNcffXVzJkzh/Xr11O2bNnjP0hEYo6ZrXDOpWZ0W9SsuD4Zzrl3nXM9ihQpEpbxp89JJJr0o8QeOXKEwYMH+x1HRHwQFyURbrt376Zo0aJ+x/BFSkoKvXv35qWXXmLu3Ll+xxGRCFNJZMGePXsStiQA7rvvPqpVq0arVq10TmyRBKOSyIJEXdyUrlixYnz88ceUKlWK+++/3+84IhJBKoksSOTFTelKlizJP//5T+bPn8+aNWv8jiMiEaKSOI4jR46wb9++hC8JgG7dupE3b15tEiuSQOKiJMwszcwmhON4Q3v37gVI6MVN6UqUKMHNN9/MSy+9pLkJkQQRFyURzk1g169fD6A5Cc+QIUMoVqwYvXr14tChQ37HEZEwi4uSCKcGDRoAcPbZZ/ucJDoUL16cUaNGsWjRIjp16qRDdojEOZVEJoL3Rq9du7aPSaJLp06dePzxx3njjTd47rnn/I4jImGkksjEzz//DMBDDz1E/vz5fU4TXfr3788ll1zCoEGD2LFjh99xRCRMVBKZSD9WU506dXxOEn3Sz2C3b98+xowZ43ccEQkTlUQm0k82FO9npDtRZ599Ns2aNeP1118nHg4UKSL/Ky5KIlybwKbPSaSkpOToeONJhw4d+O6771i5cqXfUUQkDOKiJMK1CezGjRs5/fTTyZcvX46ON560bt2aXLly8corvp1dVkTCKC5KIlw2btyYEOe1PhnFixenbdu2PPfcc3+d5lVE4odKIhMbN27U+ogsGDZsGEeOHOGee+7xO4qI5DCVRAgHDhxg69atKoksqFixIgMGDODVV1/VcZ1E4oxKIoTt27fjnKNcuXJ+R4kJDzzwAK1bt+bOO+/kiy++8DuOiOQQlUQI6TuIlS5d2ucksSE5OZmJEydSqFAhnepUJI6oJEJI39v61FNP9TlJ7ChevDh9+/blrbfeYsWKFX7HEZEcEBclkdP7Saxbt44lS5YAKonsuuOOOyhevDgPPPCA31FEJAfERUnk9H4SVatW5ZFHHgFUEtlVuHBh7r77bmbNmsWsWbP8jiMiJykuSiJcChYsqAP7nYDevXtz5plncsUVV3D77bfrkB0iMUwlkQnNRZyYAgUKsHr1anr37s2YMWPo16+f35FE5ATl8jtANNPe1ieuYMGCjBkzhiNHjjBy5EjatWv31wmcRCR2aE4iE1WrVvU7QkwzM4YPH85pp51Gnz59OHr0qN+RRCSbVBKZ0ClLT17BggUZPnw4y5YtY+zYsX7HEZFsUklkQiWRMzp27EjLli0ZMGAAq1at8juOiGSDSiITWtyUM8yMCRMmUKJECS666CI+//xzvyOJSBapJDJRtmxZvyPEjbJly7Jo0SIKFixI586dOXjwoN+RRCQLVBIZaNCgAZdddhlJSXp5clL58uUZP348X375JcOGDfM7johkQVx8Cobr9KWS86644go6duzII488oqPFisSAuCiJcJ2+VMJj1KhRFCtWjG7dunH48GG/44hIJuKiJCS2lCxZkqeffprly5fTp08fpk2bxt69e/2OJSIZUEmIL9q3b0+7du0YN24c7du354wzzmD+/Pl+xxKRY6gkxBdmxuuvv86mTZtYtGgR5cqV48orr+Trr7/2O5qIBFFJiG/MjJSUFBo2bMjs2bPJkycP119/PQcOHPA7moh4VBISFcqVK8cLL7zA6tWrGT16tN9xRMSjkpCokZaWRrNmzRgzZgyHDh3yO46IoJKQKHPnnXfy448/8thjj2mvbJEooJKQqHL55ZfTrFkzHnjgAVq3bq2z2on4TCUhUSUpKYmPPvqI4cOHM2vWLMaPH+93JJGEppKQqGNm9O3bl0suuYS7776bbdu2+R1JJGHFRUno2E3xJykpiXHjxnHw4EHq1q1Lw4YNWb9+vd+xRBJOXJSEjt0UnypVqsSIESOoVKkS69ev5x//+Afr1q3zO5ZIQomLkpD41bt3bxYsWMDixYtJTk6mVatW/PLLL37HEkkYKgmJCZUqVWLGjBls2bKFtm3bavNYkQhRSUjMuOCCC5g0aRLz5s3j1ltv1eaxIhGgkpCY0rFjR+677z4mTZpEo0aN+P777/2OJBLXVBISc4YOHcr48eP54osv6Ny5M0ePHvU7kkjcUklIzElKSqJHjx6MHj2aBQsW8Oijj/odSSRuqSQkZnXt2pXrr7+e++67j7ffftvvOCJxSSUhMcvMmDhxIvXq1aNjx46sWLHC70gicUclITEtX758zJgxg+LFi9O4cWNeeeUVvyOJxBWVhMS8008/naVLl1KvXj06derEqFGj/I4kEjdUEhIXypQpw0cffUTLli0ZPHgwv/76q9+RROKCSkLiRu7cuXn88cfZt28fTzzxhN9xROKCSkLiSs2aNenQoQNjxoxhx44dfscRiXkqCYk7Q4YM4cCBA7Rt25bPP//c7zgiMS1kSZhZ8SxcikYybCg6n4QEq1KlCuPGjWPdunVccMEFzJo1y+9IIjHLQh0kzcz+BH4CLJPHJzvnyocj2IlITU11y5cvP+nx/OMf/6Bw4cLMnj07B1KJX3bu3EnTpk05cOAA69evxyyzP2WRxGVmK5xzqRndltnipm+cc2c65yqGugA6sL9ErVKlSnHPPffw7bffMm/ePL/jiMSkzEriH1l4fFbuI+Kbtm3bUqxYMfr06cPXX3/tdxyRmBOyJJxzfwKY2ZNmVj2z+4hEq3z58vHKK6+wY8cOmjVrxubNm/2OJBJTsrJ10zfABDP7zMxuMTOdSFpiSsuWLZk/fz5//vkntWvXZvr06X5HEokZxy0J59xE51xDoDOQAqwxs1fNrGm4w4nklHPOOYclS5ZQsWJF+vTpw6FDh/yOJBITsrSfhJklA1W9yy5gNXCXmb0WxmwiOapq1ao88sgjbNu2jauvvpoHHniAbdu2+R1LJKodtyTM7ClgHXA5MMw5V9c597hzLg2oE+6AIjmpRYsW1KxZk/nz5/Pwww9z3nnnsWTJEr9jiUStrMxJrAFqOed6OueWHnPb+WHIJBI2SUlJLFu2jF9//ZU1a9aQP39+LrroIt59912/o4lEpVyhbjCz87yrq4Gzj9kR6QCwxTmnXZwl5uTJkweAGjVqsHTpUpo3b06XLl344osvOOOMM3xOJxJdQpYE8ORxHlfezMY654bncCaRiClRogRTp06ldu3a3HjjjXzwwQckJemQZiLpQpaEcy7TrZfMLA/wOaCSkJhWuXJlnnzySW699Vaef/55evbs6XckkaiR2QH+zgt1G4Bz7gBwQ44nEvFBz549ufDCC7nvvvvYvXu333FEokZm89UvmlmxzI4CC0yMVFCRcDIznnrqKXbt2sW4ceP8jiMSNTJbJ1EEWEHmR4HdmbNxRPxTp04dGjZsyGOPPcbnn39Oo0aN6NWrF8nJyX5HE/FNZuskUiKYQyQqdOzYkdtuu41p06Yxbdo0kpOT6dWrl9+xRHyjzThEgnTo0IG0tDQ++ugjLrzwQh588EF++UVHxJfEldniJpGEU6xYMd555x0AihQpQqNGjWjevDnPPvss55+vfUcl8WhOQiSE1NRUpk2bxqZNm2jcuLHORyEJKSvHbvo/M7vCzFQoknDS0tL4+uuvKVSoEF26dOH333/3O5JIRGXlg/9Z4Hpgg5k9ZmZnhzmTSFQpXbo0kyZNYsWKFXTv3t3vOCIRlZXzSXzknOsInAdsBj4ys0/N7EYzyx3ugFlhZmlmNmHPHh1KSsLjyiuv5P7772fq1KksX77c7zgiEZPV80mUALoC3QkcimM0gdKYE7Zk2eCce9c516NIEZ00T8Knb9++FC9enMGDB/sdRSRisrJO4i1gIZAfSHPOXemce90590+gYLgDikSLwoUL079/f2bNmsWHH37odxyRiMjKnMTzzrlqzrlHnXPb4K+D++GcSw1rOpEo07t3b0477TSaN29O9+7dOXLkiN+RRMIqKyXxcAbDFud0EJFYULBgQVauXMldd93FCy+8wMMPZ/TvIRI/MjvpUBngDCCfmdXhv8dwKkxg0ZNIQipTpgxPPvkk27dv55FHHqF58+Y0aNDA71giYZHZHtfNCaysLguMDBq+D7g3jJlEYsLo0aNZsmQJzZs3Z86cOdojW+JSyMVNzrmXvRMPdXXONQ26XOmc+78IZhSJSiVLlmTevHmULFmSpk2b8tBDD+Gc8zuWSI7KbHFTJ+fcFCDFzO469nbn3MgMHiaSUMqVK8eCBQvo3bs3DzzwAOXLl6dLly5+xxLJMZmtuC7g/SwIFMrgIiLAGWecwfTp02ncuDF33XUXe/fu9TuSSI7J7HwS472fQyMXRyQ2JScnM3LkSOrVq8fIkSMZMmSI31B9p7gAABu+SURBVJFEckRWdqYbbmaFzSy3mc01s51m1ikS4URiSWpqKu3bt+fRRx/l1Vdf5cCBA35HEjlpWdlP4jLn3F6gFYFjN1UC+oczlEiseuaZZyhTpgwdO3akTZs2HD582O9IIiclKyWRvkjqCmCac05H0RMJoVSpUqxdu5aRI0fy3nvv8eqrr/odSeSkZKUkZprZWqAuMNfMSgF/hjeWSOzKly8fd9xxBxUrVmTy5Ml+xxE5KVk5VPgA4AIg1Tl3CPgPcFW4g4nEMjOjU6dOzJ07ly+++MLvOCInLKtnm6sKXGtmnYG2wGXhiyQSH2677TZKly5NamoqtWrV4qOPPvI7kki2ZWXrpsnACKARUM+76OivIsdRpkwZZs+eTY8ePThw4ACtWrXixRdf5OjRo35HE8myzI7dlC4VqOZ0vAGRbDv33HN5+umn2bVrF5deeindunXj119/pW/fvn5HE8mSrCxu+hIoE+4gIvGsZMmSrFixgmbNmjFy5Ej279/vdySRLMlKSZQEvjaz2Wb2Tvol3MFE4k1SUhKDBg3ip59+olatWuic7BILsrK4aUi4Q4gkimbNmjF79mxatGjBI488wvDhw/2OJJKprGwCO5/Anta5vevLgJVhziUSty677DJuvPFGnnjiCSpWrMj06dP9jiQSUla2broZmA6M9wadAcwIZyiReDd27FjatGnD5s2bGTVqFKtWrfI7kkiGsrJOohfQENgL4JzbAJwazlAi8S5v3rxMnz6dQYMG8cknn1CnTh3tRyFRKSslccA5dzD9FzPLBWhzWJEc0Lp167+uDxs2jC1btujsdhJVslIS883sXiCfmV0KTAPeDW8skcRQt25dli5dyogRI/j444+pUKECTZs2Zffu3X5HEwGyVhIDgJ3AF0BPYBZwXzhDiSSSevXq0bdvX5YtW8awYcP45JNPuOeee/yOJQJkYRNY59xRM5sBzHDO7YxAJpGElJqaSmpqKj///DNjxozhtttuo1atWn7HkgQXck7CAoaY2S5gHbDOOyvdA5GLJ5J4HnjgAYoWLcrtt9/OkSNH/I4jCS6zxU13EtiqqZ5zrrhzrjhQH2hoZndGJJ1IAipWrBjDhw9n/vz53Hjjjfz+++9+R5IEltniphuAS51zu9IHOOc2eue3/hB4KtzhRBLVTTfdxI8//sjgwYPJmzcvEyZM8DuSJKjM5iRyBxdEOm+9RO7wRRIRCCx2uvXWW3n55ZfZvn2733EkQWVWEgdP8DYRySF33nknBw8e5KWXXvI7iiSozEqilpntzeCyD6gZqYAiiaxy5cpccMEFTJo0SScrEl+ELAnnXLJzrnAGl0LOuYgsbjKzamb2hpmNM7O2kZimSLTp0qULGzZsoGbNmmzZssXvOJJgsnqO6xxjZpPM7Gcz+/KY4S3MbJ2ZfWtmA7zBLYGnnXO3Ap0jnVUkGnTv3p2XX36ZH374gb59+/Luu+/SpEkTbfUkERHxkgBeAloEDzCzZGAsgVKoBlxnZtWAyUAHM3sCKBHhnCJRISkpic6dO9O0aVPWr1/P8OHDWbhwIU89pQ0MJfwiXhLOuQXAr8cMPh/41jm30TuY4GvAVc65n51zvQgcGuR/trQSSSQVKlTg+++/x8wAePDBB3n22Wd59dVXNVchYXPcw3KYWTXn3NfHDLvIOTcvB3OcAfwQ9PtWoL6ZpQD3AgWAJ0Lk6wH0AChfvnwORhKJLuXLl2fPnj0sXbqUq6++mm+++YZevXoB0KhRIxYuXOhzQolHWZmTeMPM7vEO05HPzJ4GHg13MADn3GbnXA/nXEfn3KIQ95ngnEt1zqWWKlUqErFEfJH+JejAgQM0bNiQNWvWMH/+fO69914WLVrE559/7nNCiUdZKYn6QDngUwKnLv2JwOE6ctKP3jTSlfWGiYinQoUKf12vWrUqp5xyCk2aNKFfv37kzZuX66+/nnnz5vkXUOJSVkriEPAHkA/IC2xyzuX0BtvLgMpmVtHMTgE6AO/k8DREYlrw4tSmTZv+db1YsWK8/fbbHDhwgEsuuYQNGzb4EU/iVFZKYhmBkqgHNCaw5dG0E52gmU0FFgNnm9lWM7vJOXcY6A3MBr4B3nDOfXWi0xCJR6VLl6Zjx47MmTOHAgUK/O22yy67jE8//ZTk5GRGjx7tU0KJR8ddcQ3c5Jxb7l3fBlxlZjec6ASdc9eFGD6LwAmNss3M0oC0SpUqnWgskaiXlJTElClTQt5epkwZbrjhBsaOHUvZsmXp3Lkzp59+egQTSjzKypzEz2ZWPvgCzA93sOxwzr3rnOtRpEgRv6OI+GrIkCHky5ePgQMHUr58eapXr87UqVP9jiUxLCtzEu8BDjAC6yQqEjgJUfUw5hKRE1C2bFlWrVrFjh07eO+995g8eTJdu3albt26VKlSxe94EoOOOyfhnKvpnDvX+1mZwI5vi8MfTURORJUqVWjcuDGPPfYYK1as4OjRo7zwwgt+x5IYle09rp1zKwlsFisiUa5MmTJcfvnl/Otf/+LgQR3hX7LvuCVhZncFXfqZ2asE9pUQkRhw6623sn37diZPnux3FIlBWVknUSjo+mEC6yjeDE8cEclpzZs3p27dutx8882YGd26dfM7ksSQ45aEc25oJIKcDG0CKxKamTF37lzatWvHLbfcwqFDh+jevTvJycl+R5MYYM65jG8we5fAVk0Zcs5dGa5QJyo1NdUtX778+Hc8jn/84x8ULlyY2bNn50Aqkejw22+/0a5dO+bOnUutWrVo3749rVu35pxzzvE7mvjMzFY451Izui2zOYkRYcojIj4oVqwYH374IY8//jj33nsvq1evZuTIkXz22WecddZZfseTKJVZSWxyzulciSJxJCkpiYEDB9KmTRsOHjzIhRdeSFpaGosXL0Y7o0pGMtu6aUb6FTPTimqROFKlShVq1KjB9OnT2bBhA82bN2fTpk1+x5IolFlJWND1M8MdREQir2nTprzyyit88803dOvWjVDrKCVxZVYSLsR1EYkj7du3Z9iwYcybN49bb73V7zgSZTJbJ1HLzPYSmKPI513H+9055wqHPV0WaRNYkZNz6623sn79esaMGUP16tXp1asXSUnZPiCDxKGQfwXOuWTnXGHnXCHnXC7vevrvUVMQoKPAipyspKQkHnvsMc466yz69OnDtddey9GjOX1uMYlF+qogIgDky5eP1atXM3ToUKZPn06rVq3YsWOH37HEZyoJEflLgQIFuP/++7nlllt4//33ueeee/yOJD5TSYjI35gZ48aN484772Ty5MmsW7fO70jiI5WEiGRowIAB5M2bl06dOmkfigSmkhCRDJ166qm8+OKLbNiwgYYNG7Jnzx6/I4kPVBIiElL79u2ZPXs227ZtY/z48X7HER9k5XwSIpLA6tevz6WXXsqoUaMoU6YMBQoUoE2bNn7HkgiJizkJM0szswmaHRYJj7vvvptt27bRpUsX2rZty9dff+13JImQuCgJ7UwnEl4XX3wxL730Ev369aNQoUJUr16d0qVL8/777/sdTcJMi5tE5LjMjC5dugDQoEEDBg0axI4dO+jZsyfffPMNBQoU8DmhhEtczEmISOS0adOGtWvX8u677/LDDz9w00038e233/odS8JEJSEiJ6RRo0YMGDCA6dOnU6NGDRYvXux3JAkDlYSInLBHH32UzZs3c/rpp9O9e3cdFDAOqSRE5KSULVuWhx56iK+//poHH3yQgwcP+h1JcpBKQkROWocOHbj22msZOnQoXbt29TuO5CCVhIictOTkZKZOnUq/fv2YOnUqn332md+RJIeoJEQkR5gZgwcPpnTp0txyyy3MnDmTDz74wO9YcpLioiS0x7VIdChYsCDjxo1j7dq1pKWl0bJlS6ZPn+53LDkJcVES2uNaJHq0bt2apUuX0rlzZypWrEi7du249957/Y4lJyguSkJEokvNmjV5+eWXWblyJQ0bNuSpp55i7969fseSE6CSEJGwKVq0KCNGjODPP/+kf//+fseRE6CSEJGwql+/Pj179mTChAls27bN7ziSTSoJEQkrM6Nz584AzJo1i3nz5nH48GGfU0lWqSREJOzq1KlDcnIy3bt3p2nTptx4441+R5Is0qHCRSTs8uXLR69evViyZAkFChRgypQp9O3bl0OHDpGSkkKpUqX8jighqCREJCJGjx4NwK5duyhTpgwNGjTgwIEDNGnShHnz5mFmPieUjGhxk4hEVMmSJWnevDkHDhwAYMGCBUydOtXnVBKKSkJEIm7y5MkMHDiQjRs30qBBA3r27Mlrr72Gc87vaHIMlYSIRFzx4sUZNmwYFStWZMqUKeTNm5frrruOfv36+R1NjhEXJaFjN4nErrPOOosffviB66+/nmeeeYaFCxf6HUmCxEVJ6NhNIrEtb968DBs2jFNPPZUmTZowefJkvyOJJy5KQkRiX4UKFVi3bh01atRg7NixfscRj0pCRKJG/vz5ufHGG/nss8+oVq2ayiIKqCREJKr07t2bK664gu3bt9O7d28WLVrkd6SEppIQkahyyimnMHPmTLZu3Urp0qXp168fw4YNY/369X5HS0gqCRGJSvnz5+ef//wnn332GYMGDeK8887j3//+t9+xEo5KQkSi1j333MP69evZtGkTKSkppKWlsWnTJr9jJRSVhIhErVy5clG5cmVSUlJ4//33ATjzzDNp0KABW7du9TldYlBJiEhMKFeuHPPmzeOWW27hiy++oF27duzcudPvWHFPJSEiMaNevXqMGzeOl156iZUrV3Leeeexe/duv2PFNZWEiMScdu3aMXfuXLZu3cr999/PkSNH/I4Ut1QSIhKTGjVqRI8ePXjmmWeoW7cuGzdu9DtSXFJJiEjMeu6553j55ZdZvXo1o0aN8jtOXFJJiEjMMjM6d+7MVVddxZtvvsnhw4f9jhR3VBIiEvO6devGTz/9xOuvv+53lLgTFyWh80mIJLZWrVpx7rnncvfdd/P++++zZcsW7Z2dQ+KiJHQ+CZHElpSUxMSJE/n999+5/PLLqVChAhdffDHvvPOO39FiXlyUhIhIvXr12Lp1KzVq1Phr2G233cYHH3zgY6rYp5IQkbhRoEABFi1axK5du1i6dCm5c+emZcuWrFy50u9oMUslISJxpUiRIpQoUYJ69eqxcuVKcufOTY8ePXjiiSe09dMJUEmISNwqVqwYd9xxBz///DN333031atX58cff/Q7VkxRSYhIXBs+fDhbtmzhyiuvZP369Vx33XU45/yOFTNUEiKSECZPnkyfPn1YuHAhr7/+uooii1QSIpIQChcuzIgRI6hWrRrXXXcdHTp04OjRo37HinoqCRFJGLlz52bRokX07t2bN954g6FDh/odKerl8juAiEgkFStWjDFjxrB//34efPBBmjZtykUXXeR3rKilOQkRSThmxrPPPkv58uW54oorePnll/2OFLVUEiKSkPLly8f48eP5/fff6dq1K0OGDOG7777zO1bUUUmISMJq0aIFa9asoVChQgwdOlTrKDKgkhCRhFazZk22bdtG69atmTdvnjaNPYZKQkQSXoECBbjkkkv44Ycf/joN6v79+31OFR1UEiIi8NcWTpUqVcLMKFSoEC+//DK///67v8F8ppIQEQHOOeccSpcu/bdhXbt25dxzz2Xfvn0+pfKfSkJEhMBmsYsXL2bKlCmsW7eOGTNm0KlTJ7777jvatm3L6tWr/Y7oC+1MJyLiqVixIhUrVgSgSpUqXHXVVZQvX55hw4YxZ84c+vfvz4MPPkiePHl8Tho5mpMQEcnEI488wq5du7j55psZPnw4efPmZdKkSX7HihiVhIjIcZQoUYLx48czePBgAG666SYWL17sc6rIUEmIiGTRkCFD2LVrF6VLl+aCCy5IiMN5xEVJmFmamU3Ys2eP31FEJM6VKFGC4cOHA4Gtn9Kvx6u4KAnn3LvOuR5FihTxO4qIJIDOnTtz4MABrr32Wu655x5uuOEGDh065HessIiLkhARibRTTjmF8ePHAzBlyhRmzpzpc6LwUEmIiJygIkWK8M477wBwzTXX8OSTT8bdsZ9UEiIiJyEtLY0XX3wRgH79+lG/fn22b9/uc6qco5IQETlJXbt25eDBg4wYMYJly5bRoEEDXnjhBb9j5QiVhIhIDsidOzd9+/Zl2LBhfP/993Tv3p3Vq1fH/AECVRIiIjlo4MCBfP/99+TJk4fatWvToEEDtmzZ4nesE6aSEBHJYeXLl2fq1Klcc801rFu3jqpVq5KWlsavv/7qd7RsU0mIiIRB69atefPNN1m3bh2XXnopM2fO5PHHH/c7VrapJEREwiglJYW3336bFi1a8N577/kdJ9tUEiIiEdCoUSO++uorLrroIn7++We/42SZSkJEJAJuuOEGqlevzvz58xk0aJDfcbJMJSEiEgHly5fnyy+/pHv37kycOPGvPbWjnUpCRCSCBg0aRPXq1Wnbti0ff/yx33GOSyUhIhJBKSkpLFq0iLJlyzJo0KCoP9aTSkJEJMKKFi1K7969Wbx4MQ8++KDfcTKlkhAR8cHtt99Ohw4dePjhh6N6JzuVhIiID5KTk7n99ts5fPgwc+bM8TtOSCoJERGf1KtXj0KFCrFw4UK/o4SkkhAR8UlycjK1atVi1apVfkcJSSUhIuKjOnXq8Mknn3DBBRdQu3btqNssViUhIuKjNm3aALB48WL2799P69at+fjjjzl8+LDPyQJUEiIiPmrSpAkTJ05kw4YNPPnkk+zZs4dmzZrRr18/v6MBKgkREV+ZGTfddBOVKlXi8ssv55FHHqFMmTKMHTuWpUuX+h1PJSEiEi1y587Nvffey+eff07p0qXp0aOH35FUEiIi0aZMmTIMHDiQ1atX88EHH/iaRSUhIhKFunXrRqVKlRg6dKivOVQSIiJRKF++fNx8880sWbKE119/3bccKgkRkSjVu3dvGjZsSPfu3f86vtOSJUswM7788suIZFBJiIhEqfz58/P000+zf/9+Xn31VQCmTJkCELGd7lQSIiJRrE6dOtSpU4fnn3+eo0ePsnfvXgAKFy4ckemrJEREotxdd93FmjVrmDhx4l8lUaRIkYhMWyUhIhLlOnbsSJUqVZg+fTp79uwBoFChQhGZtkpCRCTKmRmtWrVi3rx5/PTTTwDkyZMnItOO6pIws/JmNsPMJpnZAL/ziIj45YorruDQoUOsX78+otONeEl4H/g/m9mXxwxvYWbrzOzboEKoCUx3znUD6kQ6q4hItLjwwgspV65cxKfrx5zES0CL4AFmlgyMBVoC1YDrzKwasAS4ycz+Dfi7b7qIiI+Sk5N9OYNdxEvCObcAOPas3+cD3zrnNjrnDgKvAVcBNwKDnXPNgCsyGp+Z9TCz5Wa2fOfOneGMLiLiqwoVKnDttddGdJrRsk7iDOCHoN+3esM+APqY2XPA5owe6Jyb4JxLdc6llipVKuxBRUT81L1794hOL1dEp5ZNzrkvgbZ+5xARSVTRMifxIxC8RqasN0xERHwULSWxDKhsZhXN7BSgA/COz5lERBKeH5vATgUWA2eb2VYzu8k5dxjoDcwGvgHecM59FelsIiLydxFfJ+Gcuy7E8FnArBMZp5mlAWmVKlU6mWgiInKMaFncdFKcc+8653pE6oBXIiKJIi5KQkREwkMlISIiIakkREQkJJWEiIiEFBclYWZpZjYh/WQcIiKSM+KiJLR1k4hIeMRFSYiISHiYc87vDDnGzHYC3/udIxtKArv8DuETPffEpOcenSo45zI8jHZclUSsMbPlzrlUv3P4Qc9dzz3RxOpz1+ImEREJSSUhIiIhqST8NcHvAD7Sc09Meu4xRuskREQkJM1JiIhISCoJH5jZJDP72cy+9DtLpJlZOTP72My+NrOvzOx2vzNFipnlNbOlZrbae+5D/c4UaWaWbGafm9lMv7NEkpltNrMvzGyVmS33O092aHGTD8ysCbAf+JdzrobfeSLJzE4DTnPOrTSzQsAK4Grn3Nc+Rws7MzOggHNuv5nlBhYBtzvnlvgcLWLM7C4gFSjsnGvld55IMbPNQKpzLlr3kwhJcxI+cM4tAH71O4cfnHPbnHMrvev7CJyu9gx/U0WGC9jv/ZrbuyTMtzQzKwtcAUz0O4tknUpCfGNmKUAd4DN/k0SOt7hlFfAzMMc5lzDPHRgF3A0c9TuIDxzwoZmtMLMefofJDpWE+MLMCgJvAnc45/b6nSdSnHNHnHO1gbLA+WaWEIsbzawV8LNzboXfWXzSyDl3HtAS6OUtco4JKgmJOG95/JvAK865//M7jx+cc7uBj4EWfmeJkIbAld6y+deAZmY2xd9IkeOc+9H7+TPwFnC+v4myTiUhEeWtvH0B+MY5N9LvPJFkZqXMrKh3PR9wKbDW31SR4Zwb6Jwr65xLAToA/3bOdfI5VkSYWQFvIw3MrABwGRAzWzaqJHxgZlOBxcDZZrbVzG7yO1MENQRuIPBNcpV3udzvUBFyGvCxma0BlhFYJ5FQm4ImqNLAIjNbDSwF3nPOfeBzpizTJrAiIhKS5iRERCQklYSIiISkkhARkZBUEiIiEpJKQkREQlJJiISRmc0zsxw9r7GZFTWz24J+vyjRjqoqkaOSEIk9RYHbjnsvkRygkpCEZGb9zayPd/0pM/u3d72Zmb1iZuPMbHnweR/MrIWZTQsax1/f4M3sMjNbbGYrzWyad2yqY6eZ4X28cw0M9YZ/YWZVveGlzGyOl2GimX1vZiWBx4CzvB0Rn/BGX9DMppvZWi+/hfHlkwSikpBEtRBo7F1PJfAhm9sbtgAY5JxLBc4FLjSzc4GPgPreoRUArgVe8z647wMu8Q7ithy4K3hiWbjPLm/4OKCfN2wwgcNXVAemA+W94QOA75xztZ1z/b1hdYA7gGrAmQT2bBc5aSoJSVQrgLpmVhg4QOAwKakESmIh0N7MVgKfA9WBas65w8AHQJqZ5SJwboS3gQYEPpw/8Q4D3gWocMz0jnef9AMdrgBSvOuNCBwMD+8wDr9l8nyWOue2OueOAquCxiFyUnL5HUDED865Q2a2CegKfAqsAZoClYA/CHybr+ec+83MXgLyeg99DehN4KRRy51z+7xFO3Occ9dlMsnj3eeA9/MIJ/Z/eSDo+omOQ+R/aE5CEtlCAmWwwLt+C4E5h8LAf4A9ZlaawDkA0s0HzgNuxvuWDywBGppZJfjrqJ9VjplWVu5zrE+A9t79LwOKecP3AYWy91RFToxKQhLZQgJHZl3snNsB/AksdM6tJlAWa4FXCXxYA4GTBgEzCRTHTG/YTgJzJFO9I7wuBqoGTygr98nAUOAyM/sSaAdsB/Y5534hsNjqy6AV1yJhoaPAikQpM8sDHHHOHTazfwDjvLPaiUSMlluKRK/ywBtmlgQcJLCISySiNCchIiIhaZ2EiIiEpJIQEZGQVBIiIhKSSkJEREJSSYiISEgqCRERCen/AZx3gq8ul7F4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# pipeline 1d extraction (for comparison)\n",
    "jpipe_x1d = Table.read(x1dfile, hdu=1)\n",
    "print(jpipe_x1d.columns)\n",
    "# plot\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "ax.plot(jpipe_x1d['WAVELENGTH'], jpipe_x1d['FLUX'], 'k-', label=\"jpipe_x1d\")\n",
    "ax.set_title(\"JWST Pipeline x1d extracted spectrum\")\n",
    "ax.set_xlabel(\"wavelength\")\n",
    "ax.set_ylabel(\"Flux Density [Jy]\")\n",
    "ax.set_yscale(\"log\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boxcar Extraction\n",
    "\n",
    "Extract a 1D spectrum using a simple boxcar.  Basically collapse the spectrum in the cross-dispersion direction over a specified number of pixels.\n",
    "\n",
    "Limitation: currently it is assumed there are no bad pixels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fixed width boxcar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define extraction parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'slit[0] slice')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAACkCAYAAAAT31biAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAQIUlEQVR4nO3de6wlBX0H8O9v7+7yWFceihQQX5Vg0IbabNW29qGohdYWU22rtSk2GmIirTZ9oaaN/kGiTV82mjZUraQ+0FqJxD9USm21TUuEmrYiGBAfgCJQRFZ0YR+//nEOccXdO4c9h3vmcj+fZHPPnJkz85u5c+bc786c+VV3BwAAgOXatOwCAAAAEM4AAABGQTgDAAAYAeEMAABgBIQzAACAERDOAAAARkA4AwAAGAHhDIClq6o3VNW7p48fU1XfqqqVg0z7uKrq6TTnzjj/d1XVd6rqpkOo7V+q6hXTxy+tqo8/0HkAwCyEMwBGpbu/0t0P6+69yfeGo/s5ursvvG+gqs6oqmur6ttV9Ymqeux+83xZkrMWUNt7uvt5884HAA5EOANg3auqRyb5UJI/SnJskiuTvH+pRQHAAyScAbBmquoPq+rmqtpZVZ+vqjMOMM19ly1urqoLkvxkkrdOL2N860Fm/UtJru7uf+juXUnekOT0qnrSjHUdXlXvrqr/q6o7q+rTVXX8AaZ7WVX9237DT66qy6rqjqr6elW9bvr8pqo6v6q+MJ3nB6rq2FlqAWDjEs4AWBNVdWqS85L8aHdvT/KzSb602mu6+/VJPpXkvOmljucdZNInJ/nv/V53d5IvTJ+fxTlJjkpycpJHJHllku+s9oKq2p7kn5J8NMmJSZ6Y5PLp6N9K8oIkPz0d940kb5uxFgA2KOEMgLWyN8lhSU6rqi3d/aXu/sKC5v2wJN+833PfTLJ9xtfvziSUPbG793b3Vd1918Brnp/klu7+s+7e1d07u/uK6bhXJnl9d9/U3fdkcibvRVW1ecZ6ANiAhDMA1kR3X5/kNZkElVur6uKqOnFBs/9Wkoff77mHJ9k54+v/PsnHklxcVV+tqj+pqi0Drzk5k7NzB/LYJJdML5G8M8k1mYTT77tUEgDuI5wBsGa6+73d/cxMwksnefMsL5thmquTnH7fQFVtS/KD0+dnqWt3d7+xu09L8uOZnBX7jYGX3ZjkCauMO6u7j97v3+HdffMs9QCwMQlnAKyJqjq1qp5dVYcl2ZXJd7r2zfDSr+fgIeg+lyR5SlW9sKoOT/LHSf6nu6+dsbZnVdUPTXur3ZXJZY5DtX0kyQlV9ZqqOqyqtlfV06fj/ibJBffdzr+qjquqs2epBYCNSzgDYK0cluRNSW5PckuSRyV57Qyve0sm39f6RlX91YEm6O7bkrwwyQWZ3Hzj6Ule/ABq+4EkH8wkmF2T5F8zudTxoLp7Z5LnJvmFTNbnuiTP2q/mS5N8vKp2JvnPaU0AcFDVPcvVIgAwDtOzUZ/P5Ozb73f3387wmnck+eUkt3b3Ex/kEgHgkAhnAAAAI+CyRgAAgBEQzgAAAEZAOAMAABiBzWu5sJXt23rzcUev5SIBAABG494vfvX27j7uQOPWNJxtPu7onHTBq9ZykQAAAKPxxV97/ZcPNs5ljQAAACMgnAEAAIyAcAYAADACwhkAAMAICGcAAAAjIJwBAACMgHAGAAAwAsIZAADACAhnAAAAIyCcAQAAjIBwBgAAMALCGQAAwAgIZwAAACMgnAEAAIyAcAYAADACwhkAAMAICGcAAAAjIJwBAACMgHAGAAAwAsIZAADACAhnAAAAIyCcAQAAjMDmZRcArLHqZVcAAIxJ17IrYGrwzFlVvbOqbq2qz+733LFVdVlVXTf9ecyDWyYAAMBD2yyXNb4ryZn3e+78JJd39ylJLp8OAwAAcIgGw1l3fzLJHfd7+uwkF00fX5TkBQuuCwAAYEM51BuCHN/dX5s+viXJ8QuqBwAAYEOa+26N3d1JDnqHgao6t6qurKor9+68e97FAQAAPCQdajj7elWdkCTTn7cebMLuvrC7d3T3jpXt2w5xcQAAAA9thxrOLk1yzvTxOUk+vJhyAAAANqZZbqX/viT/keTUqrqpql6e5E1JnltV1yV5znQYAACAQzTYhLq7X3KQUWcsuBbWodLQeN3xO9uPbcGBaMYK8xk6tnqPjY5Pw/GY+4YgAAAAzE84AwAAGAHhDAAAYASEMwAAgBEQzgAAAEZAOAMAABgB4QwAAGAEBvucwWqOfvi3H/RlzNKXa9Oc/arWYhlJMtTZZfOmffMvY6DOoWVs2bR3cBlD22Jzrb6MTQPjZzHLtto6sC6LqGPL4LrOv98M1TlUwyzzGLKnVwan2Ttn76I9+4aXMWQR23vfGvRgWkSdw8uYf/9ehJWBdd000GFp3+CRc3gemweOBSuZ5T00UOcM+83ugffRllq9zlmOe0Prsoj1GLJ34P/+Z1nGPftW/xN1Xw+fX9g9MM3QPGb7nc63jD0zrMeQPfuG5zG0Ljfedcyq4zevDP9t8IgjVv+bcKjOnmF7Dx0PZvmdzbuPz1TnwDQ3rDLOmTMAAIAREM4AAABGQDgDAAAYAeEMAABgBIQzAACAERDOAAAARkA4AwAAGAHhDAAAYAQ0oWYu27beOzjNUFPkoYaYszRrHWpAujLQuHOWZSyisfJg8+aBOmdpnDi0jCNWds81Pkk2DzRKPWzTnlXHDzVaTYabYR+5aXjfO2zT6usyVMcsjWmP3HTP6vNYQLPhLbX69lwZ2P/XymCdA++RvTM0Yx2ax5DdPfyxN1TH0Homw3UO/c42zbDvbR3Yf4eOSbM08R2ax9YMv5e3D7xXh5qoz9Lc/PDBRu1Drx9expZafXutzNAse1evvr2GlnHUpiMGlzEGe3v138c9PfweWoRv9eqfAXt79ffhlhn2i90D89g1cHjeNcP78J6B5uV395bBeewamGbXiQPj920dXMaQoSbs9w6MT5JdvXod9+ybZVus/jkwNI9ZPquG/Psq45w5AwAAGAHhDAAAYASEMwAAgBEQzgAAAEZAOAMAABgB4QwAAGAEhDMAAIAR0OeMB90ddx+56vihDk3DXUaGDfX+GurFNss0mxZQ6NAyeoZ+P0PzGOr5toi+XEP92mbZ3oM94Wbod7UmPfYWsL0Gl7EGfcz2LeCdNtT/brCGBfSOGVzGDOs51E9wEb/zRfxOZ+mtuJpFbO9Zahg6Hgxtz0X0dxza3rO91+fb3sn823yoh2SS7F6D99GQofWc5Xe6CPfuW71v1iLe63uGljFwzJllW+zZt/r2nOVvgz1z7hd7B2pIhtdlaPwsR8XBecywLeadx76FfCRfetAxy38HAwAAIJwBAACMgXAGAAAwAsIZAADACAhnAAAAIyCcAQAAjIBwBgAAMAL6nDGXm247ZnCaWXpaPRRslPVMkqyTda21aaUDAAfU6+PjMlmj3nPrwSy90h5Mc4WzqvpSkp1J9ibZ0907FlEUAADARrOIM2fP6u7bFzAfAACADct3zgAAAEZg3nDWST5eVVdV1bmLKAgAAGAjmveyxmd2981V9agkl1XVtd39yf0nmIa2c5Nk5ZFHzbk4AACAh6a5zpx1983Tn7cmuSTJ0w4wzYXdvaO7d6xs3zbP4gAAAB6yDjmcVdW2qtp+3+Mkz0vy2UUVBgAAsJHMc1nj8UkuqUkjoc1J3tvdH11IVQAAABvMIYez7r4hyekLrIWHqGU381sr62Y9F9FAer2sKwDAWpvj7yS30gcAABgB4QwAAGAEhDMAAIAREM4AAABGQDgDAAAYAeEMAABgBIQzAACAEZinCTWwHulRBgAwSs6cAQAAjIBwBgAAMALCGQAAwAgIZwAAACMgnAEAAIyAcAYAADACwhkAAMAICGcAAAAjIJwBAACMgHAGAAAwAsIZAADACAhnAAAAIyCcAQAAjIBwBgAAMALCGQAAwAgIZwAAACMgnAEAAIyAcAYAADACwhkAAMAICGcAAAAjIJwBAACMgHAGAAAwAsIZAADACAhnAAAAIyCcAQAAjMDmeV5cVWcmeUuSlSRv7+43LaQqgHWie9kVALBeVS27AsbmkM+cVdVKkrclOSvJaUleUlWnLaowAACAjWSeyxqfluT67r6hu+9NcnGSsxdTFgAAwMYyTzg7KcmN+w3fNH0OAACAB+hBvyFIVZ1bVVdW1ZV7d979YC8OAABgXZonnN2c5OT9hh89fe57dPeF3b2ju3esbN82x+IAAAAeuuYJZ59OckpVPb6qtiZ5cZJLF1MWAADAxnLIt9Lv7j1VdV6Sj2VyK/13dvfVC6sMAABgA6lewyY9VXVbki/v99Qjk9y+ZgXA7OybjJn9k7GybzJm9k/G4rHdfdyBRqxpOPu+hVdd2d07llYAHIR9kzGzfzJW9k3GzP7JevCg360RAACAYcIZAADACCw7nF245OXDwdg3GTP7J2Nl32TM7J+M3lK/cwYAAMDEss+cAQAAkCWGs6o6s6o+X1XXV9X5y6oDqurkqvpEVX2uqq6uqldPnz+2qi6rquumP49Zdq1sTFW1UlWfqaqPTIcfX1VXTI+f76+qrcuukY2pqo6uqg9W1bVVdU1V/ZhjJ2NQVb8z/Uz/bFW9r6oOd+xkPVhKOKuqlSRvS3JWktOSvKSqTltGLZBkT5Lf7e7Tkjwjyaum++P5SS7v7lOSXD4dhmV4dZJr9ht+c5K/6O4nJvlGkpcvpSpI3pLko939pCSnZ7KfOnayVFV1UpLfTrKju5+SZCXJi+PYyTqwrDNnT0tyfXff0N33Jrk4ydlLqoUNrru/1t3/NX28M5M/Lk7KZJ+8aDrZRUlesJwK2ciq6tFJfj7J26fDleTZST44ncS+yVJU1VFJfirJO5Kku+/t7jvj2Mk4bE5yRFVtTnJkkq/FsZN1YFnh7KQkN+43fNP0OViqqnpckqcmuSLJ8d39temoW5Icv6Sy2Nj+MskfJNk3HX5Ekju7e8902PGTZXl8ktuS/N30stu3V9W2OHayZN19c5I/TfKVTELZN5NcFcdO1gE3BIGpqnpYkn9M8pruvmv/cT25ralbm7Kmqur5SW7t7quWXQscwOYkP5Lkr7v7qUnuzv0uYXTsZBmm33M8O5P/QDgxybYkZy61KJjRssLZzUlO3m/40dPnYCmqaksmwew93f2h6dNfr6oTpuNPSHLrsupjw/qJJL9YVV/K5PLvZ2fyHZ+jp5fqJI6fLM9NSW7q7iumwx/MJKw5drJsz0nyxe6+rbt3J/lQJsdTx05Gb1nh7NNJTpneNWdrJl/SvHRJtbDBTb/D844k13T3n+836tIk50wfn5Pkw2tdGxtbd7+2ux/d3Y/L5Dj5z9390iSfSPKi6WT2TZaiu29JcmNVnTp96owkn4tjJ8v3lSTPqKojp5/x9+2bjp2M3tKaUFfVz2XyXYqVJO/s7guWUggbXlU9M8mnkvxvvvu9ntdl8r2zDyR5TJIvJ/mV7r5jKUWy4VXVzyT5ve5+flU9IZMzaccm+UySX+/ue5ZZHxtTVf1wJjer2ZrkhiS/mcl//Dp2slRV9cYkv5rJHZk/k+QVmXzHzLGTUVtaOAMAAOC73BAEAABgBIQzAACAERDOAAAARkA4AwAAGAHhDAAAYASEMwAAgBEQzgAAAEZAOAMAABiB/wc9Ui7EMkzFaQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x1080 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15, 15))\n",
    "plt.imshow(image[20:35,0:100], norm=norm_data, origin=\"lower\")\n",
    "plt.title(\"slit[0] slice\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "ext_center = 27\n",
    "ext_width = 8\n",
    "bkg_offset = 4\n",
    "bkg_width = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot cross-disperion cut showing the extraction parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Cross-dispersion Cut at Pixel=50')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAF1CAYAAAAHsfZRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXxU5dn4/8+VPWEn7JkEUBZFTAjiQrTVKlqwWqxVwZ9ScUHrUvWx1ro9Wrfq09LH5evyVK0rKljqRl0qWqGtoIIKiEQUURICTCBAIBIgy/3745yJk0km2zkz50xyvV+v88rMmTPnXLkCc819n/vcR4wxKKWUUgBJXgeglFLKP7QoKKWUaqBFQSmlVAMtCkoppRpoUVBKKdVAi4JSSqkGWhSUZ0RkmIgYEUmxn78pIud5HVeIiPyfiPy313H4jYh8LiLHOdzHTBH5j0shKRdpUejEROT/E5HlIlIlIpvtD91jvI4rGmPMFGPM017HEWKM+aUx5o5Y7FtERonIX0Vkm4hUisgqEblGRJLb8N7ficgcB8du9QNZRBaJyF773842EXlJRAYDGGMOMcYs6ujxnRCRb0Wk2o6rSkTejnj9v0Rki4jsEpEnRCTdizgTmRaFTkpErgHuA34PDATygIeBqVG2T4lfdP7g1e8sIgcCHwKlwKHGmF7AmcAEoIcXMUVxhTGmOzAK6A3c63E8IacaY7rby0mhlSLyY+B64ARgKHAAcJtHMSYuY4wunWwBegFVwJktbPM7YD4wB9gFXAQMAV4DtgPrgFlh2x8BLLe3DQL/a6/PsPdRAewElgEDoxwzGZgNbAPWA5cDBkixX18EXGQ/HgEsBirt7eeF7ccAV9r72Ab8EUgKe/0CoBjYAfwDGBrx3suBr4BvAMH6sCu3f7fPgLH2tk8Bd4a9d5adl+12noZE7PeX9n53Ag8BEiUPc4DXW/jbHAdsjFj3LTAJmAzsB2rsv/HKKPu4Hvga2A2sAX5mrz8Y2AvU2e/fGeX9DX8L+/nlwOrwWOzHbwB/CttuLvBE2L/DvwCbgTLgTiDZfm0m8J8O/NtuOHYzrz0P/D7s+QnAFq//Pyba4nkAusTgj2p9cNSGPmyjbPM7+4PlNKwWYybwL6zWRAYwDtgKHG9vvxSYYT/uDhxlP74EWABkYX3oHwb0jHLMXwJfALlAX+A9oheFF4Cb7NgygGPC9mPs9/bFagF9Gfa+qVgf3AcDKcDNwJKI9y6035sJ/Bj4GOubsNjvG2xv+xR2UQCOxypA44F04P8B/4rY79/t/eTZuZscJQ9bgPNb+NscR5SiEPa3m9PKv4EzsYp8EjAN+C7s92r1Aznib9EP+CfwbDOxDMIqqMcD52AV6h72ay8Dfwa6AQOAj4BLmosBWIVVTJtbHo7IQ9DO79tAQdhrK4FpYc/72X+XbK//TybSot1HnVM2sM0YU9vKdkuNMa8YY+qx/gMdDfzWGLPXGLMCeBz4hb1tDTBCRPoZY6qMMR+Erc8GRhhj6owxHxtjdkU53lnAfcaYUmPMduDuFmKrweoCGGLHE9kH/j/GmO3GmBKsbrKz7fW/BO42xhTbv//vgXEiMjTsvXfb7622j9MDOAjrm32xMWZzM/Gcg/UN+BNjzD7gBmCiiAwL2+YeY8xOO6b3sAprc7Kxvj3HjDHmr8aYTcaYemPMPKwWzBHt3M0DIrIT68N2M3BNM8fZAlwKPA3cD/zCGLNbRAYCJwNXG2O+M8aUY7XIpkeJN98Y0zvKclnYpucAw7D+bbwH/ENEetuvdcdqWYaEHvupS873tCh0ThVAvzb0mZeGPR4CbDfG7A5btwHIsR9fiNW3/IWILBORU+z1z2J10cwVkU0i8gcRSRWRH4SdDPw87BilEfuP5jqsb+4f2aNdLmgh9g32vsH6sLhfRHbaH2jb7f3kNPdeY8w/gQexunvKReRREenZTDxDwuM1xlRh5Tl8v1vCHu/B+pBqTgUwOMprrhCRX4jIirA8jMUq/O1xpf2hnGOMOccYszXKdguwWolrw4r3UCAV2BwWw5+xWgwdZox53xhTbYzZY4y5G6sl8QP75Sog/G8Xehz+b1q1QotC57QU2IfVNdSS8ClyNwF9RST8W1UeVl8wxpivjDFnY/2n/h9gvoh0M8bUGGNuM8aMAYqAU7C+Lf7bfH8y8BB7f5uxuo7C9998YMZsMcbMMsYMweqielhERoRtErmfTfbjUqwuivBvmpnGmCVRfm+MMQ8YYw4DxmAVvt80E9ImrA86AESkG9Y3/rJov0ML3gF+3sLr32F1x4WOlQz0Dw+5pZ3braLHgCuwuk56A6uximOr7++Au7DO4QwWkVCLrRTr32C/sL9Dz7B/C5Exfx72JSJy+b8Wjm34/vf6HCgIe60ACBpjKhz9dl2MFoVOyBhTCdwCPCQip4lIlv3tfYqI/CHKe0qBJcDdIpIhIvlYrYM5ACJyroj0t7uadtpvqxeRH4nIofYH1y6s7pj6KKG9CFwpIgER6YN1MrRZInKmiATspzuw/vOH7/c3ItJHRHKBq4B59vr/A24QkUPs/fQSkTNbOM7hInKkiKRifRjvjRL/C8D5IjLOHub4e+BDY8y30fbdgluBIhH5o4gMsuMYISJz7K6QL4EMEfmJHdfNWOcxQoLAMBGJ9v+3G1a+ttr7Ph+rpRD+/oCIpHUg9kZE5IfA+VjdjOcB/09EcuwuuLeBP4lITxFJEpEDReTY5vZjrGGu3aMsv7SPlSciR4tImv1v9DdYrZ/37d08A1woImPsPN6MdV5ItYMWhU7KGPMnrD7gm7E+HEqxvjm+0sLbzsbqr92EdZLwVmPMO/Zrk4HPRaQKq+94ut0nPwhrFNMurG+Li7G6lJrzGFZX00rgE+ClFmI5HPjQPt5rwFXGmPVhr7+KdYJ4BfA61igXjDEvY7Vk5orILqxvyFNaOE5PO64dWN1DFVijmRqx8/DfwN+wWjwHEqV/vDXGmK+BiVi5/lxEKu39Lgd220X9MqxzOmVYxWpj2C7+av+sEJFPmtn/GuBPWC3GIHAo339wgnXS+HNgi4hs68jvAGB3sz2DNXS1zBjzb6y/w5MiIliFIg1r9NMOrH8nTrrNegCP2Psqw/o3OSXUEjDGvAX8AetcQwnW3/NWB8frksQYvcmOSiwiYoCRxph1XseiVGejLQWllFINtCgopZRqoN1HSimlGmhLQSmlVAMtCkoppRok/MyY/fr1M8OGDfM6DKWUSigff/zxNmNM/8j1CV8Uhg0bxvLly70OQymlEoqINDvNjHYfKaWUaqBFQSmlVAMtCkoppRpoUVBKKdVAi4JSSqkGWhSUUko10KKglFKqgRYFpZRSDbQoKKWUaqBFQSmlVANXioKI/Jd94+3VIvKCff/U4SLyoYisE5F5ofvBiki6/Xyd/fqwsP3cYK9fKyI/diM2pZRSbee4KIhIDnAlMMEYMxZIxrp37f8A9xpjRmDdU/VC+y0XAjvs9ffa2yEiY+z3HYJ179WH7ZvBK6WUihO3uo9SgEwRSQGysG5sfjzWjboBngZOsx9PtZ9jv36CfZPvqcBcY8w+Y8w3wDrgCJfiU0p1MccdZy2qfRwXBWNMGTAbKMEqBpXAx8BOY0ytvdlGIMd+nAOU2u+ttbfPDl/fzHsaEZGLRWS5iCzfunWr019BKeXQHYvv4I7Fd3gdhu8lQp4cT50tIn2wvuUPB3YCf8Xq/okZY8yjwKMAEyZM0PuJKuWxtRVrvQ6hiUDA6wia8mOeIrlxP4VJwDfGmK0AIvIScDTQW0RS7NZAACizty8DcoGNdndTL6AibH1I+HuUUj425/Q5XofQxBz/heTLPEVy45xCCXCUiGTZ5wZOANYA7wFn2NucB7xqP37Nfo79+j+NMcZeP90enTQcGAl85EJ8Siml2shxS8EY86GIzAc+AWqBT7G6dl4H5orInfa6v9hv+QvwrIisA7ZjjTjCGPO5iLyIVVBqgcuNMXVO41NKxd4t790CwO0/ut3jSL539dXWz/vu8zaOcH7MUyRXbsdpjLkVuDVi9XqaGT1kjNkLnBllP3cBd7kRk1Iqfkp3lba+UZytWOF1BE35MU+REv4ezUop7z059UmvQ0gIiZAnneZCKaVUAy0KSinHbnjnBm545wavw/C9RMiTdh8ppRyrqK7wOoQmRo3yOoKm/JinSGKNBk1cEyZMMMuXL/c6DKWUSigi8rExZkLkeu0+Ukop1UCLglLKsWvfvpZr377W6zAaOfLIFRx7rL+mlfBjniLpOQWllGPVNdVeh9DEqlV7SUvzV/e4H/MUSc8pKKU6peTkf5OamsLevRO9DsWX9JyCUqrL+O6776ivr2Pfvv3U1ta2/gbVQIuCUsqxq9+6mqvfutrrMBoEg0H7kWHTpk2exhLOb3lqjhYFpVSnYxWFFcAKSkpKvA4noeiJZqWUY/dN9tFUpISKwn8BUFLynLfBhPFbnpqjLQWlVKfzffcR2lJoJ20pKKUcu/z1ywF46CcPeRyJxSoKz5KamsqGDYu8DqeB3/LUHG0pKKUcy0zNJDM10+swGpSXl5OSMoz09AN91VLwW56aoy0FpZRjs0+a7XUIjQSDQVJT00hPz/BVUfBbnpqjLQWlVKcTDAZJS0slIyPdV0UhEWhRUEo5dvGCi7l4wcVeh9HAKgpWS2HXrl1UVlZ6HRLgvzw1R4uCUsqx7MxssjOzvQ6jQTAYJDe3jPz8KgA2bNjgcUQWv+WpOXpOQSnl2N2T7vY6hAZ79+6lsrKSn/3sI048sRevvGINS83Pz/c6NF/lKRptKSilOpXy8nIABg4cyNChQwG9VqE9tKWglHLs/FfPB+DJqU96HMn3F649+eRPePPNgaSlpfmmKPgpT9FoUVBKOZbbM9frEBqEisLevd3Zvl3Izc31TVHwU56i0aKglHLs9h/d7nUIDUJFIS0tDYC8vDzfnGj2U56i0XMKSqlOJXROIS0tFbCKgl9aColAi4JSyrFzXzqXc1861+swAKul0KNHD5KSkgEYOnQomzZtoqamxuPI/JWnaLT7SCnl2Ojs0V6H0CAYDDJgwABOOMF6PmRIHvX19WzatKlhNJJX/JSnaLQoKKUc++9j/9vrEBoEg0EGDhzIf9shLVyYB1jDUr0uCn7KUzTafaSU6lRCRSEkL88qCn452ex3WhSUUo5Nnz+d6fOnex0G8H1RmDIFpkyB3FxrGKgfTjb7KU/RaPeRUsqxcYPGeR0CADU1NVRUVDBw4ECKi611WVlZ9O/f3xdFwS95aokWBaWUY9cfc73XIQCwdetWgEbdR+CfYal+yVNLtPtIKdVphC5c82tRSARaFJRSjv38xZ/z8xd/7nUYLRaFDRs2YIzxIqwGfslTS7T7SCnl2MTARK9DABrPkHrKKd+vz8vLo6qqip07d9KnTx+PovNPnlriSlEQkd7A48BYwAAXAGuBecAw4FvgLGPMDhER4H7gZGAPMNMY84m9n/OAm+3d3mmMedqN+JRSsXVt0bVehwA0bilcGxZS+BTaXhYFv+SpJW51H90PvGWMOQgoAIqB64F3jTEjgXft5wBTgJH2cjHwCICI9AVuBY4EjgBuFRHv/npKqYQTDAbJyMige/fujdaHrlXQ8wqtc1wURKQX8EPgLwDGmP3GmJ3AVCD0Tf9p4DT78VTgGWP5AOgtIoOBHwMLjTHbjTE7gIXAZKfxKaVi76cv/JSfvvBTr8NouEZBRDjuODjuOGu9X4qCX/LUEje6j4YDW4EnRaQA+Bi4ChhojNlsb7MFCJ35yQFKw96/0V4Xbb1SyudOGH6C1yEATa9mDunfvz/p6emeX9Xslzy1xI2ikAKMB35ljPlQRO7n+64iAIwxRkRcO+0vIhdjdT01fANQSnnnqqOu8joEwCoKzc1vlJSU5Iub7fglTy1x45zCRmCjMeZD+/l8rCIRtLuFsH+W26+XAeG3HwrY66Ktb8IY86gxZoIxZkL//v1d+BWUUp1BtJYCWCebvS4KicBxUTDGbAFKRSQ0J+wJwBrgNeA8e915wKv249eAX4jlKKDS7mb6B3CSiPSxTzCfZK9TSvnclOemMOW5KZ7GUFdXx9atW6MWBT9cwOaHPLXGresUfgU8JyJpwHrgfKyC86KIXAhsAM6yt30DazjqOqwhqecDGGO2i8gdwDJ7u9uNMdtdik8pFUOnjjrV6xCoqKigvr6+oSicdVbj1/Py8hputpOamupBhP7IU2tcKQrGmBXAhGZeanJWxViXFF4eZT9PAE+4EZNSKn4uO/wyr0NodOEawGURIeXl5WGMYePGjQwfPjze4Vkx+SBPrdFpLpRSnULkFBd79lhLiF+GpfqdTnOhlHJs0jOTAHjnF+94FkNkUTj5ZGv9okXWz/Crmr3ihzy1RouCUsqxaYdM8zqEqJPhhQQCAcDbouCHPLVGi4JSyrFZh83yOgSCwSCpqan07t272dczMzMZMGCApxew+SFPrdFzCkqpTiEYDDJgwACsOTeb54dhqX6nRUEp5dhxTx3HcU8d52kMLV24FuL1BWx+yFNrtPtIKeXYzHEzvQ6hSVGYObPpNnl5ebz11lsYY1psUcSKH/LUGi0KSinH/PBhFwwGyc/Pb3gerSh899137Nixg759+8YvuFBMPshTa7T7SCnlWE1dDTV1NZ4d3xhDeXl5o5bCtm3WEi50rYJXJ5u9zlNbaFFQSjl24rMncuKzJ3p2/J07d1JTU9OoKJxxhrWE8/oCNq/z1BbafaSUcuyi8Rd5evzWrlEI8foCNq/z1BZaFJRSjp2bf66nx29rUejXrx8ZGRmeFQWv89QW2n2klHJsT80e9tTsaX3DGGlrURART69V8DpPbaEtBaWUYyc/Z000tGjmIk+OHyoKAwYMaHXbvLw8z040e52nttCioJRy7NIJl3p6/GAwSFJSEtnZ2Q3rLo0SUl5eHm+++WacImvM6zy1hRYFpZRj08Z6O9FbMBikf//+JCcnN6ybFiWkoUOHsnnzZvbt20d6enqcIrRj8jhPbaHnFJRSjlXuraRyb6Vnx29uiovSUmuJFBqWWlbW7C3gY8rrPLWFFgWllGNT505l6typnh2/uaIwY4a1RPLyWgWv89QW2n2klHLsyiOv9PT4wWCQUaNGtWlbL69q9jpPbaFFQSnl2OkHn+7ZsZub4qIlXt5sx8s8tZV2HymlHNu2Zxvb9mxrfcMYqKqqorq6us1FISMjg0GDBnlSFLzMU1tpS0Ep5dgZL1qTDHkx/r6tF66F8+oCNi/z1FZaFJRSjv164q89O3a0ovDrFkLKy8tj9erVsQyrWV7mqa20KCilHDt19KmeHTtaUTi1hZDy8vJ4/fXX436zHS/z1FZ6TkEp5diWqi1sqdriybGjTXGxdq21NCcvL4/q6moqKipiHV4jXuaprbSloJRybPr86YC35xT69+/faP0ll1g/Fy1q+p7wKbT79esXy/Aa8TJPbaVFQSnl2PXHXO/ZsYPBINnZ2aSmprb5PeEXsI0fPz5WoTXhZZ7aSouCUsqxySMme3bs5q5mbo1XVzV7mae20nMKSinHSitLKa1sZqKhOOhIUcjOziYzMzPuVzV7mae20qKglHJsxsszmPFyMxMNxUF7rmYO8epmO17mqa20+0gp5djNP7zZs2NHaync3EpIQ4cOjXtR8DJPbaVFQSnl2KQDJnly3Orqanbv3t1sUZjUSkh5eXmsWrUqRpE1z6s8tYd2HymlHFu/Yz3rd6yP+3FbmuJixQpriSYvL48tW7awb9++WIXXhFd5ag9tKSilHLvg1QuA+I+/b6koXH219bO56xTg+xFIpaWljBgxIhbhNeFVntpDi4JSyrHbjrvNk+NGu5q5LcKHpcarKHiVp/bQoqCUcuzYYcd6ctyOzJAaEn5Vc7x4laf20HMKSinH1m5by9ptUSYaiiEnLYWcnBxEJK5Fwas8tYdrRUFEkkXkUxH5u/18uIh8KCLrRGSeiKTZ69Pt5+vs14eF7eMGe/1aEfmxW7EppWLrkr9fwiV/vyTuxw0Gg/Tq1YuMjIx2vzc9PT3uN9vxKk/t4Wb30VVAMdDTfv4/wL3GmLki8n/AhcAj9s8dxpgRIjLd3m6aiIwBpgOHAEOAd0RklDGmzsUYlVIx8PsTfu/JcVu6cO33bQgpLy8vrlc1e5Wn9nClpSAiAeAnwOP2cwGOB+bbmzwNnGY/nmo/x379BHv7qcBcY8w+Y8w3wDrgCDfiU0rFVlFuEUW5RXE/bktTXBQVWUtL4n0Bm1d5ag+3uo/uA64D6u3n2cBOY0yt/XwjkGM/zgFKAezXK+3tG9Y3855GRORiEVkuIsu3bt3q0q+glOqo1eWrWV0e/zuZtVQUliyxlpbk5uZSWlqKMSYG0TXlVZ7aw3FREJFTgHJjzMcuxNMmxphHjTETjDETIudQV0rF3xVvXMEVb1wR9+O2VBRuvNFaWhIIBKiurmbHjh0xiK4pr/LUHm6cUzga+KmInAxkYJ1TuB/oLSIpdmsgAJTZ25cBucBGEUkBegEVYetDwt+jlPKxP574x7gfc//+/ezYsaNDw1FDAoEAABs3bqRv375uhRaVF3lqL8ctBWPMDcaYgDFmGNaJ4n8aY84B3gPOsDc7D3jVfvya/Rz79X8aq+32GjDdHp00HBgJfOQ0PqVU7B2ecziH5xwe12OWl5cDHbtGISS8KMSDF3lqr1hevPZbYK6I3Al8CvzFXv8X4FkRWQdsxyokGGM+F5EXgTVALXC5jjxSKjGs2GJNMjRu0Li4HdPJhWsh8S4KXuSpvVwtCsaYRcAi+/F6mhk9ZIzZC5wZ5f13AXe5GZNSKvaufsuaaCiec/o4uXAtZNCgQSQlJcWtKHiRp/bSaS6UUo7dN/m+uB+ztZbCfW0IKSUlhcGDB1NWFp/Tl17kqb20KCilHPOiO6S1ojCujSHl5OTEraXg526jEJ37SCnl2LKyZSwrWxbXY5aXl9OtWze6devW7OvvvGMtrQkEAnErCl7kqb20paCUcuw3C38DxP+cQksnme+80/rZ2h3YAoEA77SlerjAizy1lxYFpZRjD578YNyP2VpRaKtAIMCuXbvYtWsXPXv2bP0NDniRp/bSoqCUcmzsgLFxP2YwGHTl5jihYallZWUxLwpe5Km99JyCUsqxJaVLWFLaykRDLnOzpQDxuVbBizy1l7YUlFKO3fiuNclQvPrKa2tr2bZtW8IVhXjnqSO0KCilHPvzKX+O6/G2bduGMabFovDnNoY0ZMgQID5FId556ggtCkopx0b3Gx3X47XlaubRbQwpPT2dAQMGxKUoxDtPHaHnFJRSji3+djGLv10ct+O1Zd6jBQuspS3ida1CvPPUEdpSUEo5duuiW4H49ZW3pSj86U/Wz1NPbX1/gUCAb7/91oXIWhbvPHWEFgWllGNPTH0irsdzY9rscDk5Obz//vuu7Ksl8c5TR2hRUEo5dkCfA+J6vGAwSHp6umvXFQQCASoqKqiuriYzM9OVfTYn3nnqCD2noJRy7J317/DO+vhMFQHfX6MgIq7sL/wCtliKd546QlsKSinH7vyXNdHQpANamWjIJW5duBYSfq2CG1dJRxPvPHWEFgWllGPP/uzZuB4vGAw2fJBH82w7QorXBWzxzlNHaFFQSjmW2ys3rscLBoMcdthhLW6T246QcnJygNgXhXjnqSP0nIJSyrG31r3FW+veisux6uvr2bp1a6vdR/PmWUtbdOvWjT59+sS8KMQzTx2lLQWllGP3/OceACaPmBzzY23cuJHa2lpyW2kKPPKI9XPatLbtNx4XsMUzTx2lRUEp5djcM+bG7VjFxcUAjBkzxtX9xqMoxDNPHaVFQSnl2KDug+J2rDVr1gBw8MEHu7rfQCDAJ5984uo+I8UzTx2l5xSUUo4tWLuABWvbONGQQ8XFxfTr14/+/fu7ut9AIEAwGGT//v2u7jdcPPPUUdpSUEo59qel1kRDp45uw0RDDhUXF7veSoDvRyBt2rSJYcOGub5/iG+eOkqLglLKsflnzY/LcYwxrFmzhjPOOKPVbee3M6Twq5pjVRTilScntCgopRzrl9UvLsfZunUr27dvb9NJ5n7tDCkeF7DFK09O6DkFpZRjLxW/xEvFL8X8OKGRR23pPnrqKWtpq3gUhXjlyQltKSilHHvgwwcAOP3g02N6nPaMPAoVhJkz27bvnj170r1795gWhXjlyQktCkopx16d/mpcjlNcXEz37t1bnfeoI0Qk5tcqxCtPTmhRUEo51iujV1yOs2bNGg4++GDXpsyOFOuiEK88OaHnFJRSjs1bPY95q9s40ZADsRqOGhLrohCvPDmhLQWllGOPLLcmGpo2to0TDXVAZWUlmzZtcn16i3CBQIDNmzdTW1tLSor7H4/xyJNTWhSUUo69cc4bMT9Ge0YeAbzRgZACgQB1dXUEg8GGi9ncFI88OaXdR0opx7JSs8hKzYrpMdpbFLKyrKU9Yj0sNR55ckqLglLKsTmr5jBn1ZyYHqO4uJj09HSGDx/epu0fftha2iPWN9uJR56c0qKglHLs8U8e5/FPHo/pMdasWcOoUaPa3Nf/4ovW0h6xbinEI09O6TkFpZRjC2csjPkxiouLmTBhQkyPkZ2dTXp6OmVlZTHZfzzy5JTjloKI5IrIeyKyRkQ+F5Gr7PV9RWShiHxl/+xjrxcReUBE1onIKhEZH7av8+ztvxKR85zGppSKj9TkVFKTU2O2/+rqar755puYjjyC2F/AFus8ucGN7qNa4NfGmDHAUcDlIjIGuB541xgzEnjXfg4wBRhpLxcDj4BVRIBbgSOBI4BbQ4VEKeVvT614iqdWPBWz/a9duxZjTEyvUQiJZVGIdZ7c4LgoGGM2G2M+sR/vBoqBHGAq8LS92dPAafbjqcAzxvIB0FtEBgM/BhYaY7YbY3YACwH/3shUKdUg1h927R155ERXLwqunlMQkWFAIfAhMNAYs9l+aQsw0H6cA5SGvW2jvS7a+uaOczFWK4O8vDx3gldKddiimYtiuv/i4mKSkpIYNWpUm9+zaFHHjhUIBCgrK6O+vp6kJHfH4sQ6T25w7TcWke7A34CrjTG7wtrAscYAACAASURBVF8zxhjAuHUsY8yjxpgJxpgJbt+STynlP2vWrOHAAw8kPT095scKBALs37+fbdu2xfxYfuRKURCRVKyC8JwxJjRZeNDuFsL+WW6vLwNyw94esNdFW6+U8rnHPn6Mxz5+LGb778icR7NnW0t7xXJYaqzz5AY3Rh8J8Beg2Bjzv2EvvQaERhCdB7watv4X9iiko4BKu5vpH8BJItLHPsF8kr1OKeVz8z6fx7zPYzPRW01NDV9++WW7Rx79/e/W0l6xLAqxzJNb3DincDQwA/hMRFbY624E7gFeFJELgQ3AWfZrbwAnA+uAPcD5AMaY7SJyB7DM3u52Y8x2F+JTSsXYO794J2b7/vrrr6mtrY3LSWaIbVGIZZ7c4rgoGGP+A0Sb3PyEZrY3wOVR9vUE8ITTmJRSnUc8Rx4BDBgwgJSUlJhOoe1nOs2FUsqxh5c9zMPL2jnRUBuFbsF50EEHxWT/kZKSkhgyZEhMikIs8+QWLQpKKccWfLmABV8uiMm+i4uLyc3NpUePHu16X2amtXRErK5ViGWe3KJzHymlHHvznDdjtu+O3m3tTQchBQIBVqxY0fqG7RTLPLlFWwpKKd+qr6+nuLg45nMeRQq1FKxToF2LFgWllGP3f3A/939wv+v7LSkpobq6ukMthTvusJaOCAQC7Nmzh507d3ZsB1HEKk9u0qKglHLs3W/e5d1v3nV9v05GHr37rrV0RKyGpcYqT27ScwpKKcdeO/u1mOw3NPLIi+4jsIrCoYce6tp+Y5UnN2lLQSnlW8XFxfTv35/s7Oy4HjfWd2DzMy0KSinHZi+ZzewlHZhoqBUdHXnk1KBBg0hKSnK9KMQqT27S7iOllGNLNy51fZ/GGNasWcP06dM79H4njYvU1FQGDRrkelGIRZ7cpkVBKeXY3876m+v7DAaD7Ny5s8Mthb85DCknJ8f1ohCLPLlNu4+UUr7k1UnmkFjegc3PtCgopRy75z/3cM9/7nF1n04nwrvhBmvpqFgUhVjkyW3afaSUcmzFFvenhCguLqZHjx4MGTKkQ+9f6rD7PhAIsGvXLnbv3t3ueZeiiUWe3KZFQSnl2Nwz5rq+zzVr1jBmzBis+3jFX2hYallZmWsztMYiT27T7iOllC95NRw1pKteq6BFQSnl2B2L7+COxR2caKgZO3bsYMuWLZ2uKLidp1jQ7iOllGNrK9a6ur/QSWYnI4/sz/QOC53LcLMouJ2nWNCioJRybM7pc1zdnxu34JzjMKSMjAz69+/valFwO0+xoN1HSinfKS4uJj09nWHDhnkaR1e8VkGLglLKsVveu4Vb3rvFtf2tWbOGgw46iOTk5A7v4+qrrcUJt4uC23mKBe0+Uko5Vrqr1NX9FRcXc9RRRznahxt308zJyWHJkiXOd2RzO0+xoEVBKeXYk1OfdG1f3333HRs2bOD88893bZ8dFQgEqKiooLq6mszMTMf7czNPsaLdR0opX1m7di3GGM/mPAoXfgFbV6FFQSnl2A3v3MAN7ziYaCiMGyOP3OL2tQpu5ilWtPtIKeVYRXWFa/tas2YNycnJjBw50tF+Ro1yHovbRcHNPMWKFgWllGOPnvqoa/sqLi5mxIgRpKWlOdrPoy6ElJOTA7jXfeRmnmJFu4+UUr7i9ZxH4bp3707v3r271LUKWhSUUo5d+/a1XPv2tY73s3//fr766itXisLFF1uLU25eq+BWnmJJu4+UUo5V11S7sp9Vq1ZRV1fnysijL790ISDcLQpu5SmWtCgopRx76CcPOd7HokWLOOuss+jRowc//OEPXYjKHYFAgBVuXAmHO3mKNe0+Ukp5yhjDvffey6RJk8jOzuajjz4iLy/P67AaBAIBgsEg+/fv9zqUuNCioJRy7Oq3rubqt9o/0dCePXs499xzueaaazj11FP58MMPXbvLmVuGDh2KMYZZs2Y1XEPRUR3NUzxpUVBKeWL9+vUUFRXxwgsvcOedd/K3v/2Nnj17urb/ceOsxanp06dz+eWX89e//pUxY8YwdepUV+dD8hsxxngdgyMTJkwwy5cv9zoMpVQ7vP3220yfPh1jDM8//zxTpkzxOqRWbd26lQcffJAHH3yQ7du3c/TRR3PddddxyimnkJSUeN+vReRjY8yEyPWJ95sopRKWMYa7776byZMnEwgEWL58eUIUBID+/ftz2223UVJSwgMPPMDGjRuZOnUqY8eO5cknn2Tfvn1eh+gK3xUFEZksImtFZJ2IXO91PEqp1l3++uVc/vrlUV+vq6vjiy++4Mwzz+TGG29k2rRpLF26lAMPPDBmMZ17rrW4rVu3bvzqV79i3bp1PP/886SlpXHBBRcwfPhwLrnkEp5//vmoV0C3lic/8NWQVBFJBh4CTgQ2AstE5DVjzBpvI1NKtSQz9ftppauqqli1ahUrV65kxYoVrFixgs8++4zq6mqSkpKYPXs211xzDSIS05hifRFySkoKZ599NtOnT2fhwoU89NBDzJ07l0ft+TUOPPBAjj322IZl6NChjfLkV746pyAiE4HfGWN+bD+/AcAYc3e09+g5BaViyxhDVVUVlZWV7Ny5k507dzY8Dl/3zTffsHLlStatW0foc6VPnz6MGzeOgoICxo0bR1FRkeOJ7trquOOsn4sWxeVwgNUiWrlyJYsWLWLx4sX8+9//ZseOHYA1iunYY49l+PDhZGZmkpWVRWZmZqMltC4tLY3k5GSSk5NJSkpqeBz5fMCAAR0urtHOKfiqpQDkAOG3JtoIHBmLAxUUFLBmzfcNkPDERiZZRFpckpKSSEpKavQ4fAmtT0lJITk5mZSUlKiPs7Ky6N69Oz169Ij6Mzc3l4KCAlJS/Pbnix1jDBUVFXz11Vd888031NTUAN//bZp73LdvXwKBAIFAgF69esX8m6mf1NXVsWHDBrZt20ZFRUXDEvl8x44d7N27l5qaGvbv39/oZ+hxbW1tq8dLS0sjEAgwbtw4ZsyY0VAEcnNzu1Tek5OTGT9+POPHj+eaa66hvr6ezz77jMWLF7N48WLefPNNtm7d6trxqqurycjIcG1/4L+i0CYicjFwMdDhi1wuuOACysvLAQhvLUW2nIwxrS719fVNfja31NXVUVtb2/AztISe79+/n8rKSqqqqti9ezdVVVVUVVU1iQkgKyuLI488kqOPPpqioiImTpxI7969O5QLP9m5cydffvklX331VZNl586dHd5vVlZWQ4EIBALk5OQQCAT4wQ9+wKGHHuribxB/lZWVrFy5stGyevVq9u7d22RbEaFPnz5kZ2eTnZ3NoEGDyMzMJDU1lbS0tEY/wx+HJobr1asXvXv3bvS4V69eXLnwSiAxZgGNp6SkJAoKCigoKODKK6/k4gUXY4zh/kn3s2fPHqqrqxstoXX79++nrq6u4XMjtIQ/r6+vJzU11fWY/VYUyoDcsOcBe10jxphHgUfB6j7qyIGuuuqqjrwt7urr66murm4oFLt37+bLL7/k/fff5/333+fuu++mrq4OgEMOOaShSEyaNKlh2t9EsHnzZn7zm9/w3HPPNawTEfLy8hg5ciRnn302I0eOZOTIkRx44IFkZGQ0FGWgyeP6+noqKiooKytj48aNbNy4seHxokWL2LRpE7W1tSQnJ/O73/2OG264wdFN4uOprKyMRx99tKHPfsOGDQ2v9evXj4KCAi677DIOOeQQBgwY0FAAsrOz6d27d0x+z+zMbNf36dTEiV5H0FQoT1lZWWRlZXkcTfP8dk4hBfgSOAGrGCwD/j9jzOfR3tPVzylUVVXx0UcfsWTJEt5//32WLl1KZWUl3bp1Y+XKlTEd3eGG2tpaHnzwQW655Rb27dvHlVdeyTHHHMPIkSM54IADXG8ah9TX11NWVsZvf/tbXnjhBY455hjmzJnD0KFDY3I8N/3gBz9gyZIljBo1qqGbJvRtdPDgwV2qu0Z1XLRzCm3qHonnApyMVRi+Bm5qbfvDDjvMqO/V1dWZZcuWmR49ephJkyaZ+vp6r0OK6l//+pc59NBDDWAmT55svvzyy7jHUF9fb5599lnTo0cP07NnT/P888/HPYb2WLx4sQHMAw884HUoKsEBy01zn8HNrUykRYtC8x566CEDmGeeecbrUJrYvHmzmTFjhgFMXl6eefnllz0vXuvXrzdFRUUGMOeee67ZuXOnp/FEc9JJJ5kBAwaYPXv2eB1KIzNfmWlmvjLT6zAaOf10a/ETP+UpWlHw3cVryh2//OUvmThxIv/1X//Ftm3bvA4HsLqKHnjgAUaPHs28efO46aabKC4u5rTTTvO8y2P48OEsXryY2267jRdeeIFx48bx/vvvexpTpGXLlvH2229zzTXXkJnpr/HuuT1zye2Z2/qGcVRRYS1+4sc8NdFcpUikRVsK0a1evdqkpqaaGTNmeB2KKS4uNgUFBQYwJ510klm7dq3XIUW1ZMkSM3z4cJOUlGRuueUWU1NT43VIxhhjTjvtNNO7d29TWVnpdSgJ4dhjrUU1D20pdD2HHHIIv/3tb3n22WdZuHChp7Hcc889fP3118yfP5+33nqLUaNGeRpPSyZOnMiKFSuYMWMGt99+OyeddFLofJdnVq9ezSuvvMKVV17p6kyiSkXy1eijjujqo49as3fvXvLz86mrq+Ozzz7zbBjcj370I2pra/n3v//tyfE76q677uLmm2/miy++YPTo0Z7Fcc455/Dqq6+yYcMGsrP9N/zz3JesSYbmnD7H40i+58UVza3xU550ltQuKiMjg0cffZT169dz++23exZHaWkpubk+70ttxumnnw7g6fz569atY+7cuVx66aW+LAgAo7NHMzrbu6LZnBNOsBY/8WOemmiuTymRFj2n0DYXXnihSU5ONp9++mncj11fX2/S09PNddddF/djO1VXV2f69OljLrroIs9iuOiii0x6errZtGmTZzGozgc9p9C1/eEPfyA7O5tZs2Y1XAEdL1u3bmXfvn0J2VJISkpi4sSJLF261JPjl5aW8vTTT3PhhRcyePBgT2JQXYsWhS6ib9++3H///SxfvpwHH3wwrscuLbXmOEzEogDWiefPP//c0dxLHTV79myMMVx33XVxP3Z7TJ8/nenzp3sdRiNTpliLn/gxT5G0KHQh06ZN4+STT+amm26ipKQkbsdN9KJQVFQEwAcffBDX45aXl/PYY49x7rnn+n76jXGDxjFukAs3RHZRdbW1+Ikf8xTJbxPiqRgSER5++GHGjBnDZZddxoIFC+Jy0ViiF4UjjjiCpKQkli5dyuTJk+N23HvvvZe9e/dy/fX+vwHh9cf4P0Y/SIQ8aUuhixk6dCh33nknr7/+On/961/jcszS0lIyMjLo169fXI7ntu7du5Ofnx/XEUg7duzgoYce4swzz/R0KKzqerQodEG/+tWvOOyww7jyyiuprKyM+fFKS0sJBAKeT2XhRFFRER988EHcTtI/+OCD7N69mxtvvDEux3Pq5y/+nJ+/+HOvw/C9RMiTFoUuKCUlhdmzZxMMBnnvvfdifrySkpKE7ToKKSoqoqqqis8/jzqLu2uqqqq47777OOWUUygoKIj58dwwMTCRiQF/3cDglFOsxU/8mKdIek6hiwr1k3/66aecdtppMT1WaWkpP/rRj2J6jFgLnWxesmQJ+fn5MT3Wn//8Z7Zv385NN90U0+O46dqia70OoYlr/ReSL/MUSVsKXVRWVhYHHXQQn3zySUyPU1dXx6ZNmxK+pTBs2DAGDhwY8/MKe/fuZfbs2Rx//PEcddRRMT2WUs3RlkIXVlhYyKIYTwyzefNm6urqEr4oiAhFRUUxv4jtySefZMuWLY1uS5oIfvrCTwF47ezXPI7ke36c+8iPeYqkLYUurLCwkLKyMrZu3RqzYyT6cNRwRUVFrFu3jvLy8pgd4/HHH+fwww9PuO62E4afwAnDfTbRkA8lQp60KHRhhYWFAHz66acxO0ZnKgoT7TvBx7K18PXXX3PUUUcl3Eitq466iquOusrrMHwvEfKkRaEL06LQPocddhipqakxKwqVlZVUVlaSl5cXk/0r1RZaFLqwPn36MGzYsJgXhe7du9OrV6+YHSNeMjIyOOyww2J2sjlUQBOxKEx5bgpTnvPZREM+lAh50hPNXVxhYWHMi0Jubm7CdYdEM3HiRB555BH2799PWlqaq/vesGEDkJhF4dRRp3odQhNnneV1BE35MU+RtCh0cYWFhbz88svs3r2bHj16uL7/RL25TjRFRUXce++9rFy5ksMPP9zVfYcmKfT75HfNuezwy7wOoYnL/BeSL/MUSbuPurjQeYWVK1fGZP+dsShAbO7EVlJSQmpqKgMHDnR9313Rnj3WotpHi0IXF8uTzfv37ycYDHaqojBkyBDy8vJiUhQ2bNhAbm4uSUmJ999y0jOTmPTMJK/DaOTkk63FT/yYp0jafdTFDRkyhP79+8ekKJSVlWGM6VRFAazWwvvvv+/6fktKShLyfALAtEOmeR1CQkiEPGlR6OJEJGYnmzvTcNRwRUVFzJ071/WusZKSEo4//njX9hdPsw6b5XUICSER8pR47VTluvHjx/P555+zb98+V/fbWYtCLC5iq62tpaysLGFbCqrz0KKgKCwspKamxvVpoTtrUSgoKCAzM9PVolBWVkZ9fX3CFoXjnjqO4546zuswfC8R8qTdR6rRyebx48e7tt/S0lL69OlDt27dXNunH6SmpnLEEUe4erI5kYejAswcN9PrEJqYOdPrCJryY54iaVFQHHjggfTo0cP18wqdbThquIkTJzJ79myqq6vJzMx0vL9QUUjUloIfP+y0KHSMdh8pkpKSKCgo0KLQDkVFRdTW1vLxxx+7sr/Q1cyJmq+auhpq6mq8DqORbdusxU/8mKdIWhQUYHUhrVy50tV7EHfmohA62exWF1JJSQn9+vUjKyvLlf3F24nPnsiJz57odRiNnHGGtfiJH/MUSbuPFGAVhe+++45169YxevRox/vbs2cPFRUVnbYo9OvXj5EjR7paFBL1fALAReMv8jqEhJAIedKioAAaTjB/+umnrhSFjRs3AonbHdIWRUVFvPnmmxhjHE/4V1JSwqhRo1yKLP7OzT/X6xASQiLkSbuPFABjxowhLS3NtXs2d9bhqOGKioooLy9n/fr1jvZjjGHDhg0Je5IZYE/NHvbU6ERDrUmEPGlRUIA1zHLs2LGunWzuCkXBrfMKO3fupKqqKqG7j05+7mROfs5nEw35UCLkSbuPVIPCwkJeeeUVV7pDQkUhEAi4EZovjRkzhp49e7JkyRJmzJjR4f0k+nBUgEsnXOp1CE1c6r+QfJmnSI5aCiLyRxH5QkRWicjLItI77LUbRGSdiKwVkR+HrZ9sr1snIteHrR8uIh/a6+eJiLt3MFGtKiwspKKiouF8gBOlpaUMGDCA9PR0FyLzp+TkZI466ijHVzYn8s11QqaNnca0sf6a7G3aNGvxEz/mKZLT7qOFwFhjTD7wJXADgIiMAaYDhwCTgYdFJFlEkoGHgCnAGOBse1uA/wHuNcaMAHYAFzqMTbWTm9Nod+bhqOEmTpzIZ599xq5duzq8j87QUqjcW0nl3kqvw2iktNRa/MSPeYrkqCgYY942xtTaTz8AQn0FU4G5xph9xphvgHXAEfayzhiz3hizH5gLTBWrr+J4YL79/qeB05zEptovPz8fEdGi0A5FRUXU19fz0UcfdXgfJSUlpKenM2DAABcji6+pc6cyde5Ur8NoZMYMa/ETP+YpkpvnFC4A5tmPc7CKRMhGex1AacT6I4FsYGdYgQnfvgkRuRi4GBL725XfdO/endGjR7tWFBJ1Guj2OPLIIxERli5dyqRJHbt5Sug+Col8H+srj7zS6xASQiLkqdWiICLvAIOaeekmY8yr9jY3AbXAc+6G1zxjzKPAowATJkww8ThmV1FYWMh//vMfR/vYtWsXu3bt6hIthV69enHIIYc4GoGU6MNRAU4/+HSvQ0gIiZCnVruPjDGTjDFjm1lCBWEmcApwjjEm9AFdBoR/IgTsddHWVwC9RSQlYr2Ks8LCQkpLS6moqOjwPrrCcNRwRUVFLF26lPr6+g69P5HvuBaybc82tu3x2URDPpQIeXI6+mgycB3wU2NM+BUZrwHTRSRdRIYDI4GPgGXASHukURrWyejX7GLyHhCaqeQ84FUnsamOceNkc1csCpWVlXzxxRftfu/+/fvZvHlzQl+jAHDGi2dwxos+m2jIhxIhT07PKTwIpAML7f7QD4wxvzTGfC4iLwJrsLqVLjfG1AGIyBXAP4Bk4AljTOjOLr8F5orIncCnwF8cxqY6ILwodLSPvKsVhQkTJgDwySefMGbMmFa2bmzjxo0YYxK+pfDrib/2OoQmfu2/kHyZp0iOioI9fDTaa3cBdzWz/g3gjWbWr8canaQ8lJ2dTW5uruOWQlJSEkOGDHExMv8aPXo06enprFy5knPPbd/cNp1hOCrAqaNP9TqEJk71X0i+zFMkneZCNVFYWOi4KAwePJiUlK5xwXxKSgqHHHIIq1atavd7E/2OayFbqrawpWqL12E0snattfiJH/MUSYuCaqKwsJC1a9fy3Xffdej9XeUahXD5+fmsXLmy3e8LFYVEnw5k+vzpTJ8/3eswGrnkEmvxEz/mKVLX+Cqn2mX8+PEYY1i5ciVFRUXtfn9paSkFBQUxiMy/CgoKeOqppwgGgwwcOLDN79uwYQMDBw4kIyMjhtHF3vXHXN/6Rioh8qQtBdWEkxFIxpgu2VIIFcH2diF1huGoAJNHTGbyiMleh+F7iZAnLQqqiUAgQHZ2doeKwvbt26muru5yRSE/Px+g3V1IiX7HtZDSylJKK3020ZAPJUKetCioJkSkwyebu9pw1JDs7GxycnLaVRSMMZ2mpTDj5RnMeNlnEw35UCLkSc8pqGYVFhZy//33U1NTQ2pqapvf11WLAlhdSO3pPqqoqGDPnj2doijc/MObvQ6hiZv9F5Iv8xRJi4JqVmFhIfv372fNmjXtOmnclYtCfn4+CxcuZP/+/aSltX47kM5yjQLApAM6dqFjLHXw2suY8mOeImn3kWpW6GRze+/ZXFpaSmpqartG4HQWBQUF1NTUUFxc3KbtO8s1CgDrd6xn/Q5n96p224oV1uInfsxTJC0KqlkjR46kW7du7T6vUFpaSk5ODklJXe+fVuhkc1u7kDrDHddCLnj1Ai549QKvw2jk6qutxU/8mKdI2n2kmpWcnExBQUG7i0JJSUmX7DoCGDVqVMN0F225Z3NJSQmZmZlkZ2fHIbrYuu2427wOISEkQp60KKioCgsLefrpp6mvr2/zN//S0tIOXfDWGaSkpDB27Ng2j0AKDUdN5JvrhBw77FivQ0gIiZCnrtfGV21WWFhIVVUVX3/9dZu2r6+vp6ysrMu2FMDqQmpr91FnGY4KsHbbWtZu89lEQz6UCHnSoqCiau+VzcFgkJqami5dFAoKCigvL2fLltYnPesMd1wLueTvl3DJ33020ZAPJUKetPtIRTV27Fh69uzJggULOOuss1rdvisPRw0JDd9duXIlgwY1dxdby969ewkGg52mKPz+hN97HUITv/dfSL7MUyRtKaio0tLSOOecc5g/fz47duxodXstCm0fgbRx40agcwxHBSjKLaIo11/nkoqKrMVP/JinSFoUVItmzZrF3r17ee6551rdVosC9O3bl0Ag0OrJ5s40HBVgdflqVpev9jqMRpYssRY/8WOeImlRUC0qLCxk/PjxPPbYY1i30o6utLSUjIyMTjHE0omCgoJWi0JnupoZ4Io3ruCKN67wOoxGbrzRWvzEj3mKpOcUVKtmzZrFpZdeyvLlyzn88MOjbheaMrszDLF0Ij8/n3/84x/s27eP9PT0ZrcpKSlBRBL+5johfzzxj16HkBASIU/aUlCtOvvss8nKyuKxxx5rcbuueB+F5hQUFFBbW9vidBclJSUMHjy4TXMkJYLDcw7n8JzoXxiUJRHypEVBtapXr16cddZZvPDCC1RVVUXdTouCJXwEUjSdaTgqwIotK1ixxWcTDflQIuRJi4Jqk1mzZlFVVcW8efOafb22tpbNmzdrUQBGjBhBRkZGiyOQOtOFawBXv3U1V7/ls4mGfCgR8qTnFFSbTJw4kYMPPpjHH3+cCy+8sMnrmzZtor6+XosCrU93Ebq5zmmnnRbnyGLnvsn3eR1CE/f5LyRf5imSthRUm4gIs2bN4oMPPmD16qZD6nQ4amOhEUjNjdgqLy9n3759naqlMG7QOMYNGud1GI2MG2ctfuLHPEXSoqDabMaMGaSlpTV7wlmLQmP5+fls27at2ekuOttwVIBlZctYVrbM6zAaeecda/ETP+YpkhYF1Wb9+vXjZz/7Gc8++yx79+5t9JoWhcZaOtncGYvCbxb+ht8s/I3XYTRy553W4id+zFMkLQqqXS666CJ27NjBSy+91Gh9aWkpPXr0oFevXh5F5i+h6S5aKgqdZYoLgAdPfpAHT37Q6zB8LxHypEVBtcvxxx/P8OHDefzxxxut1+GojfXp04fc3NxmRyBt2LCB7t2707t3bw8ii42xA8YydsBYr8PwvUTIkxYF1S5JSUlceOGFvPfee6xbt65hvRaFpqJNdxEajtqZrvxeUrqEJaU+m2jIhxIhT1oUVLudf/75JCcn85e//KVhnRaFpgoKCvjiiy+anH8J3XGtM7nx3Ru58V2fTTTkQ4mQJ71OQbXbkCFD+MlPfsKTTz7J7bffTn19PeXl5VoUIuTn51NXV0dxcXHDDYvA6j6aMGGCh5G578+n/NnrEJr4s/9C8mWeImlLQXXIRRddRDAY5PXXX2+4N4AWhcaaG4G0Z88etm3b1qlGHgGM7jea0f1Gex1GI6NHW4uf+DFPkbQoqA6ZMmUKQ4YM4bHHHtPhqFGMGDGCzMzMRkUhlKvOVhQWf7uYxd8u/NoLTQAACalJREFU9jqMRhYssBY/8WOeImn3keqQlJQUzj//fO6++26OPvpoQItCpOTkZMaOHdtoBFJnHI4KcOuiWwFYNHORt4GE+dOfrJ+nnuptHOH8mKdIWhRUh1144YXcdddd/O///i+gRaE5BQUFvPzyyxhjEJFOd8e1kCemPuF1CAkhEfKk3Ueqw4YPH86kSZOoqKigb9++ZGVleR2S7+Tn51NRUcGmTZsAq6WQlJTEkCFDPI7MXQf0OYAD+hzgdRi+lwh50qKgHJk1axagrYRoQiebQ11IJSUlDBkyhNTUVC/Dct0769/hnfU+m2jIhxIhT64UBRH5tYgYEelnPxcReUBE1onIKhEZH7bteSLylb2cF7b+MBH5zH7PA9KZruzpxKZOnUq/fv044AB/f/vxSuR0F53xGgWAO/91J3f+y2cTDflQIuTJ8TkFEckFTgJKwlZPAUbay5HAI8CRItIXuBWYABjgYxF5zRizw95mFvAh8AYwGXjTaXwqttLT0/nnP/9Jjx49vA7Fl3r37k1eXl5DUdiwYQNHHnmkx1G579mfPet1CE0867+QfJmnSG60FO4FrsP6kA+ZCjxjLB8AvUVkMPBjYKExZrtdCBYCk+3XehpjPjDWBPTPAJ3nDiSd3KGHHsqwYcO8DsO3QtNd1NfXU1pa2ulOMgPk9solt5e/uhBzc63FT/yYp0iOioKITAXKjDGRE7zkAKVhzzfa61pav7GZ9dGOe7GILBeR5Vu3bnXwGygVewUFBaxdu5YNGzZQU1PTKbuP3lr3Fm+te8vrMBqZN89a/MSPeYrUaveRiLwDDGrmpZuAG7G6juLKGPMo8CjAhAkTmt7aSikfyc/Pp76+njfftHpDO2NL4Z7/3APA5BGTPY7ke488Yv2cNs3bOML5MU+RWi0KxphJza0XkUOB4cBK+5xwAPhERI4AyoDwNlLAXlcGHBexfpG9PtDM9kolvNAIpNdeew3onEVh7hlzvQ4hISRCnjrcfWSM+cwYM8AYM8wYMwyry2e8MWYL8BrwC3sU0lFApTFmM/AP4CQR6SMifbBaGf+wX9slIkfZo45+Abzq8HdTyhcOPPBAsrKyeO+994DOWRQGdR/EoO7NdSiocImQp1hdp/AGsB5YBzwGXAZgjNkO3AEss5fb7XXY2zxuv+drdOSR6iRC013s37+fXr16dcq70y1Yu4AFa3020ZAPJUKeXJvmwm4thB4b4PIo2z0BNLnW2xizHPD3LYmU6qCCggI++uijTtlKAPjTUmuioVNH+2iiIR9KhDzp3EdKxUHovEJnLQrzz5rvdQhNzPdfSL7MUyQtCkrFQejK5s5aFPpl9fM6hCb6+S8kX+Ypks59pFQc5Ofnk5mZyZgxY7wOJSZeKn6Jl4pf8jqMRp56ylr8xI95iqQtBaXioFevXqxdu5ZBg/w98qSjHvjwAQBOP/h0jyP5XqggzJzpZRSN+TFPkbQoKBUnnXkm2Ven6wjytkiEPGlRUEo51iuj8w2zjYVEyJOeU1BKOTZv9TzmrfbZREM+lAh50paCUsqxR5ZbEw1NG+ujiYZ8KBHyJNZ1ZolrwoQJZvny5V6HoVSXtqdmDwBZqf65JeseKyT8dJdYP+VJRD42xkyIXK8tBaWUY374kIvkp2IQ4sc8RdJzCkopx+asmsOcVXO8DqORhx+2Fj/xY54iaVFQSjn2+CeP8/gnj3sdRiMvvmgtfuLHPEXS7iOllGMLZyz0OoSEkAh50qKglHIsNTnV6xASQiLkSbuPlFKOPbXiKZ5a8ZTXYfheIuRJi4JSyrFE+LDzg0TIk16noJRSXVC06xS0paCUUqqBFgWllFINtCgopZRqoEVBKaVUAy0KSimlGmhRUEop1UCLglJKqQZaFJRSSjXQoqCUUqqBFgWllFINtCgopZRqoEVBKaVUAy0KSimlGmhRUEop1SDhp84Wka3Ahg6+vR+wzcVw4kXjji+NO74SNW5IrNiHGmP6R65M+KLghIgsb24+cb/TuONL446vRI0bEjv2EO0+Ukop1UCLglJKqQZdvSg86nUAHaRxx5fGHV+JGjckduxAFz+noJRSqrGu3lJQSikVpksWBRGZLCJrRWSdiFzvdTztISLfishnIrJCRJZ7HU80IvKEiJSLyOqwdX1FZKGIfGX/7ONljM2JEvfvRKTMzvkKETnZyxibIyK5IvKeiKwRkc9F5Cp7va9z3kLcvs65iGSIyEcistKO+zZ7/XAR+dD+bJknImlex9peXa77SESSgS+BE4GNwDLgbGPMGk8DayMR+RaYYIzx9VhoEfkhUAU8Y4wZa6/7A7DdGHOPXYz7GGN+62WckaLE/Tugyhgz28vYWiIig4HBxphPRKQH8DFwGjATH+e8hbjPwsc5FxEBuhljqkQkFfgPcBVwDfCSMWauiPwfsNIY84iXsbZXV2wpHAGsM8asN8bsB+YCUz2OqdMxxvwL2B6xeirwtP34aaz//L4SJW7fM8ZsNsZ8Yj/eDRQDOfg85y3E7WvGUmU/TbUXAxwPzLfX+y7fbdEVi0IOUBr2fCMJ8I8wjAHeFpGPReRir4Npp4HGmM324y3AQC+DaacrRGSV3b3kqy6YSCIyDCgEPiSBch4RN/g85yKSLCIrgHJgIfA1sNMYU2tvkmifLUDXLAqJ7hhjzHhgCnC53d2RcIzVb5kofZePAAcC44DNwJ+8DSc6EekO/A242hizK/w1P+e8mbh9n3NjTJ0xZhwQwOqBOMjjkFzRFYtCGZAb9jxgr0sIxpgy+2c58DLWP8ZEEbT7kEN9yeUex9Mmxpig/QFQDzyGT3Nu923/DXjOGPOSvdr3OW8u7kTJOYAxZifwHjAR6C0iKfZLCfXZEtIVi8IyYKQ9SiANmA685nFMbSIi3eyTcYhIN/j/27l/lAaCKI7j30ds7ESwtfAAnsAilb0ggiCk9A42gmArXkBLlVTqFTyAhUJaW49g47OYcUhhQqwykO+nWvYPPB7s/nZnhmUfeJ9/VVeegVHdHgFPS6xlYb8P1eqADnteJz5vgElmXk0d6rrns+ruvecRsRURG3V7nbJwZUIJh8N6Wnf9XsTKrT4CqMvbroEBcJuZl0suaSERsUP5OgBYA+56rT0i7oEh5a+Rn8A58AiMgW3Kn22PMrOrSd0ZdQ8pwxgJfACnU+P0XYiIPeAFeAO+6+4zyvh8tz2fU/cxHfc8InYpE8kDysv1ODMv6j36AGwCr8BJZn4tr9L/W8lQkCT9bRWHjyRJMxgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkpofBf0ecaRvJY4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "y = np.arange(image.shape[0])\n",
    "ax.plot(y, image[:,50], 'k-')\n",
    "mm = np.array([ext_center, ext_center])\n",
    "mm_y = ax.get_ylim()\n",
    "ax.plot(mm, mm_y, 'b--')\n",
    "ax.plot(mm - ext_width/2., mm_y, 'g:')\n",
    "ax.plot(mm + ext_width/2., mm_y, 'g:')\n",
    "ax.set_title(\"Cross-dispersion Cut at Pixel=50\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do the extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Fixed boxcar backgound weight image')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAABbCAYAAAC8uruvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAPSElEQVR4nO3de7RcZXnH8e+PkxuQlFtoDAQIaICmLoo0XFSK3JRLxWClNEK5VLpYSwpKK0KCFrFLrNiCWgvYlIQkXIMJlouwlAALkrZyNRACjUQShJAQuQRCwEDg6R/ve8hmPHPm7Jw5M+fM/D5r7ZW933fP3s/77n3OebLfd2YUEZiZmZlZz23W7ADMzMzMBhonUGZmZmYlOYEyMzMzK8kJlJmZmVlJTqDMzMzMSnICZWZmZlaSEyizfkzS65J2q/MxD5b0XDf1IelD9TxnfyfpfElX9nDfCyVd04tz/ZmkJZv6ejPrHwY1OwAzA0nLgVHAO4Xi3SNieHMiai8R8e16HStfy7+NiHlVzjUf2KNe5zOz5vATKLP+45iIGF5Ynm92QP2NpI5mx2BmBk6gzPq1zuE0SUMkLZR0Vi7vkPTfki7I2ztImivpt5KWSfpS4RibS5oh6RVJTwD79uDUR0t6WtKLkv5F0mb5WJtJ+rqkZyStljRL0la57nZJlxTOe4Ok6YUYLsmve1XSAkmb57ofS1qVy++T9MeFY8yQdEU+9jrgkIr+OUTSosL2nZIeLGzPl3RsD/rofcNykk7Osb4k6R8lLZd0eOHUQ3Lb10paLGlCft3VwM7ArXn49dwurun7hlDzsb8q6TFJ6yRNkzRK0h35+PMkbVPYv7v+2k7SrZJek/SgpG9JWlCo3zP30cuSlkg6vuvLb2a1OIEyGwAi4i3gr4F/kvRHwGSgA7goJze3Ao8COwKHAWdLOiK//BvAB/NyBHBKD075WWACsA8wEfhCLj81L4cAuwHDgX/PdV8ATpJ0qKQTgf2AL+e6fwX+FPgYsC1wLvBurrsDGAf8IfAIcG1FLCcAFwEjgAUVdb8AxkkaKWkwsBewg6QROUGbAMzvQR+9R9J44HLgRGA0sFV+TdFngBuArYFbOvsgIk4CfsPGp4nfrTx+FZ8DPgnsDhyT++R8YHvS7+kvFfbtrr8uA9YBHyBd5/eutaQtgTuB6/JrJwGX5/aaWVkR4cWLlyYvwHLgdWBNXv4rlwfwocJ+XwGWAK8A43LZ/sBvKo43Bbgqrz8NHFmoOx14rptYomL/M4C78vpdwBmFuj2At4FBeftzwLPAi8CBuWwz4E3gT3rQD1vn82+Vt2cAs2q8Zj7wF8ABwM+BG4EjSUneYz3sowuBa/L6BcD1hf22AN4CDi/sO69QPx54s+JaHt5NvAcX+z/vf2Jhey5wRWH7rM77obv+IiXUbwN7FOq/BSzI638FzK94/X8A32j2/e/Fy0BcPIncrP84NqpMPC6YSXoaMzcinsplu5Ceuqwp7NdBSiwAdiAlNZ2e6UEslfvvUDjWMxV1g0gT4FeQnvL8EFgSEZ1Pi0YCw4BfV54kz2m6CPhL0tOWdwuvebWLWLpyLzkpyeuvAJ8A1udtqN1HRe/rr4h4Q9JLFfusKqy/AQyTNCgiNtSItZoXCutvdrE9HGr21+aka1Hsr+L6LsD+FX0wCLh6E2M2a2sewjMbWC4HbgOOkHRgLnsWWBYRWxeWERFxdK5fCexUOMbOPThP5f6dE9qfJ/0hLtZtYOMf/IuAJ4HRkj6fy14EfkcaQqx0AmmI8HDSU5SxuVyFfaJGrJ0J1EF5/V5SAvUJNiZQtfqoaCUwpnMjDwVuVyOGolrx9kZ3/fVb0rUYU9i/eB2fBe6t6IPhEfHFPozXrGU5gTIbICSdRJpHdCppTsxMScOBB4C1ks7Lk7U7JH1YUudk8RuBKZK2kTSGNCRUy1fz/juR5jHNzuXXA38vadd87m8DsyNig6SDgL8BTibNvfmhpB0j4l1gOnBpnsjdIemjkoaS5jWtB14iDZVtyscJ/A9pKHE/4IGIWEx+2gLcl/ep1UdFc4BjJH1M0hDSkJ262K+aF0jzw/pC1f6KiHeAm4ALJW0haU/Steh0G7C7pJMkDc7LvnlOnZmV5ATKbACQtDPwfeDkiHg9Iq4DHgK+l/9wfhrYG1hGeuJzJekJBcA3SUNty0hzhHoyZHMz8DCwEPgpMC2XT8+vvy8f73fAWZL+AJgFnBkRKyJ91tE04CpJAs4BFgEPAi8DF5N+/8zKsa0AniBNCi8lItaRJlMvjjTZHuB/gWciYnXep1YfFY+3mJRk3kB6GvU6sJqUuPTEPwNfl7RG0jll21NDrf46k9SmVaTrdD057ohYC3yKNHn8+bzPxcDQOsdo1hYU0ZdPm83MBrb8pG0NadL+smbHU4aki4EPRERP3nlpZiX4CZSZWQVJx+RhsC1JH8GwiPRuuX4tf87TXkr2A04DftLsuMxakRMoM7PfN5E0zPU86TOXJsXAeFw/gjQPah1p3tolpOFYM6szD+GZmZmZleQnUGZmZmYlOYEyMzMzK6mhn0Q+RENjGFs28pRmZmZmm2Qtr7wYEdt3VdfQBGoYW7K/DmvkKc3MzMw2ybyYU/WrrzyEZ2ZmZlZSQ59AacgQBo3ZpfaOZmZmZs32dPWqhiZQ4/Zcw09/5o8kMTMzs/6vY3T1Og/hmZmZmZXkBMrMzMysJCdQZmZmZiU5gTIzMzMryQmUmZmZWUlOoMzMzMxKcgJlZmZmVpITKDMzM7OSnECZmZmZlVTzk8gl7QTMAkYBAUyNiB9I2haYDYwFlgPHR8Qr3R1rybKRHHryab2N2czMzKwBplSt6clXuWwAvhIRj0gaATws6U7gVOCuiPiOpMnAZOC87g6ktW8weN7DPQ7bzMzMrD+qOYQXESsj4pG8vhZ4EtgRmAjMzLvNBI7tqyDNzMzM+pNSc6AkjQU+AtwPjIqIlblqFWmIz8zMzKzl9TiBkjQcmAucHRGvFesiIkjzo7p63emSHpL00Nus71WwZmZmZv1BjxIoSYNJydO1EXFTLn5B0uhcPxpY3dVrI2JqREyIiAmDGVqPmM3MzMyaqmYCJUnANODJiLi0UHULcEpePwW4uf7hmZmZmfU/PXkX3seBk4BFkhbmsvOB7wA3SjoNeAY4vm9CNDMzM+tfaiZQEbEAUJXqw+objpmZmVn/508iNzMzMyvJCZSZmZlZSU6gzMzMzEpyAmVmZmZWkhMoMzMzs5KcQJmZmZmV5ATKzMzMrCQnUGZmZmYlOYEyMzMzK8kJlJmZmVlJTqDMzMzMSnICZWZmZlaSEygzMzOzkpxAmZmZmZXkBMrMzMysJCdQZmZmZiX1KoGSdKSkJZKWSppcr6DMzMzM+rNNTqAkdQCXAUcB44HPSxpfr8DMzMzM+qvePIHaD1gaEU9HxFvADcDE+oRlZmZm1n/1JoHaEXi2sP1cLjMzMzNraYP6+gSSTgdOBxjGFn19OjMzM7M+15sEagWwU2F7TC57n4iYCkwFkLR2XsxZ0otztoKRwIvNDqKJ3H633+1vX25/e7cfBl4f7FKtojcJ1IPAOEm7khKnScAJNV6zJCIm9OKcA56kh9q5D9x+t9/td/ubHUeztHv7obX6YJMTqIjYIOlM4GdABzA9IhbXLTIzMzOzfqpXc6Ai4nbg9jrFYmZmZjYgNPqTyKc2+Hz9Ubv3gdvf3tz+9ub2W8v0gSKi2TGYmZmZDSj+LjwzMzOzkhqWQLXj9+ZJWi5pkaSFkh7KZdtKulPSU/nfbZodZ71Imi5ptaTHC2VdtlfJv+X74TFJ+zQv8vqp0gcXSlqR74OFko4u1E3JfbBE0hHNibo+JO0k6R5JT0haLOnLubwt7oFu2t8W1x9A0jBJD0h6NPfBN3P5rpLuz22dLWlILh+at5fm+rHNjL+3umn/DEnLCvfA3rm8pX4GOknqkPRLSbfl7da8/hHR5wvpXXq/BnYDhgCPAuMbce5mLsByYGRF2XeByXl9MnBxs+OsY3sPAvYBHq/VXuBo4A5AwAHA/c2Ovw/74ELgnC72HZ9/FoYCu+afkY5mt6EXbR8N7JPXRwC/ym1si3ugm/a3xfXPbRIwPK8PBu7P1/ZGYFIu/xHwxbx+BvCjvD4JmN3sNvRR+2cAx3Wxf0v9DBTa9Q/AdcBtebslr3+jnkD5e/M2mgjMzOszgWObGEtdRcR9wMsVxdXaOxGYFckvgK0ljW5MpH2nSh9UMxG4ISLWR8QyYCnpZ2VAioiVEfFIXl8LPEn6eqe2uAe6aX81LXX9AfK1fD1vDs5LAIcCc3J55T3QeW/MAQ6TpAaFW3fdtL+alvoZAJA0Bvhz4Mq8LVr0+jcqgWrX780L4OeSHlb6ShuAURGxMq+vAkY1J7SGqdbedrsnzsyP6KcXhm1btg/yo/iPkP4H3nb3QEX7oY2ufx6+WQisBu4kPVlbExEb8i7Fdr7XB7n+VWC7xkZcX5Xtj4jOe+CifA98T9LQXNaK98D3gXOBd/P2drTo9fck8r51YETsAxwF/J2kg4qVkZ5bts3bINutvQVXAB8E9gZWApc0N5y+JWk4MBc4OyJeK9a1wz3QRfvb6vpHxDsRsTfp6732A/ZsckgNVdl+SR8GppD6YV9gW+C8JobYZyR9GlgdEQ83O5ZGaFQC1aPvzWs1EbEi/7sa+Anpl8kLnY9o87+rmxdhQ1Rrb9vcExHxQv6l+i7wn2wcpmm5PpA0mJQ8XBsRN+XitrkHump/O13/oohYA9wDfJQ0NNX5wc3Fdr7XB7l+K+ClBofaJwrtPzIP70ZErAeuonXvgY8Dn5G0nDRV51DgB7To9W9UAvXe9+bl2feTgFsadO6mkLSlpBGd68CngMdJ7T4l73YKcHNzImyYau29BTg5vwvlAODVwjBPS6mY0/BZ0n0AqQ8m5Xei7AqMAx5odHz1kucuTAOejIhLC1VtcQ9Ua3+7XH8ASdtL2jqvbw58kjQX7B7guLxb5T3QeW8cB9ydn1IOSFXa/3+F/0CINP+neA+0zM9AREyJiDERMZb0d/7uiDiRVr3+jZqtTnq3wa9I4+Ffa9R5m7WQ3nH4aF4Wd7aZNL57F/AUMA/Yttmx1rHN15OGKN4mjXOfVq29pHedXJbvh0XAhGbH34d9cHVu42OkXxijC/t/LffBEuCoZsffy7YfSBqeewxYmJej2+Ue6Kb9bXH9c3v2An6Z2/o4cEEu342UHC4FfgwMzeXD8vbSXL9bs9vQR+2/O98DjwPXsPGdei31M1DRFwez8V14LXn9/UnkZmZmZiV5ErmZmZlZSU6gzMzMzEpyAmVmZmZWkhMoMzMzs5KcQJmZmZmV5ATKzMzMrCQnUGZmZmYlOYEyMzMzK+n/AVifZAeD0RRNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAABbCAYAAAC8uruvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAR4klEQVR4nO3deZhcVZ3G8e+bDhAgMSEEMQlIWCJOYBQyyERBQEWWDBgQxicuEBSHmVG2GZBFGAYUfNARZBxBhjXsW1BRhBl2GFDCIiGExEggiRBCQkLCvoX85o9zityUXd1d6eqq7qr38zz3yb3n3HvPck91/XLvqSpFBGZmZmbWdf0aXQEzMzOzvsYBlJmZmVmVHECZmZmZVckBlJmZmVmVHECZmZmZVckBlJmZmVmVHECZlZH0mqQtanzO3SQ910F+SNqqlmX2JEn3SPpmjc85KvdD/1qet9Z6ou1l5/+0pNld3LfDcdXFc9R8vJu1gl79h8qsJ0maB2wMvFdI/khEDGxMjcwgIv4P2LoW55I0GXguIk7uoDyPd7M14DtQ1ur2jYiBheX5Rleot5HU1ug6mJn1Ng6gzMqUHqdJWlvSNElH5PQ2SQ9IOiVvj5B0o6QXJc2VdGThHOtKmixpmaSZwCe6UPR4Sc9IWiLpPyT1y+fqJ+lkSfMlLZZ0uaTBOe8WSWcVyr1W0iWFOpyVj3tZ0v2S1s15N0h6IaffJ2mbwjkmS/p5PvfrwGcq1HdLSQ9JekXSTZKGFs7R0fkr1qvsOhwgaZ6kbfP2wfmYpZL+LeftnvPWkXSOpOfzco6kdXLeIZLub+8aF9p7rqTfSnpV0lRJWxb2/bykP+a6/gxQe50haYCkNyUNy9snSVoh6QN5+/uSzinU98eS/ixpkaTzC9dmtcdyksZKeizX7QZJ10k6vazsY/LYWCjp6zntMOCrwHH5Md1vKtS7vC/Ok3RrPuYBSR/K/bks98P2hWNPkPR0rttMSfsX8trydV6i9Po4XIVHtJIGS7o413mBpNPlYN36EAdQZhVExDvA14DvSfor4ASgDThDKbj5DfA4MBL4HHC0pD3z4f8ObJmXPYFJXShyf2AHYCwwAfhGTj8kL58BtgAGAj/Led8ADpL0WUlfBXYEjsp5Pwb+BvgUMBQ4DliZ824FRgMfBP4AXFVWl68AZwCDgPtp38G5/OHACuCnhbyOzt9RvQDIQcAPgd0jYoakMcB5pIBgODCY1O8lJwHjgO2Aj+d+qPjYqh0TgdOADYA5pLaTg6Ff5HMNA54GdmrvBBHxFvAwsGtO2hWYX9h/V+DevH4m8JFc361yW04pP6ektYFfApNJfXUNaZwUfYhV/XEocK6kDSLiAlK//yjfXd23Sz0BXyq0923g96RrOAyYApxd2Pdp4NO5/NOAKyUNz3n/AOyd2zgW2K+snMmkcbMVsD2wB9Bjc8vMai4ivHhpyQWYB7wGLM/Lr3J6AFsV9jsGmA0sA0bntL8F/lx2vhOBS/P6M8BehbzDSHNRKtUlyvb/FnBnXr8T+FYhb2vgXaB/3j4AeBZYAuyc0/oBbwIf70I/DMnlD87bk4HLOznmHuDMwvYY4B2graPzd1QvYFTe71hgJrBJIe8U4JrC9nq5vN3z9tPA+EL+nsC8vH4IcH87/b1Vob0XFfLGA3/M6wcDDxbyBDwHfLNCv3yfFEj2B14gBbNnAgNyuzfM53gd2LJw3CeBuXl9t9JYAXYBFgAq7Hs/cHph3zdLYyGnLQbGFdp2eifXsrwvLizkHQHMKmz/NbC8g3NNAybk9buAfyzk7Z7L6k+ae/g2sG4h/8vA3T3xWvfipScW34GyVrdfRAzJS/n/kEsuAzYDbomIp3LaZsAISctLC/Bd0hsDwAhSUFMyvwt1Kd9/ROFc88vySm9CkO6EtQGzI6J0t2gY6U376fJC8qOVM/Ojl1dIgWTpmPbq0tX6rgUM6+T8FetV8B3g3Igofrpstf6MiDeApWX55X00gq57obD+BukuX3vlBh33zb2koGYs8ARwO+nO0zhgTkQsBTYiBYCPFsbO/+T0ciOABbnckvLyl0bEigr1XxOLCutvtrP9/rnzY9VphXZsy6pxVP4aKK5vRhovCwvH/jfpjqVZn+AAyqxz5wE3A3tK2jmnPUu6YzCksAyKiPE5fyGwaeEcH+5COeX7lya0P096wynmrWDVG9sZwCxguKQv57QlwFukR4jlvkJ6RLg76a7QqJxenNsTdK68vu/mcjs6f0f1KtkDOFnSAYW0hcAmpY08X2jDQn57fVTqv9dJAUvp2A913KzVrHYdJYnV213ud6Q7hPsD90bEzFyX8ax6fLeEFIhsUxg7g6P9T8MtBEbmcks6Kr9cV67jGpG0GXAhcDiwYUQMAWawahytds1Yvd7Pku5ADSv0wQciYhvM+ggHUGYdkHQQab7OIcCRwGWSBgIPAa9KOl5pUnSbpG0llSaLXw+cKGkDSZuQHoV05jt5/01Jj36uy+nXAP8iafNc9g+A6yJihaRdgK+THjVNAv5L0siIWAlcApytNNm9TdInlSZWDyK9eS0lBRY/WMPu+ZqkMZLWA74HTImI9zo6fyf1KnkS2Is0l+cLOW0KsK+kT+V5QaeyesB3DSno2ijPWzoFuDLnPQ5sI2k7SQPysV3123zsF/Pk5yNJc47ale+MPQp8m1UB0++Afypt5z64EPiJpA8CSBpZmD9X9HvS12wcLqm/pAmk+V1dtYg0b64nrE8K0F6E9+etbVvIvx44KrdtCHB8KSMiFgK3AWdJ+oDSByW2lLQrZn2EAyizCiR9GDgHODgiXouIq4FHgJ/kQGEf0gTZuaS7CheR7rhAmlA7P+fdBlzRhSJvIr35TiO9cV+c0y/Jx9+Xz/cWcITSp7suBw6PiAWRvj/oYuDSfMfiWNJjpIeBl0iTsvvlY+aT5tbMBB6sunOSK0hzZl4gPZYrfQqxs/NXqtf7IuJxUv9eKGnviHiSFIReS7qz8Rpprs/b+ZDTSddmej73H3IaEfEnUoB3B/AUlSfF/4WIWAL8PWke01LSxPgHOjnsXtLjqYcK24NI16/keNJk9QfzY847aOe7nyJ9kOGLpMnhy0kfariZVe3uzMXAmPyY7FddPKZL8t21s0hB3iLS/Khi31xIGvvTgceAW0h3Tkvfu3YwsDZpjCwjBcnDMesjtPqjdTOz3i/fiVtOmtQ/t9H1qSdJU4HzI+LSRtelGpL2JtV7s053NusDfAfKzPoESftKWk/S+qSvQniCVRPUm5akXfN3MfWXNAn4GGnSea+WH22Pz/UeSfpqj182ul5mteIAysz6igmkieHPkx6lTYzWuIW+NWke13LSV2ocmOcQ9XYiPcpeRnqEN4t2vuvKrK/yIzwzMzOzKvkOlJmZmVmVHECZmZmZVal/PQsbNrQtRm26Vj2LNDMzM1sjj05/e0lEtPcrAfUNoN5iBGtvfVTnO5qZmZk12vTvVPwZrroGUP2Wvc76U6bWs0gzMzOzmvMcKDMzM7MqOYAyMzMzq5IDKDMzM7MqOYAyMzMzq5IDKDMzM7MqOYAyMzMzq5IDKDMzM7MqOYAyMzMzq5IDKDMzM7MqOYAyMzMzq5IDKDMzM7MqdRpASdpU0t2SZkp6UtJROX2opNslPZX/3aDnq2tmZmbWeF25A7UCOCYixgDjgG9LGgOcANwZEaOBO/O2mZmZWdPr39kOEbEQWJjXX5U0CxgJTAB2y7tdBtwDHN/RubTuAPp9dEw3qmtmZmZWJ49Vzuo0gCqSNArYHpgKbJyDK4AXgI07O3706Je49ZarqynSzMzMrCHahlfO6/IkckkDgRuBoyPilWJeRAQQFY47TNIjkh55cel7XS3OzMzMrNfqUgAlaS1S8HRVRPwiJy+SNDznDwcWt3dsRFwQETtExA4bbdhWizqbmZmZNVRXPoUn4GJgVkScXcj6NTApr08Cbqp99czMzMx6n67MgdoJOAh4QtK0nPZd4EzgekmHAvOBL3V2ohWsZNl7b6xpXc3MzMx6BaXpS/UxuG1YjFtvn7qVZ2ZmZrambnvtskcjYof28qr6FF53xcqVrHz99XoWaWZmZlZz/ikXMzMzsyo5gDIzMzOrkgMoMzMzsyo5gDIzMzOrkgMoMzMzsyo5gDIzMzOrkgMoMzMzsyo5gDIzMzOrkgMoMzMzsyo5gDIzMzOrkgMoMzMzsyo5gDIzMzOrkgMoMzMzsyo5gDIzMzOrkgMoMzMzsyp1K4CStJek2ZLmSDqhVpUyMzMz683WOICS1AacC+wNjAG+LGlMrSpmZmZm1lt15w7UjsCciHgmIt4BrgUm1KZaZmZmZr1XdwKokcCzhe3ncpqZmZlZU+vf0wVIOgw4DGAA6/V0cWZmZmY9rjsB1AJg08L2JjltNRFxAXABgKRX74gps7tRZjMYBixpdCUayO13+93+1uX2t3b7oe/1wWaVMroTQD0MjJa0OSlwmgh8pZNjZkfEDt0os8+T9Egr94Hb7/a7/W5/o+vRKK3efmiuPljjACoiVkg6HPhfoA24JCKerFnNzMzMzHqpbs2BiohbgFtqVBczMzOzPqHe30R+QZ3L641avQ/c/tbm9rc2t9+apg8UEY2ug5mZmVmf4t/CMzMzM6tS3QKoVvzdPEnzJD0haZqkR3LaUEm3S3oq/7tBo+tZK5IukbRY0oxCWrvtVfLTPB6mSxrbuJrXToU+OFXSgjwOpkkaX8g7MffBbEl7NqbWtSFpU0l3S5op6UlJR+X0lhgDHbS/Ja4/gKQBkh6S9Hjug9Ny+uaSpua2Xidp7Zy+Tt6ek/NHNbL+3dVB+ydLmlsYA9vl9KZ6DZRIapP0mKSb83ZzXv+I6PGF9Cm9p4EtgLWBx4Ex9Si7kQswDxhWlvYj4IS8fgLww0bXs4bt3QUYC8zorL3AeOBWQMA4YGqj69+DfXAqcGw7+47Jr4V1gM3za6St0W3oRtuHA2Pz+iDgT7mNLTEGOmh/S1z/3CYBA/P6WsDUfG2vBybm9POBf87r3wLOz+sTgesa3YYeav9k4MB29m+q10ChXf8KXA3cnLeb8vrX6w6UfzdvlQnAZXn9MmC/BtalpiLiPuClsuRK7Z0AXB7Jg8AQScPrU9OeU6EPKpkAXBsRb0fEXGAO6bXSJ0XEwoj4Q15/FZhF+nmnlhgDHbS/kqa6/gD5Wr6WN9fKSwCfBabk9PIxUBobU4DPSVKdqltzHbS/kqZ6DQBI2gT4O+CivC2a9PrXK4Bq1d/NC+A2SY8q/aQNwMYRsTCvvwBs3Jiq1U2l9rbamDg836K/pPDYtmn7IN+K3570P/CWGwNl7YcWuv758c00YDFwO+nO2vKIWJF3Kbbz/T7I+S8DG9a3xrVV3v6IKI2BM/IY+ImkdXJaM46Bc4DjgJV5e0Oa9Pp7EnnP2jkixgJ7A9+WtEsxM9J9y5b5GGSrtbfg58CWwHbAQuCsxlanZ0kaCNwIHB0RrxTzWmEMtNP+lrr+EfFeRGxH+nmvHYGPNrhKdVXefknbAieS+uETwFDg+AZWscdI2gdYHBGPNrou9VCvAKpLv5vXbCJiQf53MfBL0h+TRaVbtPnfxY2rYV1Uam/LjImIWJT/qK4ELmTVY5qm6wNJa5GCh6si4hc5uWXGQHvtb6XrXxQRy4G7gU+SHk2Vvri52M73+yDnDwaW1rmqPaLQ/r3y492IiLeBS2neMbAT8AVJ80hTdT4L/CdNev3rFUC9/7t5efb9RODXdSq7ISStL2lQaR3YA5hBavekvNsk4KbG1LBuKrX318DB+VMo44CXC495mkrZnIb9SeMAUh9MzJ9E2RwYDTxU7/rVSp67cDEwKyLOLmS1xBio1P5Wuf4AkjaSNCSvrwt8njQX7G7gwLxb+RgojY0DgbvyXco+qUL7/1j4D4RI83+KY6BpXgMRcWJEbBIRo0jv83dFxFdp1utfr9nqpE8b/In0PPykepXbqIX0icPH8/Jkqc2k57t3Ak8BdwBDG13XGrb5GtIjindJz7kPrdRe0qdOzs3j4Qlgh0bXvwf74IrcxumkPxjDC/uflPtgNrB3o+vfzbbvTHo8Nx2YlpfxrTIGOmh/S1z/3J6PAY/lts4ATsnpW5CCwznADcA6OX1A3p6T87dodBt6qP135TEwA7iSVZ/Ua6rXQFlf7MaqT+E15fX3N5GbmZmZVcmTyM3MzMyq5ADKzMzMrEoOoMzMzMyq5ADKzMzMrEoOoMzMzMyq5ADKzMzMrEoOoMzMzMyq5ADKzMzMrEr/D7KUsLKUllSdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualize the weight images used in the fixed boxcar extraction \n",
    "wimage_fixedboxcar, wimage_fb_bkg = ap_weight_images(ext_center, ext_width, bkg_offset, \n",
    "                                                     bkg_width, image.shape, None)\n",
    "\n",
    "norm_data = simple_norm(wimage_fixedboxcar)\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.imshow(wimage_fixedboxcar, norm=norm_data, origin=\"lower\")\n",
    "plt.title(\"Fixed boxcar weight image\")\n",
    "\n",
    "norm_data = simple_norm(wimage_fb_bkg)\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.imshow(wimage_fb_bkg, norm=norm_data, origin=\"lower\")\n",
    "plt.title(\"Fixed boxcar backgound weight image\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattr__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__iter__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_ctx', '_instance', '_name', '_schema', '_validate', 'cdelt1', 'cdelt2', 'cdelt3', 'crpix1', 'crpix2', 'crpix3', 'crval1', 'crval2', 'crval3', 'ctype1', 'ctype2', 'ctype3', 'cunit1', 'cunit2', 'cunit3', 'dec_ref', 'dispersion_direction', 'hasattr', 'instance', 'items', 'mt_avdec', 'mt_avra', 'mt_dec', 'mt_ra', 'pc1_1', 'pc1_2', 'pc1_3', 'pc2_1', 'pc2_2', 'pc2_3', 'pc3_1', 'pc3_2', 'pc3_3', 'ps3_0', 'ps3_1', 'ra_ref', 'roll_ref', 's_region', 'siaf_xref_sci', 'siaf_yref_sci', 'specsys', 'spectral_order', 'v2_ref', 'v3_ref', 'v3yangle', 'velosys', 'vparity', 'waverange_end', 'waverange_start', 'wcsaxes']\n",
      "\n",
      "6e-07\n",
      "\n",
      "{'dec_ref': -70.77548664973652, 'dispersion_direction': 1, 'ra_ref': 90.75350083070477, 'roll_ref': -93.7605896, 's_region': 'POLYGON ICRS  90.752386052 -70.775854900 90.754461730 -70.775854900 90.754461730 -70.775171488 90.752386052 -70.775171488', 'specsys': 'BARYCENT', 'spectral_order': 0, 'v2_ref': 332.136078, 'v3_ref': -479.224213, 'v3yangle': 138.7605896, 'velosys': -1398.2, 'vparity': -1, 'waverange_end': 5.3e-06, 'waverange_start': 6e-07}\n",
      "\n",
      "None\n",
      "None\n",
      "None\n",
      "None\n",
      "None\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "No attribute 'wcs'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m~/Projects/software/miniconda3_2/envs/specreduce/lib/python3.8/site-packages/jwst/datamodels/properties.py\u001b[0m in \u001b[0;36m__getattr__\u001b[0;34m(self, attr)\u001b[0m\n\u001b[1;32m    289\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 290\u001b[0;31m             \u001b[0mval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_instance\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mattr\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    291\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: 'wcs'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/f2/dsq1cxxd78n9vltk6_6kwxtr0000zf/T/ipykernel_3463/357120051.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;31m# without background subtraction\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[0mwaves_boxcar\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mext1d_boxcar\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtmpval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mextract_1dspec\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mimage\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ms2d\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mext_center\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mext_width\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m,\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      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;31m# with background subtraction\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/var/folders/f2/dsq1cxxd78n9vltk6_6kwxtr0000zf/T/ipykernel_3463/353467436.py\u001b[0m in \u001b[0;36mextract_1dspec\u001b[0;34m(image, jdatamodel, center, width, bkg_offset, bkg_width, wavescale)\u001b[0m\n\u001b[1;32m    128\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    129\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 130\u001b[0;31m     \u001b[0mgrid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgrid_from_bounding_box\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mjdatamodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmeta\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwcs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbounding_box\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    131\u001b[0m     \u001b[0mra\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdec\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlam\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mjdatamodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmeta\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwcs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mgrid\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    132\u001b[0m     \u001b[0mlam_image\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtranspose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlam\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Projects/software/miniconda3_2/envs/specreduce/lib/python3.8/site-packages/jwst/datamodels/properties.py\u001b[0m in \u001b[0;36m__getattr__\u001b[0;34m(self, attr)\u001b[0m\n\u001b[1;32m    291\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    292\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mschema\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m{\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--> 293\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0mAttributeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"No attribute '{0}'\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mattr\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[0m\u001b[1;32m    294\u001b[0m             \u001b[0mval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_make_default\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mattr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mschema\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_ctx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    295\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mval\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[0;31mAttributeError\u001b[0m: No attribute 'wcs'"
     ]
    }
   ],
   "source": [
    "# extract the spectrum using the weight images\n",
    "\n",
    "# without background subtraction\n",
    "waves_boxcar, ext1d_boxcar, tmpval = extract_1dspec(image, s2d, ext_center, ext_width, None, None)\n",
    "\n",
    "# with background subtraction\n",
    "waves_boxcar_bkgsub, ext1d_boxcar_bkgsub, tmpval = extract_1dspec(image, s2d, ext_center, ext_width, \n",
    "                                                                  bkg_offset, bkg_width)\n",
    "\n",
    "# plot\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "gpts = ext1d_boxcar_bkgsub > 0.\n",
    "ax.plot(waves_boxcar[gpts], ext1d_boxcar[gpts], 'k-', label=\"boxcar\")\n",
    "ax.plot(waves_boxcar_bkgsub[gpts], ext1d_boxcar_bkgsub[gpts], 'k:', label=\"boxcar (bkgsub)\")\n",
    "ax.plot(jpipe_x1d['WAVELENGTH'], jpipe_x1d['FLUX'], 'k-', label=\"jpipe_x1d\")\n",
    "ax.set_title(\"Fixed boxcar 1D extracted spectrum\")\n",
    "ax.set_xlabel(r\"wavelength [$\\mu$m]\")\n",
    "ax.set_ylabel(\"Flux Density [Jy]\")\n",
    "ax.set_yscale(\"log\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wavelength scaled width boxcar\n",
    "\n",
    "The LRS spatial profile changes as a function of wavelength as JWST is diffraction limited at these wavelengths.  Nominally this means that the FWHM is changing linearly with wavelength.  Scaling the width of the extraction aperture with wavelength accounts for the changing diffraction limit with wavelength to first order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualize the weight images used in the scaled boxcar extraction \n",
    "wimage_scaledboxcar, wimage_sb_bkg = ap_weight_images(ext_center, ext_width, bkg_offset, bkg_width, \n",
    "                                                      data_lrs_reg.shape, waves_boxcar, wavescale=10.0)\n",
    "\n",
    "norm_data = simple_norm(wimage_scaledboxcar)\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.imshow(wimage_scaledboxcar, norm=norm_data, origin=\"lower\")\n",
    "plt.title(\"Scaled boxcar weight image\")\n",
    "\n",
    "norm_data = simple_norm(wimage_sb_bkg)\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.imshow(wimage_sb_bkg, norm=norm_data, origin=\"lower\")\n",
    "plt.title(\"Scaled boxcar backgound weight image\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract the spectrum using the weight image\n",
    "\n",
    "# with background subtraction\n",
    "waves_sboxcar_bkgsub, ext1d_sboxcar_bkgsub, sboxcar_bkgsub_image = extract_1dspec(s2d, ext_center, \n",
    "                                                                                  ext_width, bkg_offset, \n",
    "                                                                                  bkg_width, wavescale=10)\n",
    "\n",
    "# plot\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "gpts = ext1d_boxcar_bkgsub > 0.\n",
    "ax.plot(waves_boxcar_bkgsub[gpts], ext1d_boxcar_bkgsub[gpts], 'k:', label=\"fixed boxcar (bkgsub)\")\n",
    "gpts = ext1d_sboxcar_bkgsub > 0.\n",
    "ax.plot(waves_sboxcar_bkgsub[gpts], ext1d_sboxcar_bkgsub[gpts], 'k-', label=\"scaled boxcar (bkgsub)\")\n",
    "ax.set_title(\"Scaled boxcar 1D extracted spectrum\")\n",
    "ax.set_xlabel(\"wavelength [$\\mu$m]\")\n",
    "ax.set_ylabel(\"Flux Density [Jy]\")\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_ylim(1e-3, 1e-1)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the impact of the scaled boxcar is largest at shorter wavelengths.  This is the result of using the same aperature at 10 microns for both the boxcar and scaled boxcar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## PSF based Extraction\n",
    "\n",
    "While to first order the PSF FHWM changes linearly with wavelength, this is an approximation.  It is better to use the measured spatial profile as a function of wavelength to extract the spectrum.  This tracks the actual variation with wavelength and optimizes the extraction to the higher S/N measurements.  In general, PSF based extractions show the most improvements over boxcar extractions at lower the S/N.\n",
    "\n",
    "There are two PSF based extraction methods.\n",
    "\n",
    "1. PSF weighted: the spatial profile at each wavelength is used to weight the extraction.\n",
    "2. PSF fitting: the spatial profile is fit at each wavelength with the scale parameter versus wavelength giving the spectrum.\n",
    "\n",
    "Only the PSF weighted technique is currently part of this notebook.\n",
    "\n",
    "Note 1: calibration reference file for the specific LRS slit position should be used.\n",
    "\n",
    "Note 2: Small shifts in the centering of the source in the slit should be investigated to see if they impact the PSF based extractions.\n",
    "\n",
    "Limitation: currently it is assumed there are no bad pixels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PSF weighted extaction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generate the PSF profile as a function of wavelength\n",
    "For MIRI LRS slit observations, observations are made at two nod position in the slit after target acquisition.  This means that the location of the sources in the slit is very well known.  Hence, spatial profile (PSF) as a function of wavelength for the two nod positions is straightforward to measure using observations of a bright source.\n",
    "\n",
    "The next few steps generate the needed information for the nod position for which we are extracting spectra based on a simulation of a bright source at the same nod position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lrs spatial profile (PSF) as a function of wavelength\n",
    "# currently, this is just a \"high\" S/N observation of a flat spectrum source at the same slit position\n",
    "psf = datamodels.open(spatialprofilefile)\n",
    "# transpose to make it display better\n",
    "lrspsf = np.transpose(psf.data)\n",
    "norm_data = simple_norm(lrspsf, \"sqrt\")\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.imshow(lrspsf, norm=norm_data, origin=\"lower\")\n",
    "plt.title(\"The LRS Spatial Profile (PSF) Observation\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mock a LRS spectral profile reference file\n",
    "# Sum along the spatial direction and normalize to 1\n",
    "# assume there is no background (none was included in the MIRISim for the flat spectrum source observation)\n",
    "# ignore regions far from the source using a scaled boxcar weight image\n",
    "#   the aperture (psf_width) used in the scaled boxcar weight image could be varied\n",
    "psf_width = 12.0\n",
    "(wimage_scaledboxcar, tmpvar) = ap_weight_images(ext_center, psf_width, bkg_offset, bkg_width, data_lrs_reg.shape, waves_boxcar, wavescale=10.0)\n",
    "\n",
    "psf_weightimage = lrspsf*wimage_scaledboxcar\n",
    "\n",
    "# generate a 2D image of the column sums for division\n",
    "max_psf = np.max(psf_weightimage, axis=0)\n",
    "div_image = np.tile(max_psf, (psf_weightimage.shape[0], 1))\n",
    "div_image[div_image == 0.0] = 1.0  # avoid divide by zero issues\n",
    "\n",
    "# normalize \n",
    "psf_weightimage /= div_image\n",
    "\n",
    "# display\n",
    "norm_data = simple_norm(psf_weightimage, \"sqrt\")\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.imshow(psf_weightimage, norm=norm_data, origin=\"lower\")\n",
    "plt.title(\"The LRS Spatial Profile Reference Image (Normalized)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "y = np.arange(psf_weightimage.shape[0])\n",
    "ax.plot(y, psf_weightimage[:,150], label=\"pixel=150\")\n",
    "ax.plot(y, psf_weightimage[:,225], label=\"pixel=225\")\n",
    "ax.plot(y, psf_weightimage[:,300], label=\"pixel=300\")\n",
    "ax.plot(y, psf_weightimage[:,370], label=\"pixel=370\")\n",
    "ax.set_title(\"Cross-dispersion Cuts\")\n",
    "ax.set_xlim(ext_center-psf_width, ext_center+psf_width)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the spatial profile becomes narrower as the pixel values increases as this corresponds to the wavelength decreasing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extract spectrum using wavelength dependent PSF profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use the normalized PSF weight image to extract the specrum\n",
    "# use the background subtracted image from the scaled boxcar extraction\n",
    "ext1d_psfweight = np.sum(sboxcar_bkgsub_image * psf_weightimage, axis=0)\n",
    "ext1d_psfweight *= 1e6 * s2d.meta.photometry.pixelarea_steradians\n",
    "\n",
    "# plot\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "gpts = ext1d_psfweight > 0.\n",
    "ax.plot(waves_boxcar_bkgsub[gpts], ext1d_psfweight[gpts], 'k-', label=\"psf weighted (bkgsub)\")\n",
    "gpts = ext1d_sboxcar_bkgsub > 0.\n",
    "ax.plot(waves_sboxcar_bkgsub[gpts], ext1d_sboxcar_bkgsub[gpts], 'k:', label=\"scaled boxcar (bkgsub)\")\n",
    "ax.set_title(\"PSF weigthed extracted spectrum\")\n",
    "ax.set_xlabel(\"wavelength [$\\mu$m]\")\n",
    "ax.set_ylabel(\"Flux Density [Jy]\")\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_ylim(1e-3, 1e-1)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the psf weighted extraction has visabily higher S/N, especially at the longer wavelengths where the S/N is lowest overall."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting in Rayleigh-Jeans units\n",
    "\n",
    "For sources that have stellar continuum, it can be useful to plot MIR spectra in Rayleigh-Jeans units.  This just means removing the spectral shape expected for a blackbody with a peak at much shorter wavelengths than the MIR.  This is easily done by multiplying the spectrum by lambda^4 or nu^2.\n",
    "\n",
    "An example of this is given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rayleigh-Jeans plot\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "gpts = ext1d_psfweight > 0.\n",
    "ax.plot(waves_boxcar_bkgsub[gpts], (waves_boxcar_bkgsub[gpts]**4)*ext1d_psfweight[gpts], 'k-', label=\"psf weighted (bkgsub)\")\n",
    "gpts = ext1d_sboxcar_bkgsub > 0.\n",
    "ax.plot(waves_sboxcar_bkgsub[gpts], (waves_sboxcar_bkgsub[gpts]**4)*ext1d_sboxcar_bkgsub[gpts], 'k:', label=\"scaled boxcar (bkgsub)\")\n",
    "gpts = ext1d_boxcar_bkgsub > 0.\n",
    "ax.plot(waves_boxcar_bkgsub[gpts], (waves_boxcar_bkgsub[gpts]**4)*ext1d_boxcar_bkgsub[gpts], 'k--', label=\"fixed boxcar (bkgsub)\")\n",
    "ax.set_title(\"Rayleigh-Jeans plot for all extractions\")\n",
    "ax.set_xlabel(\"wavelength [$\\mu$m]\")\n",
    "ax.set_ylabel(\"Rayleigh-Jeans Flux Density [$\\mu$m$^4$ Jy]\")\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_ylim(10, 100)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional Resources\n",
    "\n",
    "- [MIRI LRS](https://jwst-docs.stsci.edu/mid-infrared-instrument/miri-observing-modes/miri-low-resolution-spectroscopy)\n",
    "- [MIRISim](http://www.stsci.edu/jwst/science-planning/proposal-planning-toolbox/mirisim)\n",
    "- [JWST pipeline](https://jwst-docs.stsci.edu/jwst-data-reduction-pipeline)\n",
    "- PSF weighted extraction [Horne 1986, PASP, 98, 609](https://ui.adsabs.harvard.edu/abs/1986PASP...98..609H/abstract)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## About this notebook\n",
    "\n",
    "**Author:** Karl Gordon, JWST\n",
    "**Updated On:** 2020-07-07"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Page](#top)\n",
    "<img style=\"float: right;\" src=\"https://raw.githubusercontent.com/spacetelescope/notebooks/master/assets/stsci_pri_combo_mark_horizonal_white_bkgd.png\" alt=\"Space Telescope Logo\" width=\"200px\"/> "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}