{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e9d3db20-05e8-44fa-88cd-779331c0a727",
   "metadata": {},
   "source": [
    "# NIRSpec/G395H Flux Calibration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c98f94a4-9c0c-46b1-aed9-831483d158d3",
   "metadata": {},
   "source": [
    "This notebook will walk you through the additional steps necessary for flux calibration of NIRSpec/G395H BOTS observations.\n",
    "Note that the process to follow will be the same for the other NIRSpec gratings (or PRISM) -- we're using G395H here as an example since that was the grating used in the JWST Transiting Exoplanet Community Early Release Science (ERS) program. \n",
    "This notebook also assumes that you've already worked through the NIRSpec/G395H reduction notebook and have produced the relevant outputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bea041c1-d654-46d2-829e-7203da4a477c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# First thing to do is set your CRDS variables so we can download the most up-to-date JWST \n",
    "# reference files. \n",
    "\n",
    "import os\n",
    "os.environ['CRDS_PATH'] = './crds_cache'\n",
    "os.environ['CRDS_SERVER_URL'] = 'https://jwst-crds.stsci.edu'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2648e65c-2641-4c6b-b516-18953a23677d",
   "metadata": {},
   "source": [
    "Flux calibration converts the extracted stellar spectra into \"real\" units (e.g., erg/s/cm^2/µm) from the pipeline default \"data numbers\" (DN/s) or electrons (e/s). This invovles correcting the wavelength-dependent instrument throughput and applying the relevant unit conversions which were calibrated during commissioning using flux standard stars. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d541685-55ac-4d5a-b379-953dca77f26d",
   "metadata": {},
   "source": [
    "We start this tutorial with the WaveCorrStep Stage 2 outputs from the NIRSpec/G395H reduction tutorial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7576df84-6b03-4430-a9c9-a0bd811ea2f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the input directory containing the WaveCorrStep outputs.\n",
    "\n",
    "indir = 'pipeline_outputs_directory/Stage2/'\n",
    "\n",
    "# Lists of the WaveCorrStep outputs for each segment.\n",
    "\n",
    "filenames1 = [indir + 'jw01366003001_04101_00001-seg001_nrs1_wavecorrstep.fits',\n",
    "              indir + 'jw01366003001_04101_00001-seg002_nrs1_wavecorrstep.fits',\n",
    "              indir + 'jw01366003001_04101_00001-seg003_nrs1_wavecorrstep.fits']\n",
    "\n",
    "filenames2 = [indir + 'jw01366003001_04101_00001-seg001_nrs2_wavecorrstep.fits',\n",
    "              indir + 'jw01366003001_04101_00001-seg002_nrs2_wavecorrstep.fits',\n",
    "              indir + 'jw01366003001_04101_00001-seg003_nrs2_wavecorrstep.fits']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a127108-dcf5-4a69-8ac7-729366f4f881",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\"> <b>NOTE:</b> This tutotial starts with the WaveCorrStep outputs in Stage 2. However, you must also ensure to have run the GainScaleStep in Stage 1, or else the flux units at this point will be DN/s instead of the expected e-/s.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e09ba036-3dcc-4350-845d-7c33ba6410de",
   "metadata": {},
   "source": [
    "Now let's create some new output directories. In this case, we'll make subdirectories in the Stage2 and Stage3 directories for the flux-calibrated products, to not confuse them with non-flux-calibrated outputs from the previous tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97f0143a-f306-48ac-bae9-68e50b7c91fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from exotedrf import utils\n",
    "\n",
    "# Generate output directories\n",
    "\n",
    "utils.verify_path('pipeline_outputs_directory')\n",
    "utils.verify_path('pipeline_outputs_directory/Stage2')\n",
    "utils.verify_path('pipeline_outputs_directory/Stage2/FluxCalibrated')\n",
    "utils.verify_path('pipeline_outputs_directory/Stage3')\n",
    "utils.verify_path('pipeline_outputs_directory/Stage3/FluxCalibrated')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ceb4e235-5bbe-4112-a2c4-ff73daa6fffe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Now define the output directory paths\n",
    "\n",
    "outdir_s2 = 'pipeline_outputs_directory/Stage2/FluxCalibrated/'\n",
    "outdir_s3 = 'pipeline_outputs_directory/Stage3/FluxCalibrated/'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8390418c-36de-42b9-b23d-5d00d34c3764",
   "metadata": {},
   "source": [
    "Now we're ready to start the flux calibration. There are a few steps that are dfferent fromt the normal workflow when producing flux calibrated outputs, including applying the flat field correction and flux calibration reference files. We skip these in the normal workflow since we care only about relative measurements when constructing atmosphere spectra."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c29e8106-4bb3-404f-9123-fe8360f15479",
   "metadata": {},
   "source": [
    "## Flux Calibration Workflow"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "480a137d-8c94-4693-9cc6-3d5a39b48ad4",
   "metadata": {},
   "source": [
    "### Flat Field Correction Step\n",
    "Applies the flat field correction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d04fe00c-6513-4b19-b0d0-b51e796dd802",
   "metadata": {},
   "outputs": [],
   "source": [
    "from exotedrf import stage2\n",
    "\n",
    "# Run flat field correction.\n",
    "\n",
    "step = stage2.FlatFieldStep(filenames1, output_dir=outdir_s2)\n",
    "\n",
    "results1 = step.run(save_results=True, force_redo=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2697f729-3f9c-4aad-a075-2948903d094a",
   "metadata": {},
   "source": [
    "As with the previous tutorial, we're only going to run the steps for NRS1, and you can just swap the filename variable for NRS2."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a98bd1f-233e-4eae-9ee6-780a59b3d49a",
   "metadata": {},
   "source": [
    "### Photom Step\n",
    "Applies the wavelength-dependent flux calibration reference file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eee144a8-6919-436d-acd8-01cbafa8c7ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run flux calibration.\n",
    "\n",
    "step = stage2.PhotomStep(results1, output_dir=outdir_s2)\n",
    "\n",
    "results1 = step.run(save_results=True, force_redo=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0523bf5-6cd3-4e73-a47f-ab3b7b78616b",
   "metadata": {},
   "source": [
    "These are the two \"new\" steps that need to be run for flux calilbration. We'll now pick up where we left off in the original tutorial and proceed to extract the spectra. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6c4e9e3-dabe-45dc-859a-e3bc3f683b65",
   "metadata": {},
   "source": [
    "### Bad Pixel Correction Step\n",
    "Interpolate any remaining bad pixels in the data to be ready for the spectral extraction.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8992c6b6-2177-44d3-b3d5-6299b40bacb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Interpolate bad pixels.\n",
    "\n",
    "step = stage2.BadPixStep(results1, [150, -100], output_dir=outdir_s2)\n",
    "\n",
    "results1 = step.run(save_results=True, force_redo=True, do_plot=False, space_thresh=15, time_thresh=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b23d6a40-d239-42fd-b8a0-799314d2bfc6",
   "metadata": {},
   "source": [
    "### PCA Reconstruction Step\n",
    "Just to produce the data deep stack."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c24a773d-4922-4f5b-a016-528d9de59a47",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Produce a deep frame. We're not going to remove any components in this tutorial.\n",
    "\n",
    "step = stage2.PCAReconstructStep(results1, baseline_ints=[150, -100], output_dir=outdir_s2)\n",
    "\n",
    "results1, deepframe1 = step.run(save_results=True, force_redo=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "700f1dae-e9eb-47d2-88e2-6e61533eb0e9",
   "metadata": {},
   "source": [
    "### Tracing Step\n",
    "Get the trace positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd7692e1-3a22-440f-939a-f0cb552f28da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Trace the 2D spectrum on the detector. \n",
    "\n",
    "step = stage2.TracingStep(results1, deepframe=deepframe1, output_dir=outdir_s2)\n",
    " \n",
    "centroids1 = step.run(save_results=True, force_redo=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "979360a7-66f8-4225-b008-b4f78bc3c87a",
   "metadata": {},
   "source": [
    "### 1D Spectral Extraction\n",
    "Now we're ready for the spectral extraction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f93d71a1-9ca4-4e1d-880b-0f5326104015",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We now need the Stage 3 pipeline.\n",
    "\n",
    "from exotedrf import stage3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a55d0d8d-fd69-4f43-b5be-278b5a24c1ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "# Open the centroids file. \n",
    "\n",
    "centroids1 = pd.read_csv(outdir_s2 + 'jw01366003001_04101_00001_nrs1_centroids.csv', comment='#')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ca41bf5-8ae8-4114-8747-7178deb9503b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the spectra.\n",
    "\n",
    "step = stage3.Extract1DStep(results1, extract_method='box', output_dir=outdir_s3)\n",
    "\n",
    "# The exact width used here is not super important. Use something similar to what you would for the \n",
    "# normal workflow.\n",
    "\n",
    "spectra1 = step.run(extract_width=10, centroids=centroids1, save_results=True, force_redo=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "148d8b85-d57a-4c28-a3e9-a9ec9f9e687f",
   "metadata": {},
   "source": [
    "After the intial extraction, there is one more step to do -- which is basically just a unit conversion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94537057-9b5e-411a-8a98-dc07d83501b4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Convert from MJy to erg/s/cm^2/µm\n",
    "\n",
    "stage3.flux_calibrate('pipeline_outputs_directory/Stage3/WASP-39_nrs1_box_spectra_fullres.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5aa70f8-cc16-4c5a-9914-cf5220bcce93",
   "metadata": {},
   "source": [
    "This step produces another spectrum file, but with the FluxCalibrated tag.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06262619-ce20-42c4-ae72-7fbbd494ed75",
   "metadata": {},
   "source": [
    "Now you can just repeat all of the above for NRS2!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3725e849-cf1f-4608-bf0a-322f20f850ce",
   "metadata": {},
   "source": [
    "## Inspect Outputs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b007673d-8403-4116-9ddf-baac5e2fba5a",
   "metadata": {},
   "source": [
    "Let's unpack the extracted spectra and inspect the results. Here we're working with both NRS1 and NRS2 results, but you can just comment out the NRS2 stuff if you haven't run those steps yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4792569c-ba73-4382-9660-95e00f4cdfae",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "import numpy as np\n",
    "\n",
    "# Open the flux calibrated spectrum files\n",
    "\n",
    "spec_1 = fits.open('pipeline_outputs_directory/Stage3/WASP-39_nrs1_box_spectra_fullres_FluxCalibrated.fits')\n",
    "spec_2 = fits.open('pipeline_outputs_directory/Stage3/WASP-39_nrs2_box_spectra_fullres_FluxCalibrated.fits')\n",
    "\n",
    "# Get the wavelength axes.\n",
    "\n",
    "wave1 = spec_1[1].data\n",
    "wave2 = spec_2[1].data\n",
    "\n",
    "# Trim the bluest wavelengths from NRS1.\n",
    "\n",
    "ii = np.where(wave1 > 2.9)[0]\n",
    "wave1 = wave1[ii]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3a47745-2722-4758-b5c5-c22887485999",
   "metadata": {},
   "source": [
    "In general, we want to work with the out-of-transit stellar spectra, so that there is no bias or imprint of the planet. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b787bde-952b-46cf-818b-db0c4e05be45",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Unpack the extracted spectra and median over the out-of-transit baseline (first 100 integrations).\n",
    "\n",
    "spec1 = np.nanmedian(spec_1[3].data[:100, ii], axis=0)\n",
    "spec2 = np.nanmedian(spec_2[3].data[:100], axis=0)\n",
    "\n",
    "# We'll use the scatter along the time axis as the error.\n",
    "\n",
    "spec_err1 = np.nanstd(spec_1[4].data[:100, ii], axis=0)\n",
    "spec_err2 = np.nanstd(spec_2[4].data[:100], axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "405a2ae1-f329-4fec-b154-243688722924",
   "metadata": {},
   "source": [
    "For a point of comparison, let's use a PHOENIX model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8bcd92a-0905-4ed7-92c7-64743936bd08",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Download a PHOENIX spectrum using the grid point closest to the parameters of WASP-39.\n",
    "\n",
    "_ = utils.download_stellar_spectra(5500, 4.5, 0, './')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e7274f1-801b-4a17-bc45-cb79b8cb8c9b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Open the PHOENIX spectra and wavelength files. \n",
    "\n",
    "mod_wave = fits.getdata('WAVE_PHOENIX-ACES-AGSS-COND-2011.fits')*1e-4  # Convert to µm\n",
    "mod_spec = fits.getdata('lte05500-4.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits')*1e-4 # Convert to µm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "424d47fa-d391-4a30-ab5c-8901ae1d6849",
   "metadata": {},
   "source": [
    "The extracted stellar spectra are measured at Earth (or close enough), whereas the PHEONIX models effectively provide flux at the surface of the star. So we need one more conversion to allow for an apples-to-apples comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa2026c1-a525-49dd-8544-72d7f9d5157a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Rescale the PHOENIX model spectrum by the inverse-square law for flux. \n",
    "\n",
    "Rs = 0.939 * 6.96e8  # Radius of WASP-39 in m.\n",
    "d = 213.982 * 3.086e16  # Distance to WASP-39 in m.\n",
    "\n",
    "# Apply the flux rescaling.\n",
    "\n",
    "mod_spec *= (Rs/d)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db3f75da-019e-4400-a20d-03fd7bedb60d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy.ndimage import gaussian_filter\n",
    "\n",
    "# And let's smooth the model spectrum for visualization purposes. \n",
    "\n",
    "mod_spec = gaussian_filter(mod_spec, 25)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "371ff55d-4b27-4997-b758-66f3e885d386",
   "metadata": {},
   "source": [
    "Now let's plot up the model and extracted stellar spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "id": "414d2ed7-25cf-408c-969e-38fc7ee9acb3",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiAAAAFTCAYAAAAA17fuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABfoElEQVR4nO3deVyU5fr48c/MAAPILrKK4oK4I4LiviJmKZktSoumtmn5tWxRyzRLxUo7lpqaaZ4WRVtUcjdMKZcUzV3MHTRAUHZkm5nfH/yc4wgqMDBs1/v1mlfM/dzPc1/zNOfM1f3ci0Kn0+kQQgghhDAhZVUHIIQQQoi6RxIQIYQQQpicJCBCCCGEMDlJQIQQQghhcpKACCGEEMLkJAERQgghhMlJAiKEEEIIk5MERAghhBAmJwmIEEIIIUxOEhAhhBBCmFy1TkCio6MZMmQIHh4eKBQKNmzY8MBz8vLyeO+992jcuDFqtRpvb29WrlxZ+cEKIYQQotTMqjqA+8nOzsbPz48xY8YwbNiwUp3z1FNPkZSUxIoVK2jevDkJCQlotdpKjlQIIYQQZVGtE5BBgwYxaNCgUtfftm0be/bs4eLFizg5OQHg7e1dSdEJIYQQoryqdQJSVpGRkQQGBvLJJ5/w3XffUa9ePUJDQ/noo4+wsrIq8Zy8vDzy8vL077VaLTdv3qR+/fooFApThS6EEELUeDqdjszMTDw8PFAq7z/Ko1YlIBcvXuTPP//E0tKS9evXk5KSwvjx47lx4wbffPNNieeEh4czc+ZME0cqhBBC1F7x8fE0bNjwvnUUOp1OZ6J4jKJQKFi/fj1Dhw69Z52QkBD++OMPEhMTsbe3B+CXX37hiSeeIDs7u8RekLt7QNLT02nUqBHx8fHY2dlV+OcQQgghaquMjAy8vLxIS0vT/w7fS63qAXF3d8fT09PgQ7dq1QqdTsfVq1fx8fEpdo5arUatVhcrt7OzkwRECCGEKIfSDGGo1tNwy6p79+78+++/ZGVl6cv++ecflErlA7uChBBCCGE61ToBycrK4ujRoxw9ehSAS5cucfToUeLi4gCYOnUqI0eO1Nd/+umnqV+/PqNHj+b06dNER0fz9ttvM2bMmHsOQhVCCCGE6VXrBCQmJgZ/f3/8/f0BmDRpEv7+/kyfPh2AhIQEfTICYGNjw86dO0lLSyMwMJBnnnmGIUOG8MUXX1RJ/EIIIYQoWY0ZhGoqGRkZ2Nvbk56eLmNAhBBCiDIoy29ouQahhoaGliuwBQsW0LRp03KdK4QQQojao1wJyKZNm/D39y91D4FOp+OPP/7ggw8+KE9zQgghhKhlyj0Nd8mSJXTu3LlUdQsLC7GwsChvU0IIIYSoZco1CHXy5Ml4enqWur5KpWLy5Mm4u7uXpzkhhBBC1DIyCPUuMghVCCGEKJ+y/IZW62m4QgghhKidKmQp9vj4eDZs2EB8fDy5ubkGxxQKBZ9//nlFNCOEEEKIWsLoBGTdunU899xzaLVaXFxcig02lQRECCGEEHczOgF59913GTp0KF999dUDd74TQgghhIAKGAOSnJzMSy+9JMmHEEIIIUrN6ATkoYce4sCBAxURixBCCCHqCKMfwSxdupThw4eTk5ND//79cXBwKFanY8eOxjYjhBBCiFrE6AQkMzOTnJwcwsPDmTt3rsExnU6HQqFAo9EY24wQQgghahGjE5CRI0cSFxfHwoULadGihSy5LoQQQogHMjoBOXjwIKtXr2bo0KEVEI4QQggh6gKjB6H6+PhQWFhYEbEIIYQQoo4wOgH57LPPmD17NrGxsRURjxBCCCHqAKMfwbz++uskJibStm1bPDw8SpwFc/z4cWObEUIIIUQtYnQCEhAQgEKhqIhYhBBCCFFHGJ2ArFq16p7Hrl69SmJiorFNCCGEEKKWMToBiY6OvuexnTt38uWXX3Ljxg1jmxFCCCFELWJ0AtKnTx8UCgU6nc6g/PZjmd69exvbRJXIysrCzs6uqsMQQgghaiWjE5C///67WFl2djZ//fUXc+fO5e233za2iSrh6+tLZmZmVYchhBBC1EpGJyB+fn4llnfr1g0rKytee+01Lly4UK5rR0dH8+mnn3L48GESEhJYv359qRc827t3L71796Zt27YcPXq0zG1nZWVx9uxZfH19y3yuEEIIIe7P6HVA7qdFixZGDULNzs7Gz8+PxYsXl+m8tLQ0Ro4cSf/+/cvdNkCvXr2MOl8IIYQQJTO6B+R+1q5dy+OPP17u8wcNGsSgQYPKfN4rr7zC008/jUqlYsOGDeVq+/r167i4uBATE0NgYGC5riGEEEKIkhmdgISGhhYrKyws5J9//uHKlSsEBwfr6ygUCjZu3Ghsk/f1zTffcPHiRb7//ntmzZr1wPp5eXnk5eXp32dkZACgVqsB6NSpU7EBtkIIIYQwjtGPYDIyMsjMzDR43bp1Cy8vL3r06EFubq6+/PaPe2U5d+4cU6ZM4fvvv8fMrHS5VXh4OPb29vqXl5eXwfWgKKkRQgghRMUpVw9Ifn4+FhYWAOzevbsi4yk3jUbD008/zcyZM2nRokWpz5s6dSqTJk3Sv8/IyNAnIc2bN+fy5cs8/fTTjB49usJjFkIIIeoqha4czxdsbGwYMGAAgwcPZsiQIbi4uFRGbAYUCsV9Z8GkpaXh6OiISqXSl2m1WnQ6HSqVih07dtCvX78HtpORkYG9vT3p6en6dUAUCgXe3t5cunSpQj6LEEIIURuV9Bt6L+V6BPPLL7/g5eXFnDlz8PDwoHPnznz44YclrgliKnZ2dpw4cYKjR4/qX6+88gq+vr4cPXqUoKAgo65/+fJl/vzzzwqKVgghhKjbypWAhISE8MUXX3DhwgWOHTvG448/zs6dO+ncuTNeXl68/PLLbN68mdzcXKOCy8rK0icTAJcuXeLo0aPExcUBRY9PRo4cWfRBlEratm1r8HJxccHS0pK2bdtSr169csfRsWNHAHr27CkDUoUQQogKYPQg1DZt2jB58mT++OMPkpKSmDNnDunp6Tz77LPUr1+fwYMHs2zZsnJdOyYmBn9/f/z9/QGYNGkS/v7+TJ8+HYCEhAR9MlKZdu3aRUZGBh07dsTBwaHS2xNCCCFqu3KNASmNwsJCoqOj2bRpE1u2bCE2NrYymqlw93t+df36dVxdXcnIyMDW1raKIhRCCCGqp7KMAam0BKSmetDNUygUuLu7c/XqVZTKSl1IVgghhKhRKn0Qal2XkJBAq1atqjoMIYQQosYqUwJSUFDAf/7zH9544w127NhhcGzy5MkVGlh15+3tjVarreowhBBCiBqpTAnIuHHjOHToEF5eXrz11lu88cYb+mM7d+6s8OCqq8GDB7Njxw4mTJhQ1aEIIYQQNVKZxoD4+flx7NgxAHJzc3n++ecxMzPjv//9L4GBgVW6DkhFKc0YkFu3bmFlZQUUTQ329vY2cZRCCCFE9VNpY0Dy8/P1f1taWrJmzRrs7Ox47LHHDI7Vdubm5vz8888ANGnSpIqjEUIIIWqeMiUgTZo0MVgNVKFQ8OWXX9KyZcsaM83WWPv27UOlUqHRaPRl//nPf6owIiGEEKLmKdMjmJs3b6JQKHB0dCx27PTp07Ru3bpCg6sKpe0+Wrt2LSNGjADAxcWFpKQkU4UohBBCVEuV9gjGycmpxOQDqBXJR1k0adKEnj17AkULlN28ebOKIxJCCCFqDrOKuEh8fDwbNmwgPj6+2P4vCoWCzz//vCKaqVY6d+5MdHQ0AJ988gknTpygd+/eVRyVEEIIUTMYnYCsW7eO5557Dq1Wi4uLCxYWFgbHa2sCcqfmzZtz4cIFLC0tjd51VwghhKgLjF4J9d1332Xo0KGkpKRw7do1Ll26ZPC6ePFiRcRZrTk6OrJ161a6dOnCSy+9JDvmCiGEEA9gdAKSnJzMSy+9hL29fUXEUyM5Ojry008/AbB8+XKDmUJCCCGEKM7oBOShhx7iwIEDFRFLjaVWqwEYNmwYjRo1olevXuzevbtqgxJCCCGqMaN3w01NTWX48OF06tSJ/v374+DgUKxOx44djWnCpMoyhei2goICli9fzvjx4/Hz8+P48eMAaLVaFApFZYYrhBBCVBtl+Q01ehBqZmYmOTk5hIeHM3fuXINjOp0OhUJhsGhXbWRubs748eMBCAgI0CcgixYtkv1ihBBCiBIY3QPSp08fLl68yOTJk2nRokWxWTBAjZqeWp4ekDvl5+ejVquJjIwkNDSUvLy8Eu+JEEIIUduYtAfk4MGDrF69mqFDhxp7qVrB3NwcgEGDBgFgY2NDTk4OZmYVsuSKEEIIUSsYPQjVx8eHwsLCioilVlAoFBw/fhwzMzMaNmpBQUEBakvLqg5LCCGEqFaMTkA+++wzZs+eXWc2oyuNdu3aAeDz8FYAtBoNAcPWM25uQlWGJYQQQlQbRj8XeP3110lMTKRt27Z4eHgUmwWjUCg4duyYsc3UTIqi21u/cX+OrB+GneulKg5ICCGEqB6MTkACAgJkquk92ForsHZohr1bADeuRLF7aRPGNf6XJVPcqzo0IYQQokoZnYCsWrWqAsIoWXR0NJ9++imHDx8mISGB9evX33ew6y+//MKSJUs4evQoeXl5tGnThg8++ICBAwdWWoz3Y22ppONjv6A0s0Rt48GZqNc5de4Gs1aaM22Mc5XEJIQQQlQHRo8BiYqK4ptvvinx2KpVq/j999/Lfe3s7Gz8/PxYvHhxqepHR0czYMAAtmzZwuHDh+nbty9Dhgzh77//LncMxggLscPc0gGVmSVuLR7D3MqZP1a2Y8vv55i1MqVKYhJCCCGqA6N7QKZNm8ajjz5a4rHk5GSWL1/O3r17y3XtQYMG6aezlsaCBQsM3s+ZM4eNGzfy66+/4u/vX64YjBHayxaAL9alotVCwa2ipGP/991Q21yifXO1vo4QQghRlxjdA3Lq1CkCAwNLPNaxY0dOnTplbBPlptVqyczMxMnJqcpiCO1ly/895YitteE4mZRL21kRmVY1QQkhhBBVzOgERKFQkJ6eXuKx1NTUKl2Gfd68eWRlZfHUU0/ds05eXh4ZGRkGr4oW2suWjfO8aNn7I33Zye2vcCMlWR7FCCGEqJOMTkCCgoJYvHgxd6/ortPp+PLLLwkKCjK2iXJZvXo1M2fOZN26dbi4uNyzXnh4OPb29vqXl5dXpcXk3vpZg/d7/xvArpicSmtPCCGEqK6MTkBmzpzJvn37aN++PfPnz2f16tXMmzcPPz8/9u3bx6xZsyoizjKJiIjghRdeYN26dQQHB9+37tSpU0lPT9e/4uPjKy2ufoE2AFjZN9GXxe6ezKNvxRMZnVlp7QohhBDVjdEJSNeuXYmKisLOzo7Jkyfz7LPPMmXKFOzt7YmKiqJLly4VEWeprVmzhtGjR7NmzRoeeeSRB9ZXq9XY2dkZvCrLtDHO/Pvvv7z4zmp9WcKZCM4d/ZklP6dWWrtCCCFEdVOu3XDz8/NL3OH11q1bpKam4uDggLW1tdHBZWVlcf78eQD8/f357LPP6Nu3L05OTjRq1IipU6dy7do1vv32W6DoscuoUaP4/PPPGTZsmP46VlZW2Nvbl6pNY3fDLS1HjyBUFjbcuLILgD4vX6R/Z1tZH0QIIUSNVZbf0HL1gDg5OfHYY4+xYsUKkpKS9OVWVlZ4eHhUSPIBEBMTg7+/v34K7aRJk/D392f69OkAJCQkEBcXp6//1VdfUVhYyKuvvoq7u7v+NXHixAqJpyL9N+I32j9ctH6Kg0dX9ixvSdShbHkUI4QQok4oVw/Ijh072LRpE5s3b+bKlSt07NiRwYMHM2TIkCpZb6MimaoHBGDc3ASWTvWg14tnORP1BjbObWje6TV2LGpUqe0KIYQQlaHSe0BCQkL44osvuHDhAseOHePxxx9n586ddO7cGS8vL15++WU2b95Mbm5uuT5AXXF7TxiVmSVaTS6XDn7K1TM/y9RcIYQQtV65ekDu5ebNm2zevJnNmzezfft28vPz9cuhv/zyyxXVTKUyZQ/Infq8cpE9y5phZd+UoLBdKBQKevhZ8uHL955CLIQQQlQnld4Dci9OTk4899xzREREkJyczK+//kqLFi34z3/+U5HN1Eru9YsG9d5Kv8jupd5oNfnsPS49SEIIIWqnCu0BuS0nJ4d///2XZs2aoVAoHnxCNVLZPSCR0Zms2ZFBWIhdsX1g+o2PI/Xafo5GjsDc0ok+L/zNjoUyHkQIIUTNYNIekHnz5jFz5kz9+z/++ANPT098fX3x8fHhwoULxjZRq6zZkUHSTQ1rdhRf8t23kTmOnl3pPOI3CnJvcv7Afwibdk1mxgghhKh1jE5Avv76axo2bKh/P2nSJNq0acPGjRtxdnZm6tSpxjZRq4SF2OHqpCIspHhmuGSKOwoF1HP0AeDmtf1cTbjJF2tTJRERQghRq5gZe4H4+HiaN28OwLVr1zh8+DB79uyhZ8+eFBYWMm7cOKODrE1Ce9kWe/Ryp9sPxDo+9gvJF7fx58r29Bx7kqSbtqyITLvvuUIIIURNYXQPiJWVlX4H2aioKGxsbOjWrRsADg4O99wpV5RMbV70T3u3AKzsijbG2/ffzgBk5uhk3xghhBC1gtEJSOfOnZk7dy6bN2/m008/ZdCgQahUKgAuXLiAp6en0UHWJeMed8TVSYWroxLPtiPp9eJZNIU5pFzaCRQlISWNHxFCCCFqkgoZhJqQkMCQIUPIyspi9uzZ+mNr167V94aI0hkQVI+kmxqSUrX4eJmjMrME4MS2F8i6EQtQ4vgRIYQQoiYp1xiQOzeja926NRcvXuTGjRvUr1/foN78+fNxc3MzPso6xEqtxNVJRdJNDefiCwyOHds8ku4jD7J1X5aMBRFCCFGjVdhmdHcnHwDt2rWjQYMGxkVYh0RGZxI27RptmqpxdVIVO56fncTvSxpz+qKMARFCCFGzlSsB+eWXX/Dy8mLOnDl4enrSuXNnPvzwQ/7++++Kjq9Oub1GyKmLeayZ5Um/QGsUCmjc8VUAWg9YhHurEcT+/hbj5iZUcbRCCCFE+Rm9EuqpU6fYtGkTmzZt4sCBA7i5ufHwww8TGhpK//79sbS0rKhYTaKq9oKBe6+S+vNvKTwxoAH9x18hI+UMh9Y9hN+QH/jwrSH3XFVVCCGEMLWy/IbKZnR3qcoE5F50Oh379u3jhqY9CyJSiT/2Nef3fUT3kQewqOeOrbWCjfO8qjpMIYQQdZxsRlfL/PpHFou2egPg6qhEU3gLgL3fdkFTmMut3ArfzkcIIYSoVBWagNzJzMyMfv368dlnnxEbG1tZzdQJd+4fk5Ono37jftRv3A+A6+d/pVALs1amVHGUQgghROkZvRT7t99+e89jSqUSe3t7OnTogJeXPCIorzZN1SSn5dCmqZpDp29h69yG9g9/Q15WIkc2PI5zkxB2xcC0MVUdqRBCCFE6Ricgzz//PAqFAigaq3DbnWUKhYKhQ4fy3XffYW1tbWyTdc6pi3lotUX/HBvqwBfrUtFqQW3jRm7mVVLj/8Cl+WDCpl2TAalCCCFqBKMfwRw5coRmzZrx4Ycfcvz4cRITEzl+/DgffPABTZs2ZefOnXz11Vf89ttvTJkypSJirnPu3EE3tJct//eUI7bWCv2+MVk3z3L6t4lciYtnQUSq7BcjhBCi2jN6FszAgQMJDg7m7bffLnbsk08+Yfv27URFRfHxxx+zcOFCrl69akxzla46zoK5H7fmg0i6sE3/PmDYRuxcO6A2h62fN6rCyIQQQtQ1Jp0F8+eff9KhQ4cSj/n7+3PgwAGgaNO65ORkY5sTd/F/ZJnB+8O/PMqtjHjyCu5xghBCCFENGJ2ANGjQgJ9//rnEYz/++KN+KfbMzEwcHByMbU7cxcK8aKxN5+E76TvuCgAHfujBrbQLVRmWEEIIcV9GJyBTp07lq6++IiQkhIULFxIREcHChQsJDg7m66+/5t133wVg165ddO7cuUzXjo6OZsiQIXh4eKBQKNiwYcMDz9m9ezcdO3ZErVbTvHlzVq1aVY5PVXOMDXUAwM7OBgD3VmF4B77OsS1jWf97WtUFJoQQQtyH0bNgXn75Zdzd3Zk9ezZvvvkmhYWFmJmZ4e/vz8aNGxkyZAgAM2bMwNzcvEzXzs7Oxs/PjzFjxjBs2LAH1r906RKPPPIIr7zyCj/88ANRUVG88MILuLu7M3DgwHJ9vuoutJctO3bsoG/fQOZ+mwbMJTP5BJdjFjCsnyMbdqcTsTNTZscIIYSoVowahJqfn8+mTZvo0KEDTZs2RavVkpycTIMGDVAqK3aNM4VCwfr16xk6dOg960yePJnNmzdz8uRJfdmIESNIS0tj27Zt9zzvTjVtEOrdbuVp6fP8bg5G9NeX9X75AmZmZvzfU46ShAghhKg0JhuEamFhwdNPP01cXFzRxZRKXF1dKzz5KK39+/cTHBxsUDZw4ED2799/z3Py8vLIyMgweNVkO//KRqlSA+DSPBTvgImkXduHVgtLfk6t4uiEEEKIIkZnCi1bttQnIFUtMTERV1dXgzJXV1cyMjK4detWieeEh4djb2+vf9X0FVvX7MjAzMIWG+fWtOr/H5ybDODYpudIiP2J3HwdIRPiZI0QIYQQVc7oBCQ8PJxZs2YRExNTEfGY3NSpU0lPT9e/4uPjqzoko4SF2KG2dqDTk1tRKs2oV78VALG/v0lm8nEKNbBonfSECCGEqFpGD0J95513uHHjBkFBQdSvXx9XV1f9MuxQNHbj2LFjxjZTKm5ubiQlJRmUJSUlYWdnh5WVVYnnqNVq1Gq1KcIzGXMV5GmL/lYqzej+/BH2rurI4Z9D6fjYeuzdOlZtgEIIIeo8oxOQgIAAAgMDKyIWo3Xt2pUtW7YYlO3cuZOuXbtWUUSmt2ZHRrFFyCys6tP12X3s/74bR9Y/RuvgL3j0raFsnFezHzcJIYSouYxOQCpznY2srCzOnz+vf3/p0iWOHj2Kk5MTjRo1YurUqVy7dk2/I+8rr7zCokWLeOeddxgzZgy7du1i3bp1bN68udJirG7CQuxYEZlGZk7R5CZbawX5BTqw9aT3yxfYs6wZ8ceWY1O/Ff3G6/BtZM6SKe5VHLUQQoi6pmqmq5RSTEwM/v7++Pv7AzBp0iT8/f2ZPn06AAkJCQYDYJs0acLmzZvZuXMnfn5+zJ8/n6+//rrWrgFSktBetmyc58XrIxxxdVIxNtSBcY8X/R3cuWhKlJ1bAAfXDiD75j+cjSvg0bfiZQM7IYQQJmX0ZnS3bdu2jUOHDhEfH8+0adNo1KgR0dHRNG/eHA8Pj4powiRq+jogD2Jh5USP0UfJSb9M/NHltOg1Sz9mx9VJxZpZnlUcoRBCiJrKpJvRJScn0717dx555BFWrFjBihUrSElJAWDlypXMnj3b2CZEBRr0atGAYCu7RtyI28Xupd7cyijqRQoLqX0JlxBCiOrJ6ATk9ddfJzk5mZMnT3L+/Hnu7FAJDg4mKirK2CZEBRob6oCrk4r+nWzIz0nBy+8lDkYEoym4JaukCiGEMBmjB6Fu3ryZ5cuX06pVKzQajcExLy8vrl69amwTogKF9rLVJxrvj83H2qEpWk0e0V+35HHbExQq7YGiREUSEiGEEJXF6B6QwsJC6tWrV+Kx1NRULCwsjG1CVCKlmSXdnvsLgOj1k0jPzCUzR8eKyLSqDUwIIUStZnQCEhQUxMqVK0s8FhERQffu3Y1tQtwlMjqTsGnXKmTWisrMErWNGy16ziLl8k5id0+hsCCbzBydzIoRQghRaYxOQGbNmsWmTZvo1asXixcvRqFQsGHDBp588kkiIyOZOXNmRcQp7rBmRwZJNzWs2WH8xnlmZKFUgGfb52jRaw5J//zCH1+3pjAvQzavE0IIUWmMTkC6du3K77//jkKh4M0330Sn0zF79mwSEhKIioqiY0dZ9ruihYXY4eqkqpBZK/mFCvoEWAPg2eYZeow+ip1bAH9+04GszFTpBRFCCFEpKmwdEIBbt26RmpqKg4MD1tbWFXVZk6rt64DcSaFQ0Krff/BoNQyt9n/lOp2O3Uu9sXFuS8Cw9QQHOTBtjHPVBSqEEKJGKMtvqNGzYO5kZWV1z03fRPXj0dgPSxsX+nQsShZ3H8mhgb2SpNSibCQr5SR/rGiLTnuUa9cLZMl2IYQQFabClmIfNWpURV1KmMi1y0c5svFp2jdXc+piHv/3lCPcsZNxrxfPotXkcXLHOP2S7fJIRgghREWosATk999/r6hLCRO7c1BrWIgdavOicpWZJR1C13Azbjd/R4aRmpbBF2tlYKoQQgjjlSkBUalUJb6USiXXrl2rrBhFJbtzUGtoL1vyCv53zNGzGwBp1/Zxds+7aHUQMiFOekKEEEIYpUxjQNzd3fn7779p0KBBsWNeXl4VFpQwrTtXRy2JnWtHlCoLrp+PRKvJo3m3aXwegayUKoQQotzK1AMyePBgzp49W+Kx4ODgCglIVL1+gUWDUn0bmWNrrUCpNMPasTkAKZe2c+CHnmTeOEvwq9ITIoQQonwqdBpubVCXpuHez/fff8+zzz4LwI0bN1gQkcnKL9/n39Pf06LnR/zzx/sEPP4rdi7t2fVloyqOVgghRHVQlt9QSUDuIgnIvUVGZ7IgIhWtJp+M68eI/f0t2j+8CkeXpowbJpvXCSFEXVeW31CZhitKLbSXLf0CrVGqLHBw78St9Mv8taYPV2O3siAilbD3ZOdjIYQQpSPTcEWZTBvjzK4vG/HqEw76spPbXybl8m/6BcyEEEKIBynTLBiVSlViuU6nQ3HHAlai9nu4uw0AgU9sIv/WDY5vHkXzbu8za+X/sftIDn06Wsvy7UIIIe6pTD0g7u7uJCYmotFoDF5arRYPD4/KilFUQ1ZqJR27h6Gu546lTdG/+/hjX7PgwzC0WtgVkyMzZIQQQtyTTMMV5Xb4z9X8uaojzw9xAqBZ13e5ceU3Lhz4GJ1Oy4rINMKmXZNERAghRDEyC+YuMgumfHq/GIvK3JpLh/7D5ZgFqG086PrsPhQKBa5OKtbM8qzqEIUQQlSyKpkFI+o2lXnR4mVNOr0BQF7Wv/y94QkKclNJuqlh1sqUqgxPCCFENWNUAvLwww9XVBz3tHjxYry9vbG0tCQoKIiDBw/et/6CBQvw9fXFysoKLy8v3njjDXJzcys9zrru9uqpt7V/5L+kJ8bw5zcdSPv3L3bF5DBubkIVRSeEEKK6MSoBiYmJqag4SrR27VomTZrEjBkzOHLkCH5+fgwcOJDr16+XWH/16tVMmTKFGTNmcObMGVasWMHatWt59913KzVOUTQ99/URjrg6qWg3cCn27p0A8Gj9DH9vfIpLBz/jbFyB9IQIIYQAqvkjmM8++4wXX3yR0aNH07p1a5YuXYq1tTUrV64ssf6+ffvo3r07Tz/9NN7e3oSEhBAWFvbAXhNRMdbsyCDppoYmbR/GzLwe9Rv3p1nXqbTsO5+UK79x9cQqog5lS0+IEEKI6puA5Ofnc/jwYYPZNUqlkuDgYPbv31/iOd26dePw4cP6hOPixYts2bLlvo+K8vLyyMjIMHiJ8gkLscPVScXYUAdeH+HIgGf/S0g3V9xbPkHgE5u5fv5XzkS9TuyVfPrLRnZCCFGnlWkhMlNKSUlBo9Hg6upqUO7q6kpsbGyJ5zz99NOkpKTQo0cPdDodhYWFvPLKK/d9BBMeHs7MmTMrNPa6KrSXrcF+MLf/3ns8nrx86PjYz1w9/g3HNo+k7cClLIiAz9em0jdAFi0TQoi6ptr2gJTH7t27mTNnDl9++SVHjhzhl19+YfPmzXz00Uf3PGfq1Kmkp6frX/Hx8SaMuG6wuCPNbdh+NF7tx3J+70fodDp0uqJFy4QQQtQtRiUgXl5eFRVHMc7OzqhUKpKSkgzKk5KScHNzK/Gc999/n+eee44XXniBdu3a8dhjjzFnzhzCw8PRakvep0StVmNnZ2fwEhWrU2srlErwbWQOQP1GfUg4s4bdS73JyykaUCzjQoQQom4xKgE5fPhwRcVRjIWFBQEBAURFRenLtFotUVFRdO3atcRzcnJyUCoNP9Lt/WtkvbWqc+piHlotnLtaoC/r+Nh6AM5ETQLgbFwBj74VLyunCiFEHWHUGJDExER27NhBbGwsN2/eBMDJyYmWLVsSEhJyz56K0po0aRKjRo0iMDCQzp07s2DBArKzsxk9ejQAI0eOxNPTk/DwcACGDBnCZ599hr+/P0FBQZw/f57333+fIUOG3HMjPVH5wkLsWLMjgzZN1Zy6mEdYiB0LIjrSd9wVfl/SmFsZ8VjaNiQzBzJzNCz6MdVgLIkQQojap1wJSEFBAW+++SZLly5Fo9Hg7u6Oo6MjAKmpqSQkJKBSqRg3bhzz58/HzKx8ec7w4cNJTk5m+vTpJCYm0qFDB7Zt26YfmBoXF2fQ4zFt2jQUCgXTpk3j2rVrNGjQgCFDhjB79uxytS8qxt2DUwG27svibFwB3p0mcSbqdbJunMHc0pGAxyPBqj6zVqbIwFQhhKjFyrUXzOTJk/nqq6/45JNPeOqpp7C3tzc4npGRwbp163jnnXd46aWXmDt3boUFXNlkLxjTmrHsKh++4oX/o+sozM8k/thy2j20HDO1Ha+PcJSeECGEqEHK8htargTE3d2d8PBwnn/++fvW++abb3j33XdJSKg5AwwlATE9hUJBnxdOoDC3I+n8ryScWUujDi/j5NUT30bmpGVpCQuxk2RECCGqubL8hpbr2UhmZiYNGzZ8YL2GDRuSmSkDCsX9hYaGsnF5WwD6jhvMmajXSb36B+ZWzhQ+vRszC1sWrZNxIUIIUZuUaxZM165dmTNnDunp6fesk5GRwZw5c+jWrVu5gxN1w8aNG/V/t2xsgVJpAUDBrRT+WtMXTWEuhVpkhowQQtQi5eoBWbRoEf369cPLy4vg4GBatmyJg4MDAOnp6cTGxvLbb79ha2trMI1WiAfxdClaK8S95XAS//kZ2wZ+RC/3xafHB+jajmLRulTW7MiQRzJCCFHDlWsMCEBaWhpLlixh27ZtxMbGkpqaCoCjoyMtW7Zk0KBBvPLKK/rEpKaQMSBVK/i1OH5f5kuzoCk0bD8aTWEu0ct9AXD1eZTWwV8AoDaHrZ83qspQhRBC3KUsv6HlXojMwcGBqVOnsmfPHpKSksjPzyc/P5+kpCT27NnDlClTalzyIapen47WoNNhZqagX6A1ZuaWOHr1wn/oT+h0Wn5f0piMpKPkFcjqqUIIUZPVqr1gRM13e+2PQg2cvJCHTgcdBn+Hg3snVGZWAFw69BnJl3ZwNq6AfuPjGDc3QcaHCCFEDVOpCUhKSgrR0dGV2YSohRRoUSgVtG2mxtZaoS9v0ulNgp7eg0+PmZzc9iKndryKTqvhbFwBSTc1rIhMq7qghRBClEmlJiB79uyhb9++ldmEqIW0Oh3olJy6mMfYUAdcnVT0C7RGbeOGtb03VvaNAbh+YRO7lzUlO/U8AJk5OukFEUKIGsKovWCEqBQ6DXY2ZvqZLrdnu1y7nsDZuAIUCiU+PT/k3B/TATgY0Z/6jfvTbtAKFkSkcvx8nn7PGZkpI4QQ1VO5EpD27duXql5GRkZ5Li/qOJ1Ox0uPORVLHpZMcWfWyhR+P5yDudoBgL7jrpByJYpLf83j4l8f06zLFHbF5ACwIjJNEhAhhKimypWAnDlzhjZt2uDv73/feleuXCE+Pr5cgYm6TaFQlFh+6mLRwFTnJgMZPH4vSisFNO5PbkYc5/78ACu7RqjruVO/cV8yc3SMm5vAkinuJo5eCCHEg5QrAWnbti0+Pj5888039633888/s2fPnnIFJuq25RvScWmRWawHIyzE7v8PNrVibGg7ABZEpGJhXbRD8tk9UwFo2ecTGjQdxNk4Ox59K56xoQ7SGyKEENVIuRKQoKAgtm7dWqq65VznTNRxGTk61uzIKJY03Dkm5E5rdgxBU5BD7O9vgkJJ7O53yEm/hGeb5wBPFkSksiIyTRIRIYSoJso1C+btt99m4cKFD6z38MMPc+nSpfI0IeqwyR9+R4u2fQgLKd1KtKG9bFkzyxOFsujr7NZiGABxfy9h//fduHpiFYV5GWTm6FgQkVppcQshhCi9ci/FXlvJUuw1V//hi9m17jUa+Y+jadA77F7ahFb9PuPMrkkAtH/kW+o36o1CAX0DrPWLngkhhKgYJlmKXYjqpkMLSwA82zyHQlH01ba0K9ovpk3IEi7+9TF/relHXs4NdsXkMGtlSpXFKoQQdV25EpA5c+aQkFC2fTjmzJlDUlJSeZoTolTiC/oQ9PQerO096RdoTdv+n2Bb3xe/R1bRoMlAfHvNISftAntXdSTlShS7YnLoNz5OEhEhhKgC5XoEo1KpOHDgAJ06dSpVfY1Gg4WFBYcOHaJjx45lDtKU5BFMzRUZncmaHRn6BcjCpl0j6aYGW2sF1pZKklM1aHWQlnCIvzc8gaVdIzoM/h4r+8aoLRSMGyYDVIUQwhhl+Q0tVwKiVCrp0aNHqXe71el0bNmyRRIQYVJ3JyS3DZoYR9qNeA780AMAnx4f4Nl2FAqFkn6BMjZECCHKq9ITkD59+txzoaj7+eqrr/Dx8SnzeaYkCUjdMGhiHKnJl0Ch4OzuyWgKcmjZdx429X3xbWQui5cJIUQ5VHoCUptJAlJ3jJtbtLdMQW4af37jh6NXL5wb9cOj7XMEd7aTnhAhhCijWjULZvHixXh7e2NpaUlQUBAHDx68b/20tDReffVV3N3dUavVtGjRgi1btpgoWlGTLJniTr9Aa8wtHeg77gp+D6/i3N4P2LOsGb/9lU6/8XGEvXe1qsMUQohaqVonIGvXrmXSpEnMmDGDI0eO4Ofnx8CBA7l+/XqJ9fPz8xkwYACXL1/mp59+4uzZsyxfvhxPT08TRy5qimljnNn1ZSP6BVqjUKro/fIFFAoVe75qzrk/Z3ItKYNBE+OqOkwhhKh1qvUjmKCgIDp16sSiRYsA0Gq1eHl5MWHCBKZMmVKs/tKlS/n000+JjY3F3Ny8XG3KI5i6KzI6kwURqVw7+S3/nl5D1o3TAPR6IRaVuZWMDRFCiAeoFY9g8vPzOXz4MMHBwfoypVJJcHAw+/fvL/GcyMhIunbtyquvvoqrqytt27Zlzpw5aDQaU4UtarDQXra8PsIRz7Yj6fTUVjoP34mX30tEf92SmJ8f5dT5VEImxBEZnVnVoQohRI1XbROQlJQUNBoNrq6uBuWurq4kJiaWeM7Fixf56aef0Gg0bNmyhffff5/58+cza9ase7aTl5dHRkaGwUvUXbeTEFcnFR07tMGjzTMAZF4/yh8r2hD7x2w++yGJ4FclERFCCGOUazfc6kqr1eLi4sJXX32FSqUiICCAa9eu8emnnzJjxowSzwkPD2fmzJkmjlRUZ3fuuBsZbcPn9lfQ6nTsXRVA/LGvuHp8BV2eiWZBBCyISMXVUcma2Q2rOGohhKhZKiQB0Wg0/PXXX1y9epXc3Nxix0eOHFnmazo7O6NSqYot356UlISbm1uJ57i7u2Nubo5KpdKXtWrVisTERPLz87GwsCh2ztSpU5k0aZL+fUZGBl5eXmWOV9ROt5ORfuPj6P78YbJunCb7RizHt76AbYO2ODcORuc9gH7j43h9hKOspCqEEKVkdAJy5MgRhg0bRnx8PCWNZ1UoFOVKQCwsLAgICCAqKoqhQ4cCRT0cUVFRvPbaayWe0717d1avXo1Wq0X5/7dm/+eff3B3dy8x+QBQq9Wo1eoyxyfqln6B1uyKycHWuQ02Ti05s2sS2TfOkBj7IwqlBQHDfmFBRDsWrUulvoOq2OqrQgghDBk9CyYoKIjc3Fy++OILWrduXeIPvb29fbmuvXbtWkaNGsWyZcvo3LkzCxYsYN26dcTGxuLq6srIkSPx9PQkPDwcgPj4eNq0acOoUaOYMGEC586dY8yYMfzf//0f7733XqnalFkw4kFuL2AGoCnMJXq5LwAuzQbTrNt7WNp4AGBrrWBsqOwvI4SoO8ryG2p0D8ipU6f48ccf6d27t7GXKmb48OEkJyczffp0EhMT6dChA9u2bdMPTI2Li9P3dAB4eXmxfft23njjDdq3b4+npycTJ05k8uTJFR6bqLuWTHEnMjqTJT+nkoclbUK+5NSO8bj6DOX45tF4th2JZ5tnyMzRsSAiFUCSECGEuIvRPSAdO3Zk8uTJDB8+vKJiqlLSAyLKatbKFDZuPYhN/VbsXtYUAKWZJc2CpuDZ7nmDfZNkszshRG1m0r1g9uzZw8SJE4mIiKBly5bGXKpakARElFdkdCbPPBaIq++TXNg/W1/eduBSnL1DUCj/NzhaFjUTQtRGlZ6AtGvXzuC/6hISEkhNTcXDwwMHBwfDBhQKjh07VtYmqowkIMJYg16PJzdPQ17WvyRf3Mb5fR9h7dCMpkHv4NxkoP5/O2pzKNBAn47SKyKEqB0qfQxIQECAQQIihPifccMcWLMjgzadWrBDbYtWW4ijRxeuX9xCQuyPePm9gKNnV/KKxrGy+0gO08ZUbcxCCGFq1XovmKogPSCiIkVGZ7IiMg0Aa7WCuGuJXDm8CK0mj8b+47GybwyA2kLBuGEyY0YIUbOZdC+YMWPGcOnSpRKPXblyhTFj5D/tRN0V2suWjfO82DjPizWzG7J3VSA9hnyETlPAgdW9OPzzUHIzr5GXXzRjRnbeFULUFUYnIKtWrSI5ObnEYykpKfz3v/81tgkhapWIOV50HvIfAHKzrrH/+27E7p5Cfk4KeQXQb3wcs1amVHGUQghRuSpkKfZ7jQc5d+4c9evXr4gmhKhVxoY6kHHzAApLd25c2cXxLaPJuH6U5l3fw8mrJ7tictgVU9QbIlN3hRC1UbnGgCxZsoQlS5YAcPLkSZo1a4aVlZVBndzcXC5fvsyTTz7JDz/8UDHRmoCMARGm1P/VOArzc8i6eZbjm0ZSmF+0G7OjZ3d8+8zFyq6RQX1JRoQQ1Vmlz4Lx8PAgICAAKEpAfH19adCggUEdCwsLWrVqxdixY8vThBB1Qt8Aa3YfARevjmgKizZydG85nITYtRz4oSfNu03Ho3UYKnNrAHbF5PD74Tj6BkgiIoSo2YyaBZOXl4ezszNHjhzBx8enIuOqMtIDIqpCZHQmkye/TeuOj5Kqa0NawiFif3+LW+mXAbBxbku7QV9jafO/xcsUCpg4XHbgFUJUHybbC0atVpOfn8+1a9dqTQIiRFUI7WVL6P6ld5Q0oveYxkR/0wkAd98n2P9dFyysGtB6wBc4enZDp4NF61IlARFC1EhGz4Jp1aoVcXEydVCIiubg5Kr/u2H70XQY8gP5t5I5GhlG9NdtuPL3UnJzs3g5PKEKoxRCiPIxOgEJDw9n1qxZxMTEVEQ8Qoj/b2yog/5v30bmWNk1wqKeK20fWo6mIAsru0bs/747u7etQqPRVF2gQghRDkavhNquXTv+/fdf0tLSqF+/Pq6urgbTcmUvGCHK7+DBgwQFBaHT6Rg3N4HYK/n4NFQQey4BM2t3Ui7tJPPKOhq75LF169Zis9GEEMKUTDYGBIr2hQkMDDT2MkKIEnh5efHSSy8BkJalRaFQkHlLhaOzB5k5Opq0CWHj1rF8/PHHWFtbs27dOmJiYnj//fexsbGp4uiFEOLeZC+Yu0gPiKiuIqMzWbMjg7CQou/l7b9vD0J9+eWX+eqrrwCwtrYmICCAiRMnMmzYMNk8UghhEibdC0YIYRqhvWxZM8uzaMbMHX/fZm1dtFbIkCFDyMnJITIykm+//RalUslbb711zy0ThBCiKlRIAvLdd9/Ro0cPXFxcsLOzK/YSQlS+UaNGMXPmTD766CMAHBwc+Pnnn8nIyKBfv364uLjQrVs3vv32W0rq+IyMziRs2jUiozNNHboQog4yOgH5/vvvefHFF2nbti0pKSk89dRTPP7441hYWODi4sJbb71VEXEKIR6gQ4cOTJ8+HT8/P7Zu3QqAmZkZtra2PPzww0RGRvLGG2+wb98+mjdvzscff8yVK1f056+ITCPppoYVkWlV9AmEEHWJ0WNA/P39eeKJJ5gyZQrm5ubExMTQsWNHMjMzCQkJ4cknn2TSpEkVFW+lkzEgoi7QaDSsW7eOzz//nMzMTA4dOsTgt/63A6+ZCl57UlZZFUKUjUnHgJw7d47u3bujUqlQqVRkZBRtpmVra8vkyZP54osvjG1CCFHBVCoVYWFh7Nu3jxYtWrBw4ULORr/HtZPfodUWUqiBBRGp9Bsfx6yVKQ++oBBClJHRCYi9vT15eXkAeHp6cvr0af0xjUbDjRs3jG1CCFFJlEolP/zwA1OmTMHOwQ2APcuacXj9MG5c+R2ttpBdMTn0Gx8nY0OEEBXK6HVAAgMDOX78OAMHDiQ0NJSZM2ei1WoxNzdn7ty5dOnSpSLiFEJUEmtra+bPn8+kSZN49K14ki/vwLXZYDKuH+P4ludx9XmUFj1nsSCiqP7d03+FEKI8jB4DcuDAAa5cucLw4cNJS0tj1KhRbN68Ga1WS6dOnVizZg1NmzYt9/UXL17Mp59+SmJiIn5+fixcuJDOnTs/8LyIiAjCwsJ49NFH2bBhQ6nbkzEgQhQZNDGOtBvxnNn1JukJf1HPqSU+PWdi18APlbkVvo3MWTLF/cEXEkLUGWX5Da2Uhcjy8vLIy8sz+gd87dq1jBw5kqVLlxIUFMSCBQv48ccfOXv2LC4uLvc87/Lly/To0YOmTZvi5OQkCYgQ5TRrZQpRh7LZu6ojWk0+moIsbOq3xrPtSOo37odzAzc2zvOq6jCFENVElScgFSUoKIhOnTqxaNEiALRaLV5eXkyYMIEpU6aUeI5Go6FXr16MGTOGP/74g7S0NElAhDBSt+f2Y27pxN7/BtJzzHGObXme1PhoGjR5iBdfGMUbLwXj7OxssFqrPKIRou6p9Fkwn3zyCYmJiQZl+/btIycnx6Ds0qVL+n0syio/P5/Dhw8THBysL1MqlQQHB7N///57nvfhhx/i4uLC2LFjS9VOXl4eGRkZBi8hhKGH+/qgMrei1wunUChVqMyscPTsjmf70cx5L4wmLXvj6hPKOzOWcPXfG/oZNP3GxxH8qgxgFUIUV64EZOrUqcTFxenfazQaevbsSWxsrEG969evs2LFinIFlpKSgkajwdXV1aDc1dW1WPJz259//smKFStYvnx5qdsJDw/H3t5e//Lyku5kIe7WvrkatTkoFODbyJymQW/j03Mmjh5d6PVCLJ2e2opL8yHkZsTz5zd+HN8yhuSL29DpdGh1RVN6+4+PI/g1mdYrhChSrgSkpKc2Vf0kJzMzk+eee47ly5fj7Oxc6vOmTp1Kenq6/hUfH1+JUQpRM62ITCOvAGysFCyZ4s6QAX7Uc/QBQGVuBUCDJgNpGvQ23UcdpnX/BST+8wu7l3pz+JehJMT+hA7QapFpvUIIoAKm4VYWZ2dnVCoVSUlJBuVJSUm4ubkVq3/hwgUuX77MkCFD9GVarRYoWo767NmzNGvWrNh5arUatVpdwdELUbtNG+NM++aZLPk5lbwCw2MW1s70C7RGx6cUdJlCYX4Wx359htjf36TtwGWozK2xc+3IgoiinhGZTSNE3VRtExALCwsCAgKIiopi6NChQFFCERUVxWuvvVasfsuWLTlx4oRB2bRp08jMzOTzzz+XRytCGGFsqIN+cOltt3flhaKN7BZEpAKgNi9KUKaNcdYfm+ewjz9WtOXk9pdxbNgDhUKFUqXGveVT6HTBDJoYR14BkowIUYeUOwFRKBSlKjPGpEmTGDVqFIGBgXTu3JkFCxaQnZ3N6NGjARg5ciSenp6Eh4djaWlJ27ZtDc53cHAAKFYuhKhYtxORu5OU28eOn3fFwvIKPg3NORtX1GVSmJdB/Ilv+OfPGTg37oeDZze02hD6jS8aX6ZUwG+LG5n2gwghTKbcCUjfvn1RKg2HkPTs2dOg7PYjkPIaPnw4ycnJTJ8+ncTERDp06MC2bdv0A1Pj4uKKxSCEqHhrdmSQdFPDmh0Z95xee2ePyN2KekSK/n70rXgyc3SYqe1oEjiRy4c+IyP5ONdOfQdA/cb9cfUZikvzwfpk5E79Aq31vStCiJqrXOuAzJw5s0z1Z8yYUdYmqoysAyJEcRW5vkdktOHYkdSre3Hw7Ers72/j3CQEB48u/BM9jZRL26nvHYyjZzfcfJ9AZWapv4attcLgsZCsOSJE9VBrFiKrCpKACGEadycitykVUFhYwJ6vmtPpqW2kXtvP+b0zcfTsjmuLoTRo+jBmFjYlXlN6R4SoWpKAGEESECGqxp29LGu2p5N4U4Obk4o2zSz5decxLG29iF7ui41zG9x9n6SwIBu3FkNR23iWOP5MBrQKYXqSgBhBEhAhqqdZK1P4aEIbnL1DSL22j2Zdp3Jy24v49g7HyasXaf/+hWPDHqjrGS5eaGutkP1qhDARSUCMIAmIENWbl98L3IzbTecRv7F7aRNcmoeSfysFW+c2xB9bjkKhoqHfC9Rv3A97t0CUyuJj7V8f4SjjRoSoBJW+F4wQQlQVXWE2KvN6mKuU7Ny5k+vnI/F0SOel/5sDQMt+87Cp35LUq3/y94YnOLL+ca4e/4bsm//or3F7r5pH35KVj4WoKtV2ITIhhCiJg0sLzCzrY2WpIDg4mO+//55Bgwbh5OTEuEdTcHJyYvzHiZyNK+DItf2kJ8bg5vsE107/QE7qBRQKJQ3bj8XRsyuZORYGU31dHZWsmd2wCj+dEHWH0Y9gxowZw/vvv0+TJk2KHbty5QozZ85k5cqVxjRhUvIIRojqrSxTgse8s4FTcVY4ubYirwBy0i/z1+re+uOebUdhZmFDg6YPY9ug+IKFMpBViLIx6RgQpVLJgQMH6Ny5c7Fjhw8fpnPnzmg0GmOaMClJQISovW4vGZ+dep78nGRy0i5iZd+YG5ejyE47j3Pj/jh4dKWeUwuDmTUKBUwcLuNGhHiQsvyGVsgjmHstwX7u3Dnq169fEU0IIYTR/rdkvC9hIZ0I7WXLrJUp7Irpwa30KySd20Ds72+SmXwChUKFU6PeeAe+ga1zGxZEpOr3u7mTDGgVonzK1QOyZMkSlixZAsCpU6do2rQpVlZWBnVyc3O5fPkyTz75JD/88EPFRGsC0gMiRN3Wb3wcR9Y/gbVjc1Kv/kFu5lUAVOb1cPYOxtGzO3ZuHbF2aF7sP77MVPDak5KQiLqr0ntAPDw8CAgIAODkyZP4+vrSoEEDgzoWFha0atWKsWPHlqcJIYSoEq+PcGSN03pupGaTm5vPHyvbAeDW8km8O04gLeEv/j29hpzU85ip7XDy6kODJiGYqe0o1FBiT4ms0CpEcUaPARk9ejTTp08vcRBqTSQ9IEIIKBovsiIyjf0bxmGudqBFr9m8EeYEwIrINDJzdBTkppEQu5a4v5dhbulAPSdfbJxb4+DRBXvXjiiUKv311Oaw9XPZ3VfUbrIQmREkARFC3GnWyhR2H8mhT8d792KEvhlHemYuWk0e6QkxZCYfJ+P6MZRKM+rVb4lL04excW4ts2pErWfyBOS7775j2bJl/PPPP+Tm5pYYUE0hCYgQorzCpl0j6eb/Zv3l5ySTmXKahNOr8WjzLE5ePR94DXlcI2oyk86C+f7773nxxRd5/vnn2bdvH2PGjEGj0fDrr7/i4ODAyJEjjW1CCCFqhLAQu7vWKGnErJWNWactIDPlJLcyrvBP9HsG59jUb027h1eiVJpjbuXErpgc2jfPlIGsotYzugfE39+fJ554gilTpmBubk5MTAwdO3YkMzOTkJAQnnzySSZNmlRR8VY66QERQlS0R19eT+RXw2jV/z9k3Ygl/uiye9Zt0vltvANe0783U4GVWsHYUAdJSkS1Z9IekHPnztG9e3dUKhUqlUr/uMXW1pbJkyfz+uuv16gERAghKtKtPC0JqUXLFLi1GAZA867votNqyM36l8zkE1y/sInkC5sBuHTwU+xcO5B2bR+aghys7L1p0Oxhwpdd5rMf3FGqLABQWygYN0ySElFzGZ2A2Nvbk5eXB4CnpyenT5+mT58+AGg0Gm7cuGFsE0IIUWNZqZUo/n/SAIYLl0VGO7Nmhzdh7z0DwOwlBzm4NoRjvz5jcI1zf84weG9TvzUtes1iQUQAx8/ncej0LW7l6tDooIWXOWlZ2lItVS9EVTI6AQkMDOT48eMMHDiQ0NBQZs6ciVarxdzcnLlz59KlS5eKiFMIIWqsIQPaY21/ln6B1gZJQWgvW4P3S35pATrtA6+XdeM0R9YPI/DJLWzf2xCdpoDM5BPYuvgRe9keTUE2CyI0fLE2lf+TJeRFNWX0GJADBw5w5coVhg8fTlpaGqNGjWLz5s1otVo6derEmjVraNq0aUXFW+lkDIgQoqo8+lY8kfMb0XPMCf49vRp1PTecGvUmIfZHHBt2J+H0GlL/PUBO6jkA6jfuz40rUSVey96tE827TUOhVGFl742jgx0b53mZ8uOIOqjK1wHJy8sjLy+vRv6ASwIihKgqkdGZfP71Dia+EHLHY5qi3X/bNFWz91gOeQWQ+u8Bjm4cXqZr9x13BZ1Oh06bDyhRKM0MlpKXMSWiIlR5AlKTSQIihKjuIqMzefPd+ZzfOxPf3nO5cGAuhXlpODbsgYNHFy4dnFeq63i0fpZ/T39P6+CFODcJQWVmCciqraL8JAExgiQgQoiawO+hLzi+fSIBj0cy9OGu/HUyE6XSjE6trfj9wDXOHfwSN98nyEw+Tuzvb5Xp2j1GH0NlXg+dToNSpQZ0aAtzUZlbA+DqqOR6mhYLc+k1EYZqVQKyePFiPv30UxITE/Hz82PhwoV07ty5xLrLly/n22+/5eTJkwAEBAQwZ86ce9YviSQgQoiaIDI6k9Xb03l6oP19E4DI6ExefW0iV098g4NHF/JzkslJu3Dfa1vaepGbGV+s3Mq+Ce0GfQ06LTptIdaOPihV5igU0DfAmvbN1XctxCbqmlqTgKxdu5aRI0eydOlSgoKCWLBgAT/++CNnz57FxcWlWP1nnnmG7t27061bNywtLfn4449Zv349p06dwtPTs1RtSgIihKhtWnSbTFrSCdqELMHSQkFWVjr/RE+jWdd3uXzoM3x6zKQwPxMLq/oc/fVZ0v7dX6rrWtl70+nJrWTf/AelmRVKMzXZN87QoOkgzFRQ305JUmrRrJ7b++DMWpnCrpgcfZlMGa5dak0CEhQURKdOnVi0aBEAWq0WLy8vJkyYwJQpUx54vkajwdHRkUWLFpV6SXhJQIQQtc3GPRms2ZHO0wOLHpfcPbA1v7CoB2PaGGeaBr3NpYPzsLJvikKpwtUnlOvnfiU79Z8Sr+3q8yhJ5zYWK7dz8cfRqwcJZyJo2edT6jfue98Yba0VWFsqJRmp4UyagCQmJuLm5nbP40eOHKFjx45lvm5+fj7W1tb89NNPDB06VF8+atQo0tLS2Lix+Bf+bpmZmbi4uPDjjz8yePDgEuvcnrFzW0ZGBl5eXpKACCHqpLD3rnIlLh5LW0/UFgq83cyIvZyDQqGgIC8LpZklGYmHSf13P9k3Ykm5vLPU13ZpHkqLnh9SkJvGxYOf4tV+DBbWrljZeaHV5KMpvIW52t7gHJmdU7OYNAFp0KABS5Ys4YknnjAo12q1fPTRR8yZM8fgB760/v33Xzw9Pdm3bx9du3bVl7/zzjvs2bOHv/7664HXGD9+PNu3b+fUqVNYWlqWWOeDDz5g5syZxcolARFCiJJFRmeyIjKNzBwdaQmHuJV2EXOr+pzYOhaAoKf3cHBNP3Q6zQOuVESpUqPVFP1OqOu54dPjA5QqS8zUdti7BdzzPDMVNPM059zVAvp0lF2EqwOT7gXz+OOP89RTT/HMM8+waNEi7O3tOXv2LCNHjuTEiRN88sknxjZRLnPnziUiIoLdu3ffM/kAmDp1qsFeNbd7QIQQQpTszhVcI6MdWbMjAwcbJc7el1ApFeh0Cvq8cpHUa/tx8OiCpYWCvAK4dGgBl2P+g8rcBp22QJ903P4nQF52Iie3v6J/X8+xhf7xj7WjD05evTBT2+PW4jHMLOw4fcECFEp2xcCumLhisd4eICvJSfVTIWNAtm7dygsvvIBKpeLpp59m4cKFtGnThm+//ZaWLVuW65rGPIKZN28es2bN4rfffiMwMLBM7coYECGEqByR0Zn8Z9lmnOulc/zEaRq1G0GvoObEnMklPnYbWalx2Ht04Z8/pmPj3Jp/T30PFCUet1d/LYljw57Uc2yOe6swNIU5mKsdKMi9gbmlExZWzpip7/3/5a6OSlAoZOxJBamSQajHjh2jS5cu5Ofn4+/vz4EDBzAzM66DJSgoiM6dO7Nw4UKg6LFOo0aNeO211+45CPWTTz5h9uzZbN++vVz70EgCIoQQVef2AFkHGyWxl3O5fnEL5mp7zu55l/qN+uDo1ZOMpL/R5GeRm3WN9MQjFOal3eNqCszU9hTmpeHg2Y2C3FTM1Q7odBpa9PwQa0cfFAoV6LQolKoSr2BrrSA7VyePeErJ5AnIDz/8wIQJE3Bzc2Po0KF89tlndOnShVWrVuHt7V3u665du5ZRo0axbNkyOnfuzIIFC1i3bh2xsbG4uroycuRIPD09CQ8PB+Djjz9m+vTprF69mu7du+uvY2Njg42NTanalARECCGqlztn7Zy6mGewLL0CHfm5aZir7clJv0Jm8glO/zaBBk0fJj/nOumJMaVqw6f7B+SkXaAgNxXf3uHw/xMSM/N6aDUFKFXm9z2/X6AkKGDiBOTJJ5/kl19+YcKECcydOxdLS0uOHz/Oc889x6VLl/jss8944YUXyn39RYsW6Rci69ChA1988QVBQUEA9OnTB29vb1atWgWAt7c3V65cKXaNGTNm8MEHH5SqPUlAhBCi5pq1MoWogzepZ6kg65aS/JwkLG0bcmbXm7i3GsHfG56gdfBCTv82Qf9Pc8v6FOTewMLahfyc6yVeV13PHTvXDjh6dsPOLYDCvHTUNh4oUKC29UShUBnsrXM3pQLqWSkYG+rA1n1ZnI0r0K+NUpuYNAHx9vbmm2++oW9fwzneBQUFvP/++8yfP5+CggJjmjApSUCEEKJ2ubMH5UBMLF0CW/LbbztQWrfg2M53aDtwKRf2z8HO1Z8zUa/TJmQJqdf2kn0jFoXSjLR/D9z3+kozS9TWrlg7+pCRdJiW/eaTkfQ3rj6PolCaY2Zhh7mlAzqdFqWy9EMTFAqYONyxRo1NMWkCkpmZia3tvW/O/v37DabRVneSgAghRN1wOzG5PQD19vL2bZtZcuj0LQDGhjrww/rDXLyazc2rf2Jp2xBNfhZ52YlcjllA44D/48rhL0rVnrP3AFIu78S391yyb56lYfsxKJUWoFCg1eRjZeeFpuAWKnOrB17L1VHJmtkNjfr8lcGkCUh0dPQD6/Tq1cuYJkxKEhAhhBCl8eHX19l9JJeC3DSUKguyU89hYeWMTqfBzMKWhNh1uDQfwv7vutK8+wzO752pT0Lu5/ZGgM2CpoACrB2ak3H9KB6twlCaqVEozfU7Fz+ImRJ2LDLdzsYmTUCUSiUKhYI7L3P3czCNpnSL0VQHkoAIIYSoKLNWprBxy594eAeQkpyIup4rBblpmKntyctOoDAvg7ysBG5e/YP8nGTs3QM598d0bF38yLx+7L7X9mw7EjO1PZY2nqQl/EVj//FoNXlY2jXCzMKWwvzMYivLArw+ovIe65g0ATl2rPgNSk1NZfv27fz8888sW7as2PiQ6kwSECGEEKZ056MgoGiV2WwtZ35/G98+c0k4E4Gz9wDSE2PISbuInYsfxzY9R9OgyVw6NB+dthC1jQd5Wf8CRZsE3kq/jJV9E+o5Nict4RDotNRv3A+1jQfmZireHtuGVq1a0bRpU65fv46Pjw+2trYUFBRgbn7/GT/3U202o/voo484deoUERERldVEhZMERAghRHVS0liVrzemolAoKCy4RXaeGWYqM/7ZO5smnd7g7J6pmFvVx6X5EE7vfI1OT27lj5XtcGvxONZOLbh4IBwbGxuysrIM2nF0dMTJyYkLFy7QvXt3XnrpJY4fP07//v1RKpV4eXmhUCho0KABOp2OBg0aFIu12iQgUVFRPPbYY2RkZFRWExVOEhAhhBA12S+7brIuKptbeZCRrcXSQkFuvhaFQgkUbfC3dYEXly9fZsuWLbzwwguo1WpOnjzJ5MmT2bx5MytWrODTTz/F1dWVPXv2lNiOj48P586dw8LCgpYtW9KtWzeaNGnC5MmTqz4BmTBhAr/++iuXL1+urCYqnCQgQgghaqO7e1JKcnvMpkr1v5Vhk5OTsbW1RaVS8fvvv+Pu7k56ejrHjx+nS5cujBs3jqysLJo1a8avv/4KlG5DV6MTkNDQ0GJl+fn5nD17lri4OD755BPefPNNY5owKUlAhBBCiPIx6W64GRkZxWa9WFpaEhwczBNPPMHAgQONbUIIIYQQtYzRCcju3bsrIAwhhBBC1CXKqg5ACCGEEHWPJCBCCCGEMLlyPYJp167dfXf9u5NCoShxsTIhhBBC1F3lSkACAgJKnYAIIYQQQtytXAnIqlWrKjgMIYQQQtQl5RoD0r59e06ePGlQtnr1atLS0ioiJiGEEELUcuVKQE6ePElOTo7+vUaj4bnnnuPixYsVFpgQQgghaq8KmwVTiSu6CyGEEKKWkWm4QgghhDC5cicgJc2CkZkxQgghhCiNci/F3rdvX5RKw/ylZ8+excoUCgXp6enlbUYIIYQQtVC5EpAZM2ZUdBxCCCGEqEMUOhk9aqAsWwkLIYQQ4n/K8hta7QehLl68GG9vbywtLQkKCuLgwYP3rf/jjz/SsmVLLC0tadeuHVu2bDFRpEIIIYQorWqdgKxdu5ZJkyYxY8YMjhw5gp+fHwMHDuT69esl1t+3bx9hYWGMHTuWv//+m6FDhzJ06NBii6YJIYQQompV60cwQUFBdOrUiUWLFgGg1Wrx8vJiwoQJTJkypVj94cOHk52dzaZNm/RlXbp0oUOHDixdurRUbcojGCGEEKJ8yvIbWu5ZMJUtPz+fw4cPM3XqVH2ZUqkkODiY/fv3l3jO/v37mTRpkkHZwIED2bBhwz3bycvLIy8vT//+9oydjIwMI6IXQggh6p7bv52l6duotglISkoKGo0GV1dXg3JXV1diY2NLPCcxMbHE+omJifdsJzw8nJkzZxYr9/LyKkfUQgghhMjMzMTe3v6+daptAmIqU6dONeg10Wq13Lx5k/r161ebhdUyMjLw8vIiPj5eHguZiNxz05N7bnpyz02vtt9znU5HZmYmHh4eD6xbbRMQZ2dnVCoVSUlJBuVJSUm4ubmVeI6bm1uZ6gOo1WrUarVBmYODQ/mCrmR2dna18gtbnck9Nz2556Yn99z0avM9f1DPx23VdhaMhYUFAQEBREVF6cu0Wi1RUVF07dq1xHO6du1qUB9g586d96wvhBBCiKpRbXtAACZNmsSoUaMIDAykc+fOLFiwgOzsbEaPHg3AyJEj8fT0JDw8HICJEyfSu3dv5s+fzyOPPEJERAQxMTF89dVXVfkxhBBCCHGXap2ADB8+nOTkZKZPn05iYiIdOnRg27Zt+oGmcXFxBnvPdOvWjdWrVzNt2jTeffddfHx82LBhA23btq2qj1Ah1Go1M2bMKPaoSFQeueemJ/fc9OSem57c8/+p1uuACCGEEKJ2qrZjQIQQQghRe0kCIoQQQgiTkwRECCGEECYnCYgQQgghTE4SkCq2ZMkS2rdvr1+UpmvXrmzduvW+5/z444+0bNkSS0tL2rVrx5YtW0wUbe1Q1nu+atUqFAqFwcvS0tKEEdc+c+fORaFQ8Prrr9+3nnzXK05p7rl8143zwQcfFLt/LVu2vO85dfk7LglIFWvYsCFz587l8OHDxMTE0K9fPx599FFOnTpVYv19+/YRFhbG2LFj+fvvvxk6dChDhw7l5MmTJo685irrPYeiVQsTEhL0rytXrpgw4trl0KFDLFu2jPbt29+3nnzXK05p7znId91Ybdq0Mbh/f/755z3r1vnvuE5UO46Ojrqvv/66xGNPPfWU7pFHHjEoCwoK0r388sumCK3Wut89/+abb3T29vamDaiWyszM1Pn4+Oh27typ6927t27ixIn3rCvf9YpRlnsu33XjzJgxQ+fn51fq+nX9Oy49INWIRqMhIiKC7Ozsey4fv3//foKDgw3KBg4cyP79+00RYq1TmnsOkJWVRePGjfHy8npgb4m4t1dffZVHHnmk2He4JPJdrxhluecg33VjnTt3Dg8PD5o2bcozzzxDXFzcPevW9e94tV4Jta44ceIEXbt2JTc3FxsbG9avX0/r1q1LrJuYmKhfCfY2V1dXEhMTTRFqrVGWe+7r68vKlStp37496enpzJs3j27dunHq1CkaNmxo4shrroiICI4cOcKhQ4dKVV++68Yr6z2X77pxgoKCWLVqFb6+viQkJDBz5kx69uzJyZMnsbW1LVa/rn/HJQGpBnx9fTl69Cjp6en89NNPjBo1ij179tzzB1EYryz3vGvXrga9I926daNVq1YsW7aMjz76yJRh11jx8fFMnDiRnTt3yqBGEynPPZfvunEGDRqk/7t9+/YEBQXRuHFj1q1bx9ixY6swsupJEpBqwMLCgubNmwMQEBDAoUOH+Pzzz1m2bFmxum5ubiQlJRmUJSUl4ebmZpJYa4uy3PO7mZub4+/vz/nz5ys7zFrj8OHDXL9+nY4dO+rLNBoN0dHRLFq0iLy8PFQqlcE58l03Tnnu+d3ku24cBwcHWrRocc/7V9e/4zIGpBrSarXk5eWVeKxr165ERUUZlO3cufO+4xfEg93vnt9No9Fw4sQJ3N3dKzmq2qN///6cOHGCo0eP6l+BgYE888wzHD16tMQfQvmuG6c89/xu8l03TlZWFhcuXLjn/avz3/GqHgVb102ZMkW3Z88e3aVLl3THjx/XTZkyRadQKHQ7duzQ6XQ63XPPPaebMmWKvv7evXt1ZmZmunnz5unOnDmjmzFjhs7c3Fx34sSJqvoINU5Z7/nMmTN127dv1124cEF3+PBh3YgRI3SWlpa6U6dOVdVHqBXunpEh3/XK96B7Lt9147z55pu63bt36y5duqTbu3evLjg4WOfs7Ky7fv26TqeT7/jd5BFMFbt+/TojR44kISEBe3t72rdvz/bt2xkwYAAAcXFxKJX/66jq1q0bq1evZtq0abz77rv4+PiwYcMG2rZtW1UfocYp6z1PTU3lxRdfJDExEUdHRwICAti3b5+M0alg8l03PfmuV6yrV68SFhbGjRs3aNCgAT169ODAgQM0aNAAkO/43RQ6nU5X1UEIIYQQom6RMSBCCCGEMDlJQIQQQghhcpKACCGEEMLkJAERQgghhMlJAiKEEEIIk5MERAghhBAmJwmIEEIIIUxOEhAhapCVK1eiUCg4d+6cQfnChQtRKBTMmDHDoPzmzZsolUo++eQTU4YJwPPPP1/lCyrt3r2bOXPmFCv/4IMPsLGxKfP1Ll++jEKh0L+uXr1aEWGWytKlS/XtVvV9FaIiSAIiRA3SvXt3APbt22dQvnfvXqytrYuV79u3D51OR48ePUwWY3VyrwTEWHPmzGH//v24uLhU+LXvZdiwYezfv5+HH37YZG0KUZkkARGiBvH19aVBgwbs3bvXoHzv3r08//zzHDhwAI1GY1BuaWlJYGCgqUOt1Xx8fOjSpQsWFhYma9PFxYUuXbrol/UWoqaTBESIGqZ79+4GCUhcXBxXr15l4sSJ5Obmcvz4cf2xvXv3EhgYiIWFBbGxsYwYMQIvLy+sra1p3bo18+fPR6vV6us3adKE1157rVibb731Fg0bNtTXzcvL491336Vx48ao1WpatWrF6tWrHxj71atXefbZZ3F2dsbKyopevXpx+PBhgzre3t689tprLF68mMaNG2Nvb8/QoUNJTk42qHfq1Cl69eqFpaUlPj4+/PDDDwwdOpQ+ffoARY9ZZs6cSXZ2tv7Rxe1jt504cYIePXpgbW1N27Zt2b59+wM/Q0luP5r56aefDMpff/11vL299e9XrVqFQqEgJiaGkJAQrK2t8fX15bfffkOr1TJt2jRcXV1xdXVl6tSpBv9uhKhtJAERoobp3r07Z86cITU1FShKMry8vGjRogV+fn765KSgoIBDhw7pH79cu3YNX19fvvzyS7Zs2cJLL73Ehx9+yEcffaS/9ogRI/jpp58MelF0Oh1r165l+PDh+o20nnrqKZYtW8abb77Jpk2beOihh3j22WfZunXrPeNOTU2lR48eHD16lIULF/Lzzz9Tr149+vXrx/Xr1w3qRkZGEhkZyeLFi/n888/Zs2cPEyZM0B+/desWISEh3Lhxg++//57w8HDmzp1rkMy88MILjB07FisrK/bv38/+/fv58ssv9ccLCgp45plneP7551m/fj0uLi48/vjj3Lhxo8z/Tspq5MiRDB48mPXr1+Ph4cGwYcOYOHEi8fHxfPvtt7z66qvMnTuXiIiISo9FiCpTtZvxCiHKav/+/TpAt3nzZp1Op9O99tpruuHDh+t0Op1uwoQJurCwMJ1Op9MdOHBAB+g2bdpU7BparVZXUFCgmz17ts7d3V1ffuzYMR2g27Fjh75sz549OkB36NAhnU6n0+3atUsH6LZv325wzeHDh+s6deqkfz9q1ChdmzZt9O+nT5+us7e31yUlJenLcnNzdY0aNdK9/fbb+rLGjRvrGjZsqMvNzdWX3d6mXKPR6HQ6nW7x4sU6lUqlu3Tpkr7OpUuXdCqVSte7d2+D8+rVq1fs88+YMcPgHt4+H9B99913xerfXefHH38sVfnEiRN1jRs31r//5ptvdIDuyy+/1JedOHFCB+i6dOlicG5AQIBu6NChxWK4+74KUVNJD4gQNUxAQABWVlb6no69e/fSrVs3ALp27WpQrlAo9Mdyc3OZMWMGzZs3R61WY25uznvvvUdCQgJZWVkAtG/fntatWxv8l3dERAQ+Pj76cSQ7duzAycmJfv36UVhYqH8NGDCAv//+26D35E47duygb9++ODk56c9RqVT07t2bQ4cOGdTt3bs3arVa/75169YUFBToe0oOHTpEu3btDB5veHt74+fnV+r7qFQqCQ4ONjjfysrKJDNbBgwYoP+7RYsWAPTv39+gTosWLYiPj6/0WISoKmZVHYAQomzMzc3p1KkTe/fuJSsri+PHj+uTjG7duunHhOzdu5fWrVvj6OgIwOTJk1m+fDkzZswgICAABwcHNm7cyKxZs8jNzdVPSw0LC2P+/PksWbIEpVLJTz/9xLhx4/Ttp6SkcPPmTczNzUuMLyEhgYYNGxYrT0lJ4cCBAyWe16xZM4P3Dg4OBu9vD/bMzc3Vt1HSYEwXFxdu3bpVYlx3s7KyKjaI1MLCQt9GZbrz892OoaTPbIpYhKgqkoAIUQP16NGDBQsW8Oeff6JWq+nQoQMAjRs3xt3dnb1797Jv3z4effRR/Tk//vgjL7/8MpMnT9aXbd68udi1R4wYwfvvv8+2bdtQq9UkJycTFhamP+7k5ESDBg3YsmVLibHda2qqk5MTDz30kMGYk9vu7O0oDXd3d44ePVqs/Pr169ja2pbpWhWpsLDQ4H12dnYVRSJE9ScJiBA1UI8ePZgzZw6LFi2iU6dOmJn973/K3bp147vvviMxMVG/bggUDdy887/4NRpNiYMcmzdvTqdOnVizZo0+uWnZsqX+eHBwMJ988gkWFha0b9++1DEHBwfz/fff06pVK+rVq1fWj2ygU6dOfPvtt1y6dIkmTZoARTNRjh07ZrDmiYWFBXl5eUa1VRaxsbEG7+9+tCSE+B9JQISogbp27YpSqWTLli1MmTKl2LG3334bwODHeMCAASxfvpzWrVvj7OzMl19+ec8f57CwMN5//33MzMx47733DI4NGDCAIUOG8NBDD/HOO+/Qvn17srOzOXXqFOfPn+frr78u8ZqTJk3ihx9+oHfv3kycOJFGjRqRnJzMX3/9hYeHB2+88UapP//o0aOZPXs2gwcPZubMmUDRtFs3Nzf9TB2AVq1aUVhYyOeff063bt2ws7PD19e31O2U1cKFC/Hw8KBp06asW7eO06dPY25uTmRkJKGhoZXWrhA1kQxCFaIGcnBwoE2bNuh0Ov34j9u6deuGTqfDw8ND3zsART+OvXv3ZsKECYwdO5Z27drx7rvvlnj94cOHc+vWLTIyMhgxYkSx4z/99BOvvPIKX375JYMGDWLs2LHs2LGD3r173zPm+vXrc+DAATp06MDkyZMJCQnhjTfe4PLlywQFBZXp81tZWekHwz7zzDO88847vPnmm/j4+GBvb6+vN2TIEMaPH094eDhBQUG8/PLLZWqnrMaPH88XX3zB4MGDuXLlCr/++iuWlpaEh4dXartC1EQKnU6nq+oghBDCWDdv3qRp06a88cYbxfbEqSiXL1+mSZMmrF27lmHDhukffd0u//HHH3niiScqpW2dTodGo2Hs2LEcPnyYkydPVko7QpiKPIIRQtRIH3/8Ma6urnh7e5OQkMC8efPQaDSMGTOm0tsePnw4APHx8SXO+KkMy5Yt089GatOmjUnaFKIySQIihKiRlEols2bN4tq1a5iZmREUFMSuXbvw8vKqtDY9PDwMBpa6urpWWlt3e/zxx/VrsVhZWZmsXSEqizyCEUIIIYTJySBUIYQQQpicJCBCCCGEMDlJQIQQQghhcpKACCGEEMLkJAERQgghhMlJAiKEEEIIk5MERAghhBAmJwmIEEIIIUxOEhAhhBBCmNz/A28w0ed+Eya7AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x350 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(6, 3.5))\n",
    "\n",
    "# Here we're trimming some edge effects from each detector.\n",
    "# We're also rescaling the flux by ~2% to get a better match to the model spectrum. This is not unexpected \n",
    "# and can be due to e.g., stellar brightness variations.\n",
    "\n",
    "plt.errorbar(wave1[:-5], spec1[:-5]/1.02*1e11, yerr=spec_err1[:-5]*1e11, \n",
    "             c='royalblue', ms=1, fmt='o')\n",
    "\n",
    "plt.errorbar(wave2[7:-20], spec2[7:-20]/1.02*1e11, yerr=spec_err2[7:-20]*1e11, \n",
    "             c='royalblue', ms=1, fmt='o')\n",
    "\n",
    "# And plot up the model spectrum\n",
    "\n",
    "plt.plot(mod_wave, mod_spec*1e11, c='black', lw=0.5)\n",
    "\n",
    "plt.xlim(2.8, 5.3)\n",
    "plt.ylim(0, 1.6)\n",
    "plt.ylabel(r'Flux at Earth [$\\rm 10^{-11} erg/s/cm^2/µm$]', fontsize=11)\n",
    "plt.xlabel('Wavelength [µm]', fontsize=11)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4fc9e6d-abe5-47a1-b3c9-0873fd193c0a",
   "metadata": {},
   "source": [
    "The match is pretty good!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c0c690e-0d00-467e-8fc5-afe8a813da87",
   "metadata": {},
   "source": [
    "And voilà, you now have a flux-calibrated stellar spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ba2b953-4479-43d3-a7e2-9b1b8d1b6b2a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "exotedrf",
   "language": "python",
   "name": "exotedrf"
  },
  "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.10.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}