{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Coronagraph Basics" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This set of exercises guides the user through a step-by-step process of simulating NIRCam coronagraphic observations of the HR 8799 exoplanetary system. The goal is to familiarize the user with basic `pynrc` classes and functions relevant to coronagraphy." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Import the usual libraries\n", "import numpy as np\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "\n", "# Enable inline plotting at lower left\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will start by first importing `pynrc` along with the `obs_hci` (High Contrast Imaging) class, which lives in the `pynrc.obs_nircam` module. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "pyNRC log messages of level WARN and above will be shown.\n", "pyNRC log outputs will be directed to the screen.\n" ] } ], "source": [ "import pynrc\n", "from pynrc import nrc_utils # Variety of useful functions and classes\n", "from pynrc.obs_nircam import obs_hci # High-contrast imaging observation class\n", "\n", "# Progress bar\n", "from tqdm.auto import tqdm, trange\n", "\n", "# Disable informational messages and only include warnings and higher\n", "pynrc.setup_logging(level='WARN')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Source Definitions\n", "\n", "The `obs_hci` class first requires two arguments describing the spectra of the science and reference sources (`sp_sci` and `sp_ref`, respectively. Each argument should be a Pysynphot spectrum already normalized to some known flux. `pynrc` includes built-in functions for generating spectra. The user may use either of these or should feel free to supply their own as long as it meets the requirements. \n", "\n", "1. The `pynrc.stellar_spectrum` function provides the simplest way to define a new spectrum:\n", "```python\n", "bp_k = pynrc.bp_2mass('k') # Define bandpass to normalize spectrum\n", "sp_sci = pynrc.stellar_spectrum('F0V', 5.24, 'vegamag', bp_k)\n", "```\n", "You can also be more specific about the stellar properties with `Teff`, `metallicity`, and `log_g` keywords.\n", "```python\n", " sp_sci = pynrc.stellar_spectrum('F0V', 5.24, 'vegamag', bp_k, \n", " Teff=7430, metallicity=-0.47, log_g=4.35)\n", "```\n", "\n", "2. Alternatively, the `pynrc.source_spectrum` class ingests spectral information of a given target and generates a model fit to the known photometric SED. Two model routines can be fit. The first is a very simple scale factor that is applied to the input spectrum, while the second takes the input spectrum and adds an IR excess modeled as a modified blackbody function. The user can find the relevant photometric data at http://vizier.u-strasbg.fr/vizier/sed/ and click download data as a VOTable." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Define 2MASS Ks bandpass and source information\n", "bp_k = pynrc.bp_2mass('k')\n", "\n", "# Science source, dist, age, sptype, Teff, [Fe/H], log_g, mag, band\n", "args_sources = [('HR 8799', 39.0, 30, 'F0V', 7430, -0.47, 4.35, 5.24, bp_k)]\n", "\n", "# References source, sptype, Teff, [Fe/H], log_g, mag, band\n", "ref_sources = [('HD 220657', 'F8III', 5888, -0.01, 3.22, 3.04, bp_k)]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "name_sci, dist_sci, age, spt_sci, Teff_sci, feh_sci, logg_sci, mag_sci, bp_sci = args_sources[0]\n", "name_ref, spt_ref, Teff_ref, feh_ref, logg_ref, mag_ref, bp_ref = ref_sources[0]\n", "\n", "# For the purposes of simplicity, we will use pynrc.stellar_spectrum()\n", "sp_sci = pynrc.stellar_spectrum(spt_sci, mag_sci, 'vegamag', bp_sci, \n", " Teff=Teff_sci, metallicity=feh_sci, log_g=logg_sci)\n", "sp_sci.name = name_sci\n", "\n", "# And the refernece source\n", "sp_ref = pynrc.stellar_spectrum(spt_ref, mag_ref, 'vegamag', bp_ref, \n", " Teff=Teff_ref, metallicity=feh_ref, log_g=logg_ref)\n", "sp_ref.name = name_ref" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAFgCAYAAAC2QAPxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAB9C0lEQVR4nO3dd3xTZRfA8d9JJy0tpey9tyBLZSlTUaEO3BNcKG5xIS7wdSAiIrgYKk5EnBRwIAKKDMEFInvvXSgFOs/7R9LSlpYmbdJ0nK+ffJp78zz3nlxievrcZ4iqYowxxhhTkjj8HYAxxhhjjLdZgmOMMcaYEscSHGOMMcaUOJbgGGOMMabEsQTHGGOMMSWOJTjGGGOMKXEC/R2AMcZ3RKQOMBYoDyQDAcA0VX3Tr4HlQUTOBiYAUapaN5cy0cDrQH0gEQgGflLVYYUUpjGmCBObB8eYkktE5uJMaN5ybXcDxqlqSx+fV4F6qrq5AMfoBkw+TYLzPrBPVR9zbTcD5qtq5fye0xhTctgtKmNKtrOBeekbqjoP+MRfwXhZ9ve2Chjtt2iMMUWKJTjGlGxbgMdEJDx9h6qOABCRgSKyWUS+EpFPRWSRiPwqIvXSy4pIOxH5RUTmi8gcEWma6bUGIvKd67VFIjLMtf87V5HPRGSeiNRx/VQRucdV54SIdBORy13nnOs6T2cP39v9IlIh+3tzxREhIu+KyAJXfI+JU3MR+VtENrvKdRSR1SIyz7V9dvrrIvKoiCx0tUjl+p5dr/V27ZsnIrEiUt21P0xEPnfVWSAiloQZUxhU1R72sEcJfQA9gQNAHPA+0DXb68Ncr1VzbQ8FFrqelwP2AT1c232ANTj/MAoA/gP6Zyq7PdNxFaib7VwKPON6fj3QFrgRiHbtqwtszVS+G7D5NO/tTGArcByYCvTFddvd9fq7OG9xAZQBlgM35XRsYAAwL9u5k4Beru1Rp3vPQD0gHmji2r4HZ38ggLuBt13PA4Cl/v5c2MMepeFhLTjGlGCqOgeoDTyCM4GYKyLjsxWbr6q7XM8/AjqKSG2cCcNRVf3ZdayZQFXgHKAD0BDX7S5VPQxc40ZI37rKf6qqfwL/AO+LyAJgMlBLRNzqQ6Oq/7hiuAUIA74BZomIQ0QcwA3Ae66y6UnQLe4c2yVBVX9y1X+E07/n64FlqrrGtf0p0FNEqgEHgXNFpIOqpgJdPYjBGJNPNorKmBJOVROAScAkEekK/CwiI1V1g6vIoUzFD7h+VgNqAtHpt25c9gEVgHDgkKqmZDrPb26Eczjb9nTgTVUdBRmdk8PcemPOcyYBn+G8HdYcWIKz9WUlEOKKN3PsNd09dg6x1iT391wTaJ7tWm0BqqjqZyISCIxx3U4bDbztQRzGmHywBMeYEkxE3lbVQenbqjpfRA4AkZmKRWd6XtH1cxewDectmG6ZjhcJnADaA1EiEpj+C981immzq7XEndgq42xV+t61HVTA9/afiKzi5K21RKASsMpVpBKw3fU8CWcClC7KjVNuI5f37Hptmar2yRRfeeCIiFQEpqrqxyLSBvhJRFar6lxP3q8xxjN2i8qYkq2Xa04ZAFwtOAqszlSmk+tWCsDNwCJV3QrMACqIyFmuuuHAXJwJxBJgA85bM+lz0nwOpLduHAXCRORGEbkyl9jS+wad49q+0MP31kxELs/03prhnBNnqaqmAR/i7FuDiJQBrsbZDwmcrSsVRaSy63ZWbzfOd7r3PAU4R5zzDqUnb/Nxfsfei7P/EsAKnLesAjx8r8YYD9k8OMaUYCJyB86OvGk4f6mmAUNUdbHr9WFAM5wddZvg/GXdX1U3ul5vB7wKiOsxUlVnuF5rALyB85aSA3hKVee7XhsBXAwcAa7EmWycjzNJGJrer0dELsN5y2YdsAxnJ+clwP04J/prCsSq6lU5vLd+wCAgyBVbMPA/VZ3ler0sMMZ1jEDga1f86SOiRgD9cHY+XgXch7MP0ts4+9A0BRYD/VT1oBvv+QLgWZwTKqa53udiEekAvOCKMRKYrapPnPYfzhhTYJbgGFOKuRKcuqo6wM+hGGOMV9ktKmOMMcaUOJbgGFNKichAnH1ULhSRJ/0cjjHGeJXdojLGGGOM14hIL5z92/YCqqrDs71eD+fkmUuB1sCnqjrd63FYgmOMMcYYbxCRMJwd91uoaqKIfAm85Zp0NL3M28BaVX3NNXXC56rayNux2Dw42YiIhoW5Pc9YFqmpqQQE5H/0Z0Hq+6suQEpKCoGB+fso+TNuu2bF69x2zTxTkOtV0HPbNSs+dd2tn5KSQkpKyin709LSUFXJtKsjsEVVE13bv+GcJmFOpjJ7cM5LhevnH/kM/fT8vVZEUXs4HA7NrzvuuCPfdQta3191VVXbtWvnl3MX1+utatcsP+yaeaYg16ug57ZrVnzqFrQ+kKqZfocC1wHfZNq+Hfg4W5lIV8IzGlgA9FYf/D63FhxjjDHGeMteICLTdqRrX2aTgUmqOkVEKgHrRKS+uuab8hYbReVFMTExfqvvr7oF5c+47ZoVr3MXhF2zwj23XbPiU9cb9bNZBNQRkfSlUDoDM0Uk2rXUC0AtnMvBgHMtvDR8kI9YJ+NsAgICNDU11d9hFCvt27dn2bJl/g6jWLFr5jm7Zp6x6+U5u2aeE5E0VQ3Itu98nDOY7wOSVXW4iIwEDqrqCBHpAjwI/AnUA/5Q1Xe8HpslOFnllOAkJyezfft2Tpw44aeoirZdu3ZRrVq1vAt6IDQ0lJo1axIU5NH6i8WGfZF6zq6ZZwrjepW070ZffJeVFLl9J+eU4BQV1gcnm5x60G/fvp2IiAjq1q2LiORQq3SrWLEilSpVyrugm1SVAwcOsH37durVq+e14xYlAwcO9HcIxY5dM88UxvUqad+N3v4uKyny+E4+dWhVEWEtONmEh4drQkJCln2rVq2iadOmJeJ/4OJCVVm9ejXNmjXzdyjGmFzYd2Ppkdt3sogcU9VwP4V1WtaC4yb7H7hw2fU2pniw/1dLh8z/zrGxscTGxqZvFsnbU2AJjjHGGGM8EBMTkzHyauLEiUV2VI4NEy8Gfv/9d7p160anTp1YuHAhAF988QV169bl+uuvZ+vWrXz++ecZ288++yw33ngjw4cPz/F4f/zxB1dddRUvv/wy/fr1Y9GiRQBMnjyZdu3a0a1bN7p160b9+vWZPHkyAFOmTOGee+7hxRdf5O67787oVLh582b69+/PyJEj6d+/P//884/vL4gxxpD1u3HYsGEMGzaMIUOG8OCDD2aUOX78OK1ateKRRx45pf7zzz9PxYoVT9m/evVqIiIimDFjBgCvv/46gYGBjBs3DoAjR44QHh7O4cOHAfj444/p3Lkzf//9d5bjjBw5EhFh6NChPPnkk1xwwQXMmjWrwO87Pj6e2267jQEDBhT4WCWaL2YPLM6PsLCwU2Zq/O+//07ZV9ieffZZffjhh7Ps69q1q8bGxua4nZqaquXLl9edO3eecqwLL7xQv/rqK1VV/eqrr7RXr16qqrpw4ULdsmVLRrlLLrlE4+Pjdc+ePVqtWjVNSkpSVdXBgwfr6NGjVVX10ksvzTjW8uXLtVWrVt56y0XiuhtjclcU/h/N/t14/PjxLN+LgwcP1ptvvvmU78+5c+fq4MGDtUKFCln2Hzt2TG+55Rbt2LFjluO0atVKFy9erKqq06ZN0wYNGujUqVNVVfXvv//Wjz76KMf4AI2Pj1dV1Xnz5mmbNm0K8G6zxt+/f3+vHMtdOf17AwlaBH535/SwW1QeGh67kv92HvHqMZtXj+TZmBZePebRo0cJDAwkp3W1qlSpwr59+wDYt28f7dq1A6Bjx44ZZf7++28aN25M2bJlWbVqFRUqVMgYHli/fn2+++47HnroIdatW0ft2rUz9i9fvpz9+/fn+FeRMabkevn3l1l9cLVXj9k0uimPn/242+VTUlIYMmQIY8aMAeCjjz6ic+fOLF++nKNHj2aU27NnD5999hlDhgzhgw8+yHKMJ598kqeffppbbrkly/6+ffsyY8YMzjnnHJYsWcLQoUOZMWMGV199NT/88AO33nprnvHt27cvY5TWxo0bGTx4MJ06dWLFihU8/PDDtG7dmldeeYXhw4czatQoFi9ezL59+5g+fToBAQH88MMPvPXWW3To0IFDhw5lHHfw4MH88MMPXHPNNcTFxbF3717effddQkJCuPzyyznrrLPYvn07nTt35oYbbmDfvn08/PDDNGvWjI0bN9K/f3+aNGlyyr4uXbq4fe2LIktwsjmRnMK3f+/g0tY1/B3KKebOnZul6XX9+vWnlPniiy9YsGABP//8M7NmzaJcuXKnlHn++ee59tprWbt2LYsWLWLs2LGnlBk3bhzPPPMMAM2aNSMuLi7jf87ff/+dI0ecSV6XLl1YvHgx7dq14/fffweczbeW4BhjCkv6d6OqZnSG/e+//1i1ahUvvvgiy5cvzyiblpbG0KFDGTVqVMYtpnQffvghnTt3znF6ir59+3L33Xfz3HPPAc5+KE888QRpaWkcOHDgtN95o0ePJjk5me+++45Ro0YBEBwczDPPPEPbtm35888/eeGFF5g2bRqPPvoob775Jt27d+euu+6iT58+/P3337Rp04b+/fvz999/U7VqVSZNmsTevc4VEO6//36++OILnnrqKRwOB4MGDWLSpEncc889DBgwgEsvvZTU1FSaNWvGDTfcwG+//cbBgwe5//77OXHiBAcOHMhxX3FnCU42oSSxdNorbDt4D/f2OHX1dm+3tHiie/fuGf9zAKfc7wW48sor6du3L8OHD2fatGm0b9/+lDKXXHIJb775Jh07dmTFihX06tWL3bt3Z3wx7Nmzh8TEROrUqQNA2bJlmT59Oq+++irVqlWjVq1apE+G+OqrrzJ69Ghee+01ypcvT4UKFahZs6YP3r0xpijzpKXF29K/G1WVdevWAfD1118TGhrKiBEjWLBgAUlJSYwZM4YuXboQFBTE+PHjOXToEMePH2fEiBFcccUVzJ07lyZNmjBixAi2bt3KF198QVJSEv369eOcc85hx44dfPvtt7Rv355KlSpRr149vv/+e6Kjo08b3+DBgylbtiyPPfYYzZs3Z8mSJQQFBfHZZ5/x3XffceTIkYxW9XSNGzcGoFKlSsTHx7N//36OHTtG1apVAWeL+YIFCzLK161bF4fD2a22YcOGrFy5kpSUFP777z/+/PNPypQpk3GOvn37sm7dOnr37k2lSpUYPXp0jvuKO0twcvB80Pv0m12HyhFXcvVZtfwdTr6kNzXef//91KiRtTVq27ZtGbN1VqtWjcTExCyvv/XWW9x1111Z9okII0aMAGDo0KHcdNNNAOzcuZNHHnmEsLAw1qxZwwUXXEBwcLCv3pYxxuRKRDISgyeffDJj/4kTJzh69GhGC3j6H36bN2/m3XffZciQIQC8//77GXW+//77jD8YARwOBxdeeCGPP/44S5YsAaBPnz48/PDDTJkyxa34IiIiCAwM5NChQ0yaNIny5cvz5JNPsnbt2oxjZn4vmVWsWJEyZcpkzLa8cePGLK9v2bKFtLQ0HA4Ha9eupWXLlsycOZPZs2fz888/A2R0kl6xYgXXX399RmvRa6+9xi233HLKvpxa98GGiRdfDucl+SpkGK/O2gRn5fwPXJiWLVvGL7/8QlJSEosXL6ZDhw58/fXXbNmyhalTp3LmmWfy+++/Z2w3btyYxo0bc88993DbbbfxzDPP0KlTp4zjTZgwgSeeeIJWrVrx33//8f7772f8z5SYmMiyZctOGYH13HPP0bBhQ6Kjo6lVqxa9e/cGYOHChcyaNYv27dtz8OBB3njjjcK7MMaYUi3zd+O0adO46qqrTinz5ZdfZpSZMmUK1113HeC8xf/OO+9w/Phxnn/+eR566CHCw53z1Y0ePTrj+zQ6Ojrj+7Nv375s2bKFqKgowJngjB8/ntatW+cYX3oryMiRIwkJCWHbtm3079+fFi1acMUVV/DEE0+QmJhIUlISW7ZsYc6cORw8eJDDhw/z3nvv0bp1a5YvX57Rl2jy5MncfvvtnHXWWezevZvly5ezYMECatasSXR0NKNGjWLHjh3Ex8dz++23c+TIEV577TXuu+8+atasSUJCAu+99x6NGjVizJgxNG/enHXr1nHnnXdy9OjRU/blprgME7eZjLMJDwvThKERkOIcBn1HwHP0v/A8zm7TiuDAIpuolkirVq2ymYyNKcLs/9GiYfPmzQwYMIB58+b59Dw5/XsX5ZmMbR6c7ERgyDaOV3D2tZmY+gwVUvaydvcRDiYk+Tk4Y4wxJqu33nqLLVu2MHPmTH+HUqRYgpOTwGCCB3ybZdcZjs3sPhRPUkqRbY0zxhhTCo0cOZJNmzbRp08ff4dSpFgfnFwERFTiWL3eBDTqkbGviWxj5e4AWtYoZ+uvGGMMZBmabUqu4tidxVpwTiOs/+eEdLqLNJz/8waIEkYie46c8HNkxhjjf6GhoRw4cKBY/vIz7lNVDhw4QGhoqL9D8Yi14HiooWMncUcPszO1OgePJVEzuixRYTYs2hhT+tSsWZPt27efMoeLKXlCQ0Mz5jgrLsPEbRRVNuHh4ZqQkJBl33+LfqR5nUo5lt+p0TjKVqFyZAgOa6b1KhuhYYwxRVtRHkVlLTjuOE3eUpWDbIwP5VhSOPUrlfXJ6X///Xcee+wxkpKSGDVqFJ06deKLL77gkUceoVOnTowYMYLFixfz2GOP0alTJxo1asT+/fs544wzGDhwIAEBWRPspUuXMmbMGNq0acOaNWs4++yzueOOO1BVbr75Zho3bkxaWhobNmzg7bffJjw8PGN687Jly7JlyxZuu+02OnToAMCaNWuYMmUKZcqUYf78+QwbNoyzzz6bYcOGZRm2+OSTT3L++eczefJkxo0bR0REBABbt27lmWeesZVxjTHGeI+/V/ssao8cVxNf+IPqjj9PeaTF7854vn/7Ot1x6Ngpdb3F09XEVVUff/xxvffee0851rfffqtLlixRVdWkpCSNiorSffv2aUpKij7zzDMZ5e666y4dNWqUqqpOmTJFBw0apKqqBw4c0EaNGmlKSoqmpKToxRdfrKmpqaqqunPnTt27d29GzDnJbdXy7IrCSsXGGGNyh60mXkJEVke/H4occK51IsFlIcm5Qm0FIJlA0oJCPb9VVbUlXDTCy8HCM888Q/ny5XnxxRczWkvAuRZVZoGBgQQFBREQEJBlBuO0tDTKlnW2Ss2cOZMLLrgAgOjoaEJDQ1m5ciXHjh1DVRk3bhzHjh2jQoUK3HHHHRnHeOGFFwgJCSE1NZX77ruPsLCwXFctN8YYY7zFEhxPhFWAYOetRnUEI9nuXQWRQmrysYwy3ubOauKZhYWFERUVxa5du7IkOJm98cYbDB069JRVxzdv3szGjRsz1iLZu3dvlmNERkayd+9eDhw4wKJFi5gyZQrlypXjxhtvJDg4mAEDBnDVVVdRt25dwsPDeeutt7jvvvt49913s5wn86rlxhhjjLdYguMRBxLzOqQmIYEhzl07/8pSIgDQSk2RoDJeP7s7q4lnduzYMeLi4jIW1szu008/JSEhgaeeeirL/u3bt/PEE08wdepUQkKc77Ny5crEx8dnlDly5AiVK1cmOTmZpk2bZiRIXbp0Yd68eQwYMIAWLU6uvN6jRw9eeeWVLOfJvmq5McYY4y02D447gsOcP0Wcj/TkJrOqrUgIrgiAHtgIyf6fK+eFF17g9ttvz7H1ZtKkSezdu5ennnqKFStWsHbtWgA2bNjAE088wfjx44mOjubLL78EnIvKLVq0CICDBw9y4sQJWrRowTnnnMOBAwdITXXO8Lxly5aM1XwfffTRjPOtW7eOhg0bZokhp1XLjTHGGG+wFhx3lImGKo2dyU1uHAEQURUO7MeRlgT7VkHVVnD8IGkh5cARhMORv2Hknq4m/tdff7F3716aNGlyyqrgAN9++y0PP/wwbdq04ZtvvuHAgQOMGzeO2rVrc95551GjRo2MfjqNGjXiiiuu4Oqrr+avv/5i+PDhbN26lQ8//JCAgACio6N5+eWXefDBB6lUqRL79u3LuOUUGBjIAw88QOXKlVmxYgVvvvlmRgy5rVpujDGmaLN5cIqpnObBOe18LInxkJYKZaJITk0jaM8/mQ5WGRL2ArBDK6BhFakRVcamNXeTzYNjjDFFW1GeB8duURVUSASUiQIgMHsLjSu5AaghB6h5fA3r9xy2ac2NMcYYH7MEx4vcaZlplLaJhMTkQojGGGOMKb2sD46b1M0Vcw9rGOXk2Kn1A0KQ1EQAkg7tgMq1IMAuf26slcsYY4onEekF9AP2Aqqqw7O9/i7QINOuVkBbVd3szTjsN6wb0lfMrVChQp5JzhatQqCm0tyxNct+Ca8IR3YAEK1xpB5IJKByU5/FXJxpMV251hhjSjsRCQPeAVqoaqKIfCkiPVV1TqZiP6rqVFf5SGCyt5MbsATHLZ6smLvn0HEERWQ/AGk4cJAGZYGje7OUTduyl5TUNI46IigTGkKZoCLbGb3QZV651hhjTLHREdiiqomu7d+APkBGgpOe3LjcBrzni0AswXFDUFAQ9erVc6vs2r938Obc9fx4+GoAPq75LDd2aQJN28ChzfD6mTnW+zSlB19UG8ykAecQHR7srdCNMcaYApswYQITJkzI6aXseURlID7T9hHXvlOIiAPoDYzxQoinHt/6OmSV0zBxT63adYRm42sB8HHTt7jx2htOvjgs65IIa2tfQ+OtJ5PZlicm8drN59GreZUCxWCMMcb4WvZh4iLSExiqqj1d24OBmqo6OIe6lwM1VPUNX8Rmo6h8IDz4ZEJbNizbkg09si6LkHjBKzBoUcb2itDb+ejjd3lvwSafxmiMMcb4wCKgjoikT/nfGZgpItGu/jaZDQAm+yqQUpXgiEiYiIwUkXtF5CpfnScs5GRfmuDgbLebznuUqa0nZ2zWrhAGVZqTOuC7jH0fBL9M/PfDuXDMLySlpPkqTGOMMcarVPUYMAgYKyLPA8tdHYyHAHenlxOR1sBaVT3qq1gKvQ+OiDQAngf+BGoCB1T1uXweq6rrWGeq6lmZ9uc2RK0fsFRVp4nIN8C0fL+R08jcgiOBp/anCQ0/uTZURIizbEBU1g61DwR+TeeDK3nmi7G8dE0Hm/3YGGNMsaCqs4HZ2fY9lm37b+BvX8bhjxacaOAzVX1FVR8ArhWRdpkLiMhZItIw03aAiFydw7G6AN8Ckqls+hC1h1R1GNDKdU8QoBaQPhTK+8t9u4QGnbys4jg1hwwKPZngZKxPFVkD2vaHO3+B6515V3vHWkasvpDH3556yjGMMcYYk7tCT3BUdamqfpsthuy9eg8BH4tIS9d9vKlAhRyO9QVZe2tD7kPUALYBlVzPj+cUX2pqKgMHDmTgwIGZFxPzSObWloDAoFNeDwkre2olRwBcMhaqnQl1O2d5aeTeO9m8+Jt8xWKMMcaURn4dJu7qQf2Dqq7OvF9V17v6yEzDmYi8q6ofu3nY0w1R+woYJiJVgE9yqhwQEJDbULh8ceRwiyokPIpjGsLk1N4nb0hmFhwOT+2D9bPhs+sBqPt9f2i70/maMcYY4yfFZTVxvyU4ItId6A48mEuReOAYEAns9ODQe4GITNuRrn3pnZ8ey6mSrzgCTm3BKVMmjDMTJ5JMYM4JDkBgMFRvA8DO4HpUT9qE/vUJcs5A3wVrjDHG5CEmJoaYmBgAJk6cmOrncHLll1FUItIH5+Q+DwBVRaRjttcrA9OB/wE9gSdF5FI3D5/jEDWvBJ4PObXglA0JJNmd3DKyOgw7zLye09ml0SSs+4X49/qx9oUO3Dvuc6Yt22ZrNhljjDE5KPQEx9WheCrQAZiLs5Nwk2zFzgAeUdW5qnoYiAFOmUpYRLoCNwHVROQpESlzmiFqfhEQeGoiEx7iWYteq5rl2JBWjbLrY4nYOofGyat448AdtJ7emzve/YU9R054K1xjjDGmRCj0W1Sq+gfOlZlOV+bnbNvHyGEqZ1WdD8zPYf8pQ9T8JacWnMzDyN1Rp0IYK7VSxvbm6C7UPbiARo4djN92KbePHsaN11xHz2Y2+7ExxhgDpWyiP38IyCHBKVcmiPObV+HDW8926xgRoUEkBTo7F/+Y2o59l3wEt/7gPL4o7/Msb374KS/OWmUTAxpjjDFYguNzAUGndjJ2OISJN7fnvMaVcqiRs5Rg5wzX8YTRtGoE1O4Aww7DXQvQ4LJ8FTKMW5dcxDXv/Ma2g8e8Fr8xxhhTHFmC42MBOYyiyo+0IOddveAABxGhmY5ZtSXSf7rzqRzi6/19GfLGh8xfuy+nwxhjjDEFEhsbmzFfHEV4mLglOD4WFOidf3tHiDPBCQrIYcmGGu3g8S0Q6lyp/JO0x1n4wdO8++tGr5zbGGOMSRcTE8OECRPS54yzYeKlVYDDO2tIhYaEuI6Xyz9ZmSgYshWu/RSAJ4KmUP/HAdw05hsSU4rs588YY4zxCUtwfCzHFpd8CAl0HicgII9/sqZ94JlDpF34Mt0D/uGjuP6Mfu9TUlKt87ExxpjSwxIcHwvMrcXFQ8GuRMmt4zkcODrchZarBcATO+9l8vvv2KSAxhhjSg1LcHykV+JI7kx60Gu3qIIDnf9Uud6iyoHEjIEKjQC4ffsTLH57kFdiMcYYY4o6S3B8ZL3W5Ie0swnK65aSm7ZWvYA5qW2YVXGA+5Ua9oL7lnGi06MAdNw7hdQfnvJKPMYYY0xRZgmOj3mrBSewTCS3JT9KXLDnsxWHdjq5QGfAonGQZp2OjTHGlGyW4PiYtzoZlwl2Djd3SD6OV7YyKff+zbOOe53bX98Fw8qx7/Vu/Ll2s/XNMcYY4zabB8cAEOilW1Tpt7ry2yAUWLEeidXPcW6s+ByASof+ou2nZzJ+9NOs2R3vjTCNMcaUcDYPjgEg0Eu3qNIbWQpytErV6+S4/674cRx/qyvvzf3XWnOMMcaUCJbg+Ji3EpxUV+KRr1tULpWjozKe35N0P3E3fAf3/QlAa8cGbp3fmeueHGVrWRljjCn2LMHxMW91Mk5vWXEU4HhVI0MznleoVJWoRp2gQgN44B9SqrQC4LPg55nw6hO8NW99wQI2xhhj/MgSHB+TArS4ZNapQUUArjmrVr6PUTkiJON5RNmyJ18oX5eAgT9nbP4vaDLXzO3GjW/OISExJd/nM8YYY/zFEpxiolZ0GJtH9OGsutH5Pka5MidXIY+MjMjymgQEORfsfNS5QGcFiefjff14bsxYVmw/nO9zGmOMMf5gCU4pEhEamPG8TJmIUwuUiYLwCnDeoxm7Xj4+nL0TLue9X9aRlmYdkI0xprSzYeKlXNmQwLwLFbKI0JMtOAEhZXIv2OMpGLQQXGtZ9XT8wbk/Xcpz4z/mYEKSr8M0xhhThBWXYeJiw4KzCg8P14SEhAIfZ/fhE+w8fJy2tct7ISovGlYOgKnd53JN17Z5l09NQX95BZk/AoDRcjMdrnuKTo09n1HZGGNMySIix1Q13N9x5MRacHykarnQopfcZBIU7ObnMSAQ6f4E9H0NgMH6IZ0+bcx702fbLStjjDFFliU4pVRA6GluUeWk/a1QuUXG5q1/XslHY5+0W1bGGGOKpKLXUcQUitDgYM8r3fwNJOxHk4+h7/amf9ybzHz1P6KueZvOTap5PUZjjDHFj4j0AvoBewFV1eHZXhfgPtdmXSBKVW/1dhyW4JQy9yTdT0zAIsoE5aPje9nKULYyAsg9SzjxQT/6xM+FKU2ZcdYH9Ln4Uq/N+2OMMab4EZEw4B2ghaomisiXItJTVedkKnYjEKeqH7rqtPJFLJbglDIz0zowM60DUwMLeHeyYkNC75yDjmuDJMbTd2l//tx6JW3unIA4iuyoQWOMMfmQadRUdtnziI7AFlVNdG3/BvQBMic4NwDfi8j9QFVgkpfDzTEwU0qE5qcFJ7uylZAntpNyaBvH3jyPtnu+YN+Y/6j00EKwlhxjjCkxMs17k4WIZJ/uvjIQn2n7iGtfZnWASFV9TkQa40x2mqmqV4ecWyfjUsorCY5LYPlaRDywiDUhZ1DpyH8wPAqSbMFOY4wphfYCmWeSjXTty+wIsARAVde6yuR/HaJcWIJTSoUGefefXiKqEnzjlJM7lr3r1eMbY4wpFhYBdUQkffHDzsBMEYkWkUjXvjlAfQDXvgBgt7cDyfctKhE5E1iuNlNgseTNFpx09WrVZuiZC7j+75totvo7AlKTIG4rnHkd1O7g9fMZY4wpWlT1mIgMAsaKyD6cecIcERkJHARGAC8DI0VkKNAA6K+qJ7wdi0cJjojcDLQGlgHzgQHA+94OyvheaKBvOgKf36wKS/9ozBlbf4Ctvzl3/jGZ36MvocntkygXFnL6AxhjjCnWVHU2MDvbvscyPT8M3OnrOPJzn+J54BDwBFDBu+GYwhIa7Ju7ky1qRLJLT654/r/kGwA4++B0yo2szPzZ3/rkvMYYY0xmnv6WO4Bz0p7vVPVeVR3li6CM7wUH+CbBqVQ2hISgk3nv2vr90cGrM7bPXdCfT199gLhjNgOyMcYY3/H0t1xv4GsR+UJEHheRs30RlPGdxlXKAvhsQj4RITAsMmO7Te3ySGQ1eHQjiTfGcjikKtfHT0ZfrseMuQuwLlzGGFO8xMbGZh42XmQnPvNoNXERuVhVZ4lIKHAW0FJV3/JZdH7grdXEi6rDx5LZefg4zapF5l04n0aMe4MhB54E4LsrVnNRy0zLOBzZCaObZWx+lNabS4d8SGRYqM/iMcYY4xslaTXxIBE5S1VPqOqvJS25KQ3KhQX5NLkBKBN28rNeKSJbp+LI6hyp2TVj8ybHD0SOrMKff//l05iMMcaULp4mON2BG0RkhohMFZF7fRGUKd7CMiU4FcueOmoq8tZv4JmDHL37HzZU7AFA22+68cX7r3LkuPXNMcYYU3CeJjhfANNUtS/OIeKLvR6RKfaCQk7ebqqYvQUHwOEARwBlK9elwV2fsa2sc521K7c8x0+v3cbewzYLsjHGmILxNMFZjnMOHIBqQJxXozElQlBIWMbz8OA8+p8FhhBY55yMzX5J06n8WjXenziWE8leXZbEGGNMKeJpgtMcmCIiXwNXA128H5Ip7gKDy2Q8d2e0VrVyp3YwvmXH0xx6viHbdu3xamzGGGNKB49mMlbVxUA/EQkCYoDDPonKFGvBoWXyLpRZ/W6wcBz0j4XqbdD43cgb7akmB2F8Y6Y3f42+V96Cw2ErlBtjjL/FxsYSGxubvllihom3BsoAS1Q1TUQuUtXvfBWcP5T0YeKFYebS1fSZ6brtNMzNHPj4IShT/uR2ajL8r2LG5kdRg+jZ/xmqlw/LobIxxhh/KEnDxPfgXBn0YxGZA9wkIu1EpMhmcKbwBQbnIwnJnNwABARlPN1bvg03xb1N9der8ctP0wsYnTHGmNLAowRHVXep6ihVvR64AHgdOB942xfBmeIpONi7C2pWvm8OCdWcq5Gft+Amfh11DXsPHvLqOYwxxhR94sE0/PlekEhVU1V1iaqOUNWB+T2OKXmCgwJ4J6UvVyY+450DOgIIH/g9Kbf8AMC5R7+n/NgGLPnpS+8c3xhjTJElIplzjBYi8r479TzqZGyMO4IDHYxIub7gB7prAaS6Jv4TIbBOB7hrAXFLP4c/P+CcBbeycvn71Bv4MWFlowp+PmOMMUWGiEQCUUBTEant2h0PJLpT3xIc43VeW6m8assc90XFtCSp423wxpm0OPIrO189m5WXfclZZ+ZQ3hhjTHF1Oc5JhesCrQEBUoDv3amcrwRHRO5X1bEico2qTs3PMUzJFRzopQTndOeoWBeu+ZgDCz+g+rbZVP+6C6t/PIua982ibGiwz89vjDGlVWENE1fVD4APROR8VZ3taX23homLyMhsu84Gfgc6qOp5np60KLNh4gW3Yd9Rer46H4DNI/r4/oTDymU83eSoTUC/8dQ+o5Pvz2uMMaWcP4aJu9u44u6f2r2BSkCC65Hs+mkrI5pTeO0Wlbv6joHKzdnVdACV0vZRedolfDdhKAmHDxZuHMYYY7xORDaJyEbXYxPwljv13P1NdA5wAPhHVYcDP7l+unUSU7qEFMItqiza3wJ3L6Lata9z/M7fCZVkLtr5Jglj2rP29x8LNxZjjDHe9qKq1lfV+kBX4EF3Krn1m0hVT6jqI0BZEXkeCHLt/yqfwZoSrDD64OSmUrXapHZ/iiMVziRUT9Bg5tX8N+FW9MQRv8VkjDEm/1R1YqbnW4HapymewdO1qD4WkWY4ezYbkyN/JjgAAV0fJbLroxzZt41V793KGTu/JO6VeZS97j0CG3bza2zGGGM8IyLvZdqMxM3GGY9/E6nqKlV90dN6pvQo9D44uYisVIsWj/7Arw0fIyr1AIEfX0rq+nn+DssYY4xnBPjA9XgRuNKdSkXjN5EpUQJdCU6taA9XFfcBcTg498Yn+bXp0wAEfHwpLHjNz1EZY4zxwF2qOh9Yqap/qmqaO5U8Wk28NLBh4t6xcP1+GlYpS+WIUH+HkuHjD97mnA3jaBiwGxmyBUIi/B2SMcYUO5nnwZk4cWKiqvr0i15ELsDZehOOcwR3f1XNcwSJJTjZWIJTcsUdS+LBl8YwOeAF6DUcfnoWqreF66dC2cr+Ds8YY4qdwpgHR0S+xtmKs0dEqgETVDUmr3oe3aISkaszPW8tIqM9D9UY/4gKCyasfgeSCXQmNwA7/4RRjdj742toWqp/AzTGGJOTZaq6B0BVdwFLAETktImVp31wmqY/UdW/AbfugxlTVDStXY31adUAWBzShUQNAqDywmHIc9Ec/Pcnf4ZnjDHmVDVE5FYR6SYit+KcsuY8YMzpKrmV4IjIA67ZAx9Kn0lQRNZhi3WaYqZOhTD2q3Nphz8TonmxzMNZXo/+4gq2v9SOEwe2+CM8Y4wp9kSkl4i8JSLDROTZHF4fICKLRWSe63FTHofsCJwL9Hf9rALcArQ5XSW3EhRVfR14XUSuUtVp7tQxpiiqHlWGeNfHPk7DuaFPT/hyBAB/XDSDdt/1pWbiehjXitlt3uD8mBvAYYMNjTHGHSISBrwDtFDVRBH5UkR6quqcbEWvVdXNbh72flX9NYdzdT5dJY++ubMnNyJiE/6ZYqVauVAUASDVEUzDFmc5X4ioRrtzzoWbv2VbmWYAnP/XvWx/qQ1xe7b6K1xjjCluOgJbVDXRtf0bkNOqy/eKyCMi8oyIRJ/ugKr6q4iUF5Harscw1/7fTlfPo1tMItITeAmoiHPinUjga0+OYYw/VYoIYY3reUhIKA6HwGObICDYubN+N2o9vpjEr+8ncPnH1EzeDG+35K8OY2h9QX/EWnOMMaXQhAkTmDBhQk4vZc8jKgPxmbaPuPZlNh+Yqar7RORiYBrQM7dzi8i7wNnAPpy5R21gWF4xe9qH5jqcK4sPBF4FHvGwvjF+FRIYQPpEy0Ehrqkbwk794yHk8rFw+VgOvtWb6L2LabP4QVj8IJuunUe9pqe97WuMMSXOwIEDGThw4Cn7RSQl2669QOZJxiJd+zKo6qZMmz8D00UkQFVzG8oarqotM52zmzsxe/rn6BpVPQQEqmoKUN7D+sb4XaCrFSYkJO+ZlqPv+o4taSf/+Kj3WTe2jjkfTbMBhMYYk4NFQB0RCXFtdwZmiki0iEQCiMhLIpLewNII2HSa5AZgqatvTzq3cg9PE5yuItIDCBWRSUAHD+sb43eBDmcfnJBQN5aScDj4tO1nPJh0d8aQ8tpxvzPtmcvY8s88H0ZpjDHFj6oeAwYBY0XkeWC5q4PxEOBuV7HdwNsiMhQYCuQ4iso1YnsjcD+w37W9CZjkTiwezWTsmlQnDec9sNuBH1R1zelrFS82k3HJt/TFXpyVtJSP6o7gpgGD3KqjquyIO07Zg/8R9dHJW8WL6t1Lh+uGIsE+ncjTGGOKJF/OZCwig1T17Rz236aq7+ZV39NRVAmqelxVj6nq2JKW3JjSwSHOFpyAgAC364gINcuHEVWlTpb9HTe9gbxYnU1LYr0aozHGlHbpyY2IdMm2P8/kBmyiPlMKue5QERQgnlcuc/LWb3y1zkTsco5SrPfdjfw7pxVVrhlLpQbWCdkYY7xonIj8lWlbgQ3AOFWNz6WOx31wsshrHQhjiqITAc4O/hrkRh+c7AICocoZAERcOjLLS2ckLafSR91Y8+4dpKVkH1hgjDElQ2xsbOZRVe43heffYpzz6Xzk+nkA2Aa8crpKbvXBca35kJObVPUOz+Is2qwPTsn36LuzqLTxa8p0f4T7ejXJ30FSEiEwBBIOQFAZDi6bxvHqndCPr6BmyhbiHFHEXf89dRs2827wxhhThBTSauJDVfXFTNtPquoLIvKYqo7MrZ67t6heA5YDNYFQYCNQvyABG+MvhwIrMy31MoYEFeAObaBrBGR4BQCiO/UHQIf+xcqpT9NizZtEfdyBDRW6U+/ur3AE2ASBxhiTT2eLSCXXxIBVODmCu+rpKrn7DX+/qv4mIo+o6qj0nSLyTD6DNcbvgn2QdIgjgBbXvUjcv92J+uJKGhyYC/8rz+5Wg6h6+Usg+ej3Y3KkqiSmJnIi5QSpmkpyWjIHThxgR/wO1sWtIywwjPCgcEICQlAUVeV4ynHik+KpHFaZOpF1aBDVgOCAYMoE5uN2pTGmsLwNLBeRMsBx4BYR6YBzuHmu3F1sM329h3rZXqrpaZTG+Ft6jhEU6LtWlagzzkeb7mP92EtodGQRVZe/zZGVU9E75lCuqjV+no6qsjthNysPrORQ4iFS01IJcgQRGhjK0aSjzNs+jz/2/MHxlONeO2fNsjVpENWA1pVbE1M/hirhVbx2bGNMwajqDyJSHecyUfv1ZN+axaer52kbfYqIzATWAY1x9mI2plhJ/18jxMe3jSQwmEYPfUfChoUc/PZJasX/Be+0YVO966h74zgkIMin5y9O0jSNxTsXM2/7PBbvWsymw5tyLVupTCUua3gZkcGRhAaGEhoQSoAjgEBHINGh0QQ5guhYrSPJackcTT7K8ZTjnEg5QXhQOGmaRpXwKuxN2MumI5tYH7eeQycOsS1+GxviNjB/+3zG/TWOjtU6ckHdC2heoTkNyjUgyP6tjPEbEamdabOWiNylqkPzrOfJRH+uE10MtAD+U9WZnoVZ9Fkn45Lv9g+W8dOqPbx2zZlc3qbwGiE3/T6TerOuByCVAPZfE0uVZp0L7fxF1cr9Kxny6xA2H9lMmcAytKjQgq41u9KuSjsqh1Um0BFIcloyx1KOERYYRqUylQhw+GbgxtYjW/lm/TfM2jSLHUd3ABASEELzCs05p9o5dKreiTMqnkGQwxIeY6DQOhlvAjbhnGS4FnBEVdvmWS8fCU4roBKwBtihnh6giLMEp+RLT3DevL4tfVpVK9Rzpx3exf53Lqby8Y0A/FPjWprf8iZBgaVzSqotR7ZwVexVRIVE8VC7h+hZuyfB6Su7+5GqsvLAStYdWse6uHX8s/cf/j3wL2maRlhgGO2rtuecqufQvEJzWlZqSUhASN4HNaaEiI2NJTbWObnpxIkTE1U11JfnE5HzVXV2pu07VXV8nvU8XKrhUeBiYCswGbhQVR/3PNyiyxKcku+2yUuZs3ovE25qxwUtTtsJ3zeSjnFo9suUXzomY9eaCz+lSYc+hR+LH6WkpdD/u/5sPrKZLy/5kqrhfvi38MDhxMP8vvt3luxawpJdS9h8ZDMA1cOrM7LrSM6sdKZ/AzTGDwqjBSeHc76uqg/kVc7TTghlVbU7zttTc4HEfEVnjB+lp/TirxFNwWGUv3gYevVHGbuafH8920adS9LWpf6JyQ8+WfUJy/cv5+mOTxf55AagXEg5zq9zPk91eIrYy2P58YofGd5pOA5xMOD7AazYtyKjbEpaCgt3LmTGxhnsP77fj1EbU/ylL7qZabHNOHfqedounn7jW7NtG1NspLda+nXAtgjS/BK4ajIp3z9JSsIhah1dDu/1YlPtK6h7/WtIaDl/RuhTCckJTFoxic7VO3Nh3Qv9HU6+VCtbjX6N+tGzdk8u+vIipqyeQstKLQEYuXQkU1ZPASDIEUTf+n25qflNNCrfyJ8hG1NcjXDnllR2nrbgpIrI90BfEfkCsHs5ptg52YLj1zCcWlxO4MP/EfLwv+ypeymLQzpRb+uXyIjaHPhhJKSWzCUfZmyYQVxiHINau7eae1FWLqQcF9W7iNlbZhOfFE+apvHj5h+pElaFDy/6kH6N+vHdpu/oN70fA38cyMIdCylhXReN8SlVHS8iTUTkChFx+68ET1cTfxYYDUwHxmeeOtmY4qJCuLNDaJVIn/aL84iEV6DKgA85+/FZbKjm7ItTYdELbH39Ak4c3ufn6Lxv5qaZNIxqSKuKrfwdildc3uhyTqSe4PvN37Ni/woOnDjAQ+0eok3lNjzV4SlmXzmbB9o+wPq49dz5051cFXsVMzfOJCWtZCawxniTiNwBfAH0B74SkdvdqedRgiMinwO/u2Yz3u3aNqZYeaZvc964vg1n1Ch6t4AcDqHBHR+R0Phy9gZWp/aRPwh9rSG7pz18cgKfYm53wm7+2vsXF9e72H/9oLysRYUWNIxqSOyGWH7b8RuC0KVGl4zXo0KjuL3l7Xx/xfc81+k5ktKSGPLrELpO7cqYP8Z4ddJCY0qgJqraUlUvUdWWOKeqyZOnt6iqA7+KSCtVXQH842mUxvhbubAg+raq7u8wcucIIPz6yVR+ahWrO48BoOrKSaQOj2b/4in+jc0LFu1cBEDXWl39HIn3iAgX1LmAv/f+zewts2kS3YRyIacm0MEBwVze6HK+ufQbxnYfS7PoZrz777tc/u3l/LbjtxyObIwBDmbbPuROJU8TnNnA9cAnInIDJ7szGGN8oOn5t3D03n/5vcq1BJBGxe/vYs9bfeBo8b1ttWT3EqJDo2kUVbI63Har1Q1FWR+3nraVTz8HmUMcdK/dnUm9J/Fe7/cIcgRx10938fgvj3Pg+IFCitiY/ImNjWXgwIEMHDgQCmewUUURGSsiD4rIOCDKnUqezoMzAngKiASmAI1VNfv6VMWazYNjiqrdq5dQ9bMLMrb3Nb2JShc/CZGFO1lhQagqPab14KyqZzHyvJH+DserVJW759zNgh0LmHTBJM6pdo7bdZNSk5i0YhITV0wkPCicIWcPoU+9PiXmFp4puQppJmMHcDvQCuedo3dVNS3Peh4mOA8AM1V1veuEj6vqS/mMuUiyBMcUZcmpafzywTCqb/mWxmwhQJSj3YZTttuD/g7NLRsPb+TSby5leKfh9GvUz9/heF1qWiq7j+2mRtka+aq/MW4jTy98muX7ltO1ZlcGtxtM/ShbnNUUXYWU4CwG7lbVPz2qZ8MVs7IExxQHCYkpfD49lltWDgDgRFA5gh9eiSM0wr+B5WHWxlk8/uvjfBHzBU2im/g7nCIpNS2VT1Z9wht/v8HxlOPUjqjNeTXPo2ftnrSp3MZn63AZkx+FlOBMUNWBmbYrqGqe93LdSnBERFRVs63oCeDWip7FiSU4prhITVPGvzWKu/c/n7Ev6bJJBLeIgaCiMwQ+s9F/jOaj/z7i9+t/txW683DoxCG+Xv81c7bMYcX+FShK+ZDydK7Rma41u9KyUkuqh1e321jGrwopwXkaWAqswtn3915VfSzPem4mOL+r6tnZVvQEqK2qDfIfdtFjCY4pbjQlEXm+csb28YBI9OZvCKvTzo9R5eyu2Xex//h+vrjkC3+HUqwcSz7GL9t/Yf72+SzYsYC4xDgAqoVXo1utbnSp0YWO1Tpa0mgKXSElOLuA1Zl2uZV7eNoHJ0ZVYzNtX6yqszyKtIizBMcUR/p8VSTlOOtDz6DhiX8B2BvVmsp3xUJopJ+jO6n7593pVL0TL3R5wd+hFFspaSmsPrialftX8tvO31i0cxEnUk8QFhhGm8ptaF+1PZc3vJwKZSr4O1RTChRSgjNAVSdn2u6lqj/lVc+jtagyJzcuNT2pb4zxDbn5W1j5NQ0vfAmGRwFQOe5vGFGLuKu/Iap5d7/GB3Dg+AH2H99P4/KN/R1KsRboCOSMimdwRsUzuKbpNSSmJrJk1xJ+2f4LC3Ys4LedvzH+n/Hc0OwGbmt5GxHBRbtflil+YmNjiY3NSAd83iksc3Lj4taH2t1bVAc5dfVOASJVtUT9mWAtOKbYWz6NpH3rCP715Yxdh6LOIOr6d5HKTf0W1u+7fue2H29j/Pnj6VS9k9/iKOk2Hd7E+OXjmblxJtGh0TzU7iEubXCp9dUxPlFILTi3AQ8AZfEg93B3or97VbV+tkc94L78h2yM8YlWVxHccyhxre/M2FU+7l+OTuxD8qFtfgtrW7zz3LUjso9VMN5Ur1w9Rpw7gs/6fkbtiNo8/dvTDJoziN0Ju/0dmjH5dT3QPVPuMdidSm4lOKr6aS4v2QIqxhRRUZdlnUgvInk/x17vyPp3bycxufAXedwWv41ACaRqeNVCP3dp1KJCCz646AOGnjOUP/f8Sb/p/ZixcYatZG6Ko8XZhoWvcKeSp52MewIvARWxW1TGFHl6YAMyri1pgaEcSg6igsSffK355Ui/CRAYXCixDJ43mDUH1zCz38xCOZ85aeuRrTy54En+3vc359c5n2GdhhEZXHQ6n5viy5e3qETkPdfTykA5YL1ru6Wqts+rvqdrUV0H9AbGA42AVzysb4wpRFK+HjTti+P6z7MkNwDy39fwfCVY+m6hxLI9fju1ImsVyrlMVrUjazP5wsk82PZB5m6dy/Uzr2dD3AZ/h2VKKBHpJSJvicgwEXn2NOVuEBEVkbK5FEkDPsCZazwFTHY9/nYnDk8TnDWqeggIVNUUoLyH9Y0xhcnhgGs/gfonV+5Ou+aTrGVmDobdbrX45puqsi1+G7XKWoLjLwGOAG5reRuTek8iPimea2dcy9TVU+2WlfEqEQkD3gEeUtVhQCvX3Z/s5ZoBzfM43N9AoqrOz/wA7nYnFk8TnK4i0gMIFZFJQAcP6xtj/MzRoDtaqRm3JT3MWSfedO5c+TWkpsDyz+HYQa+fMy4xjqPJR6kVYQmOv7Wr0o5pMdNoW6Utzy95nqtnXM3qg6vzrmhKtQkTJtC+fftTHpw63UxHYIuqJrq2fwP6ZC7gSoIeA4bncdouwClLMqhqkjsxezQPDnANkAosxrmy5wce1jfG+EvrG2H5VAgOR+5ZzLikFJ76+l+S/wsgcNUMZN2PGS05GhiK3Po9VG/jlVOnj6CyBKdoqBxWmbd7vc2MjTMYvWw0V8VeRb9G/Xik/SM2b47J0cCBAxk4cOAp+0Uk+4iFykDm++FHXPsyewH4n6om5TF9wVZVXZfDObur6ty8YvY0wRHgIpyT7BwG7gd+9fAYxhh/uOxN58MlLDiQLo0qMu/f1py//48sRSXlBEzoxuGL36Zc815QNvv3k2fShyhXL1u9QMcx3uMQB5c0uITO1TszccVEpqyewtLdS3mt22u2EKopiL1knYgv0rUPABGphbN7y9WZkpvBIjJLVZdlO9ZFIlIxh3O0ArzeyXgG0BWoB9QFoj2sb4wpQppUjWC1nmxVuT/pXl5r9wODkh4AoNysQTCqEfu/fQrSUvN9nn3H9wFQKaxSwQI2XlehTAWGnD2EDy78gMSURG6cdSPTN0z3d1im+FoE1BGRENd2Z2CmiESLSKSqblPVAao6QlVHuMqMziG5AdgN/JHDY6c7gXjagrNeVR9M3xCRRh7WN8YUIbWiw4jXMhnbjauW5d6YDmzq2JJFSxvTcck9AFT8axz8NY6UTg8Q2P0JCCqT2yFztPfYXgIdgUSFRHkzfONFrSu3ZmrMVB775TGeXPAk/+z9h8fPfpzggMKZRsCUDKp6TEQGAWNFZB+wXFXniMhI4CAwAkBEKgHps5E+JiLjVXVHtsP9papvZtuHiOS5DhV4Pg/OVTinSk4fX3iTqt7h9gGKAZsHx5Q2Lz37IE/I+wBMOHcBA3u2zPL6lnmTCV7wCtVStgOQEhJF4IN/Qxn3B1EO/XUof+z5gx+u/MFrcRvfSElLYeyfY3l/5fu0rNiSMd3HUDmsYLcoTcnl43lwfgNeUtUZ+anv6S2qW4BLXD9vAbzTA9EY4zdlwpzfTb+ltqBetVN/kdXpNoBqT/7LP43uZb62JTAxjqQ3O8OCMbB1sVvn2Htsr92eKiYCHYEMbj+Y17q9xvq49Vw38zobZWX85XpgYX4re5rg7FfVy1X1FlW9BShRrTfGlEbBoWEAKNCwci7zbYlw5g0v0OihWdzreJLA+J3w07PwXm8SZj4JqcmnPcfe43utFaCY6VWnFx9d9BGCcPN3N7Nk1xJ/h2RKGVXdoqr5nrfC0wTnHxHpLiK1RaQ2EJPfExtjigZHUCgAilCz/On71lSPKkPdcy6lb9IL9Ex8hf0aSfjSN0h4sT6Hv7gPUhJzrLfv2D5LcIqhJtFNmNJnCjXK1uCBuQ+w5uAaf4dkioDY2NjMw8YD/B1Pbjztg7MLyNxWWVtVG3g9Kj+yPjimtJk46W3u2D6EhbSi07C8Z31QVT5buo1GlctS5cQGNs/9gBO7V3O+LAUg7fZ5OGqevHt9LPkY53x6Dg+2fZDbWt7ms/dhfGd3wm5unHUjaZrGxxd/bMP9TQZf9sHJ5XxRqhrnTllPW3CeUNXu6Q9O9oA2xhRTgUFBAAQ5TjvhVgYR4bqza9O+bjS1mp7FuYPeoPUjM1kb6uyc7JjUjZTkkxON7j3mnALDWnCKr6rhVXmn1zucSDnBoJ8GcTjxsL9DMqWIa12rDiJyD/CXiIxyp56nCc5dItI2fUNV3RqqZYwpusKc+Q0SGJTvY1SKCKHO4DkccC1Pt+yLV1nx7t0cObgnYw4cS3CKt4blG/J6j9fZGr+Vh+c/THLa6ftdGeNFW1R1MXAT0ALnRMN58jTBWa6qf6ZviEgFD+sbY4qYsuLsN5PiCC3QcUKCQ4i+xzkMvMOaEbTc9gn737iAHUecLTgVQu3rorg7q+pZPNvxWZbsWsJLS16yhTpNYakoIucCG1T1mLuVPE1wtonIhSJSx9XJ+HEP6xtjipjdlc9jfmorPo4oeP8YqZR1iv/6aZtZ/r2zNbl8qPvz5pii67KGl3HrGbcybe00Pln1Sd4VjCm4TcDrwMsi0hdwa1E7T2cyvptsnYxxrghqjCmmJKQs/ZOHcE6gb1ZeqZS2HogicNsKqN0G9qyEGu0h0GbILa4eaPsAmw9v5pVlr1A7sjbn1TzP3yGZEkxV3wLeAhCRre5O/OdpgvOEqk5O3xCRXh7WN8YUMUGBnjbkuunRDSSPacPBgAAiU1Mp9+mlKILgvK2R1GM4wec+AKdfTdgUQQ5x8NK5LzHg+wE89stjfHTRRzQqbyv3lBaxsbHExsamb/p8mLiIvAV8CLQDHhGRL1T10bzqefTNpqqTRaSXiDwsIj2tk7ExxV9IgI8SnPCKBJ77IHEOB+XT0gAykhuA4J+fheFRxH9yM8Rt9U0MxmfCgsIY22MsYYFh3DvnXg4cP+DvkEwhiYmJYcKECUyYMAEg/6vwui97J+Mj7lTy6JtNRJ4GBgN1cGZRT3sapTGmaAkK9F0LitTtzKGAAMqnpmXZf2+9WUwPvRSAiHXfcmJsBxJ3rcpzRmRTtFQNr8q4HuM4eOIg98+9n4Rkm0PM+EShdDIOVtWLVfV+Vb0ICPOwvjGmiAlyteD45E5R7Q4cqnYGUfW6Mzn4uozdb/TvTO+H32N82294u8pwQtMSCBnfgWMjmpD03ywfBGJ8pUXFFow4dwQr96/k7p/utiTH+MJGYCww0pNOxp4mONmbotJyLGWMKTYyEhy8lOFc+ylc+X7GZlziEcqHVyHigif5KbUN+zo6G35DggK585LuDBr0YEbZsOQDBH9+HXEznoHkE96Jx/hczzo9efm8l/ln3z/c8eMd7Dy6098hmRJEVd9W1Taq+o+rg/EId+p52sk4RUSm48ymGgC2+poxxVywtzsZN+2T8VRVOZh4kPKh5bmiXS1oPy/nOgPnwdG9fLsxjUsXX0vUstdh2evsi/mISu0u8W58xid61+1NkCOIx355jGtmXMMbPd/gzEpn+jssUwKISFmc/W8quXadB+Q5yMnTTsbPA28A23E2F73tWZjGmKImyOGjTsZAQnICKWkplA/JYw6c6m2gcW9ienTPsrtS7E3s/u0TOHYQDm32WZzGO3rU7sHnfT8nIjiC2364jW/Wf+PvkEzJ8DYQCpwB7ALi3KnkUQuOK4sKB/YC1YGBwFWeHMNfRCQMGAZsBfao6jT/RmRM0eDmElT5cijxEABRoVHuxRIcytHqnSi7c2HGvqqz74bZro0bv4KGPb0cpfGm+lH1+fjij3ls/mM8/dvT/L33b4acPYTQwILNlG2KjsIeJg6sUNXXRCRYVSe6u4qCp3+6zQR6APWAuoBvZgZzk4hUFZFJIq5ljE/u7+VanGuYiDzr2t0PWKqqbwA3FHqwxhRRia4RTr6YD+fQCWeCk2cLTiZlB37HiSs+IiEgkt/Tss6MnDK1P8Tv8WqMxvuiQ6MZf/547mh5B1+u+5Kbv7uZbfHb/B2W8RI/DBNvIiIRQCUR6QJ0z6sCeJ7grFLV+1R1uKoOB273NEov6wJ8Cyd7R7paat4BHlLVYUArEemJs9f1PlexMoUcpzFFVlqac26a6uW8/xd2XGIc4PkyDaEtLyH86W1E3v0Ty8WZ5MxKPZvA5HiS5zwPY9vChO5wYIO3QzZeEuAI4P629/NGjzfYfnQ718Rew9ytc/0dlimepgOtgU9xdo/51p1KniY4G0Tk/ExrUfX3sL5XqeoXQHy23R1xTgqU6Nr+DegDbONkB6XjuR0zJSWF9u3bn/JwZarGlDjdmlTm8Qub8lTf5l4/9uFE56K/kcGR+arftGokxy//kI9SenGw12sABP39IRzcADv/hHFtYe6LYIs+Fllda3Xl876fUzOiJvfPvZ/X/niNlLQUf4dlihFV/VZVfwW2qmpb19INefJ0FNVgTl2LapiHx/C1ymRNeo649n0FDBORKkCuK8QFBgaybNky30ZoTBES4BAGdWvgk2MfTT4KQNngsvk+xjmtmnJOqy9RVT745Qr6p36ZtcD8lyEtFXo8Zcs+FFE1I2ry0cUfMeL3Ebz373tsOryJkeeNtH45xi0icgHwARAuIglAf1X9Ma96nrbgPKGq3dMfwJ35iNXX9gIRmbYjgb2qekxVH1PVN6yDsTGF42iSM8GJCI7Io2TeRITNze4AYETytYxOvpKD6kqcfh0Fw6NgWDmYMRiO2DwsRU1IQAjPdnyWJ85+gnnb5nHXT3cRn5S9Ad6YHA0CWqtqJNAWuM+dSh6vRZVtuyiuRbUIqCMiIa7tzjg7RxtjCtnR5KMEOgIJdnhn5fBmdWvS+MQHvJN6CWNT+9EvafiphZa9C6ObkfblHTZZYBF0fbPrMyYFvPWHWzM6ohtzGstUdQ+Aqu7CNQefiISfrpLvJsAoBCLSFefkP9VE5CkRKeNap2IQMFZEngeWq+ocvwZqTCl1NOkoEUERiJduHZ1ZM4okgqhYNphHezfhoq5deDfloixl/kxryPq06jhWfA4vVCF1yQRI2O+V8xvvuKjeRbzR4w02Hd7EHT/eQdyJOH+HZIq2GiJyq4h0E5FbgbIich4w5nSVPO2DU6So6nxgfg77Z3Ny5gxjjJ/EJ8cXqP9Ndo0ql+WBno04v3kVzqhRjqSUNN5dcPJr7OHoNxl13w3M+mUxDedeCEDAd4+S9OMw6DaE4C73WT+dIqJzjc6M7T6W+36+j4GzBzLxgomUCynn77CMG/wwD05HnKOfz8207xag5ekqidrogyzCw8M1IcEWizPGG+6Zcw/7ju3j85jPfXaOCcNvZ6BO4+vUzpS55j0uPKMqANvGXkStgwt5OnkA5zv+4LyAFQCk9h5BQPsBEGSzRRQFC3Ys4P6f76dR+UZMOH+CJTnFjIgcU9XT3irywjnOdY2iyr6/s6r+lls9j29RiUgVEakrIt65qW6MKbGOJh31agtOToJCnImK4qB3iyoZ+2vd/S1pQ3Zw60MvEnrNuxn7A34YwrHXz4FNp3xfGj/oUqMLY7qPYe2htdw5+067XWUyiEht15Q0W9Kfux4jAE6X3ICbCY6IOETkfyKyE/gHWADsEZGvXSc3xphTHE0+SniQT/+4IzjUOdQ4ONCRta9PYDCO0LLUqxjO2Wc0IfnpQ8y4fCVPhD5N2NEt8EFfksb3hDibYdffzqt5HmO6jWHdoXX0/74/uxN2+zskUzTMAybjnOblV+Aj18++7lR2twXnJeBPoL6qVlXVmqpaHhgO/E9EojyL2RhTGqR3MvYlh2sulaCA03cFCApw0PfMmjz36OCMfcG7lsGYM2DqjbD7X5/GaU6va62uvHP+O+w5toebvruJDXE2S3VxlctySZlfv0ZEPhWRx0RkmojE5HKoe1W1B84ZjOuqalecS0W5NdVLngmOiDiAN1X1a+CCzK+p6t84F9wMc+dkxpjSxdudjHMSEOScESIo0L2+jkEBDnbe8AvTzxh7cueqWHRSL7tt5WdnVT2L93u/T3JqMjfNuoklu5b4OyTjodMsl5RZGWCIqo4EXgRG53QsVZ3lelpNXR2GVTUNcGvtlzwTHFVNU9Wtrs2nRGSkiJyR6fVEVbVZtYwxWagqCckJlA3ybYITHORMbIIC3B8dVb3RmVxyZX+WXfMnE6s8TZfE1zksERB7P6ybDcveg5TEvA9kvK5ZhWZ80ucTKodV5q7Zd/HtereWHTJFR27LJWVQ1cmZ8oqGwH95HLOqiLwpIg+KyFtARXcC8XSYeH9gCzBARAYB36tqbB51jDGl0PGU46Rpmu9bcFz9bsTh+bRe7Zs1oH2zR9gxfSWvLInhheR34ZMrnS/OeAga9oJ+EyEs2pshmzzUKFuDDy/+kIfnPcxTvz3Ftvht3NP6Hq/Np2Q8l2n18Oyy5xG5LZeUhYiUwbnUUzfghjxOfxvOxb1bACuBSe7E7GmCE4BzafREoBPOGYMvAH5VVd+NAzXGFDsZ61D5uAUn0JXXOCT/85becV59Bi2qm7G9WypTVffC+p9IHtueoAf/htD8LRhq8icyOJK3er3F/xb9j/HLx7MrYRfPdXqOAEdhTLtishs4cCADBw48Zb+IZF85NcflkrLXU9XjwOMi0hCYKyL1VTU5p3OrahKQscCmiHQH8lya3tNvhI9xNiW1Aq5W1b6qeh9wxumrGWNKm/R1qHyd4KS6jn8wpFa+j1EjqgxB1ZplbHc8Ppq39QoAgk4cgBG1YFg50j68DI6e8l1tfCTIEcTwTsO5u/XdTN8wnf8t/h82d1uRl+NySSISLSKRACLyiJxsjtuO85ZTrhNTichFIjJdRH4Wkbn4qAVnDXC7qmY0P7nmw3Grw48xpvSIT3Z+Tfj6FtWmyj255989hFXq597Y0Vw0rlmVN/ZeytK0pjzV9wxu7tiHOeMj6Ll3ckYZx8a5JI1pR8D5zxBw1q1grQk+JyIMOnMQyanJTFwxkUphlbin9T3+DsvkQlWPubqwjBWRfbiWSxKRkcBBYAQQArwpIluBZsADqnrkNId9EngQ2AcIzu4yeXIrwRGRvqo6A7jO1YM585tJws2VPY0xpUdCknNGcF+34AQEBjEzrQNXFeAWFUD9iuG8kHINAHdWiyAowEHd5u3BleB81fV7fv/7Hy45+AGdvnuEuEWTKXfxs0j9rhAYcpojG2+4r8197D22l3f+eYeWFVtyXs3z/B2SyUVOyyWp6mOZnr/g4SEXquqy9A0R+cidSu624LwsIje7Dpx5vwK7cQ4jX+vmsYwxpUBhteAEB3hnzeDaFU7OdtG0qrO/TYOuN5JUpz7B0XXoV64ml3frwJz/LuPlmRN4PO5V+PQqZ4WrJkOzS6xFx4dEhKc6PMWaQ2t44tcnmNp3KjUjavo7LFM4qonIJ0B6nnEu0CuvSu5+M6wFZubwmAVsBF72NFpjTMmWkOxswfH1RH9B3kpwok8mONHhrpVoRAiu1xnK1XRtCr1aVOWRR57mt/avn6w8bQAnXmsL/3wGaaleicecKjQwlNHdRqMog+cNJjHVhvKXEnWBH4HNrkecO5XcbcEZq6q59lgWkVOGgBVXqampGT3FY2JiiInJbYJFY8zpxCc5W3DCg327VEOgB/PfnE6taPfnKw1wCJ37DuBE75v4+te/qPDr07Q4soYaX99J8g9PE9TzSWg3wCtxmaxqRdTixS4vct/P9zHuz3E8ctYj/g6p1PHDauK3qOr69A0R+c6dSnkmOK6ZjJef5vULgA7unKw4CAgIyG2svzHGA+nDxMMDfZzgOLyT4JQNCSQiJJBrznJ/NFZoUADX9WjP8S4zmPzbRtbO/4yHEyZTM/YBEn97h5BbvoWIKnkfyHikW61uXN34aj7870O61epG+6rt/R1SqZL5j/+JEyf6vMlSVdeLSDOgkmvXTcAdedVzayZjnDMYd8n+mohUB64k70l6jDGlzNEk50KbxWnekhXDe/NU3+Ye1ysTHMCg7o0Y9vgQvj17Cgu1JSEHV/HJp++zI+64DyI1D7d/mBpla/DUb09l3A41JZNrBNYI4FXgTqCdO/XcvXk9FOgvIjtFZIWI/CUiW4AvcN6+sqVfjTFZFMZK4pkVhUluy5UJ4p4+Z9Nk8PcAJOz4l8tGzeDFWas4cNT6i3hTWFAYL577IjuP7uT1P1/Pu4Ipzo6r6qXAF6p6A/CJO5Xc6oPjmnHwDhEJBxrgHMO+zRIbY0xuCmMlcYCiOO9bhXLOkWMDHbEMdMSyY0kFhi/uT4Nzr+OOrvUJC/Z0CjKTkzaV23Bt02uZumYqlzW8jOYVPG99M8WCq9c/5UUkEC+34ACgqgmqulxVl1pyY4w5naPJR33ewRicc1UUZdqgJzXkAGMdo3ngt7N4cfgjvP/TXySm2Ggrb7i3zb2UDynPC4tfIC3rNG2m5EgWkRhgGc51rpLcqeSd8ZXGGJNNYbXgpA/pblzF9+fySP3uEN0AufZT6Psa1OsKwPNB73PLgm6EPB/N+t1x/o2xBIgMjuTh9g+zfP9yvlr3lb/DMT6gqs+oaqyqfgFUVdUB7tQTW9cjq/DwcE1IsA5rxhRUzNcxNIluwqiuo3x6nrQ05Zd1++jauFLRW21aNWvnoE2/wAcnp574KKUXq9oNY1hMC4ID7e/N/FJVbvnhFtbHrSf2sljKh9rqQb6UeZj4xIkTE1U11Bfncd2OuhDnEg1/AKOAasCzqro6z/qeJDgi0jDbWPRzVfVXj6MuwizBMcY7un/ena41uzKs0zB/h1J0qMKelc6Vyce0BKDuiU/o26o6465rU/QStGJk/aH1XBV7FTENYniu83P+DqfUEJFjquqTe9Ei8jnOhTgrAH8Dq4BdQB9VvTqv+p7+yTAmfYVQEakHjPWwvjGmlDiadNTn61AVOyJQ9QyIqg3NLwNgc+gNvLG2Owv/WeXf2Iq5huUbcnOLm/l6/df8tuM3f4djvOOAqvYA2gAhqjpCVT8A/nWnsqcJzhxglGul0BnACg/rG2NKgeS0ZE6knvD5OlTFWt2sU4t1/qYj7F5B0opv+eG3pUxbuoWExBQ/BVc83d36bhqUa8Azvz3D4cTD/g7HFNwuyJiP789M+5PdqezpWMVlOJc2vx94CXBrumRjTOlSWCuJF2shkafue6cLwUBv12bLLyfRqmEtbulUj+5NKxPgpVmbS6qQgBBeOPcFbpx5IyN+H8FL577k75BMwfQWkfQvkbNFpKLreQecOchpedqC8wOwBWgOHMXZimOMMVmkL9NgLTinUaGB8+eV7/Nf3+k5FlkRejt7N/zD7R8uI2bcApZsPFCIARZPLSq0YGCrgczYOIOftvzk73BMwSQBCa7H3EzPfdKC85yqjnA9/0pEqnpY3xhTCqQnOIUxTLzYqtkeHt0A4RWpeSKZWd+ezcUBv59SbHbIY/zS5UPe/n0z102Io3HVctzVtQF9WlXz2krqJc3trW5n7ra5/G/x/2hTuQ0VylTwd0gmfx5T1aXZd4qI9yf6AxaKyHnpD5wdf4wxJovCWkm82At3trhHhgZxIqgcAB+EXAfDsvYfOW/BzUxJupeNoTeSeOI4w6f+QsthP/D1X9tJTbOpPrILcgTxYpcXOZJ0hHF/jfN3OCafckpuXPv/cKe+py04rwN/AQLUBmwqTmPMKdIXP7QWHPf9Gn4hNQ5vY0WFi5w7Wl0Dy6eeUm7uiashFD5MOZ+Hpt7Ca7PXcX/PRlzepob10cmkYfmGXNvkWj5d/SnXNb2OJtFN/B1SiZF5HhygyK6m62kLzkBVvVVVb1HVnoBNG2mMOUVGC04hLrZZ3O2Pask1Sc8g5es6d/R46rTlbw6czXdnL+cyx3yGT1tE7zG/8OUf20lJteUK0t115l2EB4Xz5t9v+juUEiUmJoYJEyYwYcIEKMINHZ6uRZXRXOTq2XyO1yMyxhR76S041snYfeXDnEtOVIoIce6Iqg3nPeZ83ndMjnWaLR/B4KOvsTRqKEGiPDztH85/7Re7deVSLqQc1zW9jnnb5rHlyBZ/h2MKmUcJjogcEpGNIrIJ+A1Y6JuwjDHFWUYn42C7ReWuQNftpYplQ07u7D4Uno2D9rfAY5vgwX+hy+BT6oae2Mus1guZdF0zLkqbz0NT/6bLyz+zevcR0kp5onNd0+sIdATyyapP/B2KKWSe3qK6S1Xrq2o9VT1TVSf6JCpjTLEWnxRPoCOQYEewv0MpNlJdy+ZEhGbqGilyci2rsGiIqgU9n4EH/oGGvbLUl/kv0+vrNjx2bDQz6n1Fx/gf6TNmHl1HzWXq0q0kl9JbVxXLVKRX7V58t+k7ktPcGl1sSghPb1Gd2uPNGGOySUhOICIowtZW8sCZNaMAN1ZFF4HydeHGL2HI1hyLnLHrS0YHv8OkoFEcOJrE41+uoNfo+Xy6pHQmOn0b9CUuMY6FO+ymQ2lSoEkURGSgtwIxxpQc8Unx1v/GQzd3rMPPD3flzFpR7lcKLQdlq+T6cveAf1jZciqfXF2TasEnGPH1YnqP+YW5q/cWPOBipGP1jpQPKc+sTbP8HYopRG4NExeRA8BhnMPDAdT1PBKY4JvQ/CM1NZWBA515W0xMDDExMX6OyJjiJyE5wZZp8FBggIP6lfJxzRpdAH995Hz+wHJ4vVWWl2XlV3Re+RWdAUKhR9o33DJ5KfUqhnPnefW5qn2tEj+8PMgRxLk1z+XX7b+Spmk4xCZILIjiMkxcVPPugCYia4GHVHVmtv3Xq+qnvgrOH8LDwzUhIcHfYRhTrPX/rj8BjgDe6/2ev0Mp+VKSYOXXULEhVGsDz5U/bfG0Gu35qsq9fPtfPEviIqlRMYqbOtTh5o51CCzBMyPHbohl6IKhTO07leYVmvs7nBJDRI6papGcD8LdT/OXwE8i8pqIjBaRhgAlLbkxxnjH0eSjNgdOYQkMhjOvgRrtwOGAp/bCkG0nX78m6+ghx45lXPnnAD46cR9rQ/vTIGk1z834j+6vzuODhZs5kVxkpzUpkI7VOwKwcKf1wykt3E1wjqtqIvAo4FDV9T6MyRhTzKV3MjZ+EBgCoZFQpaVzO7I6XDwq1+KTkh7n1Z5lucMxg2en/0vTp79n0q8bS1yiU7FMRRqVb8SSXUv8HYopJJ6OokoBMrrgi8gVXo/IGFPsWSfjIqDeuc6f4ZUgNOq0Ra/47RJuPvou713kbHV7fuYqurz8M58v21ai5tFpW7ktK/avIDWtZCVvJmfurkXV2zVzMcC5IjLS9bwDzttXxhgDgKpaJ+Oi4PznoG1/59w5Ia5/i4hqEBAMcTnP6ttj7mVsDnU+/zLxXN768hKmLj2DJy5qSvu60YUUuO+cWelMpq6ZysbDG2lUvpG/wymxRKQX0A/YC6iqDs/2+uNAVWA30A54RlVXezsOd1twkoAE12NGpuc2a5IxJovjKcdJ1VRrwfG3gCCo3NT5vEx55wrlD6+Gy8efLNP80lyrXxHwK3NCHqXytu+58p1FjJ69lsSU4t3y0aqSc4TZP/v+8XMkJZeIhAHv4ByYNAxoJSI9sxUrCwxW1ZdxNpK84otY3G3BeSynZctFpJ2X4zHGFHPpyzRYC04RVaejM9nJbPe/8E7nHIu/Hfw6L0cO5cM58SxZsYYJgy6kXJmgQgjU+2pH1CYqJIrl+5ZzZeMr/R1OsZJpcc3ssucRHYEtrn674FzWqQ8wJ72Aqj6dqbwDOOrFUHMNLEc5JTeu/X94NxxjTHFnCU4xVPUMZ9IzrFyOLz9+5EUeDwWOwKOT3+WVQcUzORARWlZsyYr9K/wdSrEzcODAjDniMhORlGy7KgPxmbaPuPblVDcY6A/c46Uwsyi5kx4YY/ziaJIrwbFbVMVPn1fzLPLKntvYur/4zhXWJLoJmw9vJik1yd+hlFR7gcxDKCNd+7JwJTdvA0+q6gZfBGIJjjHGq6wFpxg763a4/y8YMBOiG+RabOV//8DyaTD+PNjr9b6hPtWkfBNSNIUNcT75nWpgEVBHREJc252BmSISLSKRACJSBhgPjFbVP3w1ItvdPjjGGOMWa8Ep5qLrOx93/Qq/jIIFo6F6Wzi6h7RuT+CYfi895vWDtOPO8m+dg9Y8h3+SqtJ677cMqTiOSo3P4ap2tahdIcy/7yUHjaMbA7Dm0BqaVWjm52hKHlU9JiKDgLEisg9YrqpzXKOvDwIjgE+AM4B6rgV5w/HBiGxLcIwxXpXegmMT/RVzweHQ8xk49+GMYeYOVZh+LyHpyY2LbF9Ca9fzEfvv46adT/DOz83o0LAKt5/XgPMaVSwyK8vXiahDaEAoaw6u8XcoJZaqzgZmZ9v3WKbn/QojDktwjDFeld6CEx5sSzUUeyIn59BxbU9vOpJLVj+Wex3go+CXnE+2A5/CYxEj6NrzIi5sVcfvC3sGOAJoGNWQtYfW+jUO43vWB8cY41XpLTjhgZbglETJDXrnuH9/aO1c64yMH0Kvr9vRf8QHvP7TOhISsw+8KVxNopuw+uBq3Fls2hRfluAYY7wqPime8KBwAhwB/g7F+EBEWGjG825Mynj+U5ep0PtFOPtOaHbJKfVCJIWPkx5g8c9f0+ml2bwzfwOHj/lnrtjG5RtzJOkIe47t8cv5i7vY2NjMw8aL7P/odovKGONVCckJtpJ4CRaZaZK/2jVqwA7n8xqVK0LjTNOZHN4Or7U4pf6U4BcAeGv2JfT46TL6tm/C7efWp1Z04XVIbhLdBIC1h9ZSNbxqoZ23pIiJiSEmJgaAiRMnFtnpra0FxxjjVUeTj1oH4xIsIvTk38X1q5Rjn0ayXyOpFBmatWC5mtD8slyPc3fgdJYEDmTh70voNmoeL33wFevX/OujqLNqGNUQgHWH1hXK+Yx/WAuOMcar4pPiiQi2BKekigw92YJTs3wZzkscQ1lOMD2n5Ruu/gBUYe9/sGQ8lK8Lc06uuxioKcwOGsyRoIpEbtoPm2BhWA8CLn2DsxtX99nIq3Ih5agcVpn1cet9cnxTNFiCY4zxqvikeMqHlvd3GMZHMic4lSNDOY7zkev6VCJQpQVcMta53flB588ZD8KfHziPmbw/o3inYz9z9eQpJNbsyKCuDbigeRUcPhh51SiqkSU4JZzdojLGeJXdoirZygSf7FNaJSLk5P4gN/uaOhzOxyVjnetf3f7zKUU+D/kf3+67mA7T2jJw1AfMX7LU6yOeGpVvxMa4jaSk+XdEl/EdS3CMMV5lt6hKtqCAk60pUWHBGc/zfTupZju4+dscX4qSBCYde4Cu3/Xiw+dvY9K8NRz10hDzhlENSUpLYmv8Vq8czxQ9dovKGOM1qkp8Urwt01CCiQh3Jj1IeTnKoCAv/Y1cvxvc9yeMa5trkf6pX8K8L9kwrwarqnWjQZ8HiK7ZJN+nbFje2dF4/aH11C9XP9/HKY1iY2OJjY1N3yyyw8StBccY4zWJqYkkpyVbC04J90Pa2XyW2sP921LuqNAAHl4L9y6DG7/KtVgDdnDWrk+InnQ2c94ZzLaDx/J1uvrl6iOI9cPJh5iYGCZMmMCECRMAiuwwcWvBySY1NTV98qIsY/2NMXmzdahKl5CgAGbc14WyIV76VRJRxfmo2Ahu+AK2/Q6/jMy1eM/d7/LkaKFKvRb0uuhKmlePdPtUZQLLUDuytiU4JZglONkEBASkZ6XGGA8dSToCYC04pURokIMzapTzzcEbnQ8Ne0G3IfBcdK7FXgicBNuACQ/xebnbaHLelZzZrpNbp2gY1dDmwinB7BaVMcZr0hfatD44pUNwgI9/hYiAIwCG7oRBC6HHU9Dk4lyLX334Xc6MvYhHn3+JOX+uQvefPnlpGNWQrfFbOZFywtuRmyLAWnCMMV4TnxQPQGSw+7cKTPHlq4n4ThEc7pxLp0oL58SB6+fAJ1eABICe2gXklZQRMH0EALNbvUbnvjcTFnzqr7tG5RuRpmlsOryJZhWa+fxtmMJlLTjGGK+JT3YmOGWDrAXH+IgINOoFD/wDT+6CRzfA1R9CSM63ys5f/hAHXmzO/PeGsi8+MctrjaIaAVg/nBLKWnCMMV6T3oJjfXCMz5Wv6/wZGALNL4XaHSFhP7zd8ZSitdhDra1vsnfUp8yq1Z/mfe+nbqVy1IqsRZAjyPrheKi4DBMXb88OWdyFh4drQkKCv8Mwplh6/9/3Gf3HaJZcv4SwoMJbHdoUrm/+2sGBhCRu61LP36Gcas9/cCIO3r8o1yKpKvwbeR6RTc7j0dSFVAqvwtu93i68GEsQETmmquH+jiMn1oJjjPGa+KR4AiSAMoFl/B2K8aHL2tTwdwi5q9Lc+XPYYZj5MKyeBfE7sxQJEOXM+PmwbD7lqjRlddh2Uo4eJDCkDATZZ7eksD44xhivSZ/FuNA6nxpzOn1ehYdXwUP/Qa/hORbpeGIn+znO8Vfrs2XiDRxLsrWpSgpLcIwxXnM0+ah1MDZFT7ka0PkBZ4fkB/6B0KiMlxonJQGwPjiIOnvncPjFxoyevZYDRxNzOZgpLizBMcZ4TXxSvA0RN0WTCIRXdHZOvvkbaHkVAA2TkwFYF+xcOLQaB0ie9yofv3wPEz76mE37rU9mcWV9cIwxXhOXGEdkiCU4poir3gaumASBoVSL20x42ibWBwVlvPx40GfOJxumsXDsm1QI3MnxZv2o0rA9tLnBT0EbT1mCY4zxmsOJh6kaXtXfYRjjnkvfQICGX/ZhfVQwlBfY+RccP5hRpJPjP0iDyJXvwcr32Pf7VCo07oSj+xD/xW3cYgmOMcZrDiceJiokyt9hGOORhtXOYs7WOeiNvyDJx+HoHuftrJmPwPLPspSttGs+7JoP818i8ZY5hNRp76eo/ae4zINjfXCMMV6RpmkcTjpsfXBMsdOofCPiEuM4cOIABIdBdD0IiYB+4+GZQ7nWO/z+VWwcdwnxB3cXYrT+FxMTw4QJE9IXpj51rYwiwhIcY4xXHE0+SpqmWQuOKXYaRjUEyHlGY4cDbp7uHHKeTWUOUv/AfD4Z9zRj56wjMaXI/q4vlSzBMcZ4xeEThwGIyjQE15jioFH5PNakqt8Vap6da/279HPu/7U9u59vAcPKOScYNH5nCY4xxiviEuMArAXHFDvRodFEh0affk2qKi2g2xPQb5Jzu+ZZzp9hFTKK1GGX88nSSTDnf86Vz5eMd3ZcLkVEpJeIvCUiw0Tk2VzKXC0iG0Skr6/isE7GxhivOJzkbMGxPjimOGoU1ej0q4o7AqCba+RU+bpQox0cWOdc7PObu6H5ZfDdoyfL/zoKThyGpRMhKBwufBEanu+cdLAEE5Ew4B2ghaomisiXItJTVedkKlMP2Ads82Us1oJjjPEKa8ExxVmj8s4EJ03T8i5c6yxn35xKTZzJzi2zoPX1UPdcdnYcdrLc0onOn8kJEPsA2yZdR/yXDzpXPS+5OgJbVDV9KujfgD6ZC6jqJlWd6+tArAXHGOMVhxNdfXAswTHFUMOohhxPOc7OozupGVHT8wOElIUBM6gevwcWDcuxSK34f2DFP7DifRZ3m0LbckcIbn2Nc5blIi7TqKnssucRlYH4TNtHXPsKnSU4xhivSE9wIoIj/ByJMZ5rWP7kSKp8JTjpymb9XX5IIygv8acU6zDvOueTb+90/nxoJZQrwHl9bODAgQwcOPCU/SKSfXXSvUDmL4FI175CZ7eojDFecfDEQcqHlCfAUWTn/TImV+lDxU/bD8cdIkzs/mfG5sEgV8LjcC0FEVUHLV+PlKCsfwhsGdeHQ2/0IDluBxzZWbAY/GsRUEdEQlzbnYGZIhItIoXaQc9acIwxXrH/+H4qlKmQd0FjiqDwoHBqlK3BurjTjKRyU/XyYfyU2oZeAX/xU837aFBjE5xxBcTvhppnIWHRBKalwXPlM+rUSdkM+zfDmOYAHLlvDZEVit+yJ6p6TEQGAWNFZB+wXFXniMhI4CAwQkQEeBKoA1wjIsmq+oO3Y7EExxjjFfuP76dimYr+DsOYfGsY1ZC1B9cW+Dg1ypehX/LDBCancn/tFtDjFucL1c48WcjhAEcgpGW/w+OkY9vye6Xe1LhyBDXKhUBwWQgoHr+yVXU2MDvbvscyPVfgedfDZ0rVLSoRCRSRJ0Ukx55Sxpj8swTHFHctK7Zk4+GNHEk6UqDjVAgPJg0HSQRRrkxQ7gWH7oJ7lub4UjlJ4Oz9X1Hjncbwch32fzoQUhIhzWZLdlepSnCAcOB7St/7NsanVJUDxw9YgmOKtTaV26Ao/+z9p0DHKRd2MqmJPF2CExgMlRrDJeOckwieRsUNX8LzlTn6ekf0yE5ITS5QjKVBobZ3iUgT4DrgONAVGKaqv+fjOFVxNm2dqapnZdrfC+iHs8e2qurwzPVU9bCIHCjAWzDG5CAhOYETqScswTHF2hkVzyBAAvhr71+cW/PcfB8nIuTkr9bTJjjp2t7s/HniCBzcCG1ugEObYe8q+PuTLEXLHl4Do5uR7AgltWYHQpv2gk735TvW/Cguq4kXWoIjIgHAaCBGVdNE5EMgJVuZs4BDqro+U50rVPXzbIfrAnwLtM5UN9fZE0XkDuBM4FGMMV63/7hz4jLrZGyKs7CgMJpGN+WvvQVbWkEyzWtz2ltU2V34YtbtlES44HnY8Sd8ckWWl4LSThC0dR5sncf+P78lqseDBDb32aoHWcTExBATEwPAxIkTi+w9s8K8VXMWIMB9IvIEEANkn87xEPCxiLR0DTGbCpzyjamqX5B1IiE4zeyJqjpRVe9V1eN5BZmampox3j9ThmqMOY2MBCfUEhxTvLWp3IZ/9/9Lcpp3bgF5lOBkFxgCYdHQqBcMOwwPLIeQctD4wizFKu5fSuDnN7D03Qc5NvvFkj5TstsK8xZVHZxJyHWuW0UfA0nA5PQCqrpeRK4CpuG8jfWuqn7s5vHznD3RNTTtGqCJiLRV1T/JJiAgILfZGo0xuUhPcOwWlSnuzqx8Jh+v+pjVB1bTslLLAh8vMrQACU525evAkC3O21hrvz/l5bO2vQ/bIO23kThQuPpDaBrjHLFVChXmuz4CrFbVw67tBUC3HMrFA8dwzoToyWxHec6eqE4vq+q5OSU3xpj82ZXgXEW5anjxm7fDmMzaV2kPwJLdS7xyvAK14OREBCo0gIfXwBPbITgCujyUpYgDdT75/GZ4rjxp81+Bo/ucq5uXIoWZ4CwBKrj61YCzRSfLhAMiUhmYDvwP6Ak8KSKXunn8HGdPLHDUxpg87Ty6k4jgCFumwRR7FctUpHH5xizaucgrxwsO9NGv2YiqEBIBj2+GXsOgcgvocA8EljmlqGPu8zCqISnfPwkb5sKBDb6JqYgptFtUqnpQRB4HxrhmN6wEPJet2BnAI+kjq0QkBjhl8QsR6QrcBFQTkaeAV3ObPdGHb8kY47IrYRfVw6v7OwxjvKJT9U58vOpjjiUfIywoLF/HeKhXY2at2OXlyHKQPvnf3QudPxud77yFpWkw65EsRQOXvAlL3nRuPLkHgkJ9H58fiZayJqu8hIeHa0JCgr/DMKZYufzby6kZUZNxPcb5OxRjCmzhzoXcOftO3ur5VoGGi/vVkZ3w7gXQ/haYk70twUkdgdD1caTJRVDV/f5GmYeJT5w4MVFVi2SmZAlONpbgGOMZVaXjlI5c1vAyhpw9xN/hGFNgJ1JO0HlKZ65ucjWPn/24v8MpuMM7ICkB/aAvcnRPjkUSml1NeMu+0NzdXiFOInJMVcO9Eaa3lc6u1cYYrzmSdISE5ASqhVfzdyjGeEVoYCjtq7bn1x2/UiIaAcrVgEqNkUfWOpeHAFSyzs8Xvupz+Pxmtn18D7psMpw4nMOBihdLcIwxBbL96HYAapSt4edIjPGenrV7suXIFtbHrfd3KN4VHAaDVyFP7XUuEdHjKZIqNMt4udb6j5EZD8CI2iR//xT8NMx/sRaQJTjGmALZGLcRgPrl6vs5EmO8p0ftHgjCT1t+8nco3hdZ3dk5ue3NcN6jBF/3ETS5mNReWfvqBC0eBwteY81P7xfLRT4twTHGFMiGuA0EOgKpFVnL36EY4zUVy1SkTeU2/LS1BCY42VVsBNdNIaDLA1C7E1qxSZaXmyx4EJ6LZtM3/yNlyxLY/a9/4vSQJTjGmALZELeBupF1CXJ4eUIzY/ysZ+2erD20ls2HN/s7lMJzyyzkrl8zNlPPuCrjeb2/RxH4/gXwTmeSxrSFYwf9EaHbLMExxhTI+rj1NIhq4O8wjPG63nV74xAH0zdM93cohUfEuQbWEzvg6f0E9BufY7HguA0wsl4hB+eZwlyLyhhTwhxLPsaOozu4pOEl/g7FGK+rEl6FjtU7Mn3DdO5pfQ8BjoC8K5UUIWVPPr99DsTvghrtQdP49v42zPw3Lv3VIntRrAXHGJNvK/avQFFaViz4ooTGFEWXNbyMPcf2sGSXd9amKpZqtodmMRBZDcrV4NIP9jLhjyTeGvMKQJHtfWwJjjEm35bvWw5gCY4psbrX6k5kcCRfrPvC36EUOYFd7vd3CKdlCY4xJt+W7l5Kg3INKBdSzt+hGOMTIQEhXNH4CuZsncO2+G3+DqdoEfF3BKdlCY4xJl+OJR9j2Z5ldKnRxd+hGONTNzS9AYc4+Pi/j/0divGAJTjGmHxZtGsRyWnJdK3V1d+hGONTVcKrcHG9i/l6/dccOH7A3+EYN1mCY4zJlxkbZhAVEkXryq39HYoxPnd7y9tJSk3i7X/e9ncoxk02TNwY47E9CXuYu20uNze/2Sb4M6VCvXL1uKrxVUxbO41rm1xLw/IN/R2S38TGxhIbG5u+acPEi4uUlBR/h1DsTJgwwd8hFDvF/ZqNXz4eQbiqyVV5F/aS4n7NCptdL8/ldc0GtR5ERHAEQxcMJTk1uZCiKnpiYmKYMGFC+vUqsj2NLcHJxhIcz9kXqeeK8zX7c8+ffLnuS65tei21Igpv/anifM38wa6X5/K6ZtGh0QzrNIxVB1fxyrJXUNVCiqxIK7J3gizBMca4bd2hdQyeN5iaZWtyd+u7/R2OMYWuZ+2e3Nz8ZqasnsK4v8ZZklOEWYLjRZnuSRZ6fX/VLSh/xm3XzH3Jqcl8vuZzrou9Doc4GNdjHBHBEYVybm+wz1nhnrukX7OH2z9Mv0b9mLhiIvfPvZ+dR3cW+NzF9XoXZZbgeFFx/Z/avkQL/9wFURhxqyo7ju7g560/8/zi5+n9ZW/+t/h/hMSF8PHFH1M/qr7Pzu0L9jkr3HOX9GvmEAfDOg7j8bMeZ+GOhcR8HcNzi55jemz+F+Usrtc7JyLSS0TeEpFhIvJsDq+HisgbIvKEiLwnIo29GoBLkb13ZozJvwU7FrD+0Ho0/T91/txecztv/PUGyWnJpKSlkJyW7HykJnMs5RgHjh/g4ImD7D22l2MpxwAoE1iGDtU6cHWTq/lg+AdUL1vdz+/OGP8TEW5sfiO96vTi9T9fZ9raaTQv19zfYfmdiIQB7wAtVDVRRL4UkZ6qOidTsQeBrao6UkRaAu8C53o9Frt/mJWIpAHH81k9gIItPFaQ+v6qC85EOb+9s/0Zt12z4nVuu2aeKcj1Kui57ZoVn7ru1g8k50YRh6pmjKQSkZ7AUFXt6doeDNRU1cGZyvzqKvOra/uIq8yRAryHHAM2maiq3bYzxhhj8qcyEJ9p+4hrnztlvJrg2C9zY4wxxnjLXiDzCIRI1z5PyxSYJTjGGGOM8ZZFQB0RCXFtdwZmiki0iES69s0EOgK4+uD84+3bU2B9cIwxxhjjRSJyPnAlsA9IVtXhIjISOKiqI0SkDDAK2AU0BF5U1bVej6M0Jjgi0gB4HvgTqAkcUNXnspUZANwFnHDteldVPyrMOIsSEXEAscASIBhoANyqqsczlQnF+aHdATQCRvjiQ1scuHm9BmCfsVO4vvyWAD+q6iPZXrPPWA7yuGYDsM9ZFiKymJPXIzW9Q2ym1+1zlo0b12wARexzVlo7GUcDn6nqtwAi8p+IzFTVP7KVu1ZVNxd6dEXXIlV9HkBEvgX6AZ9kev1BCmHoXzGS1/UC+4zl5Hngr1xeexD7jOXkdNcM7HOW3feqOuw0rz+Ifc6yy+uaQRH7nJXKBEdVl2bb5QAScih6r4jsBsKAN1T1oM+DK6JUNQ3nlygiEoiz5WtNtmJ9gKGu8itE5EwRifTFvdWizs3rBfYZy0JEbgJ+A1oBZXMoYp+xbNy4ZmCfs+xaisjjQBlgqarOzPa6fc5Oldc1gyL2OSuVCU5mInI58IOqrs720nxgpqruE5GLgWlAz1MOUMqISG/gIWCGqi7L9nKhDP0rTvK4XvYZy0REmgPNVHWoiLTKpZh9xjJx85rZ5+xUL6vq7yISAPwiIvGq+kum1+1zdqq8rlmR+5yV6lFUItId6I7zF1AWqrpJVfe5Nn8Gurr+YUs1Vf1BVS8E6olI9tUWC2XoX3Fyuutln7FTXA6cEJEhQBfgbBF5MFsZ+4xllec1s8/ZqVT1d9fPVOBXnL8HMrPPWTZ5XbOi+DkrtQmOiPQBegMPAFVFpGPmYWwi8pLr1gI4O5ltcv3Dlkoi0tx1zdJtAur7Y+hfceDO9bLPWFaq+oKqPqeqI4AFwO+qOsY+Y7lz55rZ5ywrEWkqIrdl2tUIWG+fs9y5c82K4uesVN6iEpF2wFRgGTAXCAfexPnX0EFgBLAbeFtENgEtgZv8E22RkQjcJiJtgCCgGXA/MIST1+x1YJSIPIVz6N9tuRyrNHDnetlnLAcicgVwHhAsItcBbbDP2Gnlcc3sc5bVEaCviFTH2TKzDZiC81rZ5yxn7lyzIvc5K5XDxI0xxhhTspXaW1TGGGOMKbkswTHGGGNMiWMJjjHGGGNKHEtwjDHGGFPiWIJjjDHGmBLHEhxjjDHGlDiW4BhjjDGmxLEExxhT6hV0Snl/T0lvjDmVJTjGlGIiMlpEDonIja7t5iKyS0TEtf24iPwoInW9cK4bRORQQY+Ty7EfzPS8j4hscidmEQkTkTFAdAFDqCkiIzNNVW+M8TNLcIwp3R4FUoCfXNsXAceAs1zbq4HnVHVzQU+kqp8Ahwt6nFw8mOk8M4EtbtYbC3yVaZHAfFHVLcBiYFhBjmOM8R5bqsGYUk5EPgF+VtV3RWQ0ziQkQFWfEZGXgKeAMjjXb/sFaAJ8Csxz/WyCc92ZFOAD4AWgLc617lKBeFUd6TrXZlWtm+ncz2UuB+wDXgJeA+rjXMOrb/pChyIyFufaXjtwJmPrcK5cPBYYA6xW1c9EZB4wK9MxYrIvligiZV31q6vri9C11k4VVf1LRDoCq3Cuq3Ol6/0mARcCP+Jq9VHVD1x1ywAbgRqqmubRP4IxxuusBccYMwPnQnrlgEM4V1Lu63pNXCsCpwGvqerLwCPAS6qaAtwBlAPWApuB+cBxoIOqDlXVp4GLRaR19pOKSO/s5YC/cLYa/a2qdwIrgPNd5fsAjVR1kKo+DygwWVU/BOJUdZiqfpbpFH+p6l2Zj5FNA+CgZv0r7wYgzPX8MZyLpqa/vtPVOtQKWOi6bm3SK6rqcSAAqJLDuYwxhczuFxtjvgfeAmJwtkz8AVQVkQtxJgcAAnRztWokA5UAVPWwiHwHXAuEAJ8AvYAwERniqrstvXw2rU5Tbq3r5z4gwvW8Bc4Wl3Qb83hf610/92c6RmYhOFudMmsLvCYiwYDDlbQsEJEnVHWJq5XmgKoeEJGLgWXZ6ifjbO0yxviZJTjGlHKqekhEVgD3Ap1UVV1Jyyigq6vY7Thv5dwqIkHAXZkO8QYwGVimquNFpBLQUVVHAIhID04mG5n9c5pyOd07/w/okWm7fqbnqa6O0a1V9a/THCOzbWTqXOzqIFxVVVNE5ApgiYg0BHbi7JcEzr5JS13P+wAjRKSN65aWAGVd5Y0xfma3qIwx4LzdsiZT35GZwGFVPeDa/gFoJCKv4LxFVc6VBKCqK3H2n5nr2v4R+FNEXhKR/+HsK7NDRG5w1bvrNOWaAXWAW0WkPnAeEONKmmYCG0RkoqvVJ/Pto5k4E7L+InL+aY6RQVV3ueKq6drVHtjjuhV2BNiN8zuyBfCrq0zL9PeJsyPzBcBy1/aZwDxVPeHuRTfG+I51MjbGFBsi0klVF7qezwFuUdWtBTheK+A+nC1SjwE/qerS09fK8TihwETgaW+MODPGFJwlOMaYYkNEZgNzcPap2auqr3vhmNVw9sV5CRikqsn5OEYN4LiqHixoPMYY77AExxhjjDEljvXBMcYYY0yJYwmOMcYYY0ocS3CMMcYYU+JYgmOMMcaYEscSHGOMMcaUOJbgGGOMMabEsQTHGGOMMSXO/wEeEC5g7c79qAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot the two spectra\n", "fig, ax = plt.subplots(1,1, figsize=(8,5))\n", "\n", "xr = [2.5,5.5]\n", "\n", "for sp in [sp_sci, sp_ref]:\n", " w = sp.wave / 1e4\n", " ind = (w>=xr[0]) & (w<=xr[1])\n", " sp.convert('Jy')\n", " f = sp.flux / np.interp(4.0, w, sp.flux)\n", " ax.semilogy(w[ind], f[ind], lw=1.5, label=sp.name)\n", " ax.set_ylabel('Flux (Jy) normalized at 4 $\\mu m$')\n", " sp.convert('flam')\n", "\n", "ax.set_xlim(xr)\n", "ax.set_xlabel(r'Wavelength ($\\mu m$)')\n", "ax.set_title('Spectral Sources')\n", "\n", "# Overplot Filter Bandpass\n", "bp = pynrc.read_filter('F444W', 'CIRCLYOT', 'MASK430R')\n", "ax2 = ax.twinx()\n", "ax2.plot(bp.wave/1e4, bp.throughput, color='C2', label=bp.name+' Bandpass')\n", "ax2.set_ylim([0,0.8])\n", "ax2.set_xlim(xr)\n", "ax2.set_ylabel('Bandpass Throughput')\n", "\n", "ax.legend(loc='upper left')\n", "ax2.legend(loc='upper right')\n", "\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialize Observation\n", "\n", "Now we will initialize the high-contrast imaging class `pynrc.obs_hci` using the spectral objects and various other settings. The `obs_hci` object is a subclass of the more generalized `NIRCam` class. It implements new settings and functions specific to high-contrast imaging observations for corongraphy and direct imaging." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this tutorial, we want to observe these targets using the `MASK430R` coronagraph in the `F444W` filter. All circular coronagraphic masks such as the `430R` (R=round) should be paired with the `CIRCLYOT` pupil element, whereas wedge/bar masks are paired with `WEDGELYOT` pupil. Observations in the LW channel are most commonly observed in `WINDOW` mode with a 320x320 detector subarray size. Full detector sizes are also available.\n", "\n", "The PSF simulation size (`fov_pix` keyword) should also be of similar size as the subarray window (recommend avoiding anything above `fov_pix=1024` due to computation time and memory usage). Use odd numbers to center the PSF in the middle of the pixel. If `fov_pix` is specified as even, then PSFs get centered at the corners. This distinction really only matter for unocculted observations, (ie., where the PSF flux is concentrated in a tight central core).\n", "\n", "The `obs_hci` class also allows one to specify WFE drift values in terms of nm RMS. The `wfe_ref_drift` parameter defaults to 2nm between \n", "\n", "We also need to specify a WFE drift value (`wfe_ref_drift` parameter), which defines the anticipated drift in nm between the science and reference sources. For the moment, let's intialize with a value of 0nm. This prevents an initially long process by which `pynrc` calculates changes made to the PSF over a wide range of drift values. This process only happens once, then stores the resulting coefficient residuals to disk for future quick retrieval.\n", "\n", "Extended disk models can also be specified upon initialization using the `disk_params` keyword (which should be a dictionary).\n", "\n", "The `large_grid` keyword controls the quality of PSF variations near and under the corongraphic masks. If False, then a sparse grid is used (faster to generate during initial calculations; less disk space and memory). If True, then a higher density grid is calculated (~2.5 hrs for initial creation; ~3.5x larger sizes), which produces improved PSFs at the SGD positions. For purposes of this demo, we set it to False." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# The initial call make take some time, as it will need to generate coefficients \n", "# to calculate PSF variations across wavelength, WFE drift, and mask location\n", "filt, mask, pupil = ('F444W', 'MASK430R', 'CIRCLYOT')\n", "wind_mode, subsize = ('WINDOW', 320)\n", "fov_pix, oversample = (321, 2)\n", "\n", "obs = pynrc.obs_hci(sp_sci, dist_sci, sp_ref=sp_ref, use_ap_info=False,\n", " filter=filt, image_mask=mask, pupil_mask=pupil,\n", " wind_mode=wind_mode, xpix=subsize, ypix=subsize, \n", " fov_pix=fov_pix, oversample=oversample, large_grid=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some information for the reference observation is stored in the attribute `obs.Detector_ref`, which is a separate NIRCam `DetectorOps` class that we use to keep track of the detector a multiaccum configurations, which may differe between science and reference observations. Settings for the reference observation can be updated using the `obs.gen_ref()` function." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Set default WFE drift values between Roll1, Roll2, and Ref\n", "\n", "# WFE drift amount between rolls\n", "obs.wfe_roll_drift = 2\n", "\n", "# Drift amount between Roll 1 and Reference.\n", "obs.wfe_ref_drift = 5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exposure Settings\n", "\n", "Optimization of exposure settings are demonstrated in another tutorial, so we will not repeat that process here. We can assume the optimization process was performed elsewhere to choose the `DEEP8` pattern with 16 groups and 5 total integrations. These settings apply to each roll position of the science observation sequence as well as the for the reference observation." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "New Ramp Settings\n", " read_mode : DEEP8\n", " nf : 8\n", " nd2 : 12\n", " ngroup : 16\n", " nint : 5\n", "New Detector Settings\n", " wind_mode : WINDOW\n", " xpix : 320\n", " ypix : 320\n", " x0 : 914\n", " y0 : 1512\n", "New Ramp Times\n", " t_group : 21.381\n", " t_frame : 1.069\n", " t_int : 329.264\n", " t_int_tot1 : 330.353\n", " t_int_tot2 : 330.353\n", " t_exp : 1646.322\n", " t_acq : 1651.766\n" ] } ], "source": [ "# Update both the science and reference observations\n", "obs.update_detectors(read_mode='DEEP8', ngroup=16, nint=5, verbose=True)\n", "obs.gen_ref_det(read_mode='DEEP8', ngroup=16, nint=5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Add Planets\n", "\n", "There are four known giant planets orbiting HR 8799. Ideally, we would like to position them at their predicted locations on the anticipated observation date. For this case, we choose a plausible observation date of November 1, 2022. To convert between $(x,y)$ and $(r,\\theta)$, use the `nrc_utils.xy_to_rtheta` and `nrc_utils.rtheta_to_xy` functions.\n", "\n", "When adding the planets, it doesn't matter too much which exoplanet model spectrum we decide to use since the spectra are still fairly unconstrained at these wavelengths. We do know roughly the planets' luminosities, so we can simply choose some reasonable model and renormalize it to the appropriate filter brightness. \n", "\n", "Their are a few exoplanet models available to `pynrc` (SB12, BEX, COND), but let's choose those from Spiegel & Burrows (2012)." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Projected locations for date 11/01/2022\n", "# These are prelimary positions, but within constrained orbital parameters\n", "loc_list = [(-1.625, 0.564), (0.319, 0.886), (0.588, -0.384), (0.249, 0.294)]\n", "\n", "# Estimated magnitudes within F444W filter\n", "pmags = [16.0, 15.0, 14.6, 14.7]" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Add planet information to observation class.\n", "# These are stored in obs.planets.\n", "# Can be cleared using obs.delete_planets().\n", "obs.delete_planets()\n", "for i, loc in enumerate(loc_list):\n", " obs.add_planet(model='SB12', mass=10, entropy=13, age=age, xy=loc, runits='arcsec', \n", " renorm_args=(pmags[i], 'vegamag', obs.bandpass))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Companions: 0%| | 0/4 [00:00" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from matplotlib.patches import Circle\n", "from pynrc.nrc_utils import plotAxes\n", "from pynrc.obs_nircam import get_cen_offsets\n", "\n", "fig, ax = plt.subplots(figsize=(6,6))\n", "\n", "xasec = obs.det_info['xpix'] * obs.pixelscale\n", "yasec = obs.det_info['ypix'] * obs.pixelscale\n", "extent = [-xasec/2, xasec/2, -yasec/2, yasec/2]\n", "xylim = 4\n", "\n", "vmin = 0\n", "vmax = 0.5*im_planets.max()\n", "ax.imshow(im_planets, extent=extent, vmin=vmin, vmax=vmax)\n", "\n", "# Overlay the coronagraphic mask\n", "detid = obs.Detector.detid\n", "im_mask = obs.mask_images['DETSAMP']\n", "# Do some masked transparency overlays\n", "masked = np.ma.masked_where(im_mask>0.98*im_mask.max(), im_mask)\n", "ax.imshow(1-masked, extent=extent, alpha=0.3, cmap='Greys_r', vmin=-0.5)\n", "\n", "for loc in loc_list:\n", " xc, yc = get_cen_offsets(obs, idl_offset=loc, PA_offset=PA_offset)\n", " circle = Circle((xc,yc), radius=xylim/20., alpha=0.7, lw=1, edgecolor='red', facecolor='none')\n", " ax.add_artist(circle)\n", "\n", "xlim = ylim = np.array([-1,1])*xylim\n", "xlim = xlim + obs.bar_offset\n", "ax.set_xlim(xlim)\n", "ax.set_ylim(ylim)\n", "\n", "ax.set_xlabel('Arcsec')\n", "ax.set_ylabel('Arcsec')\n", "\n", "ax.set_title('{} planets -- {} {}'.format(sp_sci.name, obs.filter, obs.image_mask))\n", "\n", "color = 'grey'\n", "ax.tick_params(axis='both', color=color, which='both')\n", "for k in ax.spines.keys():\n", " ax.spines[k].set_color(color)\n", "\n", "plotAxes(ax, width=1, headwidth=5, alength=0.15, angle=PA_offset, \n", " position=(0.1,0.1), label1='E', label2='N')\n", " \n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see, even with \"perfect PSF subtraction\" and no noise, it's difficult to make out planet _e_ despite providing a similiar magnitude as _d_. This is primarily due to its location relative to the occulting mask reducing throughput along with confusion of bright diffraction spots from the other nearby sources.\n", "\n", "**Note**: the circled regions of the expected planet positions don't perfectly align with the PSFs, because the LW wavelengths have a slight dispersion through the Lyot mask material." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimated Performance\n", "\n", "Now we are ready to determine contrast performance and sensitivites as a function of distance from the star. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1. Roll-Subtracted Images\n", "\n", "First, we will create a quick simulated roll-subtracted image using the in `gen_roll_image` method. For the selected observation date of 11/1/2022, APT shows a PA range of 84$^{\\circ}$ to 96$^{\\circ}$. So, we'll assume Roll 1 has PA1=85, while Roll 2 has PA2=95. In this case, \"roll subtraction\" simply creates two science images observed at different parallactic angles, then subtracts the same reference observation from each. The two results are then de-rotated to a common PA=0 and averaged.\n", "\n", "There is also the option to create ADI images, where the other roll position becomes the reference star by setting `no_ref=True`. \n", "\n", "Images generated with the `gen_roll_image` method will also include random pointing offsets described in the `pointing_info` dictionary. These can be generated by calling `obs.gen_pointing_offsets()`. " ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Pointing Info\n", " sgd_type : None\n", " slew_std : 5\n", " fsm_std : None\n", " roll1 : [0.00370446 0.0007631 ]\n", " roll2 : [0.00431872 0.0145655 ]\n", " ref : [-0.00801918 0.0003205 ]\n" ] } ], "source": [ "# Create pointing offset with a random seed for reproducibility\n", "obs.gen_pointing_offsets(rand_seed=1234, verbose=True)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "91bc818f88c94eebaaf0348d4265dd72", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/3 [00:00" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from pynrc.nb_funcs import plot_hdulist\n", "from matplotlib.patches import Circle\n", "\n", "fig, axes = plt.subplots(1,3, figsize=(14,4.3))\n", "xylim = 2.5\n", "xlim = ylim = np.array([-1,1])*xylim\n", "\n", "for j, wfe_drift in enumerate(wfe_list):\n", " ax = axes[j]\n", " hdul = hdul_dict[wfe_drift]\n", " \n", " plot_hdulist(hdul, xr=xlim, yr=ylim, ax=ax, vmin=0, vmax=10)\n", "\n", " # Location of planet\n", " for loc in loc_list:\n", " circle = Circle(loc, radius=xylim/15., lw=1, edgecolor='red', facecolor='none')\n", " ax.add_artist(circle)\n", "\n", " ax.set_title('$\\Delta$WFE = {:.0f} nm'.format(wfe_drift))\n", " \n", " nrc_utils.plotAxes(ax, width=1, headwidth=5, alength=0.15, position=(0.9,0.1), label1='E', label2='N')\n", "\n", "fig.suptitle('{} -- {} {}'.format(name_sci, obs.filter, obs.image_mask), fontsize=14)\n", "fig.tight_layout()\n", "fig.subplots_adjust(top=0.85)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The majority of the speckle noise here originates from small pointing offsets between the roll positions and reference observation. These PSF centering mismatches dominate the subtraction residuals compared to the WFE drift variations. Small-grid dithers acquired during the reference observations should produce improved subtraction performance through PCA/KLIP algorithms. To get a better idea of the post-processing performance, we re-run these observations assuming perfect target acquisition." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "9f349852fbbb450bb9d59c9fd2df943f", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/3 [00:00" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from pynrc.nb_funcs import plot_hdulist\n", "from matplotlib.patches import Circle\n", "\n", "fig, axes = plt.subplots(1,3, figsize=(14,4.3))\n", "xylim = 2.5\n", "xlim = ylim = np.array([-1,1])*xylim\n", "\n", "for j, wfe_drift in enumerate(wfe_list):\n", " ax = axes[j]\n", " hdul = hdul_dict[wfe_drift]\n", " \n", " plot_hdulist(hdul, xr=xlim, yr=ylim, ax=ax, vmin=0, vmax=10)\n", "\n", " # Location of planet\n", " for loc in loc_list:\n", " circle = Circle(loc, radius=xylim/15., lw=1, edgecolor='red', facecolor='none')\n", " ax.add_artist(circle)\n", "\n", " ax.set_title('$\\Delta$WFE = {:.0f} nm'.format(wfe_drift))\n", " \n", " nrc_utils.plotAxes(ax, width=1, headwidth=5, alength=0.15, position=(0.9,0.1), label1='E', label2='N')\n", "\n", "fig.suptitle('Ideal TA ({} -- {} {})'.format(name_sci, obs.filter, obs.image_mask), fontsize=14)\n", "fig.tight_layout()\n", "fig.subplots_adjust(top=0.85)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. Contrast Curves\n", "\n", "Next, we will cycle through a few WFE drift values to get an idea of potential predicted sensitivity curves. The `calc_contrast` method returns a tuple of three arrays:\n", "1. The radius in arcsec.\n", "2. The n-sigma contrast.\n", "3. The n-sigma magnitude sensitivity limit (vega mag).\n", "\n", "In order to better understand the relative contributes of WFE drift to contrast loss, we're going to ignore telescope pointing offsets by explicitly passing `xoff_* = (0,0)` keywords for Roll 1, Roll 2, and Ref observations. " ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "1b3117adcf1c418ca31cc1e641305c6d", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/3 [00:00" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from pynrc.nb_funcs import plot_contrasts, plot_planet_patches, plot_contrasts_mjup, update_yscale\n", "import matplotlib.patches as mpatches\n", "\n", "# fig, ax = plt.subplots(figsize=(8,5))\n", "fig, axes = plt.subplots(1,2, figsize=(14,4.5))\n", "xr=[0,5]\n", "yr=[24,8]\n", "\n", "# 1a. Plot contrast curves and set x/y limits\n", "ax = axes[0]\n", "ax, ax2, ax3 = plot_contrasts(curves, nsig, wfe_list, obs=obs, \n", " xr=xr, yr=yr, ax=ax, return_axes=True)\n", "# 1b. Plot the locations of exoplanet companions\n", "label = 'Companions ({})'.format(filt)\n", "planet_dist = [np.sqrt(x**2+y**2) for x,y in loc_list]\n", "ax.plot(planet_dist, pmags, marker='o', ls='None', label=label, color='k', zorder=10) \n", "\n", "# 1c. Plot Spiegel & Burrows (2012) exoplanet fluxes (Hot Start)\n", "plot_planet_patches(ax, obs, age=age, entropy=13, av_vals=None)\n", "ax.legend(ncol=2)\n", "\n", "# 2. Plot in terms of MJup using COND models\n", "ax = axes[1]\n", "ax1, ax2, ax3 = plot_contrasts_mjup(curves, nsig, wfe_list, obs=obs, age=age,\n", " ax=ax, twin_ax=True, xr=xr, yr=None, return_axes=True)\n", "yr = [0.03,100]\n", "for xval in planet_dist:\n", " ax.plot([xval,xval],yr, lw=1, ls='--', color='k', alpha=0.7)\n", "update_yscale(ax1, 'log', ylim=yr)\n", "yr_temp = np.array(ax1.get_ylim()) * 318.0\n", "update_yscale(ax2, 'log', ylim=yr_temp)\n", "# ax.set_yscale('log')\n", "# ax.set_ylim([0.08,100])\n", "ax.legend(loc='upper right', title='BEX ({:.0f} Myr)'.format(age))\n", "\n", "fig.suptitle('{} ({} + {})'.format(name_sci, obs.filter, obs.image_mask), fontsize=16)\n", "\n", "fig.tight_layout()\n", "fig.subplots_adjust(top=0.85, bottom=0.1 , left=0.05, right=0.97)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3. Saturation Levels\n", "\n", "Create an image showing level of saturation for each pixel. For NIRCam, saturation is important to track for purposes of accurate slope fits and persistence correction. In this case, we will plot the saturation levels both at `NGROUP=2` and `NGROUP=obs.det_info['ngroup']`. Saturation is defined at 80% well level, but can be modified using the `well_frac` keyword.\n", "\n", "We want to perform this analysis for both science and reference targets." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "NGROUP=2\n", "F444W Saturation Limit assuming Flat spectrum in photlam source (point source): 1.29 vegamag\n", "F444W Saturation Limit assuming Flat spectrum in photlam source (extended): 10.33 vegamag/arcsec^2\n", "\n", "NGROUP=16\n", "F444W Saturation Limit assuming Flat spectrum in photlam source (point source): 3.89 vegamag\n", "F444W Saturation Limit assuming Flat spectrum in photlam source (extended): 12.93 vegamag/arcsec^2\n", "\n", "HR 8799 flux at F444W: 5.24 mags\n", "HD 220657 flux at F444W: 3.03 mags\n" ] } ], "source": [ "# Saturation limits\n", "ng_max = obs.det_info['ngroup']\n", "sp_flat = pynrc.stellar_spectrum('flat')\n", "\n", "print('NGROUP=2')\n", "_ = obs.sat_limits(sp=sp_flat,ngroup=2,verbose=True)\n", "\n", "print('')\n", "print(f'NGROUP={ng_max}')\n", "_ = obs.sat_limits(sp=sp_flat,ngroup=ng_max,verbose=True)\n", "\n", "mag_sci = obs.star_flux('vegamag')\n", "mag_ref = obs.star_flux('vegamag', sp=obs.sp_ref)\n", "print('')\n", "print(f'{obs.sp_sci.name} flux at {obs.filter}: {mag_sci:0.2f} mags')\n", "print(f'{obs.sp_ref.name} flux at {obs.filter}: {mag_ref:0.2f} mags')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case, we don't expect HR 8799 to saturate. However, the reference source should have some saturated pixels before the end of an integration. " ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "# Well level of each pixel for science source\n", "sci_levels1 = obs.saturation_levels(ngroup=2, exclude_planets=True)\n", "sci_levels2 = obs.saturation_levels(ngroup=ng_max, exclude_planets=True)\n", "\n", "# Which pixels are saturated? Assume sat level at 90% full well.\n", "sci_mask1 = sci_levels1 > 0.9\n", "sci_mask2 = sci_levels2 > 0.9" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "# Well level of each pixel for reference source\n", "ref_levels1 = obs.saturation_levels(ngroup=2, do_ref=True)\n", "ref_levels2 = obs.saturation_levels(ngroup=ng_max, do_ref=True)\n", "\n", "# Which pixels are saturated? Assume sat level at 90% full well.\n", "ref_mask1 = ref_levels1 > 0.9\n", "ref_mask2 = ref_levels2 > 0.9" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "HR 8799\n", "0 saturated pixel at NGROUP=2\n", "0 saturated pixel at NGROUP=16\n", "\n", "HD 220657\n", "0 saturated pixel at NGROUP=2\n", "44 saturated pixel at NGROUP=16\n" ] } ], "source": [ "# How many saturated pixels?\n", "nsat1_sci = len(sci_levels1[sci_mask1])\n", "nsat2_sci = len(sci_levels2[sci_mask2])\n", "\n", "print(obs.sp_sci.name)\n", "print('{} saturated pixel at NGROUP=2'.format(nsat1_sci))\n", "print('{} saturated pixel at NGROUP={}'.format(nsat2_sci,ng_max))\n", "\n", "# How many saturated pixels?\n", "nsat1_ref = len(ref_levels1[ref_mask1])\n", "nsat2_ref = len(ref_levels2[ref_mask2])\n", "\n", "print('')\n", "print(obs.sp_ref.name)\n", "print('{} saturated pixel at NGROUP=2'.format(nsat1_ref))\n", "print('{} saturated pixel at NGROUP={}'.format(nsat2_ref,ng_max))" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "No saturation detected.\n" ] } ], "source": [ "# Saturation Mask for science target\n", "\n", "nsat1, nsat2 = (nsat1_sci, nsat2_sci)\n", "sat_mask1, sat_mask2 = (sci_mask1, sci_mask2)\n", "sp = obs.sp_sci\n", "\n", "# Only display saturation masks if there are saturated pixels\n", "if nsat2 > 0:\n", " fig, axes = plt.subplots(1,2, figsize=(10,5))\n", "\n", " xasec = obs.det_info['xpix'] * obs.pixelscale\n", " yasec = obs.det_info['ypix'] * obs.pixelscale\n", " extent = [-xasec/2, xasec/2, -yasec/2, yasec/2]\n", "\n", " axes[0].imshow(sat_mask1, extent=extent)\n", " axes[1].imshow(sat_mask2, extent=extent)\n", "\n", " axes[0].set_title('{} Saturation (NGROUP=2)'.format(sp.name))\n", " axes[1].set_title('{} Saturation (NGROUP={})'.format(sp.name,ng_max))\n", "\n", " for ax in axes:\n", " ax.set_xlabel('Arcsec')\n", " ax.set_ylabel('Arcsec')\n", " \n", " ax.tick_params(axis='both', color='white', which='both')\n", " for k in ax.spines.keys():\n", " ax.spines[k].set_color('white')\n", "\n", " fig.tight_layout()\n", "else:\n", " print('No saturation detected.')" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAAFbCAYAAAAumLNCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwEklEQVR4nO3deZxdBXnw8d+DIJBIRpGAL+WF2AClKMqmRusCJYolUOvSat0apQFbt1RkcanGlUi1rigYtfiiVV/cSVAkKuBCDAmNy4sWgwQrFgkoCYhQwef945yBkzt35tzJ3c7M/L6fz3xmznLPee655z73mbM8NzITSZIkSePbYdgBSJIkSU1n0SxJkiTVsGiWJEmSalg0S5IkSTUsmiVJkqQaFs2SJElSDYvmPoqIJ0fEhojIiLgsIuaW418VEZsiYktEfLoc98mIuDUiro+ISyNibUR8LiIOnWD5O0bEK8v5L4uIKyLimMr0/SLioxHxjYhYExGfGo2hMs8pEbG+/Dm1ZVpExGsj4jvlsq+MiIPLafPK53Bp5edllcduaJm2PiJ+MMFzeWa5jksj4tsR8ZmIOHKSm3ygIuKiiDiqj8tfEBEnRsSu5Xa5MyL+pWWeD0bEjeXr+9jK+L+MiG+WP5eX+8dLI2KXcnrr/nZlRFwdEYtalr9LRJxZLv/S8jV6SWX6oeW0jIjvVfaPg8vxN0bEuyJibvn4jIjvl/H8KCLOiYhdJ7ldjo+Ii8v9+qqIeEVl2mMi4u8nt6U13Zh7zb1dLt/cO/622SsivhYR57WZNtF++6qI+NPJrq9xMtOfPv4ARwEJ7Ngyfhnw7ZZxlwJvrQw/A7gVOGqcZc8DrgNGyuEnA7cDf1RZx7+WfwfwCeALlcc/FfgJsEv58xNgUWX6qcBZleHXA4+vrHvZBM/7gy3D/wS8cZx59wduqsS9A/BR4NUdbuNN422jHr6Oy4DzWsbNAaJP6zsQ+A6wa2XcRuBu4NEt854HzKsMnwj8CNi7Mu7PgbuAh0+wv70W2ArMqbwOXwHOAnaoPOdLgbe07IcJ7N9m/2zdZgksLP9+APCfwNsmuW02AoeVfz8EuAV4WmX6OcDz+7k/+NP8H8y9o8Pm3smtz9w7/rb5Y+AS4FOty+9gv30wsA7Yo5/7S79/PNLcYJn5eeDfgBUR0e61ug14Q2ZuKee/BLgTeFw5/Srg7HJaAp+mSO6jTgb+PTPvzMw7gU8CLwGIiJ2AU4C3VeJ5a2Z+u8PY/7Fl1HMpPjjaOQzYlJk3lI/9A/Buig+SxsrMreV27Yd3A+/LzN9Vxv0CuAD4WETcv92DImJO+dhXZeYvK7F+A/h34J4J1vllYDeKDw0oXrOHAa8pXxMycyvwD8BrI2L+9jyxSky3AyspCojJODsz/6Ncxo3AN4GnVKa/CXhPROzcTXyaucy95l5zb1tbgRMoCu5t1O23mXkLxTZ803aG3QgWzc13HsXRgMNaJ2TmLZl5/uhwRARwf2BzOf3LmXlt5SG7ADdXhh/FtsnxamD0tNxhFB8CLyhPMX0jIp7REsLjI2J1RHwrIt4ZEbPbPYGI+BPg9y2xVF0PHBYRT6s8tx9l5sry8Q+IiPPLU4ffiogPRcSO5bR/ozja+J7yFNQRo6esyukPLU9XbiqH962c0lpcPq+7ozjleVJEfLcc9/XKaaVnA4uBp5breF1EnFqe/lpWeZ5PLR9/eRSXD+xfjj8pitOpn46Ic6O4pOCi0dN1bbbXAyk+YC9vM/llwB7AP4+zLZ9C8Tp/s3VCZr4oM388zuMAdgT+APx3Ofxs4BuZuU2yL5fxC+CvJlhWp3YCfj+ZB2Tmu1tG7UK5z5fT/5tiP38K0vY7D3Ovufc+5t7Mm8t/8trpZL+9DHhWRNxv8uE2w47DDmAG+fpoMinNo9j562wqf88H1tfM+ySKJNjuDQ9wPOXRj9JeFKcgR20B9qzE90flep8EHAR8LyL+OzOvoHhzbABeR/Fm/z/lzzPbrPd5wPltxgOQmWsj4kPAFyLiJxRHZT5S+W/9/sDFmfkJgCiupfo74KOZ+aKIOBpYmpmXltOfQ3HqlMy8LiKWUnwAkpk/r0z/Q2b+eUScQnH6LICjM/OuKK6XOxd4QmZ+JoprseZl5uLRuCPiYZW//xj4LHB4Zl4TEc8HVkbEwzPzwxGxN7AEeHi5nX8APJ3iNFerQ4B7yuKvdVvdEhEnAxdExOcyc0PLLPOBzZk5qWRYOoFiO95QDu9fxtnODeW6tlu5TZ4JvKccXkzxAdlWZh7VZhlzKAqQl7ZM2ggcDlzYTYyaFsy94zD3jmHubaNd7m1jHhPvt1Dk5T2BfSjeL1OORfPgHJOZd48OlP8lL+zgcR2dDSj/c3478KLR0zkt0x9HsTMv6Sha2Jli//jX8jTYjyNiJfAi4IrytPirK8t/M3B1ROyVmb9qWdazgCdMtLLMfEVEfIDijXsicHpE/GVmfh34DbBfRHyb4rqyecAdFNfedeNL5brfVT6Hq4ELo7g5YifgEZNY1t8CazPzmnL4U8AKitO1ox+k38vM35Tr+hHw0HGWtRfF6d+2MvNLEfEZilOFj64LLIobTE6lOCp0Tma+pzL5BRHxBIpTgb9g21PIvdDuFOq7IuJWin3sbOBfATLzPMoP2El4B/DmzGxNwLdRbEfJ3DsBc+82zL3bb8L9tpxndNvuxRQtmr08o/nmlb83jjdDRATwYeDdmbmuzfT9gOXAs1r+C74JeGBleIT7TnP/pvxdTcK/oPgPsZ2fl7/3a1n3Y4FryuuZJpSZ12TmayluNvgk8MZy0t9RXAP4l+V/vOcBs+qW18H6tlTiHKG4xuvDmfkE4DnAZO4s3odtLxG4h2IbVrfX1srfd1IcxWknaJ/wql5BkXhObxn/M2DPqFx3l5mryu12M9u+3gDnZ+aTKE4N7wOcUZl2LbD3OOv/I+Cn5d93lb9bn8/OlWlVp2TmkzJzQWa+vfUUZKci4iSKU89nt5mcmN/UnXnlb3OvubdqxufecXSy345u2ymbm6ds4DPIYoo30IYJ5nknxX/aF0TEzhGx7+iEiNiD4vTc4sz8VRTXlY3eIHUl8CeV5RxcjqOyvj0r0+cCvyyX+5zqerjvqN4v2daEpwfLZT2mPD0E3Jv0vkjxQQLw6PL5/boc3mmi5QH/Uy539Hk+sGZ+KLbDHOCrHa6j1X9RbB/Kdd8PeBCdnQZudRPFTSHjKo+anERxfd1BlUlfA34HHNPucRMsbxPFh/uJEfGAcvT/Bf689fqziDiIIqF/oRy1mSJBtx69eSjFdulIFNc5XjreT8u8TweOBpaWwwe0LG43tk3e0mQtxtxr7m0x03PvBDaUv9vut6XRbTtlc7NFc4NFcRH9YuCkdqf9ynlOp0gy55VvuPnAi8tpDwA+B7wBuKkcfgnwv8qHnwP8bRT9IHehuGP3HIDM/AVFEvv7cll7Aosorp2DIlksroSyFLi0fNxobDtS3Byxsuap7gqcFOXNLFHcrf507ju1thF4ZPmhtCNjk9JtwKyIODoiXkmR+O6guIYN4C9q1g/FqaK7gceUw613FY+uIyLiC4z1KeDIKG9AobiR43rgux2su9X/A+4fLX1dW2XmqnK9j6mM2wK8CvjXiLj3P/yIOJDiCEXb/ah0Tjl9cTn8CeDHwFvLI2pExG7ABylaWG0q13k3xbXDL4/7epHuQnGdccfXFGfmeZl51Hg/lefyRIqjPS+jeE0eQNHaqGoe418TKE3I3GvunWimmZp7a5ZRt99CkZd/w/b9Q9MM2YC+d9P1h+IapQ0UpyQuA+aW419FcZPJFuDT5bhPUtwYcj1FL8Yrgc9T3Nww3vIPLJfd+rOsnP62cabPqyzjFIreieuBU1uWvyfFG+975c/fVaY9nKJ9zGXAmvLvh7Q8/jiKU25122kuxQ0JayjuPF5Lccpzt3L6bIr/rH9crucLwI0UrX2gKJ6uLh//sMq4jRTXzp1GcUruAmD3cr4st/PBlTheUr4uKylaByXwtXLa/hQJ9QqKa9ROLWPYBJxYzvMUiv6elwMXAweU459bzncjRcugkyqPfe442+TrwF+Vf+9axnpr+fv+lflGKI4ozGt5/F+VcVxWvrbfK9e9S5v9rdpX800UpxI/Uln3WeXjLy2f/4vbxLsHxU1EG8tt8B3ghS2v8aXlNv0+8Mou3lc3MHafvrQy/cHlc5s17Bzgz3B+MPeae9PcW3mNL6U3ufd+5bJGt+mlwL6d7Lfl9KXAx4adH7r5ifKJSGqQKL6N7N3Ak7NyE5PqRcRZFL1nPzjsWCRNLebe/ijPtlxK8Q/JlD3SbNEsNVREPBXYLzPPHXYsU0UUd7QfnZnvGHYskqYmc2/vlZczfTMz1w47lm5YNEuSJEk1vBFQkiRJqjG0LzeJiIcAbwUemZmPKsftTtF65WfAAcBrc2yzdqL4xp/DKL7L/VpPoUjS5JmHJalzw/xGwMdT3F17aGXc24HVmfl/I+IEih6YL6g+qGzl8mrgsMzMiLgyIr6RmT9FkjQZ5mFJ6tDQiubM/GwU3zFftYiiVQ8UbVM+3uahxwLr876Lsa+g6AW5TbLesnlr3nnHXWzZvLX18feaPTKL3265Y2DTAEbmzhk3pn6tc6LpMyGebh47XeLpZp3TJZ5+rbOTeA48cv7NVL6AoSn6nYehPhc3bT9pWjz9WudUiqcuppkQT9306fJZ1c94dpm1MyNz58S4C+/AMI80t7Mn9303+VbgQRGxY0vbl+o8o/NVv4EGgDvvuIuRPXYb86KtWrGai1asBmDpuSfznpPbn1HsxzSAs9cu56WPPqPttH6tc6LpMyGebh47XeLpZp3TJZ5+rbM67bglC1m0ZOG90yo56LfjBtU8PcvDUJ+Lm7afNC2efq1zKsVTF9NMiKdu+nT5rOpFPK15GMpcfPNt937V5fZqWtE8+hWWt1J8reZv2vRJvImi2fmoORRNvbexZfNWtmzeOuEbX5J66aLKP+Vw3wfHJX+44OYhhjVZPcvDYC6WNFiteRiKXAyw137dnfBrWveMVcBjy7//rBwmInaIiH3L8RcDR4x+tWQ5/1e2Z2VrLlw30GnDiKebmKZLPN0+djrE0806p0s8/VpnN/E0VGPycN1091vjmSnxdBPTVPqsGkY8kzG0Ps0R8STghRTfM/8h4F0UXxv5Doqvl5wPnJGZvyq/oef8zDykfOzzgSMp7tq+pt1d29esuzaBxh3dqDvtNWjGMzHjmZjxTKxypHk9Rc5qlH7nYWhmLm7qftIUTYsHmheT8UysifEAHHjk/Kl5TXNmXkbx3exVvwOWtJl3A3BIZfgTwCf6GZ8kTXfmYUnqXNMuz+ipVS3XtDRB02IynokZz8SMZ2JNi2dYmrYdjGdiTYsHmheT8UxsusYzbb9G+5p112aTTg1ImrmaennGIJiLJTXB2WuXT93LM/pt9sgslp57MlBcHL5m5fohRyRpJllw/BEsOOHeOrnbTkdTlrlY0rBU8/DskVldL88jzZLUZx5pNhdLGq5eHGme1tc0S5IkSb1g0SxJkiTVsGiWJEmSalg0S5IkSTUsmiVJkqQaFs2SJElSDfs0S1If2Ke5YC6WNCz2ae6QvUElNYV9ms3FkobLPs2SJEnSAFg0S5IkSTUsmiVJkqQaFs2SJElSDYtmSZIkqYYt5ySpD2w5VzAXSxoWW851yDZHkprClnPmYknDZcs5SZIkaQAsmiVJkqQaFs2SJElSDYtmSZIkqYZFsyRJklTDlnOS1Ae2nCuYiyUNiy3nOmSbI0lNYcs5c7Gk4bLlnCRJkjQAFs2SJElSjcZd0xwR84CvA/9VjpoD/CAzF1fmWQy8BLizHPXRzDx/cFFK0vRlHpaksRpXNAO3ASdn5mqAiHgTcEmb+Z6TmZsGGZgkzRDmYUlq0biiOTNvAUYT9c7AkZn5xjazviwibgRmAR/IzF8PMExJmrbMw5I0VqO7Z0TEi4DfZ+YnWsY/FLg9MzdHxHHAKZl5THWeX12/Obds3jpmmatWrOaiFav7GbakGeq4JQtZtGThmPEHHjn/emDewAPqgW7yMJiLJQ3WeHl4ZO4c9tpvblfdM5peNF8EPC0zfz/BPLsAtwM7Z+Y9o+NtcySpKaZyy7lu8jCYiyU1w7RuORcRRwPfHU3UEbF7RMwp/z4zIkYvLTkAuK41UUuSumMelqT7NO6a5oqTgJdXhs8Afg0sB24EPhQR1wGHAC8YfHiSNO2ZhyWp1NiiOTP/tmX4tMrf7x18RJI0s5iHJek+jb08Q5IkSWoKi2ZJkiSpRmMvz+jW7JFZLD33ZADWXLiONSvXDzkiSTPJguOPYMEJ9zbMGBlmLMNkLpY0LNU8PHtkVtfLa3TLuW7Y5khSU0zllnPdMhdLaoJp3XJOkiRJagqLZkmSJKmGRbMkSZJUw6JZkiRJqmHRLEmSJNWwaJYkSZJq2KdZkvrAPs0Fc7GkYbFPc4fsDSqpKezTbC6WNFz2aZYkSZIGwKJZkiRJqmHRLEmSJNWwaJYkSZJqWDRLkiRJNWw5J0l9YMu5grlY0rDYcq5DtjmS1BS2nDMXSxouW85JkiRJA2DRLEmSJNWwaJYkSZJqWDRLkiRJNSyaJUmSpBq2nJOkPrDlXMFcLGlYbDnXIdscSWoKW86ZiyUNly3nJEmSpAGwaJYkSZJqNPKa5ohYA9xZDt6Tmce0TN8FeCdwA3AAsDwzrxlslJI0fZmHJWlbjSyaga9m5rIJpi8Ffp6ZZ0XEIcBHgScMIjBJmiHMw5JU0dTLMw6JiNMjYllELGozfRFwBUBm/hB4ZETMGWiEkjS9mYclqaKpR5rfkZlrI+J+wOURcVtmXl6ZvidwW2V4azlu6+iIkblzOHvt8jELXrViNRetWN2nsCXNZMctWciiJQvbTdpj0LH0QNd5GMzFkgZrvDw8Mrf7/+kb33IuIpYDv8vMN1XGfQt4bWZ+qxzeCuyTmfcma9scSWqKqd5ybnvzMJiLJTXDtGw5FxEHRcSJlVEHABsjYvfKqb9VwGPL+Q8Bvt+aqCVJ28c8LEljNfHyjK3A8RGxNzAH+C/gU8By4Nfl7/cC74yI1wP7AyeOsyxJ0uSZhyWpReOK5sz8JfD0NpNOq8zzO+ClAwtKkmYQ87AkjdW4yzMkSZKkprFoliRJkmpYNEuSJEk1GndNc6/MHpnF0nNPBmDNhetYs3L9kCOSNJMsOP4IFpxwb5e5kWHGMkzmYknDUs3Ds0dmdb28xvdp3l72BpXUFFO9T3M3zMWSmmBa9mmWJEmSmsaiWZIkSaph0SxJkiTVsGiWJEmSalg0S5IkSTVsOSdJfWDLuYK5WNKw2HKuQ7Y5ktQUtpwzF0saLlvOSZIkSQNg0SxJkiTVsGiWJEmSalg0S5IkSTUsmiVJkqQatpyTpD6w5VzBXCxpWGw51yHbHElqClvOmYslDZct5yRJkqQBsGiWJEmSalg0S5IkSTUsmiVJkqQaFs2SJElSDVvOSVIf2HKuYC6WNCy2nOuQbY4kNYUt58zFkobLlnOSJEnSAFg0S5IkSTUad01zRMwH3gpcBewD3JKZb26ZZzHwEuDOctRHM/P8QcYpSdOVeViSxmpc0QzsDnw6M78EEBFXR8SqzGy9e+Q5mblp4NFJ0vRnHpakFo0rmjPzypZROwC/bTPryyLiRmAW8IHM/HXfg5OkGcA8LEljNa5oroqIpwMXZ+ZPWiZdBqzKzM0RcRxwAXDMwAOUpGnOPCxJhca2nIuIo4GnA0sz8w8TzLcLcDuwc2beMzr+V9dvzi2bt46Zf9WK1Vy0YnUfIpY00x23ZCGLliwcM/7AI+dfD8wbeEBd6jYPg7lY0mCNl4dH5s5hr/3mdtVyrpFFc0QsAp4AvAb4X8B+wH8Cd2fm1og4E/jnzLw7Ig4BPp+ZB1SXYW9QSU0xFfs09yIPg7lYUjP0ok9z4y7PiIgjgM8A64BvArOBsymOdvwaWA7cCHwoIq4DDgFeMJxoJWn6MQ9L0liNK5rLu7MfUDPPewcUjiTNOOZhSRrLLzeRJEmSalg0S5IkSTUsmiVJkqQajbumuVdmj8xi6bknA7DmwnWsWdn6RVaS1D8Ljj+CBSfc2zBjZJixDJO5WNKwVPPw7JFZXS+vkS3nesE2R5KaYiq2nOsVc7GkJuhFyzkvz5AkSZJqWDRLkiRJNSyaJUmSpBoWzZIkSVINi2ZJkiSphi3nJKkPbDlXMBdLGhZbznXINkeSmsKWc+ZiScNlyzlJkjStXfzLDcMOQQIsmiVJUoMdu/ehww5BAiyaJUmSpFoWzZIkSVKN2qI5Il4XEYvLv0+KiKP6HJMkqYW5WJKGq5OWcw/OzLcBZOaHI+IdwKV9jaoHbHMkaZj60HLOXCxJk9DrlnOdFM03tgzf2vVaB+C3W+7gPSefO+wwJM1Qa1auv7dAXLRk4ZYeLNJcLEmTUM3DZ69d3vXyOima50fEq4GNwHxgv67XKkmaLHOxJA1RJzcCngLsDvw98CDg1X2NSJLUjrlYkoao9khzZt4eEa8HHgzckpl/6H9YkqQqc7EkDVcn3TOeAvwM+CjwtxFxct+jkiRtw1wsScPVyeUZJwAHAd/JzE9SXEsnSRosc7EkDVEnNwL+IjPvjIgsh2/tYzw9Y5sjScPUh5Zz5mJJmoRet5yLzJx4hoiPAj+lOMKxDvjTzHxp12vus2vWXZsvffQZww5DkrjkDxesB46snXEC5mJJ2n5nr13OgUfOj26W0cnlGUuBOcAewEOA07tZoSRpuyzFXCxJQ9NJ94zbgNcCRMShwJ19jkmS1MJcLG3r4l9u4Ni9Dx12GJpBOume8emIWBARy4CzgQ/0O6iIWBgRH4yIZRHxxjbTd4mID0TEayLiYxFxYL9jkqRhGnQuNg+rKS7+5YZxh1unSf3UyY2A6zJzTUS8D3gccFo/A4qIWcA5wMMy866I+FxEHJOZX6/MthT4eWaeFRGHULRgekI/45KkIRtYLjYPq19ai9ztPVI8epTZolmD1Mk1zQ+NiBcAGzLzbmDXPsf0WOD6zLyrHP4OsKhlnkXAFQCZ+UPgkRExp89xSdIwDTIXm4fVF90UyXXLsYBWv3VSNK8GngGcGRHHd/iYbuwJ3FYZ3lqOm+w8kjSdDDIXm4fVN+MVvOMVvcfufeg2jxmv8Pb6ZvVbJ5dnfA34emZujYjbMnNln2O6CditMjynHDepeUbmzuHstcvHLHzVitVctGJ1byKVpIrjlixk0ZKF7Sbt0YPFDzIX9yQPg7lY7fWqwLVQVqvx8vDI3O5PhHVSNH8COB/4PPDEiDg4M9/a9ZrHdwWwX0TsXJ4a/DPggxGxO3B3Zm4FVlGcPvxWeS3d98vx99qyeSv2BpU0SBeNUwhe8ocLbu7B4geZi3uSh8FcLGmwxsvDZ69dzl77ze1q2Z0Uzd/LzM8DZObnI+JPu1pjjcy8IyL+AXhfRGwGfpCZX4+Is4BfA8uB9wLvjIjXA/sDJ/YzJklqgIHlYvOwBs0jxpoKOimaH1wz3HOZeQlwScu40yp//w5o/DdhSVIPDTQXm4claVud3EhyTUT8ICK+GBHfB67ud1CSpDHMxZI0RJ0caf4MxfVtBwE/zMz/7G9IkqQ2zMWSNEQdHWkGyMzPmqQlaWjMxZoxJmpBJw1LR90zMvNHowMRcXhmXtXHmHpi9sgslp57MgBrLlzHmpXrhxyRpJlkwfFHsOCEI0cHR3qwSHOxpqTRb+/r9bxSnWoenj0yq+vlRWZOPEPER4CbgR8DCZyQmX/d9Zr77Jp116ZtjiQ1wSV/uGA9cGTtjBMwF2sqGj1abNGsYTt77XIOPHJ+dLOMTo40HwF8EZhXDj+omxVKkraLuVhTzmQLYAtmNVknRfM/ZOYaKE4HAk/ob0iSpDbMxZI0RJ0UzT+IiL8HTqY4snFLf0OSJLVhLpakIRq3e0ZEHBYR5wCbgCcBGzJzf+CFA4pNkmY8c7EkNcNELecuB2YDB2fmC4BfANjqSJIGylwsSQ0w0eUZewPPA86IiG/RWU/nxrDNkaRh6mHLOXOxJG2HgbecA4iIRwPPBm4FHpqZL+56zX1mmyNJTdGLlnNgLpak7dWLlnMdHbHIzLWZeQrwHmDnblYoSdo+5mLNRH4zoJpiUqf5MvM24EV9ikWS1AFzsSQN3qSvjcvM/+lHIJKkzpmLNVP4hSdqiil1Q4kkSZI0DBbNkiRJUo1OvhFwSrLNkaRh6mHLuSnNXCxpWIbScm4qss2RpKboVcu5qchcLKkJBtZyTpIkSZrJLJolSZKkGhbNkiRJUg2LZkmSJKmGRbMkSZJUw6JZkiRJqmGfZknqA/s0F8zFkobFPs0dsjeopKawT7O5WNJw2adZkiRJGoBGXZ4REe8G7gBuBx4JLM3MG9vMtwnYVA7ekJnPG1SMkjSdmYclqb1GFc3AbzPz9QARcTrwOuDlbeY7LzOXDTIwSZohzMOS1EajiubRRF3ageJIRztPjIjTgN2Ar2Tmd/senCTNAOZhSWpv4EVzRFwM7NVm0hsy88vlPA8EngI8c5zFnJGZayNiFnBVRByfmRv7ErAkTTPmYUmavMZ1z4iIEeBDwOsy87oO5v80xVGOj1fH/+r6zbll89Yx869asZqLVqzuVbiSdK/jlixk0ZKFY8YfeOT864F5Aw9oO/UqD4O5WNJgjZeHR+bOYa/95nbVPaNRRXNE7AG8Bzg9M2+IiGdm5uciYgdgn8z8eUQcA+yUmV8tH7OO4kaVb1eXZZsjSU0xlVrO9TIPg7lYUjP0ouVco65pBr5GEdMnIwLgNuBzwCOA84FDgJuAZRFxOLA38Ll2iVqStF3Mw5LURqOK5sw8fJzxGygSNZn5Q8a/xk6S1AXzsCS155ebSJIkSTUsmiVJkqQaFs2SJElSjUZd09xLs0dmsfTckwFYc+E61qxcP+SIJM0kC44/ggUn3NswY2SYsQyTuVjSsFTz8OyRWV0vr1Et53rJNkeSmmIqtZzrNXOxpCboRcs5L8+QJEmSalg0S5IkSTUsmiVJkqQaFs2SJElSDYtmSZIkqYYt5ySpD2w5VzAXSxoWW851yDZHkprClnPmYknDZcs5SZIkaQAsmiVJkqQaFs2SJElSDYtmSZIkqYZFsyRJklTDlnOS1Ae2nCuYiyUNiy3nOmSbI0lNYcs5c7Gk4bLlnCRJkjQAFs2SJElSDYtmSZIkqYZFsyRJklTDolmSJEmqYdEsSZIk1bBPsyT1gX2aC+ZiScNin+YO2RtUUlPYp9lcLGm47NMsSZIkDYBFsyRJklSjUdc0R8Qy4KjKqLdl5iVt5ns+cBhwD3BtZp47kAAlaZozD0tSe40qmgEy86iJpkfEPsCrgcMyMyPiyoj4Rmb+dCABStI0Zx6WpLEaVzRHxOuAu4D7Ae/PzDtaZjkWWJ/33cF4BfAXgMlaknrAPCxJYw28aI6Ii4G92kx6A3ABsCkzfxsR/wi8HzixZb49gdsqw1vLcdsYmTuHs9cuH7OSVStWc9GK1dsZvSSN77glC1m0ZGG7SXsMOpaJDCoPg7lY0mCNl4dH5s7petmNbTkXEQcBX8nMh7aMPxF4XGaeWA6/D9iYme+rzmebI0lNMVVbznWbh8FcLKkZpl3LuYj4l8rgAcDGcvwOEbFvOf5i4IiIGH3ijwW+MrgoJWn6Mg9LUntNu6b57oh4L3ATcAjw0nL8I4DzgUMy8xcR8U7g3RFxD/ARbz6RpJ4xD0tSG40qmjPzNeOM30CRvEeHPwF8YkBhSdKMYR6WpPYadXmGJEmS1EQWzZIkSVKNRl2e0UuzR2ax9NyTAVhz4TrWrFw/5IgkzSQLjj+CBSfc2zBjZJixDJO5WNKwVPPw7JFZXS+vsS3numWbI0lNMVVbzvWCuVhSE0y7lnOSJElSE1k0S5IkSTUsmiVJkqQaFs2SJElSDYtmSZIkqYYt5ySpD2w5VzAXSxoWW851yDZHkprClnPmYknDZcs5SZIkaQAsmiVJkqQaFs2SJElSDYtmSZIkqYZFsyRJklTDolmSJEmqYZ9mSeoD+zQXzMWShsU+zR2yN6ikprBPs7lY0nDZp1mSJEkaAItmSZIkqYZFsyRJklTDolmSJEmqYdEsSZIk1bDlnCT1gS3nCuZiScNiy7kO2eZIUlPYcs5cLGm4bDknSZIkDYBFsyRJklSjUdc0R8QqYHZl1COAvTPzzpb5NgGbysEbMvN5AwlQkqY587Aktdeoohn4P5n5GYCI+GPg9NZEXTovM5cNNDJJmhnMw5LURqOK5tFEXXoF8P5xZn1iRJwG7AZ8JTO/2/fgJGkGMA9LUnsD754RERcDe7WZ9IbM/HI5zxyKoxjPGGcZj87MtRExC7gKOD4zN1bn+dX1m3PL5q1jHrtqxWouWrG626chSWMct2Qhi5YsHDP+wCPnXw/MG3hA4xhUHgZzsaTBGi8Pj8ydw177ze2qe0YjW85FxD8BP83MlR3M+2mKoxwfr463zZGkppiKLed6kYfBXCypGaZly7mI2AE4FlhVHRcR+5Z/HxMRT608ZH/g2sFGKUnTl3lYksZq1DXNpb8EVua2h8AfAZwPHALcBCyLiMOBvYHPZea3Bx+mJE1b5mFJatG4ojkzv9hm3AaKRE1m/hB45mCjkqSZwzwsSWM17vIMSZIkqWksmiVJkqQajbs8o1dmj8xi6bknA7DmwnWsWbl+yBFJmkkWHH8EC064t2HGyDBjGSZzsaRhqebh2SOzul5eI1vO9YJtjiQ1xVRsOdcr5mJJTTAtW85JkiRJTWPRLEmSJNWwaJYkSZJqWDRLkiRJNSyaJUmSpBq2nJOkPrDlXMFcLGlYbDnXIdscSWoKW86ZiyUNly3nJEmSpAGwaJYkSZJqWDRLkiRJNSyaJUmSpBoWzZIkSVINi2ZJkiSphn2aJakP7NNcMBdLGhb7NHfI3qCSmsI+zeZiScNln2ZJkiRpACyaJUmSpBoWzZIkSVINi2ZJkiSphkWzJEmSVMOWc5LUB7acK5iLJQ2LLec6ZJsjSU1hyzlzsaThsuVcjeOWLBx2CGM0LSbjmZjxTMx4Jta0eIaladvBeCbWtHigeTEZz8SmazzTumhe1LAXDZoXk/FMzHgmZjwTa1o8w9K07WA8E2taPNC8mIxnYtM1noFf0xwROwBLgLcAf56ZP6pMOxWYAzwI+FpmfrnN4+cB/wxsBOYBp2Tm7f2PXJKmD3OxJE3OMI40PxL4HnBHdWREPAY4OjP/Gfgn4F0R8cA2jz8HODczzwR+BJy+vYEsOP6IgU4bRjzdxDRd4un2sdMhnm7WOV3i6dc6u4lnyBqRi5u2nzQtnn6t03imVjzdxDSVPquGEc9kDLxozsz/yMwNbSYdD1xRzvN74MfAE6szRMROwNHAleWo7wCLtjeWyp3tA5k2jHi6iWm6xNPtY6dDPN2sc7rE0691dhPPMDUlFzdtP2laPP1ap/FMrXi6iWkqfVYNI57J6MvlGRFxMbBXm0lvaHear7QnRXIetbUcV7UH8Lu8r+VHu3kAGJk7h5E9duPstcu3Gb9qxWouWrG65hlI0uQdt2ThNtfO7funfzSag/YYRjzmYkkzTWsehiIXb7n5tq6XPbSWcxGxCTh+9Dq6iHgL8D+Z+ZZy+MvAR6qJvTy6cTuwS2ZmRBxeznN4m1VsBn4L3DxBGCPAlgFOg+LDZryY+rXOiabPhHi6eex0iaebdU6XePq1zk7i2Q+YO0FcQ9OAXNy0/aRp8fRrnVMpnrqYZkI8ddOny2dVP+OZTZd5eOA3Ak5gJfBGgIjYETgYuLwc/t/ALzPz9xHxTeBRwFrgz4BV4yyvkR9QktRw5mJJamPgR5oj4kHAS4FTgPOBf8/MNeW0Uynu1n4Q8JXRIxsR8V3gtMz8dnnH9huAnwH7Aq/yjm1JmhxzsSRNzrT9RkBJkiSpV5p0ecZ2aXqv0YhYRXEdzahHAHtn5p0t820CNpWDN2Tm83oVQ8t6lgFHVUa9LTMvaTPf84HDgHuAazPz3D7F826Klle3U7TAWpqZN7aZbxN93D4RsRB4BnATkJn5ppbpuwDvBG4ADgCWZ+Y1vYyhsq75wFuBq4B9gFsy880t8ywGXgKM7kcfzczz+xFPub41lXXdk5nHtEwf5PaZB3wd+K9y1BzgB5m5uDLPYvq8fSLiIRSv0yMz81HluN2B5RRHXw8AXpuZv2rz2IG8vwapybnYPFwbj3l4bCyNy8PlOs3F28Yw2DycmVP6p3zCh1K8kR9eGf8Y4KLy752AnwIPbPP4rwKPLv9+OfCWHsf37Mrff0zR17TdfMsGtL1q10ORIDZw35mIK4ED+hTPWyt/nw68f9DbB5hF8UG9czn8OeCYlnnOoDgtDXAI8K0+xvMo4GmV4auBI1rmWQzMG8Q+08n2H/D2eTCwsDL8JuDxg94+wLOAE4B1lXHnAH9T/n0CcH6bxw3s/TXInybnYvNw7brMw2PjaVwe7uQ1mGm5eNB5eMp/jXY2pNfoBPF9pjL4CuD948z6xIg4LSLeEhGP62UMrSLidRHx6og4PSJmtZnlWGB9lnsSxXb8i37EkpmvrwzuQHGko51+bp/HAtdn5l3lcLv9YBH37U8/BB4ZEXN6HAfl8q/MzC9VRu1A0X2g1cvK1/EN5X/W/XRIub8si4h275FBbp9bMnM1QETsDByZmd9uM2tft09mfhZo7WF073Zg/HwysPfXIDU5F5uHJ2YeHquheRjMxa0xDDQPT4nLM5rQa7Tb+Mqddt+snLJscUZmri2T51URcXxmbpxsLHXxABcAmzLztxHxjxQfHie2zLcn2+6E27VdOomnsn0eCDwFeOY4i+nZ9mmjk+c73jxbexRDWxHxdODizPxJy6TLgFWZuTkijqN4XY8Zs4DeeUe5/e8HXB4Rt2Xm5ZXpQ9k+wHOBT7UZP+jtM6q6HbYCD4qIHTPz7nHmGZ1vu99fg9TkXGwe3v54zMMTa1AeBnNxJ/qWh6dE0ZyZx27Hw24CdqsMzynHVd0M7BoRUSbrdvP0Kr4TgY9NsIy15e87ImIDRQun7UpGk9he3wBObTP+JmD/yvCc7Y2lk3giYgT4IPDizPz1OMvo2fZpo5N9pZN5eioijqY4+ra0dVpmXlcZ/Abw5Yi4X2be049YKtv/noj4VhlXNVEPfPuU/hp4WuvIQW+fitHtcCvFNvhNS6Ienadn769BanIuNg93F495uL0m5eFynebien3Lw1P+8owJrKQ43dO212j5ov0eGO01ChP3Gt1u5Q0yx1aXHRE7RMS+5d/HRMRTKw/ZH7i213GU6/qXyuABlDtJNR7gYuCIiIhy+LHAV/oUzx7A2cCpmXldRDyzNZ4BbJ8rgP3K00tQ7gcRsXvltNYq7tufDgG+n5l9+8+9PO12LPBK4CER8dhqPBFxZrlfQ/E6XtevJBQRB0VE9SjYAcDGYW6fcj1HA98t38cMa/u0uHc7UMknw3p/NUQjcrF5eMJ4zMNtNCkPl+szF3emb3l4yreciynQazQi/grYJzM/UBl3KMXF6YeUO/YyYD2wN8VdyWf2MobKes+kuOHiJoqbBN6QmddU4ynnez5wJMVdpddk/+7avorijMfokY3bMvOEQW+fiHgyxQ0Fm4HfZ+abIuIs4NeZuTwidqW4I/m/KT4s3p79uyP5CIpTWuvKUbMpPtAOrsTzSuDhwHUUr+N7R/f7PsSzd7n+qyj+G98JeBXF3ckD3z6VuD4FvDwzby6Hq69X37dPRDwJeCHwVOBDwLuAXYF3ANcD8ylOZ/9qWO+vQWp6LjYPTxiPeXhsLI3Kw2VM5uKx6x5oHp7yRbMkSZLUb9P58gxJkiSpJyyaJUmSpBoWzZIkSVINi2ZJkiSphkWzJEmSVMOiWTNCRMyJiNvLNkKSpAEzD2uqs2jWTPE84EvAScMORJJmKPOwprQp8TXaUg8cAPwTcHVEnELRZP1DFI36bwVeABwI/A3wJxRf67sAeD7wZGAhsAk4MjP/pvy2o/cAPwX2Ab6cmRdHxG4UzdV/BuxF0eD9LYN5ipLUaOZhTWkWzZr2IuIxwLcy86aI+Crw3Mz8cER8EZiVmadFxMeB/w28IjMfXj7uWRRnY14IfCEzPx4RjysX+xpgY2aeWX4D008iYj7w2nL8WeUyXjzI5ypJTWQe1nTg5RmaCZ4NHB4Ry4C7gJMr034MkJk/AA6iODJBOe6zmbmF4mtKHx8R64Cnlt9V/wjgoIg4A3gl8ENg93L8xsoyPtbH5yVJU4V5WFOeR5o1rZWn727NzDdXxl1buRGl+j3yPwQeWpnvmcDlwMMyc0lE7ARcBnwR+D5wY2a+r5z3BcAt5fj55bgATurk++wlaboyD2u6iMysn0uagiJiZ+DfgP/JzMXluIOA84Hbytl+A5yZmevK6ScBB1NcS7dDZr45It5OkdTvAPYDXgbsDJwF3ACMAD/LzA9VrqW7Hngg8JXM/Eb/n60kNY95WNOJRbMkSZJUw2uaJUmSpBoWzZIkSVINi2ZJkiSphkWzJEmSVMOiWZIkSaph0SxJkiTVsGiWJEmSalg0S5IkSTX+P9K9fxxEMiaeAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Saturation Mask for reference\n", "\n", "nsat1, nsat2 = (nsat1_ref, nsat2_ref)\n", "sat_mask1, sat_mask2 = (ref_mask1, ref_mask2)\n", "sp = obs.sp_ref\n", "\n", "# Only display saturation masks if there are saturated pixels\n", "if nsat2 > 0:\n", " fig, axes = plt.subplots(1,2, figsize=(10,5))\n", "\n", " xasec = obs.det_info['xpix'] * obs.pixelscale\n", " yasec = obs.det_info['ypix'] * obs.pixelscale\n", " extent = [-xasec/2, xasec/2, -yasec/2, yasec/2]\n", "\n", " axes[0].imshow(sat_mask1, extent=extent)\n", " axes[1].imshow(sat_mask2, extent=extent)\n", "\n", " axes[0].set_title(f'{sp.name} Saturation (NGROUP=2)')\n", " axes[1].set_title(f'{sp.name} Saturation (NGROUP={ng_max})')\n", "\n", " for ax in axes:\n", " ax.set_xlabel('Arcsec')\n", " ax.set_ylabel('Arcsec')\n", " \n", " ax.tick_params(axis='both', color='white', which='both')\n", " for k in ax.spines.keys():\n", " ax.spines[k].set_color('white')\n", "\n", " fig.tight_layout()\n", "else:\n", " print('No saturation detected.')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" } }, "nbformat": 4, "nbformat_minor": 2 }