RXJ1131-XSPEC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OBSID 12833\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from astropy.table import Table\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from subprocess import call\n",
    "\n",
    "\n",
    "home_dir = '###'\n",
    "repro_dir = '12833/repro'\n",
    "os.chdir(home_dir+repro_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After downloading the data from the chandra website or using $\\textit{download_chandra_obsid #OBSID}$, we reduced the data with the following command: $\\textit{chandra_repro}$. For details on what the reprocessing does visit the following website:\n",
    "http://cxc.harvard.edu/ciao/ahelp/chandra_repro.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reduction of Data\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here I will discuss some important commands that need to be run in order to reprocess the data. Afterwards, I will include my reprocessing bash script.."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### acis_clear_status_bits \n",
    "This script clears out several ACIS status bits that are set by the Chandra pipeline and that need to be cleared before using acis_process_events. This is needed to support the bad-pixel/afterglow pipeline. \n",
    "\n",
    "### destreak\n",
    "issue: There is a feature in the serial readout of the ACIS-S4 CCD (CCD_ID=8) such that for a particular frame of data, spurious events (with moderate to small pulse heights) are occasionally reported along a single row of a node.\n",
    "solution: Removing the streaks before creating a new bad pixel file improves the streak detection efficiency and prevents misidentifying pixels with many streak events as hot pixels. \n",
    "\n",
    "### Reprocess badpix\n",
    "There are three steps in creating a new ACIS bad pixel file:\n",
    "\n",
    "    -identify known bad pixels and columns from the calibration database and pixels with bad bias values (acis_build_badpix with a custom bitflag )\n",
    "    -search for afterglows and hot pixels (acis_find_afterglow)\n",
    "    -mark pixels adjacent to afterglows and hot pixels and sort the list of bad pixels and afterglows (a second run of acis_build_badpix)\n",
    "    \n",
    "    \n",
    "### acis_process_events \n",
    "Creates new level 1 event file that has a number of different filters applied. Details on website...\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creation of Pi files\n",
    "We are assuming that you have already done the following  to create a region of interest:\n",
    "    - opened ds9 and picked the region of interest: \"ds9 ...evt2.fits &\" and RofI is called \"simple.reg\"\n",
    "    - Also created background image \"simple_bkgd.reg\"\n",
    "I saved them as physical units\n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specextract\n",
    "After creating the region of interest and background region, I created a bash script which creates all relevant files for creating a spectrum. First I will display my bash script and then I will describe each piece\n",
    "\n",
    "#!/bin/bash\n",
    "\n",
    "#Runs specextract on data \n",
    "\n",
    "punlearn specextract\n",
    "\n",
    "pset specextract infile='acisf12833N002_evt2.fits[sky=region(simple.reg)]'\n",
    "\n",
    "pset specextract outroot=simple\n",
    "\n",
    "pset specextract bkgfile='acisf12833N002_evt2.fits[sky=region(simple_bkgd.reg)]'\n",
    "\n",
    "pset specextract asp=pcadf415060985N002_asol1.fits\n",
    "\n",
    "#pset specextract mskfile=acisf01838_000N003_msk1.fits\n",
    "\n",
    "pset specextract badpixfile=acisf12833_000N002_bpix1.fits\n",
    "\n",
    "pset specextract grouptype=NUM_CTS\n",
    "\n",
    "pset specextract binspec=15\n",
    "\n",
    "pset specextract clobber=yes\n",
    "\n",
    "specextract mode=h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discussion of fits files used\n",
    "#### Primary Event File: evt2\n",
    "This is our primary file. It contains all necessary data with regards to counts per second and energy. It was created from the level 1 event file after filtering the GTI (Good Time Intervals give the time periods when the mission time line parameters fell within acceptable ranges).\n",
    "#### Aspect Solution: asol\n",
    "The aspect solution describes the orientation of the telescope as a function of time. The detected position of an event and the corresponding telescope aspect are combined for an accurate determination of the celestial position of that event.\n",
    "#### Mask: msk\n",
    "The mask file records the valid part of the detector element - ACIS CCD or HRC plate - used for the observation (i.e. the portion for which events can be telemetered).\n",
    "#### Bad Pixels: bpix\n",
    "A list of pixels identified as \"bad\"; criteria for flagging a pixel are listed in the Bad Pixels dictionary entry. Any tool that reads this file will exclude the bad pixels from its calculations.\n",
    "\n",
    "For more information see: http://cxc.harvard.edu/ciao/threads/intro_data/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this $\\textbf{specextract.sh}$ file, we now have our .pi file. Pi stands dor Pulse Invariant and is used to calculate the energy in each bin. We primarily use the PHA - Pulse Height Amplitude - which is another indicator of energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading of Data and Basic Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4VGX2wPHvIUAo0otKjQjYsFBsiFgQRQF7r+xPsa2IWNHVtezaFRaxrb0hWBAL6CKCiqCioDQFQQERURDpHZLz++O9E2+GmWR6PZ/nmSeZW89MJvPe+563iKpijDHGxKpSugMwxhiT3awgMcYYExcrSIwxxsTFChJjjDFxsYLEGGNMXKwgMcYYExcrSExGEJHqIvKeiKwRkTfSHY8xJnJWkKSBiHQRkc+9L82VIjJZRA5M8jkXicgxyTxHnE4HdgYaqOoZwStF5A4R2SYi60Vktff+Hepb31tEfheR+r5lJ4nIryJSx3te+h6IyK4i8q6ILBURFZGioPOd6Z1jo4h8ErSuofc3+9OL5QsROcy3vp2IjBWRFSJSYUctcW4QkfkisklEFovIfSJS6NvmBRHZ6r3+dSIyTUSO8K2vKiIPi8gSb5uFIjLYt97/2vuISLG3XWDb50WkbTkxVhWRN73jqIgcGbT+BhGZ7cW2UERuCFp/gIh85n3ml4jIP33r9haRqSKyynt8JCJ7+9ZfIyILRGSt9/caLCKVw8RZ5MUXcn2m8f4WKiJnpjuWeFhBkmIiUhsYDQwF6gNNgTuBLWmOK93/eC2Beaq6vZxtXlPVnYCGwMdA6Z2Lqr4HTAAGA4hIXeAJ4ApVXRPiWCXA/4DTwpxrJfAf4L4Q69YD/wc0AuoB9wPv+d7DbcDrwMXlvBa/R4BLgQuBWsDxwNHeMfwe8F5/Hdxre0tECrx1NwOdgIO8YxwFfFvOOb/wHesYYBMwTUTalbPPJOB84PcQ68SLvx7QA7hKRM72rX8VmIj7zB8BXCEiJ3rrluIuJOrj/rbvAiN8+74HdFDV2kA7YH/g6nLizCYX4T5rF6U7kLioqj1S+MD9s68uZ30fYDKuoFkDzAW6+dbXAZ4FfgN+Bf4NFPjW9wXmAOuA74EOwMu4L85NuC/BG4EiQHFfdotx/+RHAkuC4lkEHOP9fgfuy/sV7/izgLa4L7HlwC/AseW8tr2AT4DVwHfAid7yO4GtuC/g9cDFIfa9A3jF93xvL/5GvmUNvTiOA54Hhod7Lb5llb3jFIWJ+RLgk3JeUyWgt3eMxkHrWrt/sXI/D22AYuCgoOXNcRcXR3vPXwD+7VtfwztnE+/5aOCacs7j/zv2ASaF2GY08GYEn+ElwJEVbPMIMNT3fCOwt+/5G8DNIfarDPwd2BjmuA2Aj4DHw6wPfK4r+963x4EPvM/WZGAX3EXCKtz/V3vf/gOBn/jr/+cU37oC4GFgBbAQuCroXOX+b4aItSXu//I0YDuwc0XvfaY+7I4k9eYBxSLyoogcLyL1QmxzMLAA98V4O+7KM1Bl8yLuQ9caaA8ci/uyQ0TOwH3hXgjUBk4E/lTVC3CFRW9V3UlVH/Cd6wjcF/xxEcbfG1cw1cNd8Y7FfZk2Be4C/htqJxGpgruy/BBoDPQDhonIHqp6O3AP3h2Hqj5bXgAiUtV7jX/ivgwAUNUVQH9gGNCLJF+1ishMYDPuCvoZVV0ew2G64Qrvr/wLVfUX4Euge4jzFuBe/0Jgmbf4S+BaEblSRPYVEYkhlreAw2PYLzg+8Y7znW/xf4ALRaSKiOwBHIorEPz7rca9n0Nxnwf/unNFZC3uS3x/wnzOwjgTuBX3/7QF+AL4xnv+JjDIt+1PXux1cBc4r4jIrt66vri7xQNwF2gnB50n7P9mGBcCU1V1JO7i77woXlNmSXdJlo8P3Bf3C7gru+24L6KdvXV9cLf64tv+K+ACXA5hC1Ddt+4c4GPv97FA/zDnXITvapy/rtxa+ZYdScV3JON863rjrvIKvOe1vGPWDXH+w3FVIpV8y4YDd/iO/Uqo2H3rt+LuZopxhciRIbbbDXdnM6yi98BbFu8dSTXvb3BRiHWR3JHcCnwZZt0I4Gnv9xdwX7KBL9vNwHm+bQtwV/KTvc/IUn9MRHZH0gPYFsHnt9w7EtwX8Ayg0LesM/Aj7vOuwJ1h9q0JXAn0DLO+DfAvYJcw6wOfa/8dydO+9f2AOb7n+1J+DcF04CTv9wnAZb51xwTORQX/m2GOPR/vLhJ3Vz+jovc+Ux92R5IGqjpHVfuoajNcnW8T3BVbwK/qfbo8P3vbtASqAL95Sd7VuCuzxt52zXFXVNH4Jcrtl/l+3wSsUNVi33OAnULs1wT4RVVLfMt+xt3JROp1Va2L+6edDXQMsc1TwEvACSLSOYpjx0RVN6vqcGCgiOwfwyFWALuGWbertz7gIe/1V8dVkT4oIsd7cRSr6mOqehhQF7gbeE5E9ooilqa4+vqYichVuCvtnqq6xVtWH5ePugtX8DYHjhORK4P3V9UNwJPASyLSOMT6+bg7ncejCCv4Mxv8vPTzKiIXish03/9XO9ydC3ifYd++/t8r+t8sw2ucsRt/5YJeBfYVkQOieF0ZwwqSNFPVubirJn+Ss2lQ1UQL3BXmL7irnoaqWtd71FbVfbztfgF2D3eqCJZvwNW9A6VVKI0ifS0VWAo0FxH/Z64Fri45KuqqsC4D7vBVOyAiF+O+pK4EbgGe9qrBUqEK0CqG/Sbg3peD/AtFpDlwCDA+eAd1ZuPuPnqGWL9JVR/DVfvtHby+HKcAn0WxfRki8n+4HEM3VV3iW9UKKFbVl1R1u7duBHBCmENVwn0Ow11kVCb85zxmItISeBqX+2jgFdqzcQ0JwOU+mvl2ae77vaL/zWAXecedLiK/A1O85Rcm5tWklhUkKSYie4rIdSLSzHveHHcL/KVvs8bA1V598hm4qrD3VfU3XI7hYRGpLSKVRGR3XzPQZ4DrRaSjOK29fw5wV2EVfdHNA6qJSE8vp3ErUFjBPpGagiuobvRe15G4qrER5e4VhlcAj8U1HEBEmgAPAn29K+EncdVf/wh3DBGpxl+vr9B7HlhX4D2vDFQSkWree4KIHCKuCXdVcf1fbsLdJU3x1ou3b9XAecTXlDfodczzYh3mHbdARPYBRgIfqepHofYTkT2BLnh5CK+J7JFePJVF5CJcVWN5LbcCr3M3ERmKq9q8s5xt/e9RVe91ibfuPFxeo7uqLgjadZ73tpzrfWZ3Ac7CVX8hIt1FpL0XS21czmIVLm+AiFwSuDsR1yz4ZkIUsAlQE3dh9Yd3rr9R9gLvdaC/iDQV1yrwpsCKCP43S3nv4Zm4lnoH+B79gPMk/S0oo5fuurV8e+Cusl7HXYlv8H7+F6jtre+Du9J8FNdqax6+llD81fRzibf+W+Bs3/rLgR9wuYvZeC1SgJNwCffVwPUE1SX79u+Du/Ja7m23iLI5En/LqWOARb7ngXxDszCvfR/gUy/u4BYxZY4dYt8d1uMaJWzAFbxvE9SSB9jDO9c+3vNFlG0Bp8GPoPcheP0L3rojcF+C63BVQZ8CXX37FoXYd1E5r60S7kvpR1xVyy/AA0A13zYv4HJE673XvBj3xV3JW38ZMM17vatxebVevv39f8c+uDxT4Fg/4xLFe1Xw2V0U4nUVeesW8leru8DjSd++RwNfe/H9jrvyr+GtOwPXemo97kv8fWA/377P4y6ENngxPOh/b4JiDLz3/hyJv7VbmZwXLo+13ff8bu9vugJXoH0KXOL7fA/GXaAsBAZ4r1ki+d/0neNs3P9YlaDl1bzz9gr12jL5EXgDTIYQkT64D26XdMeSa0RkJa457fR0x2Kyn5efelJVW1a4cY6zqi2TF0SkO65l0/x0x2Kyk1dteIJXddgU1zR/VLrjygRWkJicJyIjcK25+qprFWRMLASXQ1qFq7aaA/yz3D3yhFVtGWOMiYvdkRhjjIlL9jUzi0HDhg21qKgo3WEYY0xWmTZt2gpVrbAvWV4UJEVFRUydOjXdYRhjTFYRkZ8j2c6qtowxxsTFChJjjDFxsYLEGGNMXKwgMcYYExcrSIwxxsTFCpIcMHjcvIRsY4wxsbCCJAcMGV/x8FGRbGOMMbGwgiTHRXonYncsxphY5cVYW506ddJc65A4eNy8kHcZ/bu1YUD3thFvE1A0cAyL7tthsj1jTB4TkWmq2qmi7fKiZ3suGtC9bWlhEK4QGNC9LUPGz2fRfT2toDDGJI0VJDko+E6kaOCY0uX+O5Fw24W6YzHGmHByumpLRHoDvVu3bt13/vzcTTYHFxB+gTuR8rbxb2eMMQGRVm3ldLJdVd9T1Uvr1KmT7lCSKpK7h0i2sYS7MSYWOV2QGFdNFel21kTYGBMLK0jSLNl3AZHmOoK3iyYuu5MxJr9ZQZJm4e4C4v1yjqb/SNHAMaWJ9sDv0dyd2J2MMfnNCpIMFe+Xc6T7D+jelkX39SxNtPt/N8aYSFjz3xSIttltoquKKmqx5ReIJVRcwce0psPGGMjx5r8B6e7ZXl7TWv+6aHqihxJufyCiuwx/gRNNc2BrOmxMbrLmvxmqvLuNQDVTQLhqJv8x/L9XVE0V6tzB+yeSJeGNyQ9WtZUk4ap+oOwXdqB57ln//YIpC1eWOUbRwDGlzXJL99myjq8/HgWNmvPxjJ9Y+MNmaNUNqtfnmWlruKR7eyistcM5w52/zLF9/M2GK6oaC9fEONyxjTG5xQqSJAk3Fpb/Cz2wHcCUhSvLbOOvzhoyfj6sXgyf3A+z3uDVqlvgbTgKOKoq8PKjAFwCMA3YaWfGtOjEPkedA627U3Tb+NLxtmJpDlxRgWCFhTH5zQqSFIk0iR1qnxMrTWbT4D4AjCw+nLElB/Kz7sx6rU49WUd91lFP1lFX1tOrdXW61vqVfX4cD6+NYXPBTtxfuRPn3zKLAvbeoSCrKJ7garRw8QbWWRLemPxjyfYUCJfEDpccP3i3+rx22aHuyfTh8PblfFmyF4dc9yZF986I6JzXHL0b17T6DWa9yZbZ71BYvIE/tRYNOp0GB5wHzQ4EkbCJ8mgT/6GOY0l4Y7KbJdvTLJIkdrjkeKAQeW3kG/De1Uwu3ocLtw6EOs1K9y0vKb/ovp5cc+ze0LobnPIEhQN/grOG8VnJvjDrTXi2O/y3K+NefoDqbC43tuDz2V2FMSaY3ZEkSXlX+hFd0a/5lRWDDmGt1uCUrXexhp1CnieQ6I7k7mHwuHkM6NoEZr0OXz8Ly2azuWAnqh1+NRx6ZWmSvrxmxNGcN5r+K8aYzBPpHQmqmvOPjh07aqq1vGl0VNsP+vCHv54Ub1d9vqeu/2cj1T/mlTle4OegD38ou4+3LJLzDvrwB9WSEj114MOqw89Vvb226v2tVL95WbWkZIfXETiP/9jB547m9Qbva4zJTMBUjeA71qq2EijcuFWR9KfwT497/21XwqLPuGP7RRQ9NC9kgjxUS6poOi0W3fw+03QPiqafy0lb7uLXgqbwzt/h+RNg+ZwKjxvPEC42NpcxucUKkgQKl/PwFxIV2Xndd9xU+CbsfTJvFB9B/25tWHRfT/p3a8PgcfMqHBY+sF1F/DEeefTxNB3wCZz4KPwxB57sAuNuh60by5yvvHOHW2edEo3JfZYjSZKYWjFtWs2Cew+mVd0CuGIyRXd+DhC2D0pAcC4k1HnKa4UVmNcdgA1/wrjbYPowqFcEvYdAqyMrPEYkMzTGM/yLMSb1LEeS5hxJqDxAcB6hzDab16o+31O3/rOe6sJJpev9+4T7vaLzBK/z51nCxrXgU9UhB7j8yagrVTf8GdE5IsmdRJs/MsakBxHmSOyOJMkiuRN48Z2xHDbtGnaT37h+2+WMKjk8omMHj6MV7jwQeV6izB3Ktk3w6QMweQjUaAAnPQZtjwXC3135h3UJdewB3dta/xJjsoTdkaT5jkS14qvzljeNVp09SvXuJqoP7K6vvz4s9DYhfi+v5VNFdyTh4gq0BNth/6UzVB871N2djL5OdcuGsOcvL/7APtZqy5jsQIR3JDZEShKFalnlv3PoVmkaJa8P4lttzbT9/sM9k9bssA2Ezo3Ek1cItW/5Q5vsB30nwPi74MvHYOFEBpz6VET7+gXeD8uJGJNjIiltsv2RrjuScnMia37VNf/cWfXJrqpb1pduX95dTKRX8uVtF+quwJ8rCZVDKePHCaoP7aF6ZwOd+Nwtrs9LmHiDj2G5EWOyC5Yj+UsqcyQRt0567QI2ff8/hnd6jbsmbwq7fbLyCRW1pgoIee6NK2H0NfD9O9D8EDjxEWi0R5njBlhrLWOyl+VIMvCOpNTSGaq319bJz1xf4fbJyieEa00VaU92LSnRATffqHpvC9W7GqpOfFi1uDjm3I2f5VCMyQxYjiSDTRrEWq1O53P/UeGmibxqj3SI9yHj55duFzwhVtljdGXi6v25s8rz9Bx/JyyaVCZ3EiubEMuY7GIFSRKF7O29bhl8/y7Di3twWfW6FW+fQOEm2wo+f3kTcgUfY+p954KeA9NegA9ugqePgnPfgMZ77nD+ZL8+Y0x6WI4kgSoa7XbwuHls+mQQt1QZTrctD/KTNgWC5h9J0rmDRZJ7qagnffAxXn1rFOf+dCObN2+k2oVvQMvOEccbaqph/7mMMaln85GkQUWd/gZ0b8stTaZDs4P4SZuWjsUV6gs00ecOFsndQWCcr3BjhwUf45avqkLf8fy6rRYMOxOWfhtxvIH3INy5jDGZywqSJAg7UOHyuW5QxH3PSG1AIUTyBR2uv0m5x6jbgvO23gLV68Erp8GKsgVGRYM4hltvgz8ak7msaitOUTVv/eR++OReuHYOZw1fFHdVTiqb1gaqoiJt3ttSfuetqrejNRvTcMBkqFI9ogmzAiIZiNIYk1zW/DcNzX+DJ5/awWOHqj57XNj9EnHuZItksMjSbeaPc8OqvNMv/Da+5f7OkNGc1xiTHFjz39QKVL34J7UC35X1ivmw/DvocX9KY0rEnUmkzYZ3PPcx0OVamDSIy7+sx/9KDiqz/1n//YLXLjt0h2qrwPqDd6tf5q4tkvMaY9IgktIm2x+puiMJNyy7qqp++oC7Ol+9JOx+iTi3XzKu4ss7pn/4ldLttm9TffJw1QfbqG5cVWb4leAOj+VNFWx3JMakHrk81a6ItBKRZ0XkzXTH4hf2KlkVZr4BzQ+GOk0j3y8R504C/x1EqOR7mVgKKkPvR2DDH/DRHaX7hJomONQAl6GS7JZ4NyazpLwgEZHnRGS5iMwOWt5DRH4QkR9FZGB5x1DVBap6cXIjjd0OTWt/mQIrfoD25yf93PHMGx+J4LlG/L+HPfd3NeCQK2Ha81zYdKmbMz5MfP73zt/DPni5MSaDRHLbksgH0BXoAMz2LSsAfgJaAVWBGcDewL7A6KBHY99+b0ZyzpSPtfXTx6qTH1Fds1S1pER1+LmqdzdV3bwupWEkqzookpkad1i+eZ3qoHaqQw/UNjeNiii+cMl3q+YyJjXI1GS7qk4UkaKgxQcBP6rqAgARGQGcpKr3Ar1SG2GcvngMxt7ifp9wN7TuBnNHw9G3QuFO6Y0tDuXNkRJR8n3iUmb8cTYvVH2AmyqPoGhgldLj+vcJ1UQ4VA97S7wbk0EiKW0S/QCKKHtHcjrwjO/5BcCj5ezfAHgSdxdzc5htLgWmAlNbtGiRuCK6PKsWq/5rZ9VXTlddNsfdidy/m+qoK1SLi1MTg08qRg6OdO72UmOud40Opg+vcKTgcE2C7Y7EmNQgU+9IwpAQy8L2lFTVP4HLyzugqj4FPAWuQ2Jc0UXqy8dBi+GEh6BeSzh7WEpOG046r9TDnvvYfzP5i8kc9s7fGdBrMGB3E8Zku0xptbUEaO573gxYmqZYojZ43DwoKYbZI6HNsa4Q8a8jt1oa+RPfUY/oW7mQbzs/Brt1hXf7wWvnw1r3pw5O1vv5k/L9u7XJqffTmKwXyW1Loh/sWLVVGVgA7MZfyfZ9EnW+ZCfbW940WnXBp67KZtbIHdepVcfsYPs21c8Gqf6rseq/d1UdfZ2rDvSEqj6LpErNGJM4ZGrVlogMB44EGorIEuB2VX1WRK4CxuJacD2nqt+lOra4/PA/KCiEtj3SHUl2KKgMXQbAXifCxIfgmxfh66eh/u7Q9jiOrVQIS5tCjQY0ZA0Ub093xMaYMHJ60EYR6Q30bt26dd/58xPb9yC4ddGYqjezRmsypNkgDmnVoNy+DvnW0iiioVrW/wHfvw3zxsLCiVC8pczqtVqDV4u78dD2M9gedP3jfz8TNSyMMSbyQRtzuiAJSPbEVgcMHMH0apfDUbdQ9ME+ZUapDYxam8+j10b92rducEPur1sKm1bB9i2w+AuXg9qzF5z1CkU3vw+ww3Hz+X02JtEiLUiiqtoSkULgNFyOo3RfVb0r2gBzycGV5gIKRYcD8U9Slfeq1oRmHYGOfy07qC93f1vIP+a+ClOfA3ZJV3TGmCDR5kjeAdYA04AtFWybN46v9ytb1xfQ7ollQJUyneUCrZrybb7yWEYMrkiNI65h8XfzqD76DqoxhM0UWmdFYzJBJBn5wANfS6tseiR9iJQXT1J9oouqRtaaKFkdBTNVJCMGR2rEm6+71nGfP2rDpxiTZCRp9N/PRWTfxBdnWUwVfpsBTQ6IeBcbdPAv0b4XN31dA1p0hq+fQShJUlTGmGhEW5B0AaZ5o/TOFJFZIjIzGYFljTVLYNNK2HV/IP+qsCKR8Pek099g5QIe7Lgm+ecyxlQoqlZbItIy1HJV/TlhESVQMpv/lpozGl47Dy4ZD83CN25I5fzqmcjfLDfa9yJ4+0K28kXhVaxoeAht+42M6JzGmOglrfmviNQD2gDVAstUdWLUEaZQUpv/TnwIJvwLbv414tF987GJarjXHEiMR/p+lB5n7D9gypNw7RzYqXFU5zTGRCbSgiSqqi0RuQSYiOuBfqf3845YAswZK+ZD7aZZPUR8VupwEZRsh+npHRjTGBN989/+wIHAl6p6lIjsiStQ8tef86FB66h2yZd6/PKaAAMxNQ8ufe8atYWWXWDaC9C5P1SqVOE5rZrLmCSJpGlX4AF87f2cDhQGfo/mGOl4JKv576Cxc3XTXU1UR1+bd016o1Ves9yYm+zO8JoC/zgh5Pvf8qbR9ncxJg4kqfnvEhGpC7wNjBORd8ii4d4TbdiEaVQrXg8N2liT3nTY+0SoXh+mPR/2/be/izHJF1VBoqqnqOpqVb0DuA14Fjg5GYFlg93FK0MbRle1lY/Kq86LuaqvciEccC7MHUMjVifuuMaYqOT0oI3JaP7rr4M/p2A891Z5lsM2D+FXGpVuY/XxqTF43Dzem/ApEwqv54FtZ/F48UkAHLxbfaYs3HHMM/u7GBOdpDT/zdZBG+Nt/huuP8LTt57D+QUfUf32ZRTd8oE1NU2XF3qxZOEcmt36HVSpVmZV0cAxVoAYE6OkNP/FDdp4ErAd2OB75LRw9eyt5DcW6S6lLYZMmnS9gWayAj57OORqy5MYk1zRNv9tpqo2BaCnY80/WF5zT8Dq49Oq1RHMaXQ8e00aDO1OhcZ7la7q380aQhiTbNFWbT0FDFXVWckLKfFiqdqqcBiP7Vvg7l3g8Ovg6FsTFaqJ1fo/4PFDoEYD6DuBwROX5vWQNMYkQkJzJCIyC1DcHUwbYAFuPhIBVFX3iy/c5Io3RxJyqI3lc+Hxg+HUp2G/M+OM0CTEgk/gpZOhYx/o/Z/SxTZUijGxSfQMib3ijCf3/Old7UbZq90kUasj4bCrYfIQ2LMntOme7oiMyQsRZYlV9Wd1I/zeBazxPV8L3J7MADNByPzHinnupxUkmeWofzC3pDmMHgDbNgOWvzIm2aJtbrSfqpb2/FLVVUD7xIaUeULWqa/4EXbaBarVTn1AJrzKhdy1/QJY8wtMfRYI8/czxiRMtK22KolIPa8AQUTqx3CMlPF1SEz8wf/8ERralW6mKNs4oh0Ti/el3f/uYfjqQ/n78R3SGpsxuS7aO5KHgS9E5F8ichfwOfBA4sNKDFV9T1UvrVOnTuIPvnox1A05z5dJgwHd27Lovp6lSfWulw+hvqzn73W/THNkxuS+iAoSETlURERVXwJOBZYBfwCnqurLyQwwI23fAut/h7rN0x2JCadpR2h2EHz9DJTY3O7GJFOkdyQX4eZqHwEcBLypqkNV9fvkhZbB1ixxP+u2SG8cJqTS5PpBfWHlT/DzpPQGZEyOiyi/oaqXA3gTWR0PvCAidYCPgf8Bk1W1OGlRZprVi93POnZHkolKk+t79oQqNWH2SNita3qDMiaHRTuM/FxVHewNk3I0MAk4A5iSjOAy1ppf3E+r2spsVWvCHsfD9+9A8bZ0R2NMzop5tEFV3aSq76tqv0h6PuaU1b+AVHJztZvMts8psGkV/Px5uiMxJmdFVZCIyIveDImB5/VE5LnEh5Xh1vwCtZpAQZV0R2Iq0uoIqFQFfhrP4HHzGDxuXumqcL8bY6JjHRJjsXqxVWtli8Ja0OIQ+HE8Q8bPLzOQY7jfjTHRibYgqSQi9QJPMr1DYtKs/sUS7dmkdTdYNpvGrEp3JMbkpJzukCgivUXkqTVr1iTuoMXbYe2v1vQ3Q1RUJTV43DyOH1MIQNeCmYAbDbho4Jiwv8dSzWVVYyafRdtqK6s6JCalZ/u630CLrWorQ1RUJTWge1s+uOcK2GlnulRy0+j4e8CH+j2Wsbmsaszks2iT7QJ0AOqr6lBgvYgclJTIMpX1Ick+ItDyMA6uNBc3rY4xJpGizW88DpTg+pDcBawDRgIHJjiuzFXah8TG2UqX4NkrA1VT/tkPB4+bV/bOougwdv3uLf7ZuUbpIv/w8rEMNR9JHMbkg2in2v1GVTuIyLeq2t5bNkNV909ahAkQ7wyJZXz6IHz8b/jHMqhSLTHHNDELN/vhDssDM1qe9Bi0Pz9lcRiTzSKdITHaZPs2ESnAqx8QkUa4O5T8sWYx1GxshUi2abSHm8990eR0R2JMzom2IHkEGAU0FpG7cUOk3JPwqDLZmiVQp1m6ozAef5XU4HHzwrfCEoGWnaMewDHS1ljxzMJoLb5MtosqR6Kqw0RkGtDNW3QB4DrkAAAdNUlEQVSSqs5NfFgZbO1v0GD3dEdhPP5cxIDubUufh6xqatkF5rzn+gFF2OpuyPj5EeU74smJRHoOYzJVtK22zgB+VdXHgPrAPSKSX9PPrVsKtZukOwoTi5ad3U8bd8uYhIq21dZtqvqGiHQBuuM6KD4BHJzwyDLR1o2weQ3U2iXdkZgKhKxq2nkfqFbHVW/tf1bYfctrjQXR3X3s0HosgnPY3YnJNtEWJIE5R3oCT6rqOyJyR2JDymDrfnM/a9kdSaYL+WVcqQBadK4w4V5eFVnRwDFRfdGHq7aqsBrOmCwSbbL9VxH5L3Am8L6IFMZwjOwVKEhq75reOEzsig5zsyau+z3dkRiTM6K9IzkT6AE8pKqrRWRX4IbEh5UYItIb6N26devEHHBt4I7ECpKsVZonmQztTqtw8/7d2kRdDRXt9vG0+DImE0TUIVFEDgW+1Gh6L2aQhHVInPwIjLsNBv4C1WrHfzyTesXb4f4i2O9M6DUo6t2jrYayaiuTzRLdIfEiYJqIjBCRPiKSn9nmdb9B1Z2sEMlmBZWhxcHujsQYkxARFSSqermqdgDuAOoBL4jIFyJyj4h09Xq75761S63FVi5o2Rn+mAsbVkS9a7TVUFZtZfJBtMPIz1XVwaraAzdw4yTgDGBKMoLLOOt+t/xILmjZxf30+pNE07M8OMcRvG/wc2vKa/JBhQWJiHQXkadF5ADv+aUAqrpJVd9X1X6R1KHlBOuMmBuatGeTVi2t3opnLpHgfW1eEpOPIrkjuRLXMut8ETkaOCC5IWUoVe+OxKq2sl7lqkwraWMDOBqTIJE0//1DVVcD14vIfeTT3CN+G1dC8VbrjJjF/M1yLy/Yly7LRtDl5heARlH1LA/XvDf4ufVSN/kikoKk9L9EVQeKSL8kxpO51i11P60zYtby9ybvcvNyBjKCST1XUzS6UVRNdCvq+W7NfU2+qbBqS1XfCXo+NHnhZDDrjJhTlmhjaNoJZo9MdygVsmHmTaaLdvTfTiIySkS+EZGZIjJLRGYmK7iMss4KklzSv1sb17P995nccWiV+I5TzvNEsAS+yXTRDpEyDJd4n0W+zYxYWpBYsj0XDOjeFtbWgnG30afKR8CxsR+nnOfG5INoC5I/VPXdpESS6dYuhZqNoCD2q1eTYWrvCvufDdNegMOvzaiLBBtm3mSTaEfuvV1EnhGRc0Tk1MAjKZFlGuuMmJsOvw5KtsOnD7gm3mFU1PEw0QZ0b0v/bm1KE/eL7uvJovt6xlSI5HuOJd9ffypEW5D8DdePpAfQ23v0SnRQGck6I+am+q2g099g6rPwbj/YvuWvdSUl7k500WQWfvyiu3OZ+jx8/y4TJoyFbZuTGlqiciP5nmPJ99efCtFWbe2vqvsmJZJMt/Y318rH5J7jH4Tq9WDigzDzdajZECpXg7W/wnZXWDxSFXjvr13eKwTuuR123Q9aHwN79oJd9weRhIdn43WZTBfRMPKlG4s8DQxW1e+TF1Li+OYj6Tt/fhxXJdu3wL8b83nzS+l88YMJi89kmB/Hw4JPXOfTbRuYuqomo36uymJtzFJtwHqtjqDUl3W0kOW0q7SQQyrNoUOlH6lECTTeGw6+HNqf72ZjjEFwbiSgvNxIqOl8YzlOLsn3158okQ4jj6pG/ADmAFuBH4CZuNZbM6M5RjoeHTt21LisWqx6e2298Zbr4juOyVotbxod/vmGP1W/ekb1ya6qt9dWfewQ1UWfJ/ycsW4X6XFyVb6//ngAUzWC79hocyQ9gDa4tpKB/EjvKI+RfTb+CcBKrZXmQExGqlEfDrwYLv0EzngRtq6n5PmeMOk/Ls9C5iV8A/FkWlwmO0U7jPzPoR7JCi4TDB43jwse/QCAVboTRQPHUDRwjP0D5pmIOh6KwD4nw+WTeb+4E3x0O7x+AWxZH1PCt7zcyOBx80o/i0C5n8tQxwnEkw+JaMsxJV9UyXYReRHor24QR0SkHvCwqv5fMoLLBAO6t4Vd2sBIWEUtG0cpT0XV8bBaba7adjW9TlgMH/4DnutBE/rGfc7gdeHG+4rmOPkg319/KkTbamu/QCECoKqrRKR9gmPKPFa1ZSJQNsErFL3TkiMq3cDQ3x/hncJbOfXmlXyjbdOW8A3XydE6O5p4RVuQVBKReqq6CkBE6sdwjOyzcSUgXHR0fk7FYiIT+i6hJ/xxGguH9uatGvdAj/ugY4+EnjfSqptQ8dloxSYRok22Pwx8ISL/EpG7gM+BBxIfVobZ+CdUr8s1x+4FWIIyHyT0b9xoD07e+i9ocSiMuRae6AxzxzD4wx9CnjPac/vvIpLx2czWz3ukcWfr68skERUkInKoiIiqvgScCiwD/gBOVdWXkxlgRti0EqrXL32aDwnKfBfv3zj4LqFPt/ZwwdtwxgtuSJYR53LSpBNhzPUwZzRsWJGQBHik+wbii+RuJls/70PGz4epz8Hwc+H9G2DKU1BSHHo7E5dIq6UuAh4TkXnA/4A3VfX35IWVYTb+CTUapDsKk0XCJuf3OQX27A0zR7DoradpNX0YfP00AB9X3RlGvcdZBTvB2g5JnUQtEE8u50TqsRbG/gOq1IBFW2HLWqjTFPa0qrxEi7Zn+57A8cBxQB3gY1zBMllVdyzqM0SnTp106tSpsR/gyS4s2FKHo3+7codVlqDMHanoDR18jqpsY3/5iQ6V5tOx0jzaV5pPI1kLwIySVmwo6k7nXv8HjfdMadzZ2jPcH/c1ld/kmspv0W3Lg5x4ZBf6zzwRWhwCZ72Sta8v1ZLSs93/AKoDJwBDibD3Y7oecfdsf3hv1VFXlD61nrK5a9CHP6hqav7GIXvLl5Ro94FPqE58SPXpbqq313G95Z/prvrtq6rbtkR1zGTEmYkCf7dSW9bryn82UX317NJFU/97ueqdDdxIBD7Z8PrSJdLv9mhnSDxDRAJtYK8DLgGe10hKrGy28U83qJ/JeWmvLxdhnjZ3w9tf8hFc9wMc+2/3GXz7chjaAb5+tuwoxWbHv9s3L1NP1sNh15Quum1hOyjZlhXTK2ebaFtt3aaq60SkC65660XgicSHlUG2boTtm8rkSKynbO5Lxd84XG/5Mstr7Qyd+8FVU+G8kW7yrTHXwiPt4aundxjKPhlxZ+XnfeYIft9pL2hxcOmi77UIdm4HM0aU2TQrX1+GiTZH8q2qtheRe4FZqvpqYFnyQoxfXDmSNUtg8D7Qewh07JPQuExmyKr6clX4aQJ8ej/8MgVqNYEuA6DjRVC5MN3RpVS4v9vNh9fnsq97wFG3MnjrSWW2uaRgDLdWGcYL7d+gz0mxTa+cT5I1+u9o4L/AT0BdoBCYEc0x0vGIK0eydIaro/7+3diPYbJGJteX+/MAg8bOVf1xguqzx7nP5+B2qtOHqxZvT8jxs02Zv9s3r7j3ZOn0HbdZ+5vqHXVVP7ozxRFmJ5I0+u+ZwFigh7qhUuoDN0R5jOziDY/i70diTDr4r6yHTPgRdj8K/vYBnP8WVKsLoy6Dxw91/SU2r4nr+Flt/lg3LfYu++24rtYusHs3mPFa6cjMJn5Rz9muqm+p6nwAVf0N6Jb4sDLIppXup/UjyQtZV18uAq27waWfwunPQ9Wa8MEN8PBe8O7V8NuMdEeYEqV/t+Jt8NPH0Kb7DrNVlm6zzymwdgn8MSfFUeauaHMk36hqh6BlM1U1RNGfOeLKkXz1NLx/PVz/I+zUKLGBGVOBcHkAvx1yOb9+4+agnzXSNRRp2snN2tjuVKhWJ6LjZ2R+KBILP4MXe8FZr8BeYaZKWrnANVboNRg65ezA5QmR0BwJcAVuNsQNuJkRA4+FwCuRHCOdj7hyJB/f6+pbt2+L/RjGeOLJQ/jzAJHkch57/2vVLx5XHXqg+wz/q7HOGXqG6vxxIXMpkRyzTJ4mzO/JEPHxx/7D9RXZvDb8NiUluv7fRapvXZaamLIYCc6RvIqbCfFd72fg0VFVz4+ujMsyG1e6q7iC3B/k2CRfKvMQD3y6DA65Av4+BfpOgPbns8sfk+CV0+A/+8HkR6LOpZTJ04T5PRkiPv68D6HoMCgsZ8oHESZtbgWLv0xNTHkgooJEVdeo6iJVPUfLzo64MtkBpt3GPy3RbjKCP38TVS5HBJp2hJ4Pc9CWx910wA1awbjbYNA+bjyq1b9kX34o2KpFsOIHaHNchZtOLWkLqxbC+uXJjysPRJsjKQROA4rwDfioqnclPLIEiitH8vIpsHkt9B2f2KBM3khlHiKSnErA3QcXc17JuzD7Lbdg/3Og6/VQf7eYjwmJe11Rv2+BfOZV06Bh63KP10Hm8VbhHVy2dQB7HnVuxPHmXE6pAsnqR/I/4DXgRtwQKdcB10VzjFQ+cNVvT7Vu3Tr2SsInD1d95fTY9zfGJ5X9VEKdK+T5Vy1Wff9G1bsaqd5ZX/Wdfqp/LqjwmNHmbOIR0fFfPk11yAERHa/NTaNU72rocirJjCnLkaR+JM1U9SxVfUBVHw48ojxGyqjqe6p6aZ06dSreOJyNq6zpr8ltdZvD8fdD/+muFdOM4a5V0yunww8fhJzDI1ZJm0Rq60ZY9Bm0OS6ic2ylCjRpD798lZx48ky0BcnnIrJvUiLJVJYjMQmUyjxEqHOVe/7aTeCEB6H/DDjiRvh9Fgw/G/6zr8ujLJlG/6NbhzxWpK8r1gR1hcdfOBG2b4a2x0Z0jv7d2kDzg2DptzuMV5awmPJItDmS74E2wAJgCyCAaq72I9m2Ge7eGY6+zdUdG5NPire5O5Lpw+DH8W7k3LotXIe+fU6FXfffodNfRZI2R/zoAa63+k0LKbr1o8jOMWc0vHYe/N+HZQZ3NH+JNEcSbUHSMtRyVf05ithSLq6C5JuXoFknaNqh4u2NyVWbVsHc9+G7UbDgYzddcP3dYf+zYb+zoF7IrwYgBQlqVdbetydfbmzKpduui/wc6/+Ah1rDMXdCl2tCb5PnklWQ/DPUcs3lVlvGmLI2roQ578GsN1xeAqBlF1eo7H0SVKsddtek3JEs+x6eOLR0hO6ozvHogVC3JZz/ZmJjyhGRFiTR5kg2+B7FuGl3i6KOzhiTvWrUd8PW9xkN18yCo2+F9b/Du1fBQ23gzYth/kdQvD018cwf6362iWFY+N26wuIvXDWeiVlUBYm/pZaq3g0cCTRNSmTGmMxXtwV0vcFNvHXJeGh/PpvnfgjDTnPz+Hx4Kyz7rnTzcAnquFpzzfsQdt7XNRYo5xwhFR0OW9e7pHui48oj0d6RBKsBtEpEIMaYLCbicok9H2a/9UPhzJddb/ovn4AnOsOTXeCLxxjQOXRT+piHG1k0GRZ/XmaAxqjyLkWHu58LJyY2rjwT7Zzts0Rkpvf4DvgBGJKc0Iwx2WgrVWDvE+GcV+G6eXD8g1CpCoy9BQbtCSP7unGuosjPhrRtM7x3tbsr6nxVbMeo2cBNvxvI9ZiYxNNqazuwTFVTVBEaO0u2G5NcEbXMWj4Hpj7n5kzfspYVNVrx9JqDebv4MJZRP/Q+5Rl/F3z2sJvYq3Uc0yJ9MBCmPQ8DF0PlwrwbBqU8SWm15R14f8C7H2Siqs6MIb6UsoLEmNSpsNXUlvUw+02Y/qqbdx5hYnE7up7eD/bq5Sbnqsjvs+CpI2HfM+GUJ+ILeO77MOIc6PO+GznYJ2n9XrJEUlptiUh/YBjQ2HsME5F+sYVojMlLhTtBxz4MbvEo9PsGjriR3eR3GHUpPNgGRl0OCz4JPzRL8XZ4tx9UrwfH3R1/PC07g1QKmycxFYt2ko2LgYNVdQOAiNwPfAEMTXRgxpjsFM1wKQO694SjbuHNbacyYI+Vbpyv7952P2s3hf3OhP3OhoZtYe2v8OePMHe0a2V1+nOuKXK8qtd187sv+gy4OabXku+iLUgE138koNhbZowxQJStpgL7HLun+6VlZzj+ATc0y4wRbvKtSYOhoBCKt/y1Q7vT3DAtibJbV9fCbOtGqFrjr7jyLCcSq2gLkueBKSIyynt+MvBcYkMyxuSq4ER20cAxQFAiu0p1N798u1PdxFOzR8KaJdBgd2jQGhq0gVq7RD3OV7l26wqfP+JyNrsflbjj5olYku0dgC64O5GJqhq6J08GsWS7MZknoxLZW9bB/UVwWH/oFnIkqLyUrGT7i8ACVX1EVYcAi0TE7kiMyTGJ6tGdNT3DC2tBkw7w04SYds+a15kk0fZs309VVweeqOoqoH1iQzLGpFuienSXd5yMS2Tvc7JL4i+ZFvWu+d4DPtqCpJKI1As8EZH6RJ9nMcaYzEtkd7gQCuvAF9YINVrRFgIP42ZJfBNQ4EwgAQ25jTHpFlEiPEXHGTxuXuoLmsJafN3wJA78/hVYtQjqFZW7eaLer1wQS7J9b+BoXLJ9vKp+n4zAEsmS7cZEJ1GJ8FiPk65E/MEDX2ZKjWvgwEvcPPYRyqiGAwmUrPlIUNXvVfVRVR2aDYWIMcZEahn1Yd8z4JuX3QRePtEk1PMt+R7vMPLGmByUqER4NMcZPG4eRQPHlFYRBX5P9pdy8Hl7fLU/bNvApNceLLNdNA0H8i35HnXVVjayqi1jsku6qopKz/vyqbBstpsBsnJh1DHlSlVXpFVb1uLKGGOCde4HL5/MuOFD6PvdPqWLy0uo53Py3QoSY0zGSVcfk9LztjoSmrSn+88Ps+i8x2Df0yu8yxjQvW1pgZErdySRshyJMSbjpOsKvvS8InDuG9CkPYy8GD66g0qUpCWmbGAFiTHGhLJTI7jwXej4N5g0mAm7PAp//BDRrhnXaz/JrCAxxphwKleF3v+BnoMo2jALHjsY3vw/WD633N1yPScSzAoSY0zWSXk/jQMvhmtmQpdrYN5YePwQeKMPLMucrnTp7LtiBYkxJuukpZ9GzYZwzB2uSfDh18L8j+CJQ+H1C+H32amPJ0g6+65YQWKMMdGoUd/NWXLNTOh6A/z0MTx5GAw7ExZ/me7o0sI6JBpjskJwP42AtPfT2LQKvnoGpjwBG/+EFp3dHUvrYxI7i2MIyX5PIu2QaAWJMSbrZGQ/ja0b3Bhdnw+FtUtgl32hywDY+2SoVJD00yfjPUnaoI2ZQEROFpGnReQdETk23fEYYwxVa8Ihl8PV38JJj8P2La6F19COMPV59zxHpbwgEZHnRGS5iMwOWt5DRH4QkR9FZGB5x1DVt1W1L9AHOCuJ4RpjMlBG99OoXBXanwdXToEzX4bqdWH0NfBIB5j6HGzfmpTTpvM9SXnVloh0BdYDL6lqO29ZATAP6A4sAb4GzgEKgHuDDvF/qrrc2+9hYJiqflPeOa1qyxiTNqpuLvhP7oUlX0OdFnDEDbD/OVBQJd3RlSujcyQiUgSM9hUkhwJ3qOpx3vObAVQ1uBAJ7C/AfcA4Vf0ozDaXApcCtGjRouPPP/+c4FdhjDFRUIUfx8PHd8PSb9wMjF1vhP3OgoLMHPYw23IkTYFffM+XeMvC6QccA5wuIpeH2kBVn1LVTqraqVGjRomL1BhjYiECbY6BvhPg3NehWh1450p47ECY8RqUFKc7wphlSkESqo1c2FslVX1EVTuq6uWq+mQS4zLGmMQSgbbHwaWfwtmvQpWaMOpSeLILzB3j7lyyTKYUJEuA5r7nzYClaYrFGGOSTwT27AmXTYTTn4firTDiXHjmGFjwabqji0qmFCRfA21EZDcRqQqcDbyb5piMMSb5KlWCdqe6Vl4nDoV1v8FLJ8JLJ8GSaemOLiLpaP47HPgC2ENElojIxaq6HbgKGAvMAV5X1e9SHZsxxqRNQWXocCH0+waOuxd+nwXPHA0jzoPlc9IdXblyume7iPQGerdu3brv/PnpG9DMGGOitmUdfPmE6ym/ZR3sfzYcdQvUbZGyELKt1VZSqOp7qnppnTp10h2KMSbPRT3Me2EtOOJG6D/DzSH/3SgY2gnG/RM2rU5OkDHK6YLEGGMyRczDvNeoD8f+C/pNg3anweRH4JH2MOW/SeslHy0rSIwxJhvUaQanPAGXfQq7tIMPbnQTbM15L+1NhnM6RxJgQ6QYY9IhacO8q8L8cTDuNvhjLjQ/BI67G5pVmM6ISkYPkZJqVpAYY9ItKUPfF2+Hb1+Gj++BDcthv7PdLI61d03I4S3Zjmu1JSJPrVmzJt2hGGNM4hVUhk5/g6u/gS7XwndvuWHrPxuU0mHrc7ogsVZbxphMkdRh3gtrwTG3w9+nQKsjYfyd8NjBMPf9lORPcrogMcaYTJGS6YDrt4JzXoULRkFBVRhxDrx6ZtIHhLSCxBhj0izqPiYV2f1ouGIy9LiPrzfukvSpfq0gMcaYNIu5j0l5CqrAIVdwxk89En/sIFaQGGOMiYs1/zXGmDRIWh+TBB7b+pH4WEFijMlkSeljkoBjWz8SrB+JMcakQk4XJNaPxBiTDZLZxySp/Vc8VrVljDEmJKvaMsYYkxJWkBhjjImLFSTGGGPiYgWJMcZkiYQPpZIgVpAYY0yWSMpQKgmQ0wWJ9SMxxpjks+a/xhiTwZI5lEpFbIgUHytIjDG5IJlDqYRi/UiMMcakhBUkxhiTJVIx3EksrCAxxpgskZLpemNgBYkxxpi4WEFijDEmLlaQGGOMiYsVJMYYY+JiBYkxxpi4VE53AMkkIr2B3sBaEcnMQWqSpyGwIt1BZAB7Hxx7Hxx7H5xI34eWkRwsL3q25yMRmRpJj9RcZ++DY++DY++Dk+j3waq2jDHGxMUKEmOMMXGxgiR3PZXuADKEvQ+OvQ+OvQ9OQt8Hy5EYY4yJi92RGGOMiYsVJMYYY+JiBUmOEZHnRGS5iMxOdyzpJCLNReRjEZkjIt+JSP90x5QOIlJNRL4SkRne+3BnumNKFxEpEJFvRWR0umNJFxFZJCKzRGS6iCRstj/LkeQYEekKrAdeUtV26Y4nXURkV2BXVf1GRGoB04CTVfX7NIeWUiIiQE1VXS8iVYBJQH9V/TLNoaWciFwLdAJqq2qvdMeTDiKyCOikqgntlGl3JDlGVScCK9MdR7qp6m+q+o33+zpgDtA0vVGlnjrrvadVvEfeXT2KSDOgJ/BMumPJRVaQmJwnIkVAe2BKeiNJD69KZzqwHBinqvn4PvwHuBEoSXcgaabAhyIyTUQuTdRBrSAxOU1EdgJGAteo6tp0x5MOqlqsqgcAzYCDRCSvqjxFpBewXFWnpTuWDHCYqnYAjgf+7lWFx80KEpOzvJzASGCYqr6V7njSTVVXA58APdIcSqodBpzo5QdGAEeLyCvpDSk9VHWp93M5MAo4KBHHtYLE5CQvyfwsMEdVB6U7nnQRkUYiUtf7vTpwDDA3vVGllqrerKrNVLUIOBuYoKrnpzmslBORml7DE0SkJnAskJDWnVaQ5BgRGQ58AewhIktE5OJ0x5QmhwEX4K4+p3uPE9IdVBrsCnwsIjOBr3E5krxt/prndgYmicgM4CtgjKr+LxEHtua/xhhj4mJ3JMYYY+JiBYkxxpi4WEFijDEmLlaQGGOMiYsVJMYYY+JiBYkxSSQi632/nyAi80WkRZhti7wm25WClk8XkYR0HDMmGawgMSYFRKQbMBTooaqLQ22jqouAX4DDffvtCdRS1a9SEacxsbCCxJgkE5HDgaeBnqr6k7eskYiMFJGvvcdh3ubDcb2vA872lhmTsaxDojFJJCLbgHXAkao607f8VeBxVZ3kVXWNVdW9RGQX4FuguapuF5E5wBmqmtcTlZnMVjndARiT47YBnwMXA/5ZGo8B9nZDggFQW0RqqervIvId0E1ElgHbrBAxmc7uSIxJIi/Z3hj4CBitqvd4y1fg7jo2hdjnYqALsAxYo6r3pjBkY6JmORJjkkxVNwK9gPN8g2h+CFwV2EZEDvDtMhI4ATgLN+y5MRnNqraMSQFVXSkiPYCJ3t3I1cBj3qi8lYGJwOXetqtF5EtgZ1VdmLagjYmQVW0ZY4yJi1VtGWOMiYsVJMYYY+JiBYkxxpi4WEFijDEmLlaQGGOMiYsVJMYYY+JiBYkxxpi4/D8zvzr27L8r6gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from xspec import *\n",
    "s1 = Spectrum(\"simple_grp.pi\")\n",
    "m1 = Model(\"phabs*po + ga\")\n",
    "comp1 = m1.phabs\n",
    "comp2 = m1.powerlaw\n",
    "# Parameter objects are accessible-by-name as Component object attributes:\n",
    "par4 = m1.gaussian.LineE\n",
    "# ...and we can modify their values:\n",
    "par4.values = 3.895\n",
    "m1.phabs.nH = 5.0\n",
    "comp2.PhoIndex = 1.5\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "Fit.nIterations = 100\n",
    "Fit.statMethod = \"chi\"\n",
    "Fit.perform()\n",
    "\n",
    "\n",
    "\n",
    "Plot.device = '/Null'\n",
    "Plot.xAxis = \"keV\"\n",
    "s1.ignore('0-0.3')\n",
    "s1.ignore('5.0-10.0')\n",
    "s1.ignore('bad')\n",
    "Plot('data')\n",
    "\n",
    "Fit.perform()\n",
    "\n",
    "Plot('data')\n",
    "\n",
    "# Get coordinates from plot:\n",
    "chans = Plot.x()\n",
    "rates = Plot.y()\n",
    "folded = Plot.model()\n",
    "\n",
    "# Plot using Matplotlib:\n",
    "plt.plot(chans, rates, '+', chans, folded)\n",
    "plt.xlabel('KeV')\n",
    "plt.yscale('log')\n",
    "plt.ylabel(r'counts/cm$^2$/sec/chan')\n",
    "plt.title(r'Spectrum of RXJ1131 OBSID 12833 Image A')\n",
    "plt.savefig('/home/crhea/Desktop/myplot.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now for some nice 2D plotting\n",
    "\n",
    "First lets create a nice image from our fits file\n",
    "\n",
    "dmcopy \"acisf12833_repro_evt2.fits[sky=region(simple.reg)]\" source.fits\n",
    "\n",
    "dmcopy \"source.fits[events][bin x=::.1,y=::.1][IMAGE]\" source_img.fits\n",
    "      \n",
    "      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up matplotlib and use a nicer set of plot parameters\n",
    "%config InlineBackend.rc = {}\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from astropy.io import fits\n",
    "from matplotlib.colors import LogNorm\n",
    "from astropy.wcs import WCS\n",
    "\n",
    "image_file = \"full_source_img.fits\"\n",
    "\n",
    "hdu_list = fits.open(image_file)\n",
    "image_data = hdu_list[0].data\n",
    "hdu_list.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWcAAAEVCAYAAAA8UiRHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XucXWV97/HPl3CzM1FBAsqljSCVBmqDByGRpFJKkCIKRU4TUhCsFuRAIR5bBap1PCrW1kv0JUXiLRZFYhEVMRUDak10ggSJCsELApYIShARZorBhN/5Y629s2bP3jNrz+zLWnt/36/Xes3e6/qsNTO//ezfep5nKSIwM7Ni2anbBTAzs/EcnM3MCsjB2cysgByczcwKyMHZzKyAHJzNzArIwdnMrIAcnHucpPskPSFpRNIvJK2UNJhZ/jlJK2q2+YKkDzXY39mS1mXeXyBpg6StklbWrLurpGvTMoSkY2qWL5N0j6THJD0g6f2Sds4sf7ukH0jaJmkox7k2XF/SyyStk/Roeh0+ImlmZvm/SLo/LcvPJP1jzfYrJP1I0lOSzp6kHCslvWOy8haBEvdI2tTtsthYDs794eURMQjMBQ4HLsksOx94paQ/A5C0OF3n4pz7fgB4B/DxBsvXAWcAv6iz7EvACyPi6cBhwJ8AF2aW3w28EfhyzrJMtP4z0nLuC/wRsD/wr5nlHwMOScvyYmCppFMzy78H/B/guznLUhZ/CuwNHCjpRd0ujO3g4NxHIuIXwI0kQTo77w3ARyT9PvBB4NyIGMm5z+si4gvAr+osezIilkfEOmB7neU/jYhH07cCngKel1n+yYj4T+DxnGVpuH5EXB0RX4mI/4mIXwMfAY7OLP9RRIxmNqkty+URcTPw2zxlqZA0O/3W8Oq0Zv5rSa+T9CJJ309r8h/KrH+QpK9J+pWkhyV9WtIzM8tfKOl2SY9L+g9Jq7K1dEknSdqY7vfbkl4wSRHPAr4IrE5fW0E4OPcRSfsDf0FSw6yKiJXAT0lqhV+JiK90sExLJT0GPExSc76yQ4f+U+DOmrJcLGkE2AwMAFe38HhHAQcDi4HlwD8CxwGHAn8l6SWVYgDvYkcN/wBgKC3frsDngZXAnsBngL/MlP+FJN9gzgWeRXItr5e0W70CSfo94DTg0+m0JD2GFYCDc3/4gqTHgfuBh4C31llnLck/9Kc6WbC0Rvt04A+BDwO/bPcxJS0iqSX+U01Z/hmYCbwQuAr4TQsP+/aI+G1EfBUYBT4TEQ9FxM9Jrv3haRnujog1EbE1IrYA7wMqgXsesDPwwYj4XURcB3wnc4y/Ba6MiFsiYntEfBLYmm5Xz6np8q8CN6T7flkLz9mmwcG5P5wSETOBY4BDgL2yCyUdDPw98G/AeyXt0ukCRsRPSGqy/5ZnfUl3pjc5RyQtzHscSfNIasSnRcSP65QjIuJ24AngbXn3m0P2Q+eJOu8H0/LtLekaST9Pv1F8ih2/r32Bn8fY0cruz7z+A+ANaUrjUUmPktS8921QprOAz0bEtojYClyHUxuF4eDcRyLiv0i+Er+nMk+SgI+SfNX+O5Ja3Zu6UT6SmttBeVaMiEMjYjCd1ubZRtLhwPXA36T545aUpcXeBQTwgvQbxRkkqQ6AB4H90t9ZxQGZ1/cD74yIZ2am34uIz9QeJE1xHQuckbZe+QVJiuNESXvVrm+d5+Dcf5YDiyRVbgqeR1IzuywingJeA7xR0iF5diZpZ0m7AzOAGZJ2r2kOt1u6HGDXdLnSZa+VtHf6eg5JK5KbM9vukm67E7Bzuu2MCcrScH1JhwFfAf4uIr5Us91Oks6VtEfatOxIklYs2bLsmu5bwC7pvtvx/zMTGAEelbQf8A+ZZcMkN1YvSK/7ycCRmeUfAV4n6aj0PAaUNCGcyXhnAj8Gnk9yg3guSWppM3B6y8/KmhcRnnp4Au4DjquZdwXwOZJa16PAvJrlbyXJg6rO/s4G1mbeD5HU9LLTUM3xa5fPTpd9guTr/Wi63r8Cu2e2XVln27MnONeG66fHeook8FWmO9NlO5EE7kfS+T8GLs2eP/CNOvs+ZoJyvCN9PTtdd+fM8s3ZbUlSF29OXx8K3JaWYyNJS5rNmXWPSOePAP9Bkop4S2b5CcCt6e/1wXSdmXXK+EOSD6ra+W8ENnT779ZTJH98ZnlJuhA4NiJO6XZZDCTdAnw4Ij7R7bJYazmtYbmlX+tPBjZ0uyz9StJLJD07TWucBbyApNZvPcbBuSQkXdXl4/8xSS+/x4C6XbvbeOyunns31Tn355P0VvwNScrjtIh4sOMF64B+/r0DTmuUhaTRiBjodjm6wefuc+9HrjmbmRWQa84tIskX0qzNIkKTr9XY86T4n5zrPgg3RsQJ0znedOw8+SqWlz/ozNpnbN+bqfkfkoFH8hiq6UnbaQ7OZtY3RHmCXlnKaWY2bTsBT+t2IXJycDazviGg46N6TZGDs5n1Dac1zMwKqMma816Ssr1hV0TEioZrt5iDs5n1jSZrzg9HxBFtK8wkHJzNrG8452xmVkBurWFmVkCuOZuZFVRZgl5ZymlmNm2uOZuZFZDbOZuZFVCZbgh6PGcrLGmo20WwHlNJa+SZus01ZyusiKFuF8F6jNMaZmYF5BuCZmYF5JqzmVkBlanm7BuCNoZvwlkvE0lrjTxTt7nmbGP4Jpz1MgG75I1629pZksl1rOYs6RBJw5K2Svr7zPwDJH1d0l2S7pR0UWbZvpK+JumLkgbTea+T9ANJGyWtkzQnnb9I0m3pstskHZvZz32Z1ydI+pGkuyVdnJm/UtK96X43Spqbzj9b0lA6nd3Oa2Rm7SXBzjvnm0jHc85M53SyrJ2sOT8CXAicUjN/G/CGiPiupJnAbZLWRMSmdP2/Aw4EzgA+DFwdER8GkPQK4H3ACcDDwMsj4gFJhwE3AvtlDyRpBnA5sAjYDNwq6fr0WAD/EBHXtvrEbXLSkGvt1nYS7DIj9+r9MZ5zRDwEPCTpZTXzHwQeTF8/LukukqC6CZgBPJVOStd5LLP5ABDp/Nsz8+8Edpe0W0RsBbak848E7o6IewAkXQOcnB6rkSeAkcxrawMHZuuESs25DApVTEmzgcOBW9JZHwKuAn4DLM2sdz7wf4FdgWMZ75XA7WlgJiJelM7fD7g/s95m4KjM+3dK+ifgZuDiiNgaEaumd1b9oRM1X9eubbok2GW3bpcin8K01khzyp8DllVqxxHxs4j404h4eUQ8Xlk3Ii6PiIOANwFvrtnPocC7gXPrHabOvEh/XgIcArwI2DPd92RlvkrSqKTRSU+wx3UiaDowW+X/TdJVU9sBSZU0z9RlbQ3Oks7P3GDbd4L1diEJzJ+OiOuaOMQ1ZHLYkvYHPg+8KiJ+Wmf9zcABmff7Aw9Akl6JxFbgEyQpkAlFxJkRMRARA02UuavcVM7KrPL/FhFnTmkHDs6JtIY7N50eqLeOJAEfA+6KiPdNtk9JB2fevgz4STr/mcCXgUsi4lsNNr8VOFjScyXtCiwBrk+3f06mPKcAd+Q5x7Jx7TM/f5D1qJIE544VQdKzgQ3A04GnJC0D5gAvAM4EfiBpY7r6pRGxusGuLpB0HPA74NfAWZX5wPOAt0h6Szrv+PRGJAARsU3SBSQtOWYAH4+IO9PFn5Y0i+SzdSPwummftDWlaDnlIpXFWkQk//kloIiYfC2blKTwtTRrH0lERL37Rrkd8TTFhufmPN5d3NYXTenMzLpOgFtrWFk4t2p9o0Q3BAtQBOs251atb5RozNCSFNPMrEVKckPQaY0+55SG9RWnNawsnNKwvuK0hplZAbm1htn05Em3OCVjTXNaw2x68qRbnJKxppUoreGas9k0uQZfIpXu23mmLivJZ4hZcbkGXyKuOZsVj2u4Vr0hmGfqMgdnK5R2BtBmargO5D3KNwTNpqYoKYKilMNazGkNM2sF1+DboCQ1Zwdn6wm9GsQihiY9t14997ZorrXGXpI2ZKZzOlpUDxDfGh5s36y9WjLY/rMVG87Iebz3dnewfdecre2KUrMrSjmsi0rUWqMAmRXrdUW5uVaUclgX+YagmU3EtfgucVM6M5uIa/Fd4pqzWbG4pjpe314Tj61hVhyuqY7Xl9ekRDXnkhTTzKwFSjTYvoOzmfWPEtWcnXO2ntW3OVVrrEStNRycC8zBZXrakVP176TkSjTYvoNzgfXlDZtpanfw9O+k5FxzNsuvlQHVwdMmJGD3nFOXOThbxzQKwg6o1jFOa1i/mqgW7CBsXee0hvUrB2CbSCFuqJYkOBegCGbWL7r+4V1Ja5SAa84GFKRG08N8fQvCaQ0rm67XaOjtAFaE62uUarB9B+c+UvTg5wBmbVeimnMBimCd4uBnfa9EY2uUpJhmZi1QouDstIZZCxQ9ZVRkHb92JemEUpLPELNic8po6jp67VxzNus+12ZtHLfWsF5VpoDn2qyNU6LWGg7O1hQHvLHK9GFllCo4F6AIZmYd4pyz2XidqGV2qiZbOY6/SZRPzMg3dZuDs3VMqwJZEYYlbeVxnBrpnNgJntw939RtDs5WOJMFq2YCY6N9Zed3Ozi69t05Idg2Y6dcU7d1vwRmNVoZrBrtKzs/+7qXUi82Xkhs33nnXFO3OTgb4IDRKIfcjuvimnJ3bZ8xI9fUbQ7OBvRuwJCG0KqaefPHLod8NWwrv0BsZ0auqdscnK1naX4aXJcMVd/XC8zSUMdbeVh3BGIbM3JN3ebgbD1JGiKGk9fV2u/6dN76oer8yrIytvKw5gXiSXbLNbWapAMlfUzStXnWd3C20srbpC6bupCGxgVIB8zpKdO3gVanNSR9XNJDku6omX+CpB9JulvSxQARcU9EvCZvWR2c26BMf6xFN502zdlcszRUTWlk89Can0lvzGfCJnb+vdZXtg+3FuecVwInZGdImgFcDvwFMAc4XdKcZsvp4NwG7f5j7acgMdG1nOg6SEPE4sw+5g0l6Yx5yTaVZQBck84brn/sbvcG7Kffd7s1mXPeS9KGzHTOuP1FfBN4pGb2kcDdaU35SeAa4ORmy9qXwbneVw5Jz5V0i6SfSFoladd0/pCksyWtlHRMVwueKltNpV0mug7VwDo/k3+eNwTL0vRGenMwhoHlVNett892Xe+8Qde/79ZJ0ho755qAhyPiiMy0Iudh9gPuz7zfDOwn6VmSPgwcLumSyXbSd8F5gq8c7wbeHxEHA78GcueGrLOyrStqW1rUDXjzkjTG4E1bYHma7liWTJpf8zPP/lqkkgOfiGvNrZXcENw11zQNqnfoiPhVRLwuIg6KiHdNtpO+C840/spxLFC5i/pJ4JT09QjwBPAb4MkOl9XqqG1l0bC33yqqLTNYMsTI4OVJgAauWnwagyclrwdP2pKslza5Y95QNUi3u9Y62f5da26tgE40pdsMHJB5vz/wQLM76X4fxc6r95XjKODRiNiWmbcfQES8J51X05XBiqheTfOc4Q+wYv4QgzdtYeS4WSwYXsOVnAtQnVeVyU1XctCTtfawMlElZdFOtwIHS3ou8HNgCbC02Z30Y8253leOeh+T0e6C2NQ1Sj1Ubv5lA+gK/XpHDTq1btUirhg4jy8OnMyC4TU7tl2fbBvDY4OyA3JvaENTus8Aw8DzJW2W9Jq0kncBcCNwF/DZiLiz2bL2Y8253leO/waeKWnn9MLm+hoi6Srg1LaUss9NVkOtl8qo3uBbP4Q0dtnM0S2MDF4O84Y4kdXMWbyJ+zmAM+dfy+BNW1gQ42vdUxlnwzXr9pI0mr68LiLOnMo+Wtk1OyJObzB/NbB6Ovvux5pz9StH2iJjCXA98HXgtHSds4AvTrajiDgzIgYiYqBtpe1TeW6WQU3AzPT8yy6XhpLUxTVJauOy0UuZzX2s5kQGb9rCFwdOZp2+VW1mN26/mTK1Okdc5ht+3Sh75f9tqoG5yZrzpE3p2qnvas4RsU1S5SvHDODjEXGnpDcB10h6B3A78LFultOaU+1UUpPiWMga1rIoWTYfRpZczlVxB+eNXsEVA+exdmAhf67jd7R/jrGdVerVnjsxpGkZlLHsgdiav2v2wxFxRDvLMxFFOLXaCpKiF69lN7+mZ49d76acViWdSaSh6g28yuuKGKa63kLWsG7VIgZP2sLcgY2sm79obC46cxPQikcSEVHvnlFuBx/x9Fi+4chc656km29zcO4BvRqci6ZhwJ5P9WZe1szRLTw+MKvajpnlSeuNK7lo/I1EK7RWBOfnHfGMeN+GebnWPVlf7Wpw7secs5VUtcacjpNRaZlRCcyVGnM1LTEfRm6YhVaRtMhYnqy3Yv5FY4YOzaYyrPd5yFCzFhvT8WQ40zpjGTBviAXDa8YEZqDasWTj6NzkfdqFu/qascOLVtQOgmS9ocnu213l4GxdVS8AZkeQm2i7akBdMgTLSFpckHQsOWf4A8nrkfNheVKDBpLUx2Kq3bdrUyQVMdyeVIcDfneVqbWGg7N1VaOBhsbVZDNjaUASdKuDGsVQkrK4Jqk9VwNxqtJlO5v6YMkQLNmRHmkUpJvVyieHW+slrTV2zTUx9YGPWsLB2VqqXTXDbEpDSsbJAKrjYLB+iAWLk9YYAFePLoVlMHdgY7JdZfjQmp5/rB+qBnmtmn7w7MQDYm3qeiatIel5ko6uM3+hpIPaVywrukZBpxNjWWePMXjTFpg3xODI+WwcncuCxWu4efGLWTpwNVctPo0TWc3cgY1oflrbznb7XkV1LI2IobFjPFO/i3izwXaq41Fb+/TKA16XA4/Xmf8E1VFwrR916ut5bYCstG+upCdGBi+H9UPMHdjI0oGr+eYjx7MPv+TK1cv460c+x2WjlwLsSHfMG9qR2khbb9S21qg3uH47Bj0qU4qjVz5Ieunp27Mj4vu1MyNiAzC7LSUyyxgTwCopjCU7ctJXxR2cE3swh00AfH3P+cx55B64HPQSePzVe/PNR45Pli+Hy4Zfn9wsnJd05a4G3eHxAagVY21M11SP0eqylemDZCK9FJx3n2DZ01pZEJu+otZuWlGuSgCttm1O93nGTz7HbO5jIWu5cvUyjl09jN4EDMBtd0DcDJoHC1nLU19O+i9sYg4xDJcOXJbUmDM9CmFsfnuic5huwKq0SmnUYmU6x+iVYNpqle7beaZumyw43yrpb2tnSnoNcFt7imRTVdR/yKmUa8KAVWnJsQz0quCy0UtZy0IWnvhVYh7JeIJ/njzY7aZH4N6fwBm3fI5P7/lKIGnzPHN0C5fqGWP2l3cUvEblm+gcKgMw1fZKbNRdvBW/y6J+WHdTmWrOk92SXAZ8XtJfsyMYHwHsCvxlOwtm/a1erTXb3K36eskQIzcMcfVJS3n8gb1hA/AS4BOwaF/45QOwz3rgHpJHKqQqo9TF4h2dWaqdWiaRN/eczVO3qqleM4r6Yd1tTQTevSRtyLxf0cnmdBPWnCPilxHxYuBtwH3p9LaImB8Rv2h/8azf5AlelaCzkDXVlMQVA+fxroOXwYEkDxR7BJgL+5yYbvNSWMtCVnNicmNwGdXHVFXyzXkHPWoU9Bo+v3CSXHbtvInSHTY9TT59u6vtnCcc+EjSsRHxtfT1cyPi3syyUyPiug6UsRQ88NHUTaUVRCVwnRN7cPXo0mp75rW3HM+Wowa5kZdyJefyzUeOZ9OeB/JL9mEVi9nEnGpPwrzHrC3fVPLB9QZsqrcPD9bfWCsGPtr3iOfEuRvOyrXukN5d3FHpJH03Il5Y+7re+37X7uDcq/+0ec4rO+DRmCBZGWmO5IGta1kIwGzu4yhu4QH2ZS0L2cQcgLFDhNakNPKWpdn1KuqWv8nfaa/+DeTVquD82g1/k2vdt+uyQo9Kpwav6723NurmmMrtNFknjWxAGvdkbeCyxa+H5XDmqmurQRjgFo6qBut1qxYlPQczY2mwZGy6oZnA18x6lane/rMj7GXlSYPY1ATiSXbNNXXbZME5Gryu9956UKc/FLKtGuqlEqqBLO1AcqmewYLhNZyz+ANJzRi4dP77uY/ZXMgHmc19LFi8ptplqtJTsBo4G7TUqDfw0rSCZPYRWJXu4ukxKg8EqGg4GFMf15pbpcmcc1dNltZ4FPgmSS15Yfqa9P2CiNij7SUsiSLnnHvh6/C4Fhvzk15/G0fnJi0vSAbR38Qc1rKIOzmIw1b9lKsWn8aZ86+tDno0bmwNGqdNWl3+au/G5WNbhUy3TXPtMXpVK9Ia+xyxfyzdcEGudZfrkkKnNU4G3gu8J/O68v6U9hbNWqXX/mEreeJKTbmSR16x6iLW6VssZA2Hzf8pCxav4cz51yYD7ZOMw5F9jFVtSqHV16n6LSDt6CINjanBj6md17TqqG0TnUev/Z7bpSztnCdrSvdflQnYBGyqmWd9rhPdi8cE1PSpJ5qfBNuR42ZVBy8aPCkZBGkOm7hs+PVsHJ1LDCedTs6JPao17DGpjDqPtppS+eqo7nf9jt6NMZy8r8zLntuYgZnm9XYtuFvK1AllslHpJOmtkh4Gfgj8WNIWSf/UmeJZ0bWre3HDnnSVEeSG044kaaCD5H0ltXHp/PezdOBqAJYOXM2VXJTsI/uUbQ3BNeN77uVRrzlcvRt52TTMQtagVUkQPif2YCFrkocBLEs+MB4fmDXhWB82fU3mnAs92P4yYAHwooh4VppjPgo4WtLr214661sNB+GPscG0MtZGJchtYg5z2MTNwy9Ox9xYw9WjS/kUpzFzdMvYjibzkuZ0jY6VlefmXKPWGJX561YtSp7KctysasuSKwbOg+U7AvfM0aRjTGXkPQfo1kpaa+yWa6Lgg+2/Cjg92/kkIu4BzkiXmTVlOjXUcTfOliVBbfCmLawfmMfgSVuYwyY+/MgyHmBfLhu9lDlsYu7ARg7n9mQf88fmc8c8TSVzjGZaSkw01nN1ONI07RLDmSezpAZv2lJ9SMDIcbOqY0vn7bHYqDw2Xs+kNYBdIuLh2pkRsQXYpT1Fsl41nU4X1adtsyP4DJ60pXpTcN7oer44cDJX3rKM1+25nJdyI+sH5rGYpI3aB7mQkcHLqzcHgTEPh60N/M0MdJS3pcfIcbOqTebO5Urexls5c9W1jAxePqYr+XQ5Vz2xXgnOT05xmdk4zQaNMYPd1wlaIzfM4qrh0xi5YRaPbd2bB9gXNsKVq5cxa9kIh/7jPfzZs4Z5G29lE3NYEEcnNdTMzcVYzJibgpPe3Gvi3LKj6FXz5kuSR2xdybkcy7dhOSyIo3l8YNaYx1e0qnmdjVWmds6TBec/kfRYnelx4I87UUBrjTJ+1Z30Rt2SIc5cdW3SSxA4Y/XnuPPcA2EZrPkA3HsZ6GA4dvUwa285vrrZ4E1bWBBHJ0OO1jw3sF3N6YAxY0evm7+oWotey6Jx3cizP611euYZghExIyKeXmeaGRFOa5RIEWpgeYNNttbYKL0QMQTXDMFyOID72bTngXzqxFdyO4fDAbDoWPgWwABsOXEQPpOkEs5Z/AFGbpjFuvmLkp6DS6bWWiOP2vJXmtGNOadljBsRrza1Ui/vblPTS923zVpmuoPuZ8dGrnToWDC8hrUs5HYO534OSG78vRZ4A5zxWthy8yAHjt7DluWD3M8BrFh1UfVBruv0rbaPtVwvb14NxkvGthbJ1rAb1eaL8CFbZmVKa3S/7m5WozYY1cu/VrpCr5u/iHUsYvZwMhLdobfcQ7wU9CvY8pFBHmJv1g/M41SS0W3PWfwBZo4uJYZnIbW/Jlovb15twUHaZG7J2DJMdFOw17tod0IRUhZ5lKOU1tMmCzgNh9lckrRvfnxgFlr1fgZP2sIVR52XjEa3Z9LmeePoXC4duKw6hvM6IGLWuOZz2Q+Bestaea6Vc9L8oaT23MQwog7M01NpSlcGTmtY1zUTcGqD5sjg5dVlIzck3bPfwZtZyNrqoEirOTG5AXjN2G3rafVIcPX2V80pD4/9ZlBZ7rxy+zTZzrnQPQTNuqbRjbq6o8gtGYLlcN7oFczi8aT9cDqWxloW7bgBSOMAmB1/uXKcycoH1B1eNLvP2nPJ/mw0znPt8R2wW6cnHlNl+RV5yNBeVH0oa2agoAXDa6pdpCuygwvVG6+im2mCej0T2z1kaZm1YsjQgSMOicM25Iux39FLCj1kqFnLTKf2Ny6oVp6WnY6twbKk9QVLhnZ0kV6/48Zb9fFU88buZ6plbsW5ZG8M1ms22EplD8yt1Cs9BM1aZjoBYsJUw/qkSVolpww7RnirBOXaPO+YJnnTKPNUgvRkPR+nu39rrExjazit0SJOa0xfK756T5S2yPbCa3TMXvj636takdZ42hF/FM/b8O+51r1DRzqtYe1RtlpXs0FxwnEw5tXJ59appTa6GdcqZfsd9L7ydN/ufgmsbXq5Bjiu5122g0e1NcSOefW2r3Cut3+UqZ2zg7OV0pigt34IGGq8nPHpik4HTadLiiEQWwswbkYeTmtYYeVJCWQfNVWR5yklndbt41uiZ0alM+umPAGt+tSQHOs6/2vgpnTWYg4sO+TtOdeuDif+XZRXmZrSOTiXhL8W75B3UKJmHjVVz1SfimLFFYjtT83INXWb2zm3iNs5d4dvtPWPVrRznvHCufF7627Kte7IwKyfAdlnqK7o5PgarjlbqeUNzPUGHbL+EyG2b5uRa6LLAx85OFtfyI5iV2/AoXbr5IeBP3gmEDQTnLvKwdk6qptBqvbZfLWv26mTqZeypHm68SESIbb9bkauqdscnK2jpho4pjI6XFmCVL/qzu9HPLV951xTtzk4d0g/f9Vsxbk3+kdu9UBJ1uMC2DYj39Rlbq3RIm6tYdZerWitocOOCK7bkG/l58uj0pm1QqvaMVuP25Zz6jIHZyu9iR7x5Jy0jRE4OJt1ymS9BF1DtioHZ6vHQWJ6pnr9WjWIv/WAAH6Xc+oyB+cO8tfo6Zmo40grA6qfhtJ67XhQ7pQEsDXn1GUOzlZ4eUaXK/IHX5HK1q0PinY1hWya0xrtI2l3Sd+R9D1Jd0p6Wzr/Akl3SwpJe2XWP1vpX2T6eoukjen02sx62zPzr8/M/4ak2ZLu69xZWlaRglsrdLMm3WvXsmklCs7d7wbTvK3AsRExImkXYJ2k/wS+BdxjBZFZAAAKy0lEQVQAfGOS7VdFxAV15j8REXNbW1Sz8fo+QHZTJTiXQOlqzpEYSd/ukk4REbdHxH11NnkCGKkzP69HgO3Almnsw0rOeeMeUaKac+mCM4CkGZI2Ag8BayLilkbrRsSqiHhPZtYrJX1f0rWSDsjM313SBknrJZ2S2f7UiLg/Il7U+jOxsmhU23XQLiEH5/aJiO1pCmJ/4EhJh+Xc9EvA7Ih4AXAT8MnMst9Pu2ouBZZLOmiynUm6StKopNEmT8F6xEQpCgfu1qv8v0m6ako7eAr4bc4J9korbJXpnOmfQX6lDM4VEfEoSY75hJzr/yoiKo1kPgL8r8yyB9Kf96T7PDzH/s6MiIGIGGiu5NYp/XTzrR8+DCr/bxFx5tR2QDM1Zw+23wxJsyQ9M339NOA44Ic5t31O5u0rgLvS+XtI2i19vRdwNLCpleW27uinm2/9dK5T5pxzWz0H+Lqk7wO3kuScb5B0oaTNJKmO70v6aJ1tL0yb330PuBA4O53/R8CGdP7XgX+OCAfngumHmqG1WYmCs4cMbREPGWqd0q8PtW3JkKF/cETwppxDhp7f3SFDy9jO2ayv9WNgbpkStXN2cDaz/vEUSc+HEihjztmsLznn3gJB0qUsz9RlrjmblYTTGS1SkrSGa85WCkWuNRa5bFajRK01HJytFIpcayxy2dqh1B9GJQrOTmuYWVNK/WFU6b5dAg7OZtZfClArzsPB2cz6h9s5m5kVUOUBryXg4Gxm/aPSzrkE3FrDrGBK3Rqi6Nxaw8ymqtStIYoucPdtKw/X1KxvuPu2lYlratY33FrDzKyAHJzNzArITenMzAqqAPnkPByczax/NDe2xl6Sss+0WtHJJ3A7OJtZ/2gurfFwN58h6KZ0Zj3IzSMbKFFTOgdnKwUHm+a4eeQEStJD0MG5T5Ut2DnYWEuUqPu2g3OfcrCzqSrbB/sYlRuCeaYuc3A2o+QBp8NK/cFeopqzW2uYUfKAY80pQODNw8HZzPqHewiamRVQiQbbd3A2s/7hgY/MzAroKUoz2L6Ds5n1F6c1zMwKKLpdgHzcztnMrIAcnM0Kwh1hLMvB2awg3BHGspxzNrM+Up7mGg7OZtZHytNF0MHZzPpIeXqhODibWR9xzdnMrIDKE5zdWsOsRNzcbrqC5IZgnqm7HJyt9PopYLm53XSVZ7R9pzWs9BywLL/ypDUcnM2sj7i1hplZAbnmbGZWQK45m5kVkLtvm5kVkNMaZmYF5bSGmVnBuOZsZlZADs5mZgXk1hpmZgXk1hpmZgXktIaZWQE5rWFmVkDdqzlLGgD+DXgS+EZEfHqi9T1kqJn1kdYOGSrp45IeknRHzfwTJP1I0t2SLk5nnwpcGxF/C7xisn07OJtZH6ncEGzZYPsrgROyMyTNAC4H/gKYA5wuaQ6wP3B/utr2yXbstEYLSep2EcxsQg/eCEN75Vx5d0kbMu9XRMSK7AoR8U1Js2u2OxK4OyLuAZB0DXAysJkkQG8kR8XYwblFIqKtkVnSaEQMtPMYReVz97m3SkScMPla07YfO2rIkATlo4APAh+S9DLgS5PtxMHZzKy16lXUIiJGgVfn3YlzzuVxXbcL0EU+9/5U1nPfDByQeb8/8ECzO1FEtKxEZmb9Js053xARh6XvdwZ+DPw58HPgVmBpRNzZzH5dczYzmyJJnwGGgedL2izpNRGxDbgAuBG4C/hss4EZHJynrF77RkmrJG1Mp/skbUznL5J0m6QfpD+PzWxzjKQNkv4lM2+lpHsz+5qbzn+GpC9J+p6kOyW9Op0/W9I3Mttfkrav/JGkl2bm12t7iaRvpPu4rz1XC9Lr8YP0fDZMcr2OkbQy8/o3mfX+aaJ9pvNXptt9o86d9MJpcG32lLRG0k/Sn3uk88+WNJROZ3e14C0iaYak2yXdkL5v9PdfuHOPiNMj4jkRsUtE7B8RH0vnr46IP4yIgyLinVPZt28ITt1K4EPAv1dmRMTiymtJ7wV+k759GHh5RDwg6TCST9T90mXnAQuBd0g6JCJ+mM7/h4i4tuaY5wObIuLlkmYBP5I0ppdR2p5yCXAosC9wk6Q/TBdfDiwiyYndKun6iNg05SvQvD+LiIcrbya4XrXWRsRJefZZYrXncTFwc0T8c/pBejHwpu4Ure0uIqlhPj0zr97ff19xzXmKIuKbwCP1lilp8PxXwGfSdW+PiMoNgTtJ2k/ulr7fiaTb0lPUv8s75rDAzHT/g+nxt5E0aK+U5WTgmojYGhH3AneTtLustr2MiCeBSttL0m23A1tynn5L1V4vku6tjQJ1Hr9J91E5rzI6Gfhk+vqTwCnp6yeAkXQqx/BqE5C0P/Ay4KM5Vu+pc5+Ma87tsRD4ZUT8pM6yVwK3R8TW9P1HgW8DX4+IuzLrvTP9Cn8zcHG6/oeA60nu/M4EFkfEUyRtKk9Nt9sPWJ/Zz2Z21NLrtb0kIirbvqjZE21CAF+VFMCVNY35x1yviPg2yTWpmC/peyTn/feZ/F3dfUbERenyUymHeuexT0Q8CBARD0raO329qovlbIflwBtJ/p6zxv399+C5T8jBuT1OZ0ctsErSocC7geMr8yLiRpI0R9YlwC+AXYEVJF9n/x/wUpLeRccCBwFrJK2NiMeyh6lTnqD+t6RONtU5Ok3r7E1S7h+m3z6gwfVKfRf4g4gYkXQi8AXg4Bz7LJNx59HtAnWCpJOAhyLiNknHZBY1+vvvK05rtFjajOZUYFXN/P2BzwOvioifTrSPiHgwEluBT5CkJCBpwH5duuxu4F7gkJrNG7WxbEnby6mqpHUi4iGS63AkNL5eme0ei4iR9PVqYBdJe020z7JpcB6/lPQcgPTnQ90rYdscDbwivRF9DXCspE9N8PffVxycW+844IcRsbkyQ9IzgS8Dl0TEtybbQeafUiS5xkqLkP8maTuJpH2A5wP31Gx+PbBE0m6SnktSy/wOSVvLgyU9V9KuJDcNr5/yWTZB0oCkmZXXJN8cKuc07nrVbPvs9Dog6UiSv9lfTbLP0pjgPK4HzkpXOwv4YndK2D4RcUnawmE2yd/j1yLijAn+/vuK0xpTpKR94zHAXpI2A29Nm9EsYfxX9AuA5wFvkfSWdN7xaU2pnk+nrTFEksZ4XTr/7cBKST9Il72ptqVCRNwp6bPAJpKbhedHxPa0zJW2lzOAj0+l7eUU7QN8Po2xOwNXR8RX0mX1rlfWacB5kraR3ARaEhGRfjg12meZ1D0PSbcCn5X0GpIP5f/dxTJ2WqO//77iHoJmZgXktIaZWQE5OJuZFZCDs5lZATk4m5kVkIOzmVkBOTibmRWQg7P1NEnb02En71Ay3Ooza5a/XtJvJT2jW2U0q8fB2XrdExExN31KxSMkw65mnU7Se/IvO14yswk4OFs/GWbHCH1IOohk6NU3kwRps8JwcLa+IGkGybgk2fFEKqPhrSV5zNDe3SibWT0Oztbrnqbk8Ve/AvYE1mSWLSF5MMFTJE967qfxK6zgPLaG9TRJIxExmN7wuwH4j4j4oKQXkOSaH0xX3RW4JyIWdKusZlkOztbTKsE5fX04ydCbB5EM3v5YRLwrs+69wDER8bOuFNYsw2kN6xsRcTvwPZJ0xhKSge2zPp/ON+s615zNzArINWczswJycDYzKyAHZzOzAnJwNjMrIAdnM7MCcnA2MysgB2czswL6/xDi3YezzQuVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hdu = fits.open(image_file)[0]\n",
    "wcs = WCS(hdu.header)\n",
    "\n",
    "plt.subplot(projection=wcs)\n",
    "plt.imshow(image_data, cmap='jet', norm=LogNorm())\n",
    "plt.colorbar()\n",
    "plt.title(r'RX J1131-1231 Image A')\n",
    "plt.xlabel('RA')\n",
    "plt.ylabel('DEC')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "from astropy import cosmology\n",
    "from astropy.coordinates import *\n",
    "from astropy.cosmology import Planck13, default_cosmology\n",
    "default_cosmology.set(Planck13)\n",
    "z_RXJ = 0.658\n",
    "d_RXJ = Distance(z=z_RXJ, unit=u.kpc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We now have the distance to our object: 4064301.700848627 kpc\n",
      "We are looking back 6.239697929569168 Gyr\n"
     ]
    }
   ],
   "source": [
    "print(\"We now have the distance to our object: \"+str(d_RXJ))\n",
    "print(\"We are looking back \"+str(Planck13.lookback_time(z_RXJ)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coordinates\n",
    "The International Celestial Reference System (ICRS) is the current standard celestial reference system adopted by the International Astronomical Union (IAU). Its origin is at the barycenter of the Solar System, with axes that are intended to be \"fixed\" with respect to space. ICRS coordinates are approximately the same as equatorial coordinates: the mean pole at J2000.0 in the ICRS lies at 17.3±0.2 mas in the direction 12 h and 5.1±0.2 mas in the direction 18 h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.coordinates import SkyCoord\n",
    "c1 = SkyCoord(ra=51.60*u.arcsec, dec=59*u.arcsec, distance=d_RXJ, frame='icrs') \n",
    "c2 = SkyCoord(ra=51.58*u.arcsec, dec=59*u.arcsec, distance=d_RXJ, frame='icrs') \n",
    "sep = c1.separation_3d(c2)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/crhea/anaconda3/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: The mpl_toolkits.axes_grid module was deprecated in version 2.1. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist provies the same functionality instead.\n",
      "  warnings.warn(message, mplDeprecation, stacklevel=1)\n",
      "WARNING: Cannot determine equinox. Assuming J2000. [aplpy.wcs_util]\n",
      "WARNING: Cannot determine equinox. Assuming J2000. [aplpy.wcs_util]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Auto-setting vmin to -8.000e-01 [aplpy.core]\n",
      "INFO: Auto-setting vmax to  8.880e+00 [aplpy.core]\n",
      "INFO: Auto-setting resolution to 64.5892 dpi [aplpy.core]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnsAAAH+CAYAAADkochhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmcXGWZ9//PZcKSgCzKJgkj3cCMSUdQ0uCSxJUlqIDjhOngjy3PQBRtxe0JqAzp6KNCnkfUMeqYQYKgTLfJuARBkHHUScct3WHLAgIdHQKKIKBIIpDk+v1xzqk+Xam961SdOvV9v17nlapTZ7mrujt99XXf132buyMiIiIi2fSCZjdARERERJKjYE9EREQkwxTsiYiIiGSYgj0RERGRDFOwJyIiIpJhCvZEREREMkzBnoiIiEiGKdgTERERyTAFeyIiIiIZpmBPREREJMMU7ImIiIhkmII9ERERkQxTsCciIiKSYQr2RERERDJMwZ6IiIhIhinYExEREckwBXsiIiIiGaZgT0RERCTDFOyJiIiIZJiCPREREZEMU7AnIiIikmEK9kREREQyTMGeiIiISIYp2BMRERHJMAV7IiIiIhmmYE9EREQkwxTsiYiIiGSYgj0RERGRDFOwJyIiIpJhCvZEREREMkzBnoiIiEiGKdgTERERyTAFeyIiIiIZpmBPREREJMMU7ImIiIhkmII9ERERkQxTsCciIiKSYQr2RERERDJMwZ6IiIhIhinYExEREckwBXsiIiIiGaZgT0RERCTDFOyJiIiIZJiCPREREZEMU7AnIiIikmEK9kREREQyTMGeiIiISIYp2BMRERHJMAV7IiIiIhmmYE9EREQkwxTsiYiIiGSYgj0RERGRDFOwJyIiIpJhCvZEREREMkzBnoiIiEiGKdgTERERyTAFeyIiIiIZpmBPREREJMMU7ImIiIhkmII9ERERkQxTsCciIiKSYQr2RERERDJMwZ6IiIhIhinYExEREckwBXsiIiIiGaZgT0RERCTDFOyJiIiIZJiCPREREZEMU7AnIiIikmEK9kREREQyTMGeiIiISIYp2BMRERHJMAV7IiIiIhmmYE9EREQkwxTsiYiIiGSYgj0RERGRDFOwJyIiIpJhCvZEREREMkzBnoiIiEiGKdgTERERyTAFeyIiIiIZpmBPREREJMMU7ImIiIhkmII9ERERkQxTsCciIiKSYQr2RERERDJMwZ6IiIhIhinYExEREckwBXsiIiIiGaZgT0RERCTDFOyJiIiIZJiCPREREZEMU7DX3hY2uwEJyOJ7Ar2vVqP3JSKpoWCvvWXxP+4svifQ+2o1el8ikhoK9kREREQyTMGeiIiISIYp2Gue01NwbLPvr/eVzfeVxHuq9rpJtLXZ76vZX1cRaVEK9pqnlf7jTkNb9b6af38Fe8lodlsV7IlknII9ERERkQwzd292G9qKmekDF5FMcXeLHh9t5ttquMbv4DZ3n1vHZolIaGKzG9CeRprdABGROukc82wb8K4artIHB9WjNSKyOwV7IiJSN4Z+sYikjX4mRUSkbgzYo9mNEJExFOw1xumo4k3a3ewOGNzS7FZIMpZHD5TZE0kfFWg0WFCgoTF7IpIVnWMKNP7GzP93DVd5Pwy7e3cdGyYiIU29IiIiIpJhyraLiEjdqBtXJH30MykiInWjAg2R9FE3bhtYvPgA3Dty2+LFBzS7Sdk0u6PZLRBpuiizV+0mIslRgUaDqUBDJEYVuhkwtkCj08w/WcNVzlGBhkhi9AeViDSPAr3M0Zg9kfRpSDeumb3MzH5uZs+a2Udi+48wsx+b2WYz22hml8ReO9zM/svMvmdm+4b73m1m95jZnWY2aGbTw/0vDq/zFzNblnfv38Qfx84fiu0/K7z/LjMr+pelmc01s/vM7AEzuyy2/ydmdmT8XiIiIiJp0Kg/wJ4A3g+8PW//DuDD7r7ezF4IDJvZ7e6+KTz+fQQLL54D/Ctwo7v/K4CZnQFcDcwF/gr8MzAj3Ep5o7s/nrdvA/AO4KvFTjKzCcCXgJOBrcA6M1sdtlVERFCBhkgaNSSz5+5/cPd1wPN5+3/n7uvDx08Dm4Ep4csTgF3hZuExf46dvg/g4f5n3H2QIOjL91gF7dvs7veVOexE4AF3H3H354B+4MzwtSeAnZXcS0Qky6Jgr9pNRJKTmmpcMzsSeCXwy3DXMoJM27uBb8SOe6+ZPQgsJcj+leTuJ8SfAj80s2EzW1hlE6cAD8Webw334e7vcPeHYvdaCAwV2Spy6qmTuPfeKdx//1QuvXT/gse8610v5O67p3DHHYezZs1LmDYt+C9zjz3g2msP4u67p3DnnYfz+tfvXfD8H//4MGbO3BOAl750Ir/+9VROOWVSpU0UKU6Vye1mzP9xqsYVSZdU/IyFY/L+A/hAlL1z998Cr8s/1t2/BHzJzN4JXA6cX8WtZrn7I2Z2CHC7md3r7v9daTML7CtWyryc2FqRFZ6T84IXwJe+9GJOPvn3bN26g3XrDmf16m1s3jwmMcqNN/6Fr371aQBOP30yV1/9Ik477VEuuuiFABx77MMcfPAL+MEPDuOEEx6hWOH1lCkTuO22Q/nwh5/ghz/cXq55IuWp8KLd5MY6G7gydSLpklhmL8zA3Rluh5c4bg+CQO+b7v7tKm7Rz+5jAEty90fCf/8AfIega7ZSW4EjYs+nAo9Uc/9KnXjiXjzwwPNs2bKD55+H/v5nOPPMybsd9/TTo9HbPvtYLpibPn1PfvSjIGh77LFdPPXULrq79yx4r8MOm8APf3gYl1/+JDfdtA2A88/fl+9+9xB+8INDuffeKVxxxei8fOeeuy933RVkDK+//uB6vWURyQjNsyeSPon9jEUZuFLHmJkBXwM2u/vV5a5pZse4+/3h07cC95c6Pu/cfYAXuPvT4eNTgE9Uej6wDjjGzDqAh4H5wDurOL9iU6ZM4KGHduaeb926k1e9aq+Cx77nPS/kQx/anz33NN70pt8BcNddz3HmmfvQ3/8MRxwxkZkz9+SIIyaybt1zu51//fUHc/nlT7Jq1bYx+088cS9mzHiYbducdesO5+abt7F9u/Pxj+/PrFm/449/3MWBB6ZmFICIiIgU0ZA/qMzsMIKxHPsBu8zsA8B04FjgXOAeM7szPPxj7n5LkUv1mtlJBIUeTxLrwg2nPdkP2NPM3g6cklcpeyjwnSC+ZCJBZe+t4bl/D3wROBi42czudPdTw4zkNe7+FnffYWa9wG0ExSPXuvvG8X0yhVmBDuNik19/+ctP8+UvP83ZZ+/D5ZcfwAUXPM611z7NtGl7MDR0OL/97Q5+9rNn2bGj8L3+8z+3c+65+3LddX9h+/bRe9x++3aeeGIXAN/+9jPMnr03O3fCqlXb+OMfg/1PPrlrfG9URDJH1bgi6dOQYM/df0/Q7ZlvkMJj4Ypd55ISrx1Z5twR4Lgir32HoFs3f/8jwFtiz28BigWidbN1606OOGJC7vnUqRN45JGdJc4Iunq/8pWDgMfZuRM+9KEncq+tXfsS7r//+YLnLV36J845Z19WrjyEM898lJ3hbfJjS/cgCNWCKyJSiiZVFkkf9cOl0Lp1z3LMMXtw5JET2WMPmD9/H1av3rbbcUcfPfpf6lvfOikX0E2aZEyeHMTQJ520Nzt2sFtxR9wHP/gEf/7zLr72tYNy+04+eRIHHvgC9t7bePvb92Ht2r/yox9t5x//cR9e9KLg20bduCKST1OviKSP/gBLoZ07obf3j9x222FMmADXXvs0mzYFwdqSJQcwNPQcN920jd7e/TjppEk8/7zz5JO7OP/8YJq/Qw4Jqmt37YKHH97JueeWn/7v/PMf4/vfP4ylSw9k48bnGRz8KzfccDBHHz2RG298huHhYLzfpz71FD/96UvYudO5447nWLAgf35qqZnWiZUMUGZPJH2s2FgwSYaZOYw0uxklnX/+vnR378X73vfHZjdFpHYKnhukE3fPDcfpMvN/r+Eqx8GwuxddrlJEaqc/wEQkmxToNYUKNETSR4OuZDdf//pflNWT5Gh1jUxLcp49M/ugmW00sw1m9u9mVnh5IBEZQ8GeiDSWMm6ZllSBhplNIVgis9vdZxBMgTW/3u0XySJ14zbG6eEmIknTWL1myC0PmXA37kRgkpk9D0wmoVWMRLJGBRoN1goFGiIilesEeBewEGAqzByq4SqHwW+BeHn/cncfs8a4mV0CfArYDvzQ3f+/mpos0maU2RMRkXEJg7LlAMeZ1ZpBeLxUNa6ZHQicCXQATwErzewcd/9GjfcTaRsK9kREpG4M2KOW3yxFlnSMOQnY4u6PAZjZt4HXAgr2RMpQsCciInVjBhOTCfb+B3i1mU0m6MZ9M8Ga6yJShoI9ERGpGzPYY0L546rl7r80s1XAeoLQ8A5ihSEiUpyCPRFpDaqybQk1Z/Yq4O6LgcXJXF0kuxTsiUhrUKDXEmoesyciidGkyiIiIiIZpr+/RESkfoxgbQsRSQ1l9kS0Vmth7fy5tPN7H68kF8cVkZroR6wxtFxammksWGHt/Lm083uvzWhVbBTsiUhqaLm0BtNyaSKSLZ24u0XPuvcyH5pa/VVshOFSK2iISO3095eIiNSPxuyJpI7G7ImIiIhkmDJ7IiJSPxqzJ5I6+pEUEZH6UbAnkjr6kRQRkfrSmD2RVFGwJyIi9aPMnkjq6EdSRETqR8GeSOqoGldEREQkwxTsiYikTasv1zahhk1EEqNke2NouTQRqVzrLdem5dJEUkzLpTWYlksTkWzJWy5tP/OhE6u/iv1Iy6WJJEXduCJJa/UuOZFqRJm9ajcRSYyCPZGktV6XnIiIZIj+nhIRkfpSwYVIqiizJ82hrk2RbFI3rkjq6EdMmkNdmyLZpGpckdTRj6SIiNSPoW5ckZRRN640lrpvRbJN3bgiqaNgTxpL3bciIiINpb+nRESkvvSbRSRV9CPZGFouTUSybOxyaRqzJ5Iq6sZtjJuAheGWbRqTJ9Vo5e+XVm57/S3MbRqzJ5I6+hGT+tKYPKlGK3+/tHLbk6SpV0RSRz+SIiJSX/rNIpIq6sYVkdai7lMRkaro7y8RaS3qPk03FWiIpI6CPRERqR+N2RNJnbbrxjWz35jZPWZ2p5kN5b32ETNzMzuoyLnnm9n94XZ+/Jrxf0VE2paqcUVSp11/xN7o7o/Hd5jZEcDJwP8UOsHMXgQsBroBB4bNbLW7P5l0Y0VEWoq6cUVSpe0yeyV8DlhEEMgVcipwu7s/EQZ4twNzw9cey/tXRKQ9KbMnkjrtGOw58EMzGzazhQBmdgbwsLvfVeK8KcBDsedbw324+wnxf6UNqUJURERSqh3/nprl7o+Y2SHA7WZ2L/Bx4JQy51mBfcWygO2xWoaMUoWotLfR8c8q0BBJnbb7kXT3R8J//2Bm3wFeD3QAd5kZwFRgvZmd6O6/j526FXhD7PlU4CdFbrOc+FqReU2oufEiIunUnXtkePv9ZhFJt7bqxjWzfczshdFjgmzeOnc/xN2PdPcjCYK64/MCPYDbgFPM7EAzOzA897YGNl9E2kWrDwuYUMMmIolpq2APOBQYNLO7gF8BN7v7rcUONrNuM7sGwN2fAD4JrAu3T4T7REpr9V/c5WT9/TVDKw8LUIGGSOqYu3oVG8nMHEaa3QxJwuyO1v4l3U70taqjTtw9N6a5+2/Mhz5S/VXsEobdvbv8kSJSrXbL7IkkR8FD69DXSkTaiII9kWZQ16dkmcbsiaSKRkqINIMyS5JVmnpFJHWU2RPJAmUKJS1UoCGSOgr2JHvaMfBRplDSwlA3rkjKKNiT7FHgI5F2DPybTZk9kdTRj1hjnB5ukjXtPIVHFEil+f2nuW3ZUmzFIBFJAc2z12CaZ6+NZT0wzPr7kyLy5tnrNB/6dPVXsbM1z55IUpTZE2mUrAdCWX9/UhlV44qkjsbsidRK48FEdqcCDZHUUbAnUqtmZbIUZNaXPs/6UoGGSOoo2BNpNdUEmQpkylP3s4hknII9kSxLQyCjgLP9KLMnkioK9kSyoBkBVaX3TEPAKY2T4Jg9MzvAzFaZ2b1mttnMXlP39otkkII9EWj97FMzAioFcVJIsmP2vgDc6u4vA44DNlfUJLOfmtl+4eN3m9kHzGzPiu8q0uIU7ImAAheRekko2AuDtdcBXwNw9+fc/akKW3WAu//ZzGYCFwEHAv9W8XsSaXEK9kTSoNUzi9LWzGyhmQ2Z2dBjT1FrN+5B0TXCbWHebTqBx4AVZnaHmV1jZvtU2MTnzWwicB5wlbsvBrrG8ZZFWoqCvcY4nWA5IS0pJIUDu6Qyi60SRLZKO6Ugd+929/Xuvv7gA2q+zOPhdaIt///LicDxwFfc/ZXAM8BlFV77i8BdwNuAm8J9+9bcUpEWo+XSGkzLpbUZLSGWTfq6xuQtlzbNfOi66q9iry69XJqZHQb8wt2PDJ/PAS5z97eWOOc14TluZvsCO919u5kdDXzc3RdU31KR1qPMnkiSsjbxcnTdds/EKdArLqExe+7+e+AhM/u7cNebgU1lTjsfGDazfmAesH94rQcU6Ek7UbAn0gylgqV6BFJJBSPRdWu9frsHic1Q7jOv99ck2Wrc9wHfNLO7gVcAny51sLu/292PB/oIijKuM7Ofm9mnzex1ZqaF2qQtqBu3wdSNK6mTRJdkM7s51cXaYHnduF3mQwPVX8VeXrobt17MbBLwRuA04DWNuKdIsymz1wqUDcmGNH4dkwqMmhlsKdBrrpSvjevu2939Fnd/nwI9aRcK9lqBfnllQxLZs/Gqtk1pDFhFCjCzk83s38zsFeHz/KlcRNqGgj2RVtVKq2YoSGwf6cnsvQf438A5ZvYmgjF+Im1JwZ5Io1US+NQjOEpTgKXsdHtJaG3cKj3m7k+5+0eAU4ATErmLSAtQsCcS14gAqZLApx7BUaFrpCkAlMq02tcsPZm9m6MH7n4ZcH0idxFpAQr2ROIqCbIa+cu33vdShq31tNrXLCXBnrt/L2/Xv9b/LiKtoYE1UG3t9HCrnqaRSJ9Gfj3qcS99D0nyRpc2i4K9FDGza4B3mNkzwCPA3cDd7v7F5rZMpDE0z16DaZ49EcmWvHn2XmE+dHv1V7FDkptnz8zuA2a4+/NmNgU4DjjW3a9M4n4iaaNu3FbVauN42lmjCjJEUsInVL8l7BcEK2jg7g+H8+wp0JO2oWCvVeV3yylYaL5iX4NGFWQ0kr7fpAg32Dmx+i1hy4GfmtlHzGyOme2f+B1FUkTBXla0WrCQZrUGMo34GjQ7yIrur++3bEji+ymdwd43gG8RjCZ8D/AzM3sw8buKpISCPZF8aQ5kyrWt0C/vev5Cr9dnU22bGhXkNjuYbrQEvtfdYMeEF1S9JWyruy929yvd/Wx37wKmJ31TkbRQsCeSJuMJguJVt9H+tFbiVtumRr2HNH5WLcbN2DlxYtVbwu40s0vGtNP92aRvKpIWqsZtMFXjSkFpDcpEyhpbjfvK7hf4T3+5R9VX2X/ic0lW434LeDmwP7AeuAu4091XJnE/kbRRZk+kkEZ354030OuvY3vT2sUqLWPnhAlVb0ly939092lAB3AF8GvgVYneVCRFFOyJFJJkli2/67Ue5tdp8mVIbxertATH2MmEqrckmVm3me3j7s+6+3p3/3q4Zq5IW1Cw166UjWmewS27B1aVfj2S/LrVUvwhkscxdjCh6i1h1wM7oydmdpCZvS3pm4qkhYK9xjidYJ6n5eUObBhlY8anmsAnXiwRyf/8K/16xAPFWttTq7R+zxR77wpOG2l5bGMnE6veEvZXd/9r9MTdHwc+kfRNRdJCBRoNpgINEcmWsQUax3VP9FuGDqj6KlPtj0kWaKwCvubuP4jt2+DuM5K4n0jaKLMnUm+1ZN5md9Dp20ofq0yVtIA0jtkD3g98xsxuNLP3m9lXAU2qLG1DwZ5IvRXq7ozvK1Q5O7iFkTldux9baL68UgUeszsqCiwr2idSgzQFe2b2GjMzd38EmAn8B3AwcDfwzkRuKpJCCvYkG9IYrBRrU6nK2fxzojF68QCw1Ni5wS3lx9aVC0arUevnnsavl9RNWoI94Hxg2Mz6gXOBte7+z+7+JXd/JqmbiqSNgj1pPYUChTQUDxQK1Co5Pvq3l8JZvfxjCxV6lMrwJdktXOvnnoavlyQiTdW47v5udz8e6AMOBK4zs5+b2afN7HVmlnj/sUgaqECjwVSgIWPkZ+3yu2t7CTKB/R1jM4K1zonXTFolJKPGFmh0de/t3xqaWvVVZtiDiRVoxJnZJOCNwGnAaxpxT5Fma7vMnpn9xszuMbM7zWwo3DcQPr8zfP3OIufONbP7zOwBM7sstv8nZnakmf2mQW9DGiHJrsZ4Nm92RxDM9ecFeoNbRgO8/K7feHdtdK3oGuUyes1SLNArlJUstF+kRmb2UzPbL3x6PvC3wIcV6Em7aLvMXhiQdYfzLBV6/bPAn9z9E3n7JxAssXMysBVYB5zt7pvM7CfABcBP3P3IMvdXZi8rkshUFcrsLcs7Jv+e+Vm/YtfMv7aybFIXu2f2/n3opVVf5Tj7dZJTr9zl7seZ2UyCuQC/Dxzp7ucncT+RtGm7zF4pZmbAPwL/XuDlE4EH3H3E3Z8D+oEzw9eeIJid/bGGNFTSoZZgKT9blZ+J6w0edvo2OtdsDAK9KIsXjenLz3zlB4Ol2jve7t9KxgI2Q5ra0ubSVI0b87yZTQTOA65y98VAV9I3FUmLdgz2HPihmQ2b2cK81+YAj7r7/QXOmwI8FHu+NdyHu7/D3R9y9xMSabG0jnJBR/74vGWFXxuxyYzM6WLRmiUs8usLBoS7raZRLAirtJI3rtj7iHcfVxswNnOpN2kYh9QUaMT8C3AX8DbgpnDfvknfVCQt2jHYmxVWZ50GvNfMXhd77WwKZ/UArMC+Yn3gC4GhIpuk2XgDkkJBR6nAKV6B2x8Ef52+LXfO0oHFLLXzWLRmye4ZvPwsXTT+LwoCC825V+kcfKWCp1o/o2quWWguQkmz2P9xlrrl0tz9euBVwAx3325mRwM/T/SmIilSdsyemR0CzAIOB7YDG4Ahd9+VfPOSZWZ9wF/c/f+FKf6HgZnuvrXAsa8B+tz91PD5RwHc/TNV3lNj9iQQH2sXD3Z6YbhnOjPt5tHjgM6ejYwMdI09J6rWzTs/Fxi2Wsar0FjCQmMOm6HZ90+tsWP2Xta9j39taFrVV5ltww2pxhVpR0Uze2b2RjO7DbiZIAv2EmA6cDlwj5ktiVU3tQQz28fMXhg9Bk4hCF4BTgLuLRTohdYBx5hZh5ntCcwHVifdZmkx1WS9ooAsCtqibRmcxcrceL5FPUvwqcaDa2eMntPfwfCa6bsHepFCkyvXkpGrV9freDKKtXYb11uz7y81M7M3mdnXzOyzZrbAzGaa2V7NbpdIo5Tqxn0LcJG7n+DuC939cnf/iLufARwH3EFQmdpKDgUGzewu4FfAze5+a/jafPK6cM3scDO7BcDddxD8Kr4N2Ax8y903Nqzl0hriXail5Ger5gdbZ0/wLTUy0EVnz0aG10xnqZ3HUbPCv0l6gf4OOns2Bpm//rwgKbpWoe7bSid5zn8/5Y6tpMu1UOBZatm2SoNMLf2WSiks0PgGQQXuL4BO4ApA/39L22i7qVeaTd24Um4qlE7fBsC5NiN3ysXhv1d7H0sHFgNBVy8EWcDcurpQfNxg/rx8jewaraUKuNg56k5NmbHduH/Xva//69CxVV/lTfbzJKde+W93f135I0WyqeSoWDM7FXg7QdWpA48A34tlw0Tay3gDjfj5sWBr0ZolLLXzwgBnI/NYlTulb3D09B4GoAcmWR/H3wr2gONHG9bruWlahv2tzJyzqfwaukm8v2rvF90zOqa/SIVy/udW6loKBJsqWi4tZX5qZh8EPu/KcEgbKhrsmdnnCWYZv55gmhGAqcD7zew0d7+kAe0TSZfxzK2XP44ur8vSBzs5atYGfmYzeK1vYJ4HQV3fbDjDp3H82s2sts3g0HdBeI1esHhR+OyOINALdfq2IOtXbPqV/PdTSVdvLVOuVBJcFpocOnq9XHav2HHSFElX19agC5gBXGpmw8CdwJ3uvrK5zRJpjJJj9tz9Le7e7+6D4dYPvJVgPJ+IlFNqTNnglmDiZIIpVmyrMzKni0MvCLpwexjInbLaNtM3O3g8yfp49Drouy7oyvVbLSjiuNKCKVrCaxecg69YmypVSzCV35ZS1wuPi7qy4/vGXKOSIFVj95oijZMqh3Oh/i3QASwG7gdenehNRVKkVLD3VzM7scD+E4C/JtQekeYZz2D/YvPZRdmmQqtXzO5gZCCYOHm4Zzo+1ehcszEI5C6A4xdsLnq7rwDbvS/I9F0YdO8eNWtD0BUcdYdC0B2cf/9q31s9FLp/f14AGsvqjdjk0fPyz610Dj5l+SRkZi8ys08CnyeYSux77v7hJjdLpGFKBXsXAF80s01m9sNw2wx8MXxNJFvGUwhQKNNUSfA4fwtL7TzOYiWXzupjJWdx6GCQtbt0RV/JW06yPvoGyWX8VnIW9HewyK8Pqnp7YemcoJiD/gJFIfkFG0kolY2LgrsoeOsPpqCJJpXu9G1jA7vo/EqXh5OmSGNmj2B5y6cJVs+YTDArQ6FkhkgmVTKp8mEEBRoGbHX33zeiYVmlatyUq3WAf6Xn5VXidvq2YKJk4Pfzjau9j0nWx8UE2btK9A3C+lnBmL6jZm0IsmL9HcHEzLHxezlNnpi4c83GXOZuzJjCcEqZyMhA1+jawAWus1vxRrFxfJKwsdW4R3fv70uHZlV9lX+wHyRZjXuPu7889rwTuNHd1ZUrbaFcNe7+wOuJVeOa2W3u/lQjGpchp4dbe2jlX7K1trtQ4UWpCYLD40bmBA8712zk0FvhqrV99FF5oAdBZm+793DWrHnMYxWrfB4rmR5MzBxeO9ctWuvEyvWoQA7vHUwTEwR3I3MI1/odbdfInK7RVUCiLun5BT7f/LZVW2wiY43v67w8euDhcmkp84SZHevudwO4+4hZ9EMhkn2lVtA4D1gPvIEg7b0P8EZgOHxNKncTwXq5C5vdkIbQL9kx4+ZKHtMAHOeGAAAgAElEQVQbPFzJWdhcp292MKdertq2An2DsIp5nGszcnPwxefeGxnoCn6R99fwy7zWyZjjYsFYVJASya3520suq9e5ZiMsAz86TBbFVxqJ3kOh4LqStkhx4/u5jf5/W5jSbtx3ATea2VfM7D1mtgx4MOmbiqRFqTF7HydYJ/Zid/8/4fZuoJtgyTQRgeKZpmglizIrUwyvCSZH9luNvkE4dBC4cHQi5VLO8GAN0pGBruDxMnLdwmMCq2IraxR7HwXaWVQVQcKITc4FuBBUIecqiAnavpKzWDxoXLqij0V+PX6ljVbn5q8lnB/4lapCzpKUv68UBnvbgeOBHwOHAHcBZyd9U5G0KBXsGcQn8MrZFb4mkj21/BItNHYs/lpvgXPCIHBRzxIG6GHmwCYuXdHH+lnTuHRWMLfeoYNhlm+wwPmMzr136ay+3YoWRuZ0BYHVYLgEW/6yaqXeR7yNlSp3bLwIIxZ0LupZwirmsahnCT+zGSyeb6y2zdzgG1g6sJilcxZz6aw+RuZ00blm42jQV6qtpeYRrFXagitlz6v1HXd/zt2/5e597v5vwCua3SiRRikV7H0KWB+mvT8Wbv9K0LX7qcY0T6TBSs0JV8m5ha6Rn+GLZf9WMQ8Yu/TZJOtju/fxaFhlG+/a3e5BBe4NviH3WnSNmQN5xRj9QbdnrtAhunf0b7l1aCsNKAotxVag6nhMuxgN9KLxhA/7NM7waWz3Pn5mM+js2ciiNcExUZZyzLJw/WXaX4vxfhaSW0Gj2i0JZvaPZnYl8EIzm2Zm8RstL3aeSNaUrMY1swOBU4lV4wK3ufuTjWle9qgat4XUq9CkUJVoOEZtZKCLRT1LWDpnMZ1rNrKSs3JB23DPdFbbZi4GDr0AHr0u+DeyfkXQhRs/fkz1bbzIodJ21RLgVfNalNmLBXwQduf2LOGqBX30XTeatQSwrT62C1oVtykzthr3yO4X+xVDc6u+yj/ZjXWvxjWzKcCbgauBdcDfAU8RLP15sLu/qp73E0mrSqZeOZRYNa67P9qIhmWVgr0WkXQQEQUu+VWny8hls1ZyVrByRjiXXlS0sX7FNAboAeBD1sdrfcNotW287VHmK7+StR4BXS3i3dyxYDcSZTejAPdq72PpwOJcEJurKo6mb4nOjQJBBX1Nsnuw9/Ght1Z9lYV2Q5JTr8xy97Xh4xcRlH/f6+7PJHE/kbQpVY37CjP7BfAT4Crg/xIsJv0LMzu+Qe0TaY5aA4di3YD9HWO7TuOFBb3heWHQt3TOYkYGupg5ZxPbPRjHt6Tf6bsOuBBm2s2sYh5Xre3jMEaC7KBfv3v1b9h1G01SXLRt1XZdVlvkEd07r3t7JWfB/GBM4aKeJblM3g2+gau9jx4G6OzZyEy7meE108cEtCM2OcgKRt3TCvRSI6XVuJvN7JNm9lXgXGBEgZ60k1Jj9q4DLnH3ae5+sruf5O4vAz4ArGhI60SSUusyaOVEAU2hMWuFlv6aPzZI6fRt0BtmuXrJTaXC/C1s9z4A/IJO5rGK9bOmjQaJEGQF/fogeIxl9UYGukpPV1Js0uJaxq8Vm1swby68RT1Lgu7n2R2M2ORgmbdrgjGLDy6YAcTGIPZ35OYMjGfwls5Z3JhVQCQLtIKGtLVSwd4+7v7L/J3u/guCOfdE0qmSX/zVLINWS8AXD0LihQR5jxf59cHz3mCLiikG6BmT8fILOgG4dFYfdl0wDGCAHvxKg97RIo2ldl4Q/MXH6sUrdQsVkRSbcqWa7t5KjwuPXWrnBYFu7+hr9oAzYpM5asUGls4Jum+jrtoH184Y/axg9zkDFfSlRkozey9x96Xu/n13/wzBJPf/kvRNRdKiVLD3AzO72cx6zOy14dZjZjcDtzaqgSJVq3eX3niuN7glCLbi4+ei6Vjmbwkyd8sY2yUZBm9R5en6WdOwuc4q5tHDAPR3sHRgMT0MYJd5MLde1MUZri8LjB0XWGoOvfF+XsUyg/HH8bZE7Qynnsm1OZ6RHNwSjEsMPxObPTK6zi+MfqbxALoe70XqIi3VuDFPmNmx0RN3HyHI8Im0hXLVuKcBZzK2Gne1u9/SmOZlRm65NDO7SAUabSI+Pq9Yt2m4Pyo+iLJXue5bgjVzD5vtDK+ZzsyBTUHBwsCm3IoTIwNdu0+YHN4zXtSwW/ZrPFWt5Y4vtHZtJF6RmzfmLvf+41200XFRoDhfAV26dBLOWwfAEd2HXPTBobOqvsqH7ctJFmi8DFgFrAHuAbqAw93975O4n0jalK3GlfpSNW4bi1fJLmNMIBYFchBk+XoYYObAptzjAXpYOrCYzp6NzGNVbqqWMfPODW6h07fluj6jKV1yFb+x48a0KX9foTbX8j7j4uvb5gd7hSZIzg9Mw0yoqm7TaGw17tTuQ/2SoZ6qr7LIvphYsAdgZnsCbwemAw8DN7j7X5O6n0ialKrGnWBm7wormF6b95qWSxOJFJqguFA3ZhTohV250Zq1i9YsYeacTbmpR65a25fL4EVB31Vr++js2ci5NoOrFvRBb1jN2hsst7ZoTdAdOo9VwT3nbwnGxUVB0mBsiyu0wsd4x70VCsbiweb8LWODP9i9Szb2mUVr5RZUyeTQWdFC7y9tY/bM7E3AV4BXAb8lWBxAmQ5pG0Uze2Z2DcGYhl8RlKr/1N0/FL623t01/UoNlNmrQTOyOZXes9Bx5bJloXgWzm81bK7nVpWIJlf2W3MJE45asYEH1waVqlwTPIcg6Js5Z9Nod3F+t3EU8NX6HitVqus2/ziobGqX/LaX6oJW1q9Jds/s9Q69s+qrfNQ+n2Q37iPAe4GJwLHh1uXuRydxP5G0KVWgcaK7v9PdP0/w19C+ZvZtM9sLrY0rjZTmX+CFqnULTbFS4HFurN38LdgDwR9d0Vi91bY5CPQuDE559LqgItVmj8A1wb6ROV2s5KzcBMtjrh+2oXPNxvJj3OqVMYq/53Jfs3Jr2xZre7HsZP41WygLljVpWi4t5gF3/467r3T3f3b3MxXoSTspFeztGT1w9x3uvhC4E/gvYN+kGyYJa+Vfho1oe7kK03LHRvsLFWlEx4fdmYv8ejrXbMxNtfIzm0HfIHBhsHIGBP1P62dN4/d0woVgc4OCDQgqdxetWZK7vl9pwXx9hGvJ5k+2XCg4rUStU9rEu1mLBcDxjGAUDFuZYslW/h7OsGDqlYlVb5UKhxjdYWbfr+DY683sA8DPzezD43lfIq2sVLA3ZGZjFjh0908QTKh8ZJKNkgZIc7askHKTADdCrXPuQfHu3mXBnHMPrp3BPFbRwwAP+7RgwuQwg9d3HVxMkO172KfRNxt8qrHaNrPaNvPg2hmsYl4u+LPLfGzhRpgdiwLAcWfCav0M8ucfzL9//nXHM/F1NcG61F3CY/YuATZXeOzXCX7PHQaca2b/Y2arw7Ho1ZcMi7SoosGeu5/j7rvNp+fu17j7Hsk2SyRPvQO8Wn/xF2tHfvas3D3DaVEitjXoxo1Wijh+weZcF+7FBJk9CAI+GM343eDBuL0H185ggJ4gw9fLaDFIrL25TFkUWFUbQNdzXr78SZ2jfVHb8+9T7OtVydjA+H0lcUlOqmxmU4G3kvtTqExb3H/k7le7+/nu/gqgE7gc+DXB8CSRtlA0d25m7yhx3rMEawtW+teVSLrUuzChWAFEofuFAcqITWbYp+emVVnFPB5cO4NHZwMXjB7+FYo712bQB/RdAD2zBjiLlfitxlG+gXksYVXPvN27Q5MMekoFX/lTrRQqsihUjFHo2PhULMUouGsYM1sILATY72/2r7W69iAzG4o9X+7uy/OO+TywCHhhle07DVgMHAjcBXzO3W+opZEirahUN+7pJbb5wA1mpuVmRCqt2o0fHz6faTcH+8KpRdbPmsbV3kffdaPZu0qsXzGN4xdsZh6rsLnOyEBXsApHWO1Lf2xak0JZyEoyneNdhq7QnHqFzssf21hoUur5RV6TpnD35e7e7e7dkw+ueWGKx6NrhNuYQM/M3gb8wd2Ha7j2l4EPA68GlgP/18zOrrWhIq2m5kmVzewFwD3u3lX2YMnR1CttoNAUJPnLhcVEK2IM90zn+LWbsa3O4vmVF7yf4UGgd+mKPiZZX65rF2IVv8XaV8n+pFTSjkqre8tdUxI0duqVw7qn+DlDF1d9lc/aP5ecesXMPkMwDdgOYG9gP+Db7n5OuWub2S/c/dWx5/sAv3T3GVU3VKQFlZpn7xzgRnffVeT1owiWm1mTYPuyQsulydhApL8Dn2rY7BF8sBPb6vitFmT0LgimWinVfVtI32Aw9q+zZyMPLpiBzfXaV85ICwVvLWDscmmHdk+96J1DBf6qKePz9tGK59kzszcAH3H3t1V4fD/wAPAJd38uXE3jp+7+mqobKtKCStW7vxi4w8yGgWHgMYK/po4GXg88DlyWeAuz4aZwA7iomQ0ZF/3irU28ACF6TlA1S3/QdbvYDC4IArZHZ1ce6EX5k4d9GizYzPCKMDt43QiOYfhoNnGQ4kuUlWp3uX31lL/SSCVt0/dlGiyMPb4o6RUxauDAO4CLzOx+4G+Ab5rZMe5+f3ObJpK8UtW4XwCOB/4dOBh4c/j8YeBcd/8H/ZC0Gf1CrU7+HHvRkmCDW8YEXDPnbAqCtnBevUMHg6Bvu/eVvcWhg+SuAQSB3mUO/R1cuqJvtJAhyvDlLztWbN67qN358rtX6y0+B2GxNuTvT2LcnsYB1izJatzcPdx/UmlWLzz+bHefDrwU+ACwBNgHuMbMHqrq5iItqORMlu6+E7g93ESkGvkBySC7T7/SsxF64CsGN8zawLnMgGuCSZN9rdEXHhqffiXfoRfA8IrpHEU4Vq8XOns2BmvjEguEomxZfK3e/DamVbnsXb3bnubPIuWiFTTSyN3/CgyFm0jbqLlAQ2qjAo02UWjSYMjti6+LC+Qybn60BZW1azcX7c7d7kEhxhk+DYAptpnD+kd/jjt7NgaTKveOXrfkJMOlAptmdpG2U/dsS7/XsQUaB3f/jb9jqPrFKpbbBxJbG1ek3ZWaekUk3eKT8qZNsS7QsDt3xCaP6VKNJkNev2IaA/Rw6aw+Xhurqu0bDKpuLwauWtAXLKcGzBzYFGT2eqYH11s2OnlyZ8/G0YmKi60yUWri4kJVsUlLw0opzZCx95rkcmkiUj0Fe9K6iq3o0Ozgr9A4uHhbl5Fb3SJaD3fpwGI6ezYyxTazdGAxEE6YfEEQ6K2fFUyv8lrfwFErNnDUrA0M0INPDTKBZ7EyWJEjV4yxZXQy5TCozC2XVqzN+Z9jI7t461Ed3OyvuwCNGbMnItUpG+yZ2afN7IDY8wPN7P8k2yyRcWj2uqjlqkSjgG92ByMDXeHYOhiZ08Vhs4PpUnoYYMlsp+86coHd+hVBt+1KzmIlwbKeNnuEmQObGBnoGpstzH+/veTWyy0Y9DVindykV7vIWHasVaUx2DOz6QX2vSHRm4qkSCWZvdPc/anoibs/CbwluSaJJKDaQKAewWGhLtBYBmvRmiVAGHzFp0PpJZepWxKO77xqbR8D9DAy0MUAPRy/dnPuXJ9qY7KFY+7VH1taLDRik2sPvCqp2C12bQVjbSNtwR7wLTO71AKTzOyLwGeSvqlIWlQS7E0ws72iJ2Y2CdirxPEira9cYFJNMFhompD+DpbOWcyiniWMzOli0ZolQdAXFlVEmTof7GRkThdHzdpADwNBle3AYmyr08MAKzmLo2aNju3LVdpGllF+Ddlq3lMlAZuCurYWVeNWuyXsVcARwM+AdcAjwKykbyqSFpUEe98AfmRm/2Rm/4tgGpavJ9sskRQbb+VktLZrLyyds5jONcE0KWOWNpsfTCx81KwN0Bt03UaZvagQY6bdzEy7Odc9CwTTrqzZOPZepZYlK/ee8oPHNEpT+9LUFol7HtgOTCJYHGBLsdWhRLKobAmUuy81s7uBkwADPunutyXesmzJLZfWllp6WokCxvteos8jDOji1bMjy4IpU4bXvJUBenLFGjMHNgXVtZFe6FyzjXmsAoKgMXeNaNoVAIpULMffQ7SqRrkijbRKUxvT1JbGWh49CMbspa66dh3wPeAEgtWhvmpm89x9XnObJdIYFc2zZ2YvBY5x9/80s8nABHd/OvHWZZDm2Wtj+VOZRGPposeRZUAvLOpZwirmjc349XeMCfqijN9SO6/wdSKlMnwi4zJ2nr0Duzv9DUOfrvoq37WzE5tnz8y63X0ob9+57n5DEvcTSZtKqnEvAlYBXw13TQG+m2SjRFpOpePd4sflj6VbxmggN39Lrmu3s2djEOSFVbQjNpl5rMoVawD4BZ1j7xUFer3BfYtOu1Kq3WnrkkyyPUldO22fYQOksRoXeIuZXRHfyKW9RbKvklz7e4ETgV8CuPv9ZnZIoq0SaTXVZszyjw+rcUe7YDty+0fmTwa2MEKQ4Vvk1wfdu/O3sJSgyMPmBhl6n2ocNWtDcE5s0umROUXa2MglyMYryfYkde20fYYNkNLl0p6JPd4beBuwuUltEWm4Sgo0nnX356InZjYR0Bpr0jitlh0ptBJF9G+U3Qv35Yop8rtZ88fUReva9nfk5uXr9G1BwDcwOl7v0ll9rOQshv2tuTVyd2tLf4FxfJV+xvU+TqQB3P2zse1TwBsIeqlE2kIlwd5PzexjwCQzOxlYCdyUbLNEYuqZHWlEEFKsyCH+b/h4ZE5XEPBF4+1md4wZa9e5ZiPM35LrygVY5NcDBMUZ/aPHj8zpYunA4tF5+JaRK+BgcMvoPHzzCwSUxT7j/CXpKl0+rQ0zWjKqBZZLmwx0lj1KJCMqCfYuAx4D7gHeBdwCXJ5ko6QNNCvz04ggpNx7y1v/NbesWRiE5YK/6LXZHUH3bjhnXg8DDPdMDzJ84RQui3qCCZp/Pz8YJ3/prD5+P2hctbYPv9Jgdkewfm48swfB41JrDOePM4zvFykgjWP2zOweM7s73DYC9wFfSPSmIilSydQru8zsu8B33f2xBrRJ2kEjg4VCVahJVKZWes0y4+SCAC/v+FhANnPOpvAxuYrepSxm0ZolXE0fq5jHgwtmsN6n8VpWMo9VLB40jl+b11YYWyRSrF3NCOxUOTyqxT6LKNhLmbfFHu8AHnX3Hc1qjEijFc3shcvK9JnZ48C9wH1m9lhYxdQ0Zra3mf3KzO4ys41mtiTc/00zu8/MNpjZtWa2R7j/AjPrCx/vZWYDZvaAmf3SzI6MXfej4f77zOzUIvfuCM+7P7zOnuH+vvA+12m9xYSMJxNYbWHCeO5Tyy/mQuf05r0e3xfrDl7UsySXqVs6sJirFvTx4NoZcCEcv2Az59oMlg4sZrv3jV5v2e7XKdm2+L+Vvp/xaqHgJnEt+FmkbQUNd/9tbHtYgZ60m1LduB8gWE7mBHd/sbu/iGDJmVlm9sGGtK6wZ4E3uftxwCuAuWb2auCbwMuAlxPMkn5hgXP/CXjS3Y8GPgdcBblFsucDXcBc4MtmVuh/n6uAz7n7McCT4fWkEZL+hVevbuVa2lnonNi8emOOiQeUYbFGZ8/GXCHG+hXTAOibDX3XwcWMztdnl3nQRRwFjVGhyOy8rtzxLq3WgsFJQ2W8eCWaVDkNY/bM7Gkz+3O47fY4kZuKpFCpYO884Gx3z/3P7e4jwDnha03hgb+ET/cIN3f3W8LXHPgVMDU8ZjsQHX8mo0u9rQLebGYW7u9392fD9/sAwXQzOeFxbwrPI7zO28PHfwnv8yfgOaT1NKurstzzcEweMLrKBWOfj8zpCrJ5ob7ZcINv4OLw+SQLKnRz4gUa4RYfJ1hR4YYCutrps2ukGe6+X7i9MP9xsxsn0iil/pzaw90fz9/p7o9FXaTNEmbdhoGjgS+5+y9jr+0BnAtcAuDuA7FTpwAPhft3mNmfCJbOmQL8InbcVnYvy38x8FQs/Z87xt3/X7gvutfCcBMprtKlyeJLmcUycH6lYVuDWZDsMs8Fhds9CO6u9j4mWR/bvY+ZA4uD8YADXbmxfp1rtuXm9RuxLUADxjVKVuVWp0jZmL3vAMcDmNl/uPs/NLk9Ik1RKrNXKkPV1OyVu+9091cQZO9ONLMZsZe/DPy3u68pcKoV2Ocl9ldybiHLge4iW3ZkvDuqaoU+j3LZu1IKjamLZeJsqweBYLi82nDPdCDI5EGwXm7fYHheOEdfZ8/G3Nx+IwNdo0u2lbp//uN60fdPloz5Py5F1bjx/7c11Yq0rVLB3nGx8Q3x7WmCcXFN5+5PAT8hGGeHmS0GDgY+VOSUrcAR4bETgf2BJ+L7Q1OBR/LOfRw4IDyv2DHtRVmfsSopBKn1MwvP6/Rto1Ox5L02QA9nsZIl/c4U28yiNUu4dFYfSwcWB3Pzxebji6ZxAcoXYSQVlOn7J5NSNvWKF3ks0laKBnvuPiE2vmG/vLEOTevGNbODzeyA8PEk4CTgXjO7EDiVYJzhriKnrwbODx/PA/4rHOO3GpgfVut2AMcQjPvLCY/7cXge4XW+V793JplWKMOXv4JFBfPzjdjkMZW6ubF2kFtJY7hnOq/1DazKfauSezwy0DW6akc0IXN+8Ud/XjvSMJ6x3sdLYpxUVeMeF0tSHKsCDWlXFsQwrcPMjiUojphAEKx+y90/YWY7gN8CT4eHftvdP5F37t7ADcArCTJ688OiE8zs48D/IpiD6QPu/oNw/y3Ahe7+iJl1Av3Ai4A7gHPc/dkq2+8wUsM7r5DGWaVHvb8W0SobkbBrNjeVSryQIwrcesl18+bOqWXeQX1fCRT5PujE3XPdpZO6p3vH0DervvRmO37Y3bM11EUkJVou2Gt1iQd7ki3xX66x9XGBsVm52R10rtkYZP7imbn5eedGAWCxcXrNpqCyBSnYE0m7SpZLk3ahrrDGyO++LSW+Fm2UwZu/ZTTQiwV2QaVtbGWM+VuCMX5RYYZvG30trV/r8QR6aX1PbSZlY/ZEBAV7EqeMSv1UOulwpZ95vEAjeh4PAAe3jD5eNnpctO7uyJyuoDAjL0As2+ZWCqD0/ZsaCvZE0kXBnkgSigUe41z2bUwVbv59Co3DiwLC+GvLKKzUOL60BX1pa4/kOJamAg0RQcGeVEu/ZCtX6LMqFFBVc71SwVfs9RGbHDzujd0z/u9gma7c/o7aMpDxtiSpXbN4LfDzl6bl0iLhTAvvNLOPmdkV0ZboTUVSJNmfMMmedv0lm5RKP8+osCKaJLnYnH7FChwK7S9172LZv0rp+yQZLfK5prBb9nsEy1kOE6yvLtJWFOw1xunhJu2knr+YB7cAFWR14gEhjF0HNxJlh0q1r0WCCkmN5c1uQBlT3X1usxsh0izqxm2Mm9B6uTJe+cUV+V16+YUX8ZUx4sdW0o3bSGlpR1Ky/v4C0f9vC1NajfszM0vFyk8izaBgT6QRqplupZTYtCm5lTCia+Yvf1ZuzF2hbF8zZD2LmPX3l8cxdu6aUPWWsNnAsJndZ2Z3m9k9ZnZ30jcVSQsFe9IaSgUjjQxUKljSrOC/4yl2yL9+vAgj//rzt1R3/f7Y+e2RgZKkOezYMaHqLWGnESyDeQrBkJq3oaE10ka0gkaDaQWNAtp91YRK339+NW6p6V1q/Tzb/WshNRi7gsaE41/hkwf/s+qr/GWfgxNdQcPMjgPmhE/XuPtdSd1LJG2U2ZPma7fgIj+DVuXEygXPyR+TV20bqm1LvSmrmBnuxs4dE6rekmRmlwDfBA4Jt2+Y2fsSvalIiijYE6mHaoKV8QZUtczfV2uA2Shpa49kzT8Br3L3K9z9CuDVwEVNbpNIwyjYE6mHRk46XMm9qplTL82U8Ws9Tuoye4ABO2PPd4b7RNqC5tkTqZdqxru1avBVTFJj/bL2OY1Hi4yndDd2PJ+6SZVXAL80s++Ez98OXNvE9og0lDJ70l4qraatRSN/ETcz41VoGpkWCEJaXst8xsaunROr3pLk7lcDC4AngCeBBe7+uURvKpIiCvakvZT7hdkqv1Ab0c5icwNWMo2Mul/blwM7JlS/JcjMrnL39e7+L+7+BXe/w8yuSvSmIimibtzG0HJp0nrKBXWluhVbJWiWehldLs0t8eCtBicDl+btO63APpFMUrDXGDeFG6gCTESyJ74UZGr+jzOzi4H3AJ2xFTMM2BdY27SGiTSYgj1pLy0yyL0lFJruRZ+tOLAjNYWuNwI/AD4DXBbb/7S7P9GcJok0nsbsSWuqZkxYtRMOS2302UpkRw1bAtz9T+7+G3c/GziA0SE1RyRzR5F0UrAn9TWegflJTUysIKQwFVE0T5Y/eyc1wV7EzN6PVtCQNqa1cRtMa+OKSLaMXRvXpnc73xiq/jIzLbG1ccPxeq9x92fC5/sAP3f3Y5O4n0jaKLMnIuOX5UyVZIFW0JC2pgINkVaVpoKItLRDms+B55vdiN3EV9Aw4Ey0goa0EQV7km5pCmjSRp+LpJEzNoeWAu5+tZn9BJgd7lrg7nc0sUkiDaVuXEm3QtN7ZFFW35e0p/QVaOwFvIxgfr0DgNPN7Ipk7yqSHsrsSWvJajYrq+9L2k9UjZsu3wP+BAwDzza5LSINp2CvMdpzuTR1wY6fPkNpDbHl0khjsDfV3ec2uxEizaJu3Ma4iWA5oYXlDsyUdg5S6tUt286foYxPY4cGLCTd/8f9zMxe3uxGiDSLMnsitSqVdVOQJs3WrO/BhDJ7ZnYEcD1wGLALWO7uXyhzzj1hiyYCC8xshKAb1wDXPHvSLhTsidRKAZ3I7pLrxt0BfNjd15vZC4FhM7vd3TeVOOdtibREpMUo2JPs07g3kcZKINhz998BvwsfP21mm4EpQKlg71Hg3cDRwD3A19w9fSMKRRKmMXuSfWkP9DTtijRCo77PokmVq92qYGZHAq8Eflnm0K8D3QSB3mnAZ6u7k0g2KNiTbGjlgCntwWiraeXvhSQl+H1mZgvNbMjMhvjzY8Gkyn5I7xAAABXFSURBVNVucFB0jXArWOxhZvsC/wF8wN3/XKZp0939HHf/KjAPmDPe9yrSitSNK9lQ7heZunLbh77ODefuywmnX7HObq/xMo+7e3epA8xsD4JA75vu/u0KrpnLGbr7DjMthyvtScGetAcFAK1HAXprSq4a14CvAZvd/eoKTzvOzKLsnwGTwudRNe5+9W+pSPoo2BOJU4CRHvo6tKbkqnFnAecC95jZneG+j7n7LUWb4j4hkZaItBgFeyJx6g4WGZ+Egj13HyTIyIlIlRTsNUZ7LpeWRQr0RApZPuaZJjcRSRUFe41xU7gBXNTMhohIgto38ztaOetcpGBPJF009YpIFmi6kUCzP4f2DPREJOUU7Im0kmLBjIKMgD6H5ovG7FW7iUhiFOyJVKuZ2SMFM7VrdtavXTRgBQ0RqY7G7IlUSwFXa9LXrTGcaEUMEUkJBXsiIlJf6pYVSRUFeyIiUj/JTaosIjVquzF7ZjbXzO4zswfM7LJw39fM7C4zu9vMVoULbRc696PhefeZ2amx/b8xsyPN7CcNehsiIiIiFWmrYM/MJgBfAk4DpgNnm9l04IPufpy7Hwv8D9Bb4NzpwHygC5gLfDm8nqRJuw3Cb7f3K+mnalyR1GmrYA84EXjA3Ufc/TmgHzjT3f8MuYW2JxH8d5XvTKDf3Z919y3AA+H1AB4jGJL8RNJvIBOSDFCaOQi/GYGXig4kbVSNK5I67TZmbwrwUOz5VuBVAGa2AngLsAn4cJFzf5F37hQAdz8h3PeO8N+FxGeUl7GyGqC0wvtK8woPaW6blDOUe6RqXJHUabdgr9Ai2g7g7gvCbtkvAj3AikrPLWA5+WtFlj9HsiLNQUta2wXpbpuU0x177OqWFUmXduvG3QocEXs+FXgkeuLuO4EB4B+qPVckR0GLiIikSLsFe+uAY8ysw8z2JCi4WG1mR0NuzN7pwL0Fzl0NzDezvcysAzgG+FWD2i0i0hpUoCGSOm3VjevuO8ysF7gNmABcC2wG1pjZfgRdtXcBFwOY2RlAt7tf4e4bzexbBGP6dgDvDTOBIiISiQo0RCQ1zF1DyBrJzBxGmt0MEZE66cTdc2Oa7aBu54yhUicUtsKG3b27/IEiUq1268YVaT+ai08aSd24IqmjYE8kbeodnLVLwUi1n5uC4OQo2BNJFQV7ImnTLsFZPdUy3Y0+ZxFpEwr2RKQ56plZU+CWHlpBQyR1FOyJjFcUtLRqt2Cz2q0ALZuiFTSq3UQkMW019UoTnR5uUomoS65Y11zaVqiI2tIKbS0k7e2TVjC6YlBUoCEiqaGpVxpMU6+ISLbkTb2yf7czq4apV36gqVdEkqLMnoiI1I8mVRZJHY3Zk+Zo1jixVh1XJyIiUiMFe9IczRonpvFp9acAWvKpQEMkVRTsiaRNqwVPCqAlTitoiKSOxuyJpI2CJ2llqsYVSR1l9qQ6rZZ1ajQt2SXtTpMqi6SOgj2pThJZp3IBTysFRFqyS9qdJlUWSR0Fe9J85QIeBUQiIiI105g9ERGpL43ZE0kVBXuNoeXS4lphCTERqYaWSxNJMS2X1mBaLk1EsiVvubS9u50jalgu7QEtlyaSFI3ZExGpRT0Kh1qp+KhSKtAQSR1144qI1KIeQxGyOJxB3bgiqaPMnqRLFjMdkjx934iIFKVgT5qn0C/oLGY6JHn6vkkXLZcmkirqxpXm0S9okeyJVtAQkdRQsCciIvUTFWiISGqoG1ekVhonJrK7qEBD3bgiqaFgT6RW6oYWEZEWoGBPRALKVEo9KLMnkjoas9cY7btcmpZGax36Okntxi6XpgINkVTRcmkN1jbLpSnIq5w+K2lpeculWbdjNSyX5louTSQpyuxJMhS8VE6flWSNcggiqaIxeyKNojFxIiLSBAr2RBpFGTwREWkCBXsiIiIiGaZgT6QW6pIVEZEWoQINkVqoS1akCM29IpI2yuxJ+1FWTiRBmlVZJG2U2ZP2o6ycSIKU2RNJG2X2JF2UdRMREakrZfYao32XS6uWsm4irWh0ubRcN279mdlc4AvABOAad78ykRuJZIyWS2uwtlkurZxoiTAtFSbS4vKXSzvO4Yc1XOewksulmdkE4NfAycBWYB1wtrtvquFmIm1FmT1pjijAU6AnkjGJjdk7EXjA3UcAzKwfOBNQsCdShoI9ERGps0S6cacAD8WebwVelcSNRLJGwZ6IiIyLmS0EFgbPDqfGzN5BZjYUe77c3WNjAbH8EwjSiCJShoI9EREZlzAoWw5g9vJaA7DHS43ZI8jkHRF7PhV4pMZ7ibQVTb0iUm+aPkbaWmKTKq8DjjGzDjPbE5gPrK5360WySJk9kXpT0Ym0tWQKNNx9h5n1ArcRTL1yrbtvrPuNRDKo7TJ7ZjbXzO4zswfM7LJwn5nZp8zs12a22czeX+Tc883s/nA7P7b/N/F/RUTaV3LLpbn7Le7+t+5+lLt/qv5tF8mmtsrshfM0fYnYPE1mtpqgousI4GXuvsvMDilw7ouAxUA3wf9mw2a22t2fbNgbEBFJPS2XJpI27ZbZy83T5O7PAdE8TRcDn3D3XQDu/ocC554K3O7uT4QB3u3A3PC1x/L+FZFiNKYx45LL7IlIbdoqs0fxeZqOAnrM7O8JArb3u/v9FZw7BcDdT4j/SzAFwcK6t14kCzSmMYuGyh8iIs3SbsFesXma9gL+6u7dZvYO4FpgToXnFpKbhqCKc0REWlVsyhR3deOKpEu7deMWm6dpK/Af4b7vAMdWca6IiOSoG1ckbdot2Cs2T9N3gTeFx7yeYLHtfLcBp5jZgWZ2IHBKuE9EqqVxexkWFWhUu4lIUtqqG7fYPE1mdiXwTTP7IPAX4EIAM+sG3u3uF7r7E2b2SYKAEYKCjiea8DZEWp/G7WWcMnUiaWLuGkLWSGbmMNLsZoiI1Ekn7p4b02x2jMPVNVznjOEyy6WJSI3aKrMnIiJJ0zx7ImmjYE9EROpIwZ5I2ijYExGROoqqcUUkLRTsiYhIHSmzJ5I2CvZERKSOlNkTSRsFe41xeriJiGRRsRWDRCQFNPVKg2nqFRHJlvypV17q/P/t3XuQJWV9xvHvE1YQwYis0Vp2CYu60VALrCSVwigokIpALFYSUVaNmGCsVKAgIZRgaShSqVSFXMREEywExRjLSwgoEOVSSIJRIWBYliUILEJwYcNFgSByEfjlj35n9zB1ZnZ2dmbPTPP9VHXNOW/3292/7pk5z+k+fZpTpzGfP/CrV6RZ4pE9SdIM8jSuNNcY9iRJM8gLNKS5xrAnSZpBHtmT5hrDniRpBnlkT5prfmbUKyBJkqTZ49W421h3Na4k9cdzr8bNpcDLpjGbB6vq0JlbK0kbVZXDaIaz58C018+jdZ3qtFOtaS6sa1/rmo3fwbmwrqOua9T7dUvrcnBwmCODp3ElSZJ6zLAnSZLUY4a90bl4Dkw76uVbVz/rmo2atnS+s7Guo65r1PtV0jzlBRrPb9cDfbs9UR9rAuuab6xL0pzhkT1JkqQeM+w9v5096hWYBX2sCaxrvrEuSXOGYW+eS3JXkpuSrE5yfWv7qyTfS7ImyYVJdpmg791Jbk2yLsmpA+17Jrk2ye1JvpRk+21VT1v+sJqOSnJzkmeTTHgaKclObbq1Sb6Q5IWtfaQ1tXXYmrq+nOT8tl9vSfL61r5rkitaXVckeem2qmdg3aZdF3B2ku2S3JDkkoF5ztv9lWT3JKvafro5yYkD46a8v5J8Osn9SdYOtA1dfpL3JfnEBPO5NMmNrd8nk2w31VqG2Bj2kpye5J62fVYnOby1L0xyVZIfT7ROk22LdP6u/V9ak2S/Ka6bpAkY9vrhoKpaUVVj/7CvAJZX1T7AbcCHxndo//D/HjgM2AtYlWSvNvoM4MyqWgY8BBw72wUMMb6mtcBvAldP1CHJYuAE4JerajmwHXB0Gz0XaoJp1NX8LXBpVb0W2Be4pbWfClzZ6rqyPR+F6dYFcCKb6hkzn/fX08AfV9UvAvsDxw38bW3J/joPGP8lw1uyXce8o6r2BZYDPwcctRXzGu/Mtn1WVNXXWtsTwJ8AJ2+m70Tb4jBgWRs+AJy1FesnCcNeL1XV5VU1difya4AlQyb7FWBdVX2/qp4CvgisTBLgYOD8Nt1ngbfBxiMBa9tRgq15gdhiVXVLVd06hUkXADsmWQC8CLh3rtYEU6sryc8CBwLntj5PVdXDbfRKunrguXW9aeCIyw1JXjw7FQw31f2VZAnwG8A5A23zen9V1Yaq+q/2+FG6ILu4jZ7y/qqqq4EfbcHyd2tH8W5P8pcDff6vPVwAbE9389oJ59WOEn4lycVJ7kxyfJKT2npdk2TXzdT/WFX9B13om8zQbdHa/7E61wC7JFnUhqvbNlqb5IDNzF9SY9ib/wq4PMl3k3xgyPjfBb4+pH0x8IOB5+tb20Lg4YGwONYOcBrwlnaU4IiZWPkJbK6m4Z2q7gH+Grgb2AA8UlWXMzdqgmnWBbwSeAD4THvBPSfJTm3cK6pqA3QhA3h5az8ZOK6qVgAHAI/PTAlDTbcugI8BHwSeHWib7/troyRLgdcB17am2dxfK4B3AnsD70yy+8B6XAbcDzzKphA9meXAu+jeFP458JOqeh3wHeC9A9Md3061fnqyU9ITmGhbTPS/6V3AZW0b7Qus3sLlSc9bhr357w1VtR/dqY/jkhw4NiLJh+lOKX1+SL8MaatJ2gG+BZyX5PfoTpHOlglrmkx7sVkJ7AnsBuyU5D3MjZpgmnXRHZHZDzirveA+xuZP134L+GiSE4BdBoLTbJju/norcH9VfXf8qCGTz6f9BUCSnYF/Af5w4OjaRGZif11ZVY9U1RPAfwN7jI2oqrcAi4Ad6I6abs5VVfVoVT0APMKm7+K7CVjaHp8FvIouZG4A/mYa6zzMRPv/OuB3kpwO7N2OmkqaAsPePFdV97af9wMX0r0TJ8kxwFuBd9fwL1NcD+w+8HwJcC/wIN1pkwXj2qmq3wc+0vqtTrJwxgti4pqm4NeAO6vqgar6KXAB8KvMgZrasqZb13pgfVWNHR06ny78AdyXZBFA+3l/W8ZfAO8HdgSuSfLaGSliiK2o6w3AEUnuovsYwcFJ/on5v79I8gK6oPf5qrpgYNRs7q8nBx4/Q/cmYaMWAi+ie0O0JfN6duD5s2Pzrar7quqZqnoW+BRbsH2aoduCCf43tdPaBwL3AJ9LMniEUdIkDHvzWLorT1889hj4dWBtkkOBU4AjquonE3S/DliW7qrH7ekuZLioBcOrgLe36Y4BvtqW8aqquraqTqN7Qd59yHxnpaYpdr8b2D/Ji9rnvg4Bbhl1TW05066rqv4X+EGS17SmQ+iO3ED34n1Mezy+rpuq6gy6L8KdlbC3lXV9qKqWVNVSut+/b1TVe+b7/mq/e+fS/e59dNzobbq/kuw8EKgWAIcD35uheS8aeHokU/87HTN0W7T296azP93HMTYk2YPuSPCn6LavV+lKU1VVDvN0oPss141tuBn4cGtfR/eZl9Vt+GRr3w342kD/w+mu1r1jrO/AfP+zzeefgR1a+wV0p3HW0l0dmm1Y05F07/ifBO6j++zOsJr+lO7FbC3wuYF1H1lNM1TXCroQsAb4CvDS1r6Q7krG29vPXVv7x1tNNwJfGKt3rtU1MJ83A5fMhd/Bra0LeCPdacc1bPobPHxL91f7uQH4aVvmsZMs/33AJwbW/5K2TV9B98ZuTavj48CCzdQyfl53AS8bP47u7+umNu+LgEXj+vwI+HFbxl6t/Ry6q+Un2xah+6aAO9r8x6Y/pm2jG4BvAnvOxr53cOjj4O3SJEmSeszTuJIkST1m2JMkSeoxw54kSVKPGfYkSZJ6zLAnSZLUY4Y9aQSSPDNwj8+Lk+wybvwfJXkiyUsmmceiJJe0x28eePzudgurNUm+nWTfgT6HJrk1ybokpw6075nk2nT3Vf1S++5FkuzQnq9r45e29r2TnDeDm0SSNEsMe9JoPF5VK6pqOd33kR03bvwquu9HO3KSeZxEd+eC8e4E3lRV+wB/BpwNkGQ7uu8vOwzYC1iVZK/W5wzgzKpaBjxE951utJ8PVdWrgTPbdFTVTcCSJD8/9ZIlSaNg2JNG7zt0N3oHurspADvT3RZs1ST9fgu4dHxjVX27qh5qT6+hu90UdLezWldV36+qp+huUbay3fHhYLrbsAF8Fnhbe7yyPaeNP6RND939Uo+eapGSpNEw7Ekj1I62HUJ3B4Ixq+junvBN4DVJXj6k3550R9yeHD9unGOBr7fHi+nurDJmfWtbCDxcVU+Pa39Onzb+kTY9dHf0OGAzy5ckjZhhTxqNHZOsBn4I7ApcMTDuaOCL1d1g/gLgqCH9FwEPTLaAJAfRhb1TxpqGTFaTtE/WB7ob1+822TpIkkbPsCeNxuNVtQLYA9ie9pm9JPsAy4ArktxFF/yGncp9HHjhRDNv8zkHWFlVP2zN64HdByZbAtwLPAjskmTBuPbn9GnjX0L3GUPa8h+fWrmSpFEx7EkjVFWPACcAJyd5AV2wO72qlrZhN2Bxkj3Gdb0NWDpsnu2iiQuA366q2wZGXQcsa1febk8XJC+q7gbZVwFvb9MdA3y1Pb6oPaeN/0ZtuqH2L9DdmF6SNIcZ9qQRq6obgBvpwtfRwIXjJrmQcRdCVNVjwB1JXt2aFgBjn987je5zdf/Qvt7l+tbnaeB44DLgFuDLVXVz63MKcFKSda3vua39XGBhaz8J2Ph1LcBBwL9Ot25J0raRTW/SJc0nSY4EfqmqPpLkRGBxVX1wGy17B+DfgTcOXNghSZqDDHvSPJbk/cDrgeXAO6rqf7bRcpfRhct/2xbLkyRNn2FPkiSpx/zMniRJUo8Z9iRJknrMsCdJktRjhj1JkqQeM+xJkiT1mGFPkiSpxwx7kiRJPWbYkyRJ6jHDniRJUo8Z9iRJknrMsCdJktRjhj1JkqQeM+xJkiT1mGFPkiSpxwx7kiRJPWbYkyRJ6jHDniRJUo8Z9iRJknrMsCdJktRjhj1JkqQeM+xJkiT1mGFPkiSpxwx7kiRJPWbYkyRJ6jHDniRJUo8Z9iRJknrMsCdJktRjhj1JkqQeM+xJkiT1mGFPkiSpxwx7kiRJPWbYkyRJ6jHDniRJUo8Z9iRJknrMsCdJktRjhj1JkqQeM+xJkiT1mGFPkiSpxwx7kiRJPWbYkyRJ6jHDniRJUo8Z9iRJknrMsCdJktRjhj1JkqQeM+xJkiT1mGFPkiSpxwx7kiRJPWbYkyRJ6jHDniRJUo8Z9iRJknrMsCdJktRjhj1JkqQeM+xJkiT1mGFPkiSpxwx7kiRJPWbYkyRJ6jHDniRJUo8Z9iRJknrMsCdJktRjhj1JkqQeM+xJkiT1mGFPkiSpxwx7kiRJPWbYkyRJ6jHDniRJUo8Z9iRJknrMsCdJktRjhj1JkqQe+39r07tm/3y98gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x648 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import aplpy\n",
    "fig = aplpy.FITSFigure(image_file)\n",
    "fig.show_colorscale(cmap='jet')\n",
    "fig.add_colorbar()\n",
    "fig.colorbar.set_axis_label_text(r'Photon Flux $photons/cm^2/s$')\n",
    "fig.add_scalebar(0.2*u.arcsec)\n",
    "#fig.scalebar.set_length(0.02 * u.arcsecond)\n",
    "fig.scalebar.set_color('white')\n",
    "fig.scalebar.set_corner('top left')\n",
    "fig.scalebar.set_label('%.2f Kpc' %(sep/u.kpc))\n",
    "fig.save('../../RXJ1131.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## True Color"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "#!/bin/bash\n",
    "\n",
    "punlearn dmcopy \n",
    "pset dmcopy infile=\"full_source.fits[energy=200:1500][bin x=::.1,y=::.1]\"\n",
    "pset dmcopy outfile=soft_img.fits\n",
    "pset dmcopy clobber=yes\n",
    "dmcopy\n",
    "\n",
    "punlearn dmcopy \n",
    "pset dmcopy infile=\"full_source.fits[energy=1500:2500][bin x=::.1,y=::.1]\"\n",
    "pset dmcopy outfile=med_img.fits\n",
    "pset dmcopy clobber=yes\n",
    "dmcopy\n",
    "\n",
    "punlearn dmcopy \n",
    "pset dmcopy infile=\"full_source.fits[energy=2500:8000][bin x=::.1,y=::.1]\"\n",
    "pset dmcopy outfile=hard_img.fits\n",
    "pset dmcopy clobber=yes\n",
    "dmcopy\n",
    "\n",
    "punlearn dmimg2jpg\n",
    "pset dmimg2jpg infile=soft_img.fits\n",
    "pset dmimg2jpg greenfile=med_img.fits\n",
    "pset dmimg2jpg bluefile=hard_img.fits\n",
    "pset dmimg2jpg outfile=truecolor.jpg\n",
    "pset dmimg2jpg clobber=yes\n",
    "dmimg2jpg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wBDAAIBAQEBAQIBAQECAgICAgQDAgICAgUEBAMEBgUGBgYFBgYGBwkIBgcJBwYGCAsICQoKCgoKBggLDAsKDAkKCgr/2wBDAQICAgICAgUDAwUKBwYHCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgr/wAARCAFWAVYDASIAAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEAAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSExBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElKU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3uLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD+f+iiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKcTD5KqsbCQMd7lxtK4GABjgjnJyc5HAxy2gbVuoUUUUCCiiigAooooAKKKKAJnspo7CPUWeHy5ZnjVVuEMgZAhJKA7lX5xhiAGIYAkq2IaKKSv1Er9QooopjCiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKVVDBiXAwMgHPzc9B/P8KSm4tJPv8A1/VwCiiptP0+/wBWv4NL0uymubq5mWK2treMvJLIxAVFUcsxJAAHJJqW1FXewm1FNt2SIaKmuXsHht1s7aaORYSLppJw6ySb2IZAFGxdhQbSWOVZs4YKsNCd0Cd1sFFFFMYUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRToJEimSWSBZVVgWjcna4B6HBBwfYg+9D2Gkm0m7DaKKKBBRRRQAUUUUAFFFFABRRTo5ERHVoFcuuFZicocg5GCOcAjnIwTxnBANJN6uw2iiigQUUUUAOELtC1wCu1WCkFxuyQSMDOSODkgYHGeoy2iigbtpYKKKKBBRTpJEdEVYFQouGZScuck5OSecEDjAwBxnJLaBtJPR3CiiigQUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFACqjuGZEJCjLEDoM4yfxI/OkoopvlsrfMAooopAFFFFABRRRQAUUU6SZ5URGCgRrtXagBIyTyQPmOSeTk4wOgAoGrWdxtFFFAgoop0EiRTJLJAsqqwLRuTtcA9Dgg4PsQfeh7DSTaTdgnEKzOttIzxhj5bugVmXPBIBOD7ZP1NNoooWwN3bdrBRRRQIKKKKACiiigAop0881zM9zczNJJIxaSR2JZmJySSepNNoV7Ddru2wUUUUCCiiigCaysptQmaCB4VZYZJCZ7hIl2ohcgFyAWwpAXqxwqgsQDDRRS1uLW4UUUUxhRRRQAUUUUAFFFFABRRRQAqqGDEuBgZAOfm56D+f4UlFFNtNLT/g/15AFKgQn94xAweQM844/Wkooi+WV7XAVHeJxJG5VlOVZTgg+tDKFCkODkZIGfl56H+f40lFUpvkcXt6vR97bXa09ACiiioAKKKKACiiigAopWQqFYkfMMjDA98c+nSkpyjKLswCiiikAUUUUAFFFFABRRRQAU6cQrM620jPGGPlu6BWZc8EgE4Ptk/U02ijqO+lrBRToJ5raZLm2maOSNg0ciMQysDkEEdCKbRrcNLeYUUUUCClATYSWO7IwMcEc55/Kkopp26XAKKKKQBRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUU6SZ5URGCgRrtXagBIyTyQPmOSeTk4wOgAptA3ZPQKKKm0/T7/Vr+DS9Lsprm6uZlitra3jLySyMQFRVHLMSQABySaTairvYltRTbdkiMiHyVZZGMhY70KDaFwMEHPJPORgYwOTniZks7r7JbWFu8cxTZcvcXSbJJC7YZcqojXYUBDFuVZtwDBVhnEKzOttIzxhj5bugVmXPBIBOD7ZP1NSJbQtYSXjahCskcyItqyv5kisHJcELt2rtAOWBzIuAw3FZeyf8AWv8AX/BKnyxqLm6O2j36dL799tmQ0UUVYgop0M81u5kt5mRirKWRiCVYEEcdiCQfUGm0aj0t5hRRRQIKm1B7CS/nk0u2mhtWmY20NxOJZEjydqs4VQ7AYBYKoJ5wOlJZJZSXKpqFxLFCc75IYRIw44wpZQecdx/SmSRoiIyzq5dcsqg5Q5Iwcgc4APGRgjnOQJ05v6sV7J8vtL+Vr+nT8n6+Y2inTyJLM8scCxKzErGhO1AT0GSTge5J96bVLYTSTaTuFFFFAhXR4nMciFWU4ZWGCD6UswhVwLeRmXauS6BTuwMjAJ4ByAe4GcDOAEw+SqrGwkDHe5cbSuBgAY4I5ycnORwMctpy5XJ8t7a7lbK2j2/r/Pf1CiiikSFFFFADoZEics8CyAqw2uTgEggHgjkE5HbI5BGRTaKKB3bVgooooEFFFFABRRRQAUUU6GNJXKvOsYCsdzg4JAJA4B5JGB2yeSBk0bDScnZDaKKKBBRRRQAUUUUAFFFFABRRRQAUUUUAOmEKuBbyMy7VyXQKd2BkYBPAOQD3AzgZwG0UUDbu72CiirVsVS3fUv7TSKRHSH7LHvSWWN0cOylV2BQFCtuIJ80YDDeQ0ovd2d0uttXu2r2S9H+GqvFNc23l/X/A8yuIJmha5WFjGjBXkCnarEEgE9idpx9D6U2nQGFZka5jZ4ww8xEcKzLnkAkHB98H6Gm1Otymo8qa/r/gfje5JJO9zcPc3skkryFmd2fLMxz8xJznnk+v600TOsLW4C7WYMSUG7IBAwcZA5OQDg8Z6DBJC8SI7FSJF3LtcEgZI5APynIPBwcYPQg097uQ2wtIxsTgyKjtiVgWw7AnGQGIGMce5JLbvFRSVl/X6adi/eUpObd2vVu/+d7t3v8AMipVUMGJcDAyAc/Nz0H8/wAKmsItPllcaleSQRrBIyNDAJGaQKdi4LKAC+0Fs/KpJAYgK0FUvdd2r/0107P06dGY3u2hUKA/vFJGDwDjnHH60lWruE29hbb7W0BuIfMSWG53ybRJIpDqHPlsSPusqnaqMBhtz1lCENvYg4+UAZyc9/TjNKEvaJJW69vz+Wi+XUp8m8ZX+/7te33CzSJK4ZIFjAVRtQnBIABPJPJIye2TwAMCm1YaO2t7KSC7sbhbtpkMMplCosY3h1KFcsxbZhgwC7GBDbgVhknmlRI5ZmZYl2xqzEhFyTgegySfqT60ShUpTcJxcWtGmmn9z1+9IFNTu3+ln/wP16dRtFFOgELTItzIyRlh5jogZlXPJAJGT7ZH1FJ6IEuZpDaVXdAyo5AYYYA9RnOD+IH5UlFOMpRd07CCiipnvZpLCPTmSHy4pnkVlt0EhZwgILgbmX5BhSSFJYgAs2Zd+gnfSxDRRRTGFFFFABRRRQAUUpKbAAp3ZOTngjjHH50lNq3W4BRRRSAKdGISjmWRlYL+7CoCGbI4JyMDGTnnkAY5yG0UDTs9iSYiZ5LmOGOJDJxEjHC5yQBuJYgdMkntk81HRRTfLZJL+v69fXUG3JtvqFFFFIQUUUUAFFFFABRRRQAUUUoYBCuwZJBDc5HXj0//AFU0k93YCbUBp8d5NBpU0k1sk7/Z7i4hEUkkeflLoHZUOBkqGbBJG5sA1FGYQjiWNmYr+7KuAFbI5IwcjGRjjkg54wXJLL9nktxclYyVdotxw7DIBwOMgM3J7E+vLZJnlREYKBGu1dqAEjJPJA+Y5J5OTjA6AChqfKubt5dNNbdfkr76lQ5YR3u0uqvd7flre2/RaEj6ffx2EeqSWUy2s0zxQ3LRkRvIgQugboWUSISByA656ikXzbTZMPKPmxNgHY+AdyHI52nrjOCOGHY1FVjStK1HXNRh0fR7N7i6uHCW9vEMvK56Io6sxPAUcsSAASQKzb5Ytyatr939bihJxu3v0e1vv8r31S6kcUlxFDJ5UxVJAElVZMbxncAR3GVB9Mge1E8wuHad0AkeRmfYqqvPYKoAXvwOOmAMVZhsbNNKk1G8Nwz70WFIFQx4ZZRl5NxMbBkUiMrl13ncmAWh0yS6iv4msIGkuC2LdYy28SHhSu0g7gSCPcDg9Kr2kLNxXSzv37rt09Ne5tGNSUo05tpN7WfW2qXW99LbvzTIXd5XMkjlmY5ZmOST61OLzzdOGnTlFWF3lgKWib3d/LBDycOVCpkAlgpzgDexp+lJBIZFuGjWNYy1wZJFBdAVISMlWKuSMAgH73OFDGlvL3TriAPBYSR3Ui/6VIZIzGzGR2JSNY18sbTGAATgo5zhwqOb563vu/VvfVpvW9ne+jd3vfUl0oxpwm+t7J+S30fnZab+SuVZpElcMkCxgKo2oTgkAAnknkkZPbJ4AGBTamtrmGCG4il0+GZpoQkckjOGgberb02sAWwpT5gww7cbtrK25tzbSmB5Y2ZSQ3luGAIJH3hwemcgkEEc0RTd/Lz/AKe/+e2pm3OTcmuu+na/Tb7iW20bVbvTLnWbawka0s2jW6uQvyRs5IRSem5trEL1IRjjCkitU1wt4AsN88i+VCphjm3cIx3gKOwO8t2ByT35jnjSKZ4o51lVWIWRAdrgHqMgHB9wD7UWaim+t/zt/V/lpqVKLjfTbTXR7X2+/wAtr2eg2iilITYCGO7JyMcAcY5/Omle/kQJRSlyUEeBgEkfKM8479e3+c0lDUVswCiiikAUU6eZ7iZ7iQKGdizBECqCTngAAAewGBTaFsN2TdtgooooEFWNV1O81rUZtW1Gd5bm5cyXM8sru80h5eRmckszNliSeSx6dKr0UuVc1+oaf1b89/x9b2VptQ0+/wBJv59L1Symtrq2maK5triMpJFIpIZGU8qwIIIPIIqGiihXtruJXsr7hRU1zqF/eQ29teXs0sdpCYrWOSQssMZdnKID91d7u2Bxl2PUmoaFe2oK9tQooopjCiiigAooooAKKKKAHQxpK5V51jAVjucHBIBIHAPJIwO2TyQMmpLnT7+zht7m8spoo7uEy2skkZVZow7IXQn7y70dcjjKMOoNQ0Utbg29LfPf/hgooopgFOmheBwjlSSqt8jhhggEcgnnB5HUHg4IIptFA9LCooY4ZwvBOTn06cevSkooptrlSt/X/A/4cQU6Gea3cyW8zIxVlLIxBKsCCOOxBIPqDT45fs8YltrqVJjvRwo2gIVx94HJyCwIwBj1yQF0+9m02/g1G2SFpLeZZI1uLdJYyykEBkcFXXjlWBBHBBFQ7tPQqTdNJwfvb729NVt59vPUtWGhahcaVPr0UNpLb28UhmSS+jWRAGij3+XvEhw08e3ghsN94RybW295pIuTFqsN5eWkMFxHaeVciCQMyP5LHcsihFlKuyAZYb1DqWDq/Q9EuPEETaXo+kXV5qk11EllDaPveQEPujWBULyMTtIYEBQjAhi64q2JuFuoRpgla6Z18nygQ6ybvl2bTkngY9z06GiEW+fmmlo+17fJt7q17Rdr9LsdN1NKbT1fRatXS0s/eV9r/aVkl16nR9Gj02118+MPDllqt1Po63WnXi60kjwyC4t5HlURXCi4DRmSJ1G9kDSvgNA5Ro+GWt2ej/2hd+F7h4jhxrbXX/Ev8iUWwjmyqgjyzcR7/mO1riFXVGBST0/4N+AV8XNb/Da88VnStJufDjxaydL1j+0vtgaVLmW8Nml2nmxWcTGdhFG2Hs2XaWiluI/oP/hlmDx7pviK71fwdp2iy+CtJeXR5Y9RXU9L1KwguIjdSWMFzZ3t4YTJLJftP5zpGupWzSROlzKw+fxWZU8utWnNOElq7VOVNyjSivdUlF88l0cfdbkmmpL6WplfNh/bUm5zi17tlJzk7c0OZO8m1JuNoyd2uZOfK5/EE2mzpPKiXc9pqN0t7NcPqc0YiktBGJEMc5I82SQrMmVAEmUVS3mFRmz6bbzyJd2bxW1pLyoub+OWSMBlRiwQBupzgJkrkgEAmvbx8FdS8La5odv4WmvdGiutZtdT1DULWBo9S8P6ablJNPvLhmkaG1uZI72227LhQjNbpIUaYPLwfjXU9Xg12XXtR8cL4w1vUZRL9p1dZNQKCP7PIbwyzExyFvKMTKwkVUjlR/uqx7cLVxVeb9nHZRdnpJ88uWFk9GnZtNOVoptNwfMRjcDVpYd17e0pu9le0lNubUfhTcZKN5Pkg7yUXzVIM4eG0ge0XUmiBhgmjju4/t0ayyMxdgY0I3bdqYLBWCtjcRvVTG9tCthHeLqELSSTOjWqq/mRqoQhySu3a24gYYnMbZCjaWn1yeA3Qs7WztoUgREcW43BpVjRJG373LhmQtw2zLEoqBtoSLTHlnm03yvLuYld3+1zpAEEcbtIhDkZf5cKMhmI2hSzAD0YtuHO3by0vZ/1q1toj52pGpG/PZWeui0aVuW6untut3d6asqwCFpkW5kZIyw8x0QMyrnkgEjJ9sj6im1bEkR0+OVdHjPkl45p8SYYuDsLHfgMMMVAAB287uar3FvcWlw9rdQPFLE5SSORSrIwOCCDyCD2q9pWf5p/l+PYudPkpprVd7SW62u9Hs9t7Nq62tppljFDb3l9qbeRcW0zbrW1ZzFOocLC2/YpJIiLMjMFSYH5mBjqjV2w0TUdYkg0/RdJuLy8lSWURWn75jGikn5EBZdoSRiT/Dg4AGTr6B4Xn1e90eabW79m1Jlhnt9M0+S6vUtzMltmKMMqS7gxRIzIm4xsh2Bo2fnnXhRTc5d/1ey8k9Xa9vUIYd+8221FJvayTUm77Jaa6u1ldvvzs0aROFSdZAVU7kBwCQCRyByCcHtkcEjBptWriFbm0k1e41yGS5eZS9u/mmZyxkLOWKbTjauctk+amA2H2sjtbp0jtU01nkumU2zBG3uMsuEA4YFuOhOVwCOQdlJW1/y9d7C5PaVZKCenk3vay0vq7pdVd230IpJ5pUSOWZmWJdsasxIRck4HoMkn6k+tNooq7JGbbk7thRRRQIKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAv6VFqVpNDd6dqxtBcwzQvdLK8SohRkljZsDdmNvmVd25ZAuCW21SEiCFojApZmBEhJ3KADkDnGDkZyM/KMEc5ltbCW8hkktjloULzBuAqAfe3Hjr8uDgksoGScVBSdNwlq07pPS3na9uvk9u2uu0lGMVOMLX6vrZJO11snfXVq7V1YnNiUNqZ7hY0ul3CR0cKi72Qk/LyAVJ+Xd6dQQOph0xLdU0rwnqN/e6z9o8nUraLy5w98XdIXsxDKwuME7BNkkGfci5wawtIv9OsJVubK8vLO8OEiulmGy3ygVnO1dzBsvwuCgxy54rqtPtNE8R+L/DniTwP4gg0vW7u6hlv7S/u4LREv/PmYzW7xLBb21v8AJFtR3iaNpAAfLUS1nKTUVzylyvmvFX1dmldK7a1Sso3kk1dcyOqnUrYWvD6rGLlNNKTd4wk1FRcpJxlHW7dk+101p93/APBOv9n74T+O303VLHxp450ue58aadJ4WstB1Ke0t5obS6lkubq1We2ltry1Biillh3yuk2nxJJNK8iC4/U/9n79kL/hZf7NUvhC18A6TpmqRRpb2FnpGhtaHTo7OTzoLOS5lW4eVortZpjcoHimM5m+zyCUmT8//wBlXW7v9kvVPCdnqWq6e3j3xGbO3eVfEEtsL+9t9cuorv7XbpNMJw8d6qO5H2iF/JWGCJoJVg/Tj4F/tT3PhT4fajqem/Dt9bibUbq1v0h1XddvKrXNw5nS+CsEe3a2MYd/K33KKqLbqs7/ABuHp4WdWkk3UjOU6k+aStF01CkkrOKfKpO858jcpVqbbaml9zgU6eQUa2Em3Go6coqVS8ocsmrPmtJKVo8ynJzm2lJU+bkPiX9of/gms3w8+ENz8PvH/iDQ9Pt7W6ivLbU7PwnbxW1tNp8iJaqZHdLaQyx3T4huLlYoX8qK2UQqS35jfGH9ny3g8P3fijSLHVLXRtV1Et4ZefU7a7sESGzkTVb7U0t7m5W21AO1gTDE4YPPsSHd5dsv7E/8FDP2vNH8K/EPR/itqn23R54L22sbu41NDNaNpb3dq15HbvJIlnbbjJYI9xckIRthEhlXyo/yU+PPxG17xJ8Hb/4kXXimMaSNPg8MLdW1hHqEGrXsIsri2gdZ7VPJhNvapKredJIHh3OgScKfocswOVUMfUcsM6c60kk5t0oJOMW5QjKUrzivgtL33BwqQTU3LLG1Z5rl88ROsoTp2mlyxVNOa5Eow0jZKDUORqPNzQlGLnI+YvEdr/ZyL4dF5aTtaW7oongaKVAZXl3ZY4DbFX5CRgzFQu/eawJpElcMkCxgKo2oTgkAAnknkkZPbJ4AGBW74v0/TU1OeGwae3htPLiQ6wXW8lAt12lo+QgYxnai7hH5ioZHUK5p2F5f6Jq02p6HPPbT6bIZbS4tisrQNvVVfzlxgrkFZF6PtK4JyPVmpQjaUr+7eN7O97WvytpX67tW1V3r8fjZTxc+eFN009FGV+aNk2lJNcycUuV6+nw6ZrFCF2KQcfMSc5Oe3pxihEMjiNSMscDcwA/M8ClELtC1wCu1WCkFxuyQSMDOSODkgYHGeozoanpunm71Ga38TWMy27ExyR28kS3jF1U+Snljap3M43iP5IzkKxVDcZc9Ta6Su7aaLlT/AD6XerdnZ28uXLCF27PS2nRu1/v/AM9tTbh0i98baRp/9j3TW+laRbxwa1dTabLImjQyTxK11dzW1uS1u9xP8mQ8in90FYiMyeleD/gfrPxA8NyfDqw+DOo2XiSHVEW0hOrsLeZykaSWs8FzPCE1FVTzhscBjOITCpkh2YvwjeX4oa/cWur+OZdGS7udMis5L+5nnAlgkhVpyxguXl8uCLc6KIly6beEitpf0w+AH7J198UtO1jxtputTt458TX2lxeIvCyQRR3ltHCk1wz/AG+PzH0+32T21/bFUbyP7MgtDFNKkkUvzeYVMbQpRjFxjKLXI5TfLGUpXU5JJuUIybUlZqVRxp6TleX22SZHSqxpYitKShJNOTaVmlJfzP2lLldNz9ooxWkbxlKi4fl7478DaL4L8RjSrrxHpWj66t5ZWd9ZvNFeW9oZYA1xOJ7GSaGS2QeWMbXkYXMqMrtCWbiJtONhYXF7c+Fr9Va9eG1nuXYRRbFIljbCLvlUy27ZDKF4DIwkGPsz9o/4C+DNF8MeJfit8PtAg0G4gs5ZblLTSjq9lbwubq3RrS4uI3t7f7RfwXDJLFdymN5WjilEsNvb3XzB8U9I1vRpLbUPFGn39nbFY5NC0q7spWZbC4e5aYG5uYlYyJOs0edhVpRO8eFUFuzA1a86adROetpTva8uXm5Fa8b6qVueTcJNpuNpvy86y/6nmVXCRrSl7OUXZwvyJq7U3GMVGS05k+aLSkleVNqHB3Fv5MUEm2UedEXzJFtB+dl+U5+YfL145yO2TFUzvYGwjjjtphdCZzNM04MbRkJsUJtyrAhyWLEMGUALtJaGvUje2p885KTulb+lr89woqSSKGO6aE3IeNZCvnRKSGGfvANg+4BwfXFR1dvd5v8Ah/uBpxk0+gUUUUhBRRToTCrk3EbMu1sBHCndg4OSDwDgkdwMZGcgGld2uNooooEFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRUtrceR5iO0uyWIpIkUu3d3XPByAwVsd9vY4IZNIkrhkgWMBVG1CcEgAE8k8kjJ7ZPAAwKV3ctxhyJp69iS5uYZ4beKLT4YWhhKSSRs5adt7NvfcxAbDBPlCjCLxu3M0QchDHgYJBPyjPGe/Xv/AJxSU4CHyWZpGEgYbECDaVwcknPBHGBg5yeRjlxbp/D/AFfciMFay9f16/h9yCCZ7eZLiMKWRgyh0DKSDnkEEEexGDUupXFteXst5aWMdqk0jOLSDeY4AWJEaF2ZioGACzFvUk8lbLVL3T7a7tLV0CX1uILkPCrEoJEkABYEod0a8rg4BGcMQWebbGy8g2n74S7lnEh5XHKleh5AIIxjLZzkbZTtJtx8r/j5LfTrp2KilLd2sn6Pbt19e26W920PinTI49PhtrtYruaGdLOSAtFcvh1hfy2BWQ4eQKSD95sdTXo/h2HR9R8BXKeCtLkltrSa2ttWt576a0i1G1lunH2y5WSdoLYhhbRPLuaGN2tMYkcsPORLqTagNV0XWZ5b1bUz3U8RMbo3O/ac7pCFOWbAOd55Vd573T/ht4q1O4tPBmjeHtO1DUtRvGeC90TSV1BCs9mo8lPKDESyb4yqSKFikZGje3/ek8uNnh6bjzWbi1JapK6d5Nxad48trr3pJxvok2vdweHzKjRrVMsp80ndNThCceV8l5Kyac00lLRxUe8otr6s+BfwJ8AeIPC2k20txp2ivqvjK0PhLSbezluvKN5DLdQwLBqiGO/bP9jQ5S5NtCdQEjy5lDN9oW1946+A2q7bSz1m/wBIkvm1zwhP4ftpxY2FrHql5a3GlMgkjhWCSG+EiXNxDHapBO8gkc2Vnv8ACP2BdB8San8P7jx/8TNW8U6wdQ0bTvEGgyfEBtTVLe90xWe7eWOS6hgaNbU2SwagkkUbpcCOVwZJoV/Vf4UfsZ/E74u/B/V/D+qa3f6He65bXjXs1trnnpbTwyLGj2jWkdqsQnWBbl3iNvOr3kuxoJWV4/n6ksVjcbRi5pqEt24UpuPtdeVOnOUYUnGNJtRxGsJJziqra+oxOYSx2Eg1Ug6S9ySkpKKb9m17G1NrljFS5pcvvStePvcsfyQ/aA0vwn4N8H3PiH4T/swX9rrNkLvQ59fT4c20iXdta3V9DJfRWka/ZrWY2dhfyeSLW0eOaZiZUjLQyfL3xh8F6/8ACx9Z8FT+JdGtriHS59DuI9O1rSLWW2u7acXBNxsRJBJ5ls0MTEI0sEqeW4VZYB+zX7bX/BMPV5dE8ceGZvEkzaZ40lN7fwQ6lJHO14I79biVIIvIt5MRTWcC7pERnAmuC7I4n/LH9s/9mLVPgzFb2U+t6JfaX488Tw22u6PfXV34ZvXurCP7RcX5k1TEGn3E3m3EUkUhktzI0csULlo47XsjKp9WjCrUcZRtKzTm+WV7xdRuVo04xSacmqk0pKMbSOXFZzPBTj9Wnz0HGUY2qTc7xgprmulOLprnjzQkoOEeRx5Wmvke8FxrWprbaf4jt9UAtLq8mupoYra7ithEYyJG8xfNl+zxJJ5AlkALFQS7TA85YWmtTwyWFhpLzfabdpyFshI/lRB3eRW2llVQjlmUgYRt3ANbfiy+tP7GW/bX73UNR1y8kvr+a909AJdstxGswl892Z33OWDRoQQCCfvNgXcdpbWiWcllcxahFcyi7MrgJswgRAm0MrhhJuJJzlQAu0lvdlTcW4K3u6W5eX4ftfE7ttt6drrRpR+Mq1MRTpJtyc3dyctfenfmu3ZyUoxTUrJyUozsk2RBkNqUeeTcsgMce3K4IO45zweF7c+owM2oH02azv5723kT5R/Z8NtcqEimaRSdyvuZ08pZF4OQ3lkt2axrOq6ajNDoegw2S3GnW0NyjTpdBisUBeRGZS0TvLG8h2kFRKYuFDBs9ru4Nmlt9oGxSwCBcMAdpwTjlcjIGSAcnAzk6UpTrRbk5LTur2VtHrtotFfS+mmnPFvl5akbabdb26tbNX1s7PTuz3r4K+PZfhZ8SLOBbCfQtYNz9psNd0rXba6k09YLlo/JhS6ea2jit5opJVaZlZo4m3XEUUxlX2Dw9+1Z4r+BVtdeAvDnjrxvLp7HWmJ0uCPUtF1NsDT7YJa3qxxfYYxFPHHEkIBexijjlDOfsfyzpF7putPd3lvpSK/iLXVsreCxuprnWPLZlaQCGNo4plbcoIZYjLI+IiqrIENDuvD/AIWiC20r6veP4eu/s0ul3UtqLaSezzIJHBTe0Qkfco3BxEY23LIvleT/AGZTxThOdSUXGE+eEWknFuDbTtJQeqXx8zlGOsWfU5bn/saVSlhbTlGMrNxtUbvaMU5zqpSkkuaKi4c9NTUuePM/fP2j/jrD8bfDVv4n8R/DGey1nXrnW/EmhJfG4vp7SFLV7Fyt79ohLmJNLsVJktpZvNspZWmCuY6+fvix4ln8b38Xi2z8F63YWFzahbBNX1SS/igsoGNtbwW00kSyC3gTZAAXfmNQxyMVn+KbrX9Q1tPEfi/xXNd3+oW0P2fUTcMRHGEiRAeN4jSIGMKqhV8sIudhUM8S/wBhfYbWy8FS3sitA329BeF42dIofNKxbFZVZ4mm3twVdFwphYnvqUp4KtHm0qSbckmpQjrKVtFFJ3k9JRvpdt3VuHE+2qSmuVqEeZaWU5c9nFyg5tWtGDvCMmlOST25edqa2toZ4biWXUIYWhhDxxyK5adt6rsTapAbDF/mKjCNzu2q2r418UW3jLxHJ4svINtzfRF763t4jGkc43KpV3eR5sqscjySHzHdpNxJ/eNiU6cp1KaclyvS63s+q8/Xtqj5+PPOjFvRvdb2el90tN0n1WumgUUUVqWFFFFABTjBMsK3LQsI3YqkhU7WYAEgHuRuGfqPWm0UajVtbhRRRQIKKKKACiiigAooooAKKKKACiiigCZLmFbCSzbT4WkkmR1umZ/MjVQ4KABtu1twJypOY1wVG4NDRRSSSEklcKKme2hWwjvF1CFpJJnRrVVfzI1UIQ5JXbtbcQMMTmNshRtLQ0JpgmnsFWLDTLnUfM+zy26+WjO3n3ccWQqO5xvYbjhCABkliqjLMoNeih3a0Y+j/r+vQdJBNEiSSwsqyrujZlIDrkjI9RkEfUH0ptFOnnmuZnubmZpJJGLSSOxLMxOSST1Jo1KfJbS/9LX8fw8x/wBrkjvft1kPszLL5kQgdh5RzkbSSTx2OSeOtRU6AwrMjXMbPGGHmIjhWZc8gEg4Pvg/Q1JbPZojtPEWcRsI1IJVicAZwQVxlmzyCQoIxmmkrPyX3+X66+eve1eo9ZWu36LS97LTy08la1i14f0bVdflk0+zaVLWPZPqNx5MskNpEGCfaJhErMEQyYLBSfmwASwB7/4P2OgeLvEWj+GPEdsiR213u8vW9RaSGDTwWnaTyAqSSGNJriVYoXzcOY1W3kJbfw+q3PhmOxkXR/DmpWsl0wa3lvtQinXyBLNkAC3QlsCBd6sOYpuCJQkWv8N7+41bxXHd6z4WtNchs7CcywXaSwpJIzOYGnntlE2TdSwqJC6klkiaWNDleKvd0ZTbcU1o07Sg3ZN2vyydrrkb5WrJ32dYbFvDuNWdKUpLRQtDm1VoqPNze9KUocq0hezkpK8an6k/8Evv2lbP4L+AtQ1zX/h1camLzR57+21Tw5p006aRpv2W+knt1na4uPOtUltLaIBmTZPcahJJMGMrv+n/AOz5+2d8Ofg58DNF+KPxvbSfBMbefcXnh6bRZbM6ZcTO87KkU4Dyz+WXLSqv7/55IwA+xvwB0jxv4w8XeErG+8T/AA6e7j8DWlroWtagmpnT7+ICYRLYhlkKzzeV9uLlopICqbzFFM128v1n8Df2g9J0vTofGF9GviCbVbm60y28WeDonls017VZrq886FprXykYWxSMSssMitJbwX7MIT9hwngp5DVdejh+edWPvKhPXmSdSVOEKkpN8seSpKfwzi3yKpJ86/RIrLcZhp1q9NQlKpedRVE3VjGEIqTk4uEee6aa56kUpc/LG3P+hX7Zn7cPh/4l+HL8+B7O31+4tlUa3okukyTtLZgTG+SJJJF+0u0BWKIQ7wXuLSTyzFcRyD8aP2jvid8KvilNL4R+D2kHxRd7LmPQoHu7galpsNzPBGbK102OYLJdQXca3Jk3oJpp7gxWpt0gkg9r/aW0P4vfGP4c6v8AFn4o6H4R1i2is7CLT7rRdZF9rWr2aWzzR3Zt7QC2eyKnULgK3KQ3VzKZljtRHJ80+L/CeieAtdiHwm1TStPOpad9l1dZ0uHu1huUXMgmVpJk3Hzbma2inlkSJhDK7x/agni4fHuvi418NUq1KPuy5G3yX9m58s24qVuSpdKkoyqJ1MMpThHmgYLA1f7Jl9Vb5YyV5KNRxlOSp8vJbWz55qM2oOioThUk6k4J/PvjfxXp+u297BfWkFvfbx9rgsLT7LFNeJKwE4TYQsYiLxi3RbdIspsT5WVue+z6VdapfxWpis7XbM9qdQuWmaJVyyR74kAkkbAjDbAhL7iEXLL13xJutC8S+K7vQ28W63LBpOr3cMV9rd9Hqf2az8wu7/abUssxluZJZAsY8oPMWEknmtIOcgTw1qurQQanc6ja6dawyC+u9Og+2iMebIUaGKVoSkZLxrteTOSzZJYIPpKM6UcLCfJbmjf3bpLRWXJrZ/Dp113T5l+eYirWxWHjjbSXMrKOspRsnJrllKT0Tauk4uV7cr5lGlqafYo5NEu5F8+yuWjAtblZoHOWEkgdXZGJ2xgMnysqg5PBNGrVhDfyrLqFlcxQtp0SzbzcxwyY81VBQEhpHDODhdzBQWxtRiI0vZo7CTTlSHy5ZkkZmt0MgZA4ADkblX5zlQQGIUkEquOltL3Y2drX0t5vv3fV308jCVWVRK7vJWT6aW369Oi0t1vdt8mtaxNeLqE2rXLXCQJCk7TsXWJYxEqBs5CiMBAOgUAdBirVvrOrWtoyXsbzxyWx+yG6uJVESlkjd4grqGLJF5JB3DZuGMqpXNLuUEZc7QSQueATjJ/QflSzGFnBt42Vdq5DuGO7AycgDgnJA7A4ycZOkZRptKN+Xtsr6dFdWulo0+ZKzVioOUKLjGTS0sk2ut9lb9Ne+pc17UZNUulubrUGuZoo47YSfOUMMMSRRFWc7sbUxtIAUAAAD5VRYH1Kyu57HQH/ANFb7RPNbb2S3gZljw+d2FDvGqsSOXwSxZcVBBM0LXKwsY0YK8gU7VYgkAnsTtOPofSgiHyVZZGMhY70KDaFwMEHPJPORgYwOTnjFq/XX+u1i3KaT5IximtraW2009bWtqm/5rtoop0ME1w5jt4WdgrMVRSSFUEk8dgASfQCrukZJOTsldkt1NZvbW0NpDOjJGftJluA6vIXb5kUKNg2eWpBLElSc4IVYKKKd5Pdt+rv8vRdF0WiJSUVYKKdDBNcOY7eFnYKzFUUkhVBJPHYAEn0AptK6Ks0r20LFyv9nyNb2eqpOk1vGZXtt6qdyq5jYOqklG+U8FdyZUsMMa9FFJKyEnLlSbvb+vRfIKKKKYBRRRQAUUUUAFFFFABRRTjGghWUTqWZiDGAdygAYJ4xg5OMHPynIHGQaTdxtFFFAgoopXYO5ZUCgnIVc4HtzzTsuW9wEp0xhZwbeNlXauQ7hjuwMnIA4JyQOwOMnGS2ikO+lhVcqGUAfMMHKg988enSkoopuUmkn0EOgELTItzIyRlh5jogZlXPJAJGT7ZH1FOd4XhUCII6DaSoP7zljuOTweQMAAYHrkmOihWTvv8A1v6/1Yrmajy/1/X3ettB00E1u4juIWRiqsFdSCVYAg89iCCPUGnx28UkYdr6JCd+UYPkbVyOikfMflHuOcDmoqKmzsNSgpN2uvN/5F/Rv7Jtw97qYmkQQzxGKIxBhI8DiFhv3blEnLYUEADDBmUibS/F+raXB9jgaJYGtDbyxxWkK+em9pVWYlD5wErK+H3Z8uMZHlxlMwyIYViEChlYkyAncwIGAecYGDjAz8xyTxi/app4kmsxardGRWjaaGGRvIjRkc3EQ3ruYqkgKyDaFYnAOCmU6VOUZOeu2n6roreer0tqbR53UhKkleNra9XbpK/q7Jr3b7pHd+GPjjc+DNH0CS30zTnXS4rn7JZ6TLJB+9drwq94GyLiSM3CbSyyLJbsYHZhvRe9+GHiS78VeGXtbPwf4XuGspL/AExrl2kGoTE2EjG9nnMqI9skUUhm8q4I2NJ58ctuLeOPxOfUbFrLU7qbxhrrXd7Z2yIj2Q2X43I00cz+fkRxyRgodsm8xqSsRGB0Pwbn1298VWt5HcNLYafCNR1tw1vE0NvYqXIEkiShN0ZEKu8Th2n8oRTHEbRXrV6eBkpSUJQ+Fzg9pJppyXLPl5Z2jeUtZJxcWlboyupRnjZygpK/Ktoe/bm5eaFR+xu5TinNRbioxSdlY+4fgJ4jtvGXiLxdoHwo1W9sPA/iO2WwtfC9vfXGoXsNxPFHqQS7eNlMs6iK4tHmuVQtHALP7VKkVxdWP0n+2Z+w/wDDDVvhjefEL4maLdaVqPi7ULPVLmbT9RsNJtpLyGIw3NtcSXMzRCe9SS58yaKV4kFhG3zywA3Pyr8BPjR8PNT+J3g/wz4a8I28Umnfu5bDTYLCFdan+2Q3dtdpekr5FzPcTwQXIjRj5LghIbZFCeieK/25tL8feC2+FPgDwFrYuptCt9LtLrS4J5Y9HzaXs1ukJjk1C7gkv7RrayktjLDOosvLikTyUa3+Tiq+NzOtUxb5K8l7rfsYSlJxSc27RpxTjT5nUlGUp3jyKVRyg/0rC0s1wmFivZSnUoSfPCTjF1JXvFKak4SWrjUnebcfbSjVcXCmvir4vfAseBta1bTNX8Na3Jr1pETDpen6XpskK3AZroXCvYrKk9obaaFxlYVdZiVdlthG/lE+i6dY+MLi10+xhv7AxXc1tFdXskWbYJKUkWWRIC7+WokjOwB22AxksYq9u/aE+NHwj1Eya94S+H1zYzan4ejk8OalO1vdXz3e6OGS8vpZJ5xBO8KSSCKNUdBJaMrLsEsng2o3Nqdl7Ebm4tbqZzcm6t0E2/cjOgmIbccbPn4PzN8gDEN9JhvaV61SrKUnCSi1eUW1ouaN1GEbwacbRW6snJRcn8ZnmDwUZqFKouaMYtpQcIx57NL2fNN2hFU4aVJKMXLnvUbQ3UzKYZ7K70owLp7/AGeFhbLvjPmSN5csiBQ7Hc/zsGb90qjCrgVLW4ubayuvJW3KTosMvmxxu4BYOCm4FkOYwC6YOCVJw5DO13Ur/WNcvNX1TUze3V1dyTXN4c5nkZizScgH5iSeQDz0FVa7YxTjta+ttXbru7N69bK+9k20fOVfZt8vKrK+1lrzN307Pb87WCnQCFpkW5kZIyw8x0QMyrnkgEjJ9sj6ikRgjhmQMAclWzg+3HNT6ZbWdzeINSungtVdDdSwojyJGXCsyIzoJGAOQm4ZxyQMkaS5Y0+a/wDn/XRLe5mm072vbW297a7adF39NSvRRUsVv51tLKiys0W1mCRZVUJwWZs/L8xQDjB3dRwCm0hxhKbsvP8ADUiooopkhRRRQAU4QTNC1ysLGNGCvIFO1WIJAJ7E7Tj6H0ptFGo1bqFK4QH92xIwOSMc45/Wkopp+61b+v6/pCCiiikAUUU6CNJZkiknWJWYBpHB2oCepwCcD2BPtQ9ENJyaSG0UUUCHRzPEjooUiRdrbkBIGQeCR8pyByMHGR0JFNoooG22kn0CiiigQUUUUAFFFFABRRUsE8aRvBPFuRuflChgwVgp3FSQMtkqPvY9QCE7pFQUZSs3Yioqay1C/wBNma5069mt5GhkiaSCQoxjkQo6Ej+FkZlI6EMQeDUNGtyNb+QUUU6SCaJEklhZVlXdGzKQHXJGR6jII+oPpTuikm02lsNoqZ9Pv47CPVJLKZbWaZ4obloyI3kQIXQN0LKJEJA5Adc9RQmoX8dhJpcd7MtrNMks1sshEbyIHCOV6FlEjgE8gO2Opqb3+H+u/wA/6ZHNde7Z6/8AD7dfL77ESOY3EigZU5G5QR+R4NXbuC3nvLqbQ9GufsU0kv2EXLGSSKNCHOXQKruseNx2gYYnaOMQedJb20f2S9lAdg8qBtoWRSdpADHOFbhiByzAdMl+nsuoapBDqd/DGkpWFrq+MrRwLgIHbywXKoMHChjhQArfdN1HKNGz9flby727bJW3aNoqhTnzTb0XSzXVvSzvZbWV7trohk1/NIgCXE+54FiuC85YSKpG1cYGFAVAFOeUz6AbXwx8QWPhXxXF4ivHtneyAuLWw1HS47yyv5o3V0trqORgpgcr82Q3IAwM71xHvZpLCPTmSHy4pnkVlt0EhZwgILgbmX5BhSSFJYgAs2Vtpr/SpodQgDROyl7eUp2yV3qSOoIOGHIZcggjjNKNrNddr7r1tpdX6O261SLhUqVKilUbst7X2ej2ate7v7yu3vd3Oxgn+IGntaaGltc2E2iadJ9mg1G/jintT9tCzmHzNjRESrIGQBiiiZmG3zHHdeM/FviTwX/bvw40zStE01/C7xwXH9l28YsNdtoLi6kF7PBdShbhne+iSBhEWNrM6jIVt/kOp+fZqtjqOlXEbx26SxxXZ8lCkkY2TLEApy0ZgcMCd23cd4bhovtGbThEwvY7qSB0uZYWjWFgPL8pPLCgtzGSzlskup25jJkic8dFR1Vrv4UnGT+ByavyzXKnFNxdrqSTneoepDMWqcadOreFlsopRvBpqLt7qlJ1OblfNyzS5k43dy5uDc+GrXQorK2a+tL1oIH0+UeY4kL71kUKftBY7Ajo2FVSrA70IzoNSsktr57rS4pru6wsEpXalsC252REwu44CjIKqrPhd2xldf6/Jfwyn7DBBcXV1NLeXFonkrOjmNhF5SYiSNHQsqoi4L46KgWDVdV1HXNRm1jWLx7i6uHL3FxKcvK56ux6sxPJY8sSSSSSa3VSvzpzd9NPek+V2jtd6bbLRbWa0XlVa/1iN5Jp2SabTWmis0o7RS13bu7e9JkMYhKOZZGVgv7sKgIZsjgnIwMZOeeQBjnIRSgDb1JOPlIOMHPf14zU0FnNPCixWE8kk84jtnQEq7AfMgGPmb5k6HjPQ7hiKaF4HCOVJKq3yOGGCARyCecHkdQeDggilGaUv6fl/Xn6CnTmqak1p3166q9+tu3T7xtOmheBwjlSSqt8jhhggEcgnnB5HUHg4IIptKzBgoCAYGCRn5uep/l+FNJNN3/4P9bmd1bYSldg7llQKCchVzge3PNLNGkThUnWQFVO5AcAkAkcgcgnB7ZHBIwabQpPlt0/r5/1sDi4yswooopCCiiigB0EL3EyW8ZUM7BVLuFUEnHJJAA9ycCm0UUdR6W8wpQU2EFTuyMHPAHOePypKmeymjsI9RZ4fLlmeNVW4QyBkCEkoDuVfnGGIAYhgCSrYakoP10/rzJbStchooopDCiiigAooooAKKKKACiiigAooooAKmsrKbUJmggeFWWGSQme4SJdqIXIBcgFsKQF6scKoLEAw0Une2gndrQKdDC87lEKghWb53CjABJ5JHOBwOpPAySBRNGkThUnWQFVO5AcAkAkcgcgnB7ZHBIwabRutC7csrSX4/rqFOIh8lWWRjIWO9Cg2hcDBBzyTzkYGMDk54BIghaIwKWZgRISdygA5A5xg5GcjPyjBHOScwtM7W0bJGWPlo7hmVc8AkAZPvgfQUdR2Sje6f36f1b8duqIBC0yLcyMkZYeY6IGZVzyQCRk+2R9RSMUIXYpBx8xJzk57enGKSiqTsn5/wBaf1toTfSwUUU6CCa5mS2toWkkkYLHGikszE4AAHUmk2kgScmkldsbRToIXuJkt4yoZ2CqXcKoJOOSSAB7k4FNo6hZ2v0ClRQ7hWcKCcFmzge/HNKDD5LK0bGQsNjhxtC4OQRjknjByMYPBzxLLZxRafDejUIGklkdWtU3eZGqhcO3y7cMWYABi3yMSFBUs46O7V0rX/paq+1/+AKTUba7/wDB8vL0216EFFOkmeVERgoEa7V2oASMk8kD5jknk5OMDoAKbSG7J6EspubvzdQuJ/MZpf3rySgu7Nk5wTluhyfUjPUZjV3QMqOQGGGAPUZzg/iB+VJRTi5Qd4u3pp5dPLT002HOXO7vfrfW77/1fvfUtQXepySXV2GlmMsTG9dizblZh8znv85U5P8AFtPXFNi1G4tLaS2065uIFurfyb9EnIW4USCQKQAMqGSNtpz8yBvTFenTTzXDiS4mZ2CqoZ2JIVQABz2AAA9AKTblBQfwrZWXW3lrst9rK1ipyjUg+e7b769b31vbXXzd231G0UUUGYUUUUAFFFFABRRRQAUUUUAFFFFABRRU2oWU2m38+nXLwtJbzNHI1vcJLGWUkEq6Eq68cMpII5BIpXV7Cur26kNFFFMYUUUUAFFFFABRRRQAUVNevYSTK2nW00UfkxhlnnEjGQIA7AhVwpfcQuMqCFJYjcYaSd1cSd1ewUUUUxhThPMsLWyzMI3YM8YY7WYAgEjuRuOPqfWm0UWTGm47MKdAIWmRbmRkjLDzHRAzKueSASMn2yPqKbRQ9gTs07XJtPuYbO/gvLnT4buOKZXktbhnEcygglGKMrBT0O1lODwQeajnjSKZ4o51lVWIWRAdrgHqMgHB9wD7U2ilb3ritHe2vz/Lb57/ACFKOEEhQ7SSA2OCRjI/UfnSUVJEk0kMgSQBEAd1aQDPO0YBPzH5ugycZPQGr/d+e3lv/lf5+o0pSdkR0UUVIhxMPkqqxsJAx3uXG0rgYAGOCOcnJzkcDHLaKkLRzGKLy44to2vJ8x3fMTubr0BxwOgHGckuMXKSiuv9dbf8D0KbTV3pby/r18/Ujop0IhZyLiRlXa2CiBjuwcDBI4JwCewOcHGC2kK2lwpxhdYVuCV2sxUAON2QATkZyByMEjB5x0OG0UAra3CilcoXJjUhc/KGOSB9eM0lOS5ZNXv/AF8hDoZ5rdzJbzMjFWUsjEEqwII47EEg+oNNoopWQ7tq19B08L28z28hUsjFWKOGUkHHBBII9wcGkQIXAkYhc/MVGSB9OM0lOhELORcSMq7WwUQMd2DgYJHBOAT2Bzg4wXF8rTavb8f+H9fn1HZSlZaX8/1/UbRToYXncohUEKzfO4UYAJPJI5wOB1J4GSQKbSFZpXCiiigQUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAU6cQrM620jPGGPlu6BWZc8EgE4Ptk/U02ijqO+lrBRRRQIsXg06X9/pyvCqpEpgnm8x3fZ+8cEIoC7wSFPKh1XL4LVXoopJWVg0SSXT5hRTp5nuJnuJAoZ2LMEQKoJOeAAAB7AYFNprYbsm7bE1lp9/qUzW2nWU1xIsMkrRwRl2EcaF3cgfwqisxPQBSTwKhoopa3J1v5BRU2n20N5fwWdzqENpHLMqSXVwrmOFSQC7BFZio6narHA4BPFQ0XV7BdXsTO9gbCOOO2mF0JnM0zTgxtGQmxQm3KsCHJYsQwZQAu0loaVXKhlAHzDByoPfPHp0pYxCUcyyMrBf3YVAQzZHBORgYyc88gDHOQ7RjHR3/AOH/AK9Nhwjrb838/wDhl8iS9sptPmWCd4WZoY5AYLhJV2ugcAlCQGwwBXqpyrAMCBDRRSV7aiV0tQooopjCnCF2ha4BXarBSC43ZIJGBnJHByQMDjPUZIIJrmZLa2haSSRgscaKSzMTgAAdSabSKSsrtafr/WtgooqWK1822lujcRIItoCO/wAzknooHJ4BJPQY65KgjaW4RhKbtHzf3asioopWUKFIcHIyQM/Lz0P8/wAapRbTfb+v6sSJRRRSAdOYWmdraNkjLHy0dwzKueASAMn3wPoKbRRQtENvmbYUUUUCCiiigAooooAKKKKACiiigAooooAKKKKACilVyoZQB8wwcqD3zx6dKSm1Gys/6/UAp00aROFSdZAVU7kBwCQCRyByCcHtkcEjBptOnhe3me3kKlkYqxRwykg44IJBHuDg1PUpL3G7fP8Arv8AoWdQ1KzvbaGC20C0tGixvmt3mLTfu40+bzJGA+ZGk+UD5pnH3QipWmkSVwyQLGAqjahOCQACeSeSRk9sngAYFEciIjq0CuXXCsxOUOQcjBHOARzkYJ4zghtKMVFWX9fiyYxjThyxe/q+vd3f3O3T0KdJGiIjLOrl1yyqDlDkjByBzgA8ZGCOc5AbRVDTVnoFFFFAiWe8luI0jkSIBOhjgRT91V5KgE8KOvfJ6sSWTSJK4ZIFjAVRtQnBIABPJPJIye2TwAMCpEvZo7CTTlSHy5ZkkZmt0MgZA4ADkblX5zlQQGIUkEquIalK3Qp1Kk2+b+rbdNPT8+joIJrmZLa2haSSRgscaKSzMTgAAdSabTpJnlREYKBGu1dqAEjJPJA+Y5J5OTjA6ACm01cHyJaf1p/ncdGYQjiWNmYr+7KuAFbI5IwcjGRjjkg54wW0UUxN3SQUUUUCCiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKkilhililNsJAhBkjkY7ZOeny4IBHHBz15HZbuW2lkH2S08pFyFBkLMw3EgseASAQMgKPlHGc1FQpe61bf+vl8vQtr2dT3Xe3Xp+K/Na9tbBRRRQQFFOjjR0dmnVCi5VWBy5yBgYB5wSecDAPOcAtoG00k+4VIXmuTFAsYJUbIxHGAWyxPOBljk9Tk4wOgAqOinHl5lzbf1e3yC8kml1CiiikIKKUlNgAU7snJzwRxjj86Sm1brcAooopAFFFOnme4me4kChnYswRAqgk54AAAHsBgUdR6W8xtFFFAgooooAKUMAhXYMkghucjrx6f/qpKKak47AFFFFIAooooAKVHeJxJG5VlOVZTgg+tJRTjKUJKUXZoAooopAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABU99pmpaY6R6lp89uzxpIizxFCyOiyIwyOQyOjA91ZSOCKgopq19f6en4b/h53T5rq2wUUUUhhUzvYGwjjjtphdCZzNM04MbRkJsUJtyrAhyWLEMGUALtJaGik1cTV7BSo5Q7lA6EcqD1GO9LPM9xM9xIFDOxZgiBVBJzwAAAPYDAptVGUotNaMqSV2lqgop0MaSuVedYwFY7nBwSASBwDySMDtk8kDJptILNK4UUUUCCiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigCSB5AJIoxH+8jIYuq8AENwT0PHbk9O5Bjoop+7bRajbbVn0CiiikIKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigD//Z\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(\"truecolor.jpg\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}