diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 12215c1f..947bfb8d 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -8,7 +8,7 @@ on: jobs: build: - uses: OpenAstronomy/github-actions-workflows/.github/workflows/publish_pure_python.yml@v1 + uses: OpenAstronomy/github-actions-workflows/.github/workflows/publish_pure_python.yml@924441154cf3053034c6513d5e06c69d262fb9a6 # v1.13.0 with: upload_to_pypi: ${{ (github.event_name == 'release') && (github.event.action == 'released') }} secrets: diff --git a/.github/workflows/ci_workflows.yml b/.github/workflows/ci_workflows.yml index ae073270..4a0aed1b 100644 --- a/.github/workflows/ci_workflows.yml +++ b/.github/workflows/ci_workflows.yml @@ -4,7 +4,7 @@ on: [push, pull_request] jobs: retrieve_cache: - uses: spacetelescope/webbpsf/.github/workflows/retrieve_cache.yml@develop + uses: spacetelescope/webbpsf/.github/workflows/retrieve_cache.yml@beda656c80a0254e6f80649d9c9c49235634522f # v1.4.0 with: minimal: true tests: @@ -49,12 +49,12 @@ jobs: steps: - name: Checkout code - uses: actions/checkout@v4 + uses: actions/checkout@d632683dd7b4114ad314bca15554477dd762a938 # v4.2.0 with: fetch-depth: 0 - name: Set up Python ${{ matrix.python }} - uses: actions/setup-python@v5 + uses: actions/setup-python@f677139bbe7f9c59b41e40162b753c062f5d49a3 # v5.2.0 with: python-version: ${{ matrix.python }} @@ -62,7 +62,7 @@ jobs: run: pip install tox tox-conda>=0.2 - name: Get WebbPSF Data - uses: actions/cache/restore@v4 + uses: actions/cache/restore@0c45773b623bea8c8e75f6c82b208c3cf94ea4f9 # v4.0.2 with: path: ${{ needs.retrieve_cache.outputs.cache_path }} key: ${{ needs.retrieve_cache.outputs.cache_key }} @@ -90,6 +90,6 @@ jobs: - name: Upload coverage to codecov if: ${{ contains(matrix.toxenv,'-cov') }} - uses: codecov/codecov-action@v2 + uses: codecov/codecov-action@e28ff129e5465c2c0dcc6f003fc735cb6ae0c673 # v4.5.0 with: file: ./coverage.xml diff --git a/.github/workflows/download_data.yml b/.github/workflows/download_data.yml index 5392f15d..9802d8ce 100644 --- a/.github/workflows/download_data.yml +++ b/.github/workflows/download_data.yml @@ -18,7 +18,7 @@ # - cron: "0 0 * * 0" # jobs: # download_webbpsf: -# uses: spacetelescope/webbpsf/.github/workflows/download_data.yml@develop +# uses: spacetelescope/webbpsf/.github/workflows/download_data.yml@beda656c80a0254e6f80649d9c9c49235634522f # v1.4.0 # with: # minimal: true @@ -27,8 +27,13 @@ name: download and cache data on: workflow_call: inputs: + url: + description: URL to gzip file + type: string + required: false + default: https://stsci.box.com/shared/static/0dt9z6b927iqgtify2a4cvls9hvapi6k.gz minimal: - description: minimal dataset + description: dataset is minimal (as opposed to full) type: boolean required: false default: true @@ -41,17 +46,25 @@ on: value: ${{ jobs.download.outputs.cache_key }} workflow_dispatch: inputs: + url: + description: URL to gzip file + type: string + required: false + default: https://stsci.box.com/shared/static/0dt9z6b927iqgtify2a4cvls9hvapi6k.gz minimal: - description: minimal dataset + description: dataset is minimal (as opposed to full) type: boolean required: false default: true schedule: - cron: "0 0 * * 0" release: + push: + branches: + - develop env: - DATA_URL: https://stsci.box.com/shared/static/qxpiaxsjwo15ml6m4pkhtk36c9jgj70k.gz + FULL_DATA_URL: https://stsci.box.com/shared/static/qxpiaxsjwo15ml6m4pkhtk36c9jgj70k.gz MINIMAL_DATA_URL: https://stsci.box.com/shared/static/0dt9z6b927iqgtify2a4cvls9hvapi6k.gz jobs: @@ -59,7 +72,7 @@ jobs: name: download and cache WebbPSF data runs-on: ubuntu-latest steps: - - run: wget ${{ (github.event_name == 'schedule' || github.event_name == 'release') && env.MINIMAL_DATA_URL || inputs.minimal && env.MINIMAL_DATA_URL || env.DATA_URL }} -O ${{ runner.temp }}/webbpsf-data.tar.gz + - run: wget ${{ (github.event_name == 'workflow_dispatch' || github.event_name == 'workflow_call') && inputs.url || env.MINIMAL_DATA_URL }} -O ${{ runner.temp }}/webbpsf-data.tar.gz - run: mkdir ${{ runner.temp }}/data/ - run: tar -xzvf ${{ runner.temp }}/webbpsf-data.tar.gz -C ${{ runner.temp }}/data/ - id: cache_path @@ -68,7 +81,7 @@ jobs: run: echo "version=$(cat ${{ steps.cache_path.outputs.cache_path }}/webbpsf-data/version.txt)" >> $GITHUB_OUTPUT - id: cache_key run: echo "cache_key=webbpsf-data-${{ (github.event_name == 'schedule' || github.event_name == 'release') && 'mini' || inputs.minimal && 'mini' || 'full' }}-${{ steps.version.outputs.version }}" >> $GITHUB_OUTPUT - - uses: actions/cache/save@v4 + - uses: actions/cache/save@0c45773b623bea8c8e75f6c82b208c3cf94ea4f9 # v4.0.2 with: path: ${{ runner.temp }}/data/ key: ${{ steps.cache_key.outputs.cache_key }} diff --git a/.github/workflows/retrieve_cache.yml b/.github/workflows/retrieve_cache.yml index aab0fc77..1dc0efee 100644 --- a/.github/workflows/retrieve_cache.yml +++ b/.github/workflows/retrieve_cache.yml @@ -13,7 +13,7 @@ # ... # jobs: # webbpsf_data_cache_key: -# uses: spacetelescope/webbpsf/.github/workflows/retrieve_cache.yml@develop +# uses: spacetelescope/webbpsf/.github/workflows/retrieve_cache.yml@beda656c80a0254e6f80649d9c9c49235634522f # v1.4.0 # with: # minimal: true # tests: @@ -21,7 +21,7 @@ # steps: # ... # - name: retrieve WebbPSF data cache -# uses: actions/cache/restore@v4 +# uses: actions/cache/restore@0c45773b623bea8c8e75f6c82b208c3cf94ea4f9 # v4.0.2 # with: # path: ${{ runner.temp }}/webbpsf-data # key: ${{ needs.webbpsf_data_cache_key.outputs.cache_key }} diff --git a/CITATION.cff b/CITATION.cff index 52d12674..4aa58f20 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -56,4 +56,4 @@ abstract: >- including its Wide Field Instrument and a preliminary version of the Coronagraph Instrument. license: BSD-3-Clause -version: 1.4.0 +version: 1.5.0 diff --git a/README.rst b/README.rst index 41881436..fb08f447 100644 --- a/README.rst +++ b/README.rst @@ -21,6 +21,14 @@ WebbPSF: Simulated Point Spread Functions for the James Webb and Nancy Grace Rom .. image:: https://img.shields.io/badge/ascl-1504.007-blue.svg?colorB=262255 :target: http://ascl.net/1504.007 + +**ADVISORY: WebbPSF IS BEING MIGRATED TO A NEW REPOSITORY: STPSF (Space Telescope PSF)** + + **To reflect its broader support for Roman as well as James Webb, WebbPSF is being migrated to a new repository: STPSF (Space Telescope PSF).** + **This transition is being done in such a way as to maintain back-compatibility for existing code, and existing installations will continue to run as-is.** + **This transitional period is ongoing now. Please do not submit pull requests to this webbpsf repo at this time.** + **The STPSF github repository will soon be available for use instead.** + WebbPSF produces simulated PSFs for the James Webb Space Telescope, NASA's flagship infrared space telescope. WebbPSF can simulate images for any of the four science instruments plus the fine guidance sensor, including both direct @@ -29,8 +37,12 @@ imaging, coronagraphic, and spectroscopic modes. WebbPSF also supports simulating PSFs for the upcoming Nancy Grace Roman Space Telescope (formerly WFIRST), including its Wide Field Instrument and a preliminary version of the Coronagraph Instrument. -Developed by Marshall Perrin, Joseph Long, Shannon Osborne, Robel Geda, Bradley Sappington, Marcio Meléndez, -Charles-Phillipe Lajoie, Jarron Leisenring, Neil Zimmerman, Keira Brooks, +.. note:: + + The current Roman WFI optical model was provided by Goddard Space Flight Center circa 2021 (the Cycle 9 reference data); a new optical model is currently being implemented in WebbPSF. + +Developed by Marshall Perrin, Joseph Long, Shannon Osborne, Robel Geda, Bradley Sappington, Marcio Meléndez, +Charles-Philippe Lajoie, Jarron Leisenring, Neil Zimmerman, Keira Brooks, Justin Otor, Trey Kulp, Lauren Chambers, Alden Jurling, and collaborators, 2010-2024. Documentation can be found online at https://webbpsf.readthedocs.io diff --git a/docs/index.rst b/docs/index.rst index 8795cea0..1f3086bd 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -61,6 +61,7 @@ Contents jwst_measured_opds.ipynb jwst_detector_effects.ipynb jwst_matching_psfs_to_data.ipynb + jwst_ifu_datacubes.ipynb jwst_large_psf.ipynb jwst_optical_budgets.ipynb jwst_psf_subtraction.ipynb diff --git a/docs/jwst_ifu_datacubes.ipynb b/docs/jwst_ifu_datacubes.ipynb new file mode 100644 index 00000000..cfb01d95 --- /dev/null +++ b/docs/jwst_ifu_datacubes.ipynb @@ -0,0 +1,616 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "74185321-c269-45be-96c5-6e39e8d317ee", + "metadata": {}, + "source": [ + "# Simulating IFU mode and Datacubes" + ] + }, + { + "cell_type": "markdown", + "id": "54d3884f-b190-43b3-9dff-f169c4833678", + "metadata": {}, + "source": [ + "This page describes how to simulate IFU mode data and spectral datacubes for JWST's NIRSpec IFU and MIRI MRS.\n", + "\n", + "
\n", + "IFU support is a relatively recent addition, and efforts are still in progress to refine and improve the fidelity of IFU mode simulations. \n", + "
\n", + "\n", + "\n", + "## Selecting IFU mode simulations\n", + "\n", + "These instruments have a `mode` attribute, which can be set to either 'IFU' or 'imaging'. Selecting IFU mode has the following effects:\n", + "\n", + " - The normal `filter` attribute for selecting spectral bandpass is superceded by a `band` attribute for selected IFU bands, the specific details of which differ for NIRSpec and MIRI. A `get_IFU_wavelengths()` function is added, which allows looking up the wavelength range for each band. \n", + " - The PSF simulation gets an added step for adding \"IFU broadening\" effects, which are empirical models for slightly broadening/blurring the simulated PSF to better match the observed PSF FWHMs. Physically this is a simplified model for optical blurring effects due to imperfect wavefront quality in the IFU image slicer optics, for example.\n", + " - For NIRSpec IFU simulations only, the PSF output is rotated by an additional 90 degrees to match the orientation of JWST pipeline output dataproducts created with the \"orient='IFUalign'\" option in the Cube Build step." + ] + }, + { + "cell_type": "markdown", + "id": "4eac3034-b525-4576-b3c8-ddfd02d7a7b7", + "metadata": {}, + "source": [ + "Note there are *three options* for computing PSFs in IFU mode. \n", + "\n", + "1. If you only need one wavelength, use regular `calc_psf()` with the `monochromatic` keyword to specify a wavelength.\n", + "2. If you want a datacube at all wavelengths, use `calc_datacube()` with a list or array of wavelengths. This is the recommended path for many typical use cases. \n", + "3. If you want a datacube at all wavelengths, you can also use `calc_datacube_fast()` which achieves a much-faster calculation runtime by making a simplifying assumption in the optical calculation, and currently by not including the detector effects or distortion modeling steps.\n", + " \n", + " * Specifically, it assumes that the wavefront optical path difference in the IFU exit pupil is independent of wavelength. This assumption is reasonably true for both MIRI and NIRSpec IFU modes within the current level of fidelity." + ] + }, + { + "cell_type": "markdown", + "id": "ef5ecc18-2412-4ca3-bb1f-7f1c8ee8f9e6", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## Datacube Coordinate Frames for IFU PSFs: SkyAlign vs IFUAlign\n", + "\n", + "
\n", + "WebbPSF is intended to simulate IFU PSFs in IFUAlign orientation only.\n", + "
\n", + "\n", + "For IFU observations, the JWST pipeline ``cube_build`` step (see [pipeline docs](https://jwst-pipeline.readthedocs.io/en/latest/jwst/cube_build/main.html)) can produce output datacubes in two different reference frames or orientations: \n", + "\n", + " - ``coord_system='skyalign'``. In these cubes, north is up and east is left. This is the default orientation.\n", + " - ``coord_system='ifualign'``. These cubes are build on the local IFU slicer plane's native coordinate system. The X and Y coordinates corresponds to the along-slice and across-slice directions, respectively (which are often labeled as the $\\alpha$ and $\\beta$ axes for MIRI MRS; see [JDox](https://jwst-docs.stsci.edu/jwst-mid-infrared-instrument/miri-instrumentation/miri-mrs-field-and-coordinates#MIRIMRSFieldandCoordinates-CoordinatesMRScoordinateframes).)\n", + "\n", + "PSF simulations are always produced in the ifualign orientation only. For PSF modeling work, we recommend reducing your science data using the ifualign frame. This allows PSF models and orientations to be consistent between multiple datasets, independent of position angles on the sky. \n", + "\n", + "If your science data has been reduced in skyalign frame, you can rotate the PSF based on the aperture position angle to match your data. However this will introduce numerical interpolation noise, particularly for spatially undersampled IFU data. That can be avoided by working in IFUalign frame, as we recommend here. \n", + "\n", + "\n", + "**Pixelscales:** By default, IFU PSFs are made using the same default pixel scales as the pipeline defaults (0.1 arcsec/pixel for NIRSpec, and variously 0.13 - 0.35 arcsec/pixel for MIRI depending on which MRS channel). If you have reduced your science data using a different pixelscale in the JWST pipeline, e.g. 0.123 arcsec/pix, simply set e.g. `miri.pixelscale = 0.123` to the same value for creating the PSF model. " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "7ce1366f-dcc8-47b7-9e73-f0ebbee3f1ba", + "metadata": { + "editable": true, + "nbsphinx": "hidden", + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "**WARNING**: LOCAL JWST PRD VERSION PRDOPSSOC-065 DOESN'T MATCH THE CURRENT ONLINE VERSION PRDOPSSOC-067\n", + "Please consider updating pysiaf, e.g. pip install --upgrade pysiaf or conda update pysiaf\n" + ] + } + ], + "source": [ + "%matplotlib inline\n", + "import webbpsf\n", + "import astropy.units as u" + ] + }, + { + "cell_type": "markdown", + "id": "e25e70a8-7a8b-456f-9719-e52951207bac", + "metadata": {}, + "source": [ + "## NIRSpec IFU example\n", + "\n", + "For NIRSpec, the `band` attribute is derived from the `prism` and `disperser` elements. " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "b3ec32d8-5ede-46fb-8f51-da8c06039fbe", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Band is PRISM/CLEAR\n" + ] + } + ], + "source": [ + "nrs = webbpsf.NIRSpec()\n", + "nrs.mode = 'IFU'\n", + "\n", + "nrs.disperser = 'PRISM'\n", + "nrs.filter = 'CLEAR'\n", + "print(\"Band is\", nrs.band)" + ] + }, + { + "cell_type": "markdown", + "id": "d950963b-b4df-42ff-8729-04ad44e17889", + "metadata": {}, + "source": [ + "The wavelength sampling can be obtained using the `get_IFU_wavelengths()` function. By default this returns the same wavelength sampling as the pipeline uses. But if desired you can also specify some other number of wavelengths `nlambda`, for instance to reduce simulation runtimes when the full spectral resolution is not needed. " + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "e68eb5e3-229f-4bcf-a8a0-5133788ba767", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "PRISM/CLEAR default wavelength sampling uses 940 wavelengths\n" + ] + }, + { + "data": { + "text/latex": [ + "$[0.6,~0.605,~0.61,~\\dots,~5.285,~5.29,~5.295] \\; \\mathrm{\\mu m}$" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "allwaves = nrs.get_IFU_wavelengths()\n", + "print(f\"{nrs.band} default wavelength sampling uses {len(allwaves)} wavelengths\")\n", + "allwaves" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "7417534b-6550-4830-9f23-a590aad7357f", + "metadata": {}, + "outputs": [], + "source": [ + "# let's get a subset for a faster PSF sim runtime\n", + "waves = nrs.get_IFU_wavelengths(nlambda=50)\n", + "\n", + "# convert waves from microns to meters\n", + "# (this is a work around for a current issue with the poppy library upstream)\n", + "waves = waves.to_value(u.meter)" + ] + }, + { + "cell_type": "markdown", + "id": "7f5180f4-2890-47db-aea7-b89c6d176130", + "metadata": {}, + "source": [ + "Compute a datacube:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "212dcbf1-4fef-43a3-8675-a7dbad84ace1", + "metadata": {}, + "outputs": [], + "source": [ + "cube = nrs.calc_datacube(waves)" + ] + }, + { + "cell_type": "markdown", + "id": "8662db59-3a15-4a6d-b3de-341a4189a05a", + "metadata": {}, + "source": [ + "The output FITS file has the same extensions as a regular PSF calculation, but each extension contains a 3D datacube rather than a 2D image: " + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "58429d9f-621b-4bcd-96c1-917abf609318", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Filename: (No file associated with this HDUList)\n", + "No. Name Ver Type Cards Dimensions Format\n", + " 0 OVERSAMP 1 PrimaryHDU 1288 (192, 192, 50) float64 \n", + " 1 DET_SAMP 1 ImageHDU 1290 (48, 48, 50) float64 \n", + " 2 OVERDIST 1 ImageHDU 1343 (192, 192, 50) float64 \n", + " 3 DET_DIST 1 ImageHDU 1344 (48, 48, 50) float64 \n" + ] + } + ], + "source": [ + "cube.info()" + ] + }, + { + "cell_type": "markdown", + "id": "dd3b4280-6ec4-4f62-b6ea-8108780dfe20", + "metadata": {}, + "source": [ + "The `calc_datacube_fast` routine does a simplified calculation, much faster: " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "228264f6-0f16-489d-99a7-4dadbea39bc7", + "metadata": {}, + "outputs": [], + "source": [ + "quickcube = nrs.calc_datacube_fast(waves)" + ] + }, + { + "cell_type": "markdown", + "id": "7b12112c-d610-4579-8072-f644c01d4a83", + "metadata": {}, + "source": [ + "Note that in this case, the output FITS file only contains the first oversampled extension:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "de797cb5-3b4f-4a4a-8152-59e68f7bc78f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Filename: (No file associated with this HDUList)\n", + "No. Name Ver Type Cards Dimensions Format\n", + " 0 OVERSAMP 1 PrimaryHDU 159 (192, 192, 50) float64 \n" + ] + } + ], + "source": [ + "quickcube.info()" + ] + }, + { + "cell_type": "markdown", + "id": "e550c2b5-20fa-493d-ae39-317f7c5a4cdf", + "metadata": {}, + "source": [ + "The display_psf function works with datacubes, but you have to specify which slice of the cube to display. " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "2ffd3c45-4d4c-4a0f-9a41-6d28e7618eac", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh0AAAHHCAYAAAAbLeozAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABwb0lEQVR4nO3deVzU1f4/8NewgyD7IiqKoibhFpRbptbNNVMz0zQ1u9nPxFzIFivXa5nllmWaZnq7mdpNLfteMy1N09AbKJo3KzONFBBRAUVlm8/vD2J0hDlvmI1hfD0fDx4PmDOf89mHM+dz3u+j0zRNAxEREZGNudT0BhAREdGtgY0OIiIisgs2OoiIiMgu2OggIiIiu2Cjg4iIiOyCjQ4iIiKyCzY6iIiIyC7Y6CAiIiK7YKODiIiI7IKNDiIiIrILNjqIiIjILm75RseaNWug0+ng5eWFP/74o0J5t27dEBcXZ/Te8h83NzfUq1cPQ4cOxfHjxyut/8CBAxg4cCCioqLg6emJ8PBwdOzYEc8++2yVty0lJaXCa5X9TJkyBQAwc+ZM6HQ65OTkVFpvXFwcunXrJq7f2qTtspfyY3jq1Cnla7a2c+dOPPHEE7jttttQp04d1K9fH/3790dqamql7798+TImTZqEyMhIeHl5oW3btli/fr1F23Dt2jW4ubnBz88PzzzzjEV13ay6+3ezb7/91uS1vn//fsP7Ll26hOeffx49evRAaGgodDodZs6cabLeQ4cOYcCAAYiMjISPjw9uu+02zJ49G1euXLF0l0U1cZ0RORK3mt4AR1FYWIhXXnkF//rXv8T3rl69GrfddhuuXbuGffv24dVXX8WuXbvw888/IzAw0PC+//znP3jwwQfRrVs3vPHGG6hXrx4yMzORkpKC9evXY8GCBWZvb/k23CgyMtLs+m51ffv2RXJyMurVq2e3dS5btgznz5/HxIkTERsbi3PnzmHBggXo0KEDvvrqK9x7771G73/ooYfwww8/4PXXX0fz5s3x8ccf49FHH4Ver8ewYcPM2gadToedO3di1qxZeOedd/DMM8+gefPm1ti9au+fKa+99hq6d+9u9Fr5FwEAOH/+PFasWIE2bdpgwIABeP/9903W9dNPP6FTp05o0aIFFi9ejJCQEOzZswezZ89GamoqPv/8c/N2topq4jojcijaLW716tUaAK1Xr16ai4uLlpaWZlTetWtX7fbbbzd67w8//GD0nlmzZmkAtA8++MDo9XvuuUdr2rSpVlxcXGG9paWlVd62G9dnahtuNGPGDA2Adu7cuUrLb7/9dq1r167i+q1N2i57KT+GJ0+erNHtOHv2bIXXLl26pIWHh2v33Xef0ev/+c9/NADaxx9/bPT6/fffr0VGRmolJSUWbcu2bdsqrd8S1dm/yuzatUsDoP373/9Wvk+v12t6vV7TNE07d+6cBkCbMWNGpe99+eWXNQDab7/9ZvT6U089pQHQLly4IG6XoygoKKjpTSCqtlv+8Uq5559/HsHBwXjhhReqvWxCQgIA4OzZs0avnz9/HiEhIXBzq9ih5OJSOw79zz//jEcffRTh4eHw9PREVFQURo4cicLCQgDA448/jsaNG1dYrvxRSmX+/PNPPPTQQ6hbty78/f3x2GOP4dy5c0bvOX78OIYNG4awsDB4enqiZcuWWLp0aZW2+dy5c3jqqafQsGFDeHp6IjQ0FJ07d8bXX39tchlT3d7S/luyrWFhYRVe8/X1RWxsLP7880+j1zdv3gxfX18MHjzY6PXRo0cjIyMDBw4cENenUt5rlpaWZlE9N6rO/lmi/JFLVbi7uwMA/P39jV4PCAiAi4sLPDw8lMuXX9dHjhzB4MGD4e/vj6CgICQlJaGkpAS//PILevXqBT8/PzRu3BhvvPGG0fLmXmfl6z148CAefvhhBAYGomnTpobl9+7di/vuuw9+fn7w8fFBp06d8J///KfSbf/f//6HRx99FP7+/ggPD8cTTzyBvLy8Kh0/IkvVjv98duDn54dXXnkFX331FXbu3FmtZU+ePAkAFbqlO3bsiAMHDmDChAk4cOAAiouLrba9paWlKCkpMfqxtsOHD+POO+/E/v37MXv2bHz55ZeYO3cuCgsLUVRUZHa9AwcORExMDD799FPMnDkTn332GXr27Gk4Pj/99BPuvPNOHD16FAsWLMD//d//oW/fvpgwYQJmzZol1j9ixAh89tlnmD59OrZv3473338ff/vb33D+/PlqbWdV9t/Sbb1ZXl4eDh48iNtvv93o9aNHj6Jly5YVGrCtW7c2lFti3rx5ACpvdGiaVuFaM/UjMbV/KomJiXBzc0PdunXRs2dP7N27t8rL3mzUqFEICAjA008/jd9//x2XLl3C//3f/+G9995DYmIi6tSpU6V6HnnkEbRp0wYbN27EmDFjsGjRIkyePBkDBgxA3759sXnzZtx777144YUXsGnTJmVd1bnPHnroIcTExODf//43li9fDgDYvXs37r33XuTl5WHVqlVYt24d/Pz80K9fP2zYsKHC+gYNGoTmzZtj48aNePHFF/Hxxx9j8uTJVTyCRBaq6a6Wmnbj44rCwkKtSZMmWkJCgqG7trLHK/v379eKi4u1S5cuadu2bdMiIiK0e+65p8JjlJycHO3uu+/WAGgANHd3d61Tp07a3LlztUuXLlVr225+rbKf8vVb6/HKvffeqwUEBGjZ2dkm3zNq1CitUaNGFV4v34bKXps8ebLR62vXrtUAaB999JGmaZrWs2dPrUGDBlpeXp7R+8aPH695eXmJXeC+vr7apEmTTJZX9nilsteqsv+WbuvNhg8frrm5uWkpKSlGrzdr1kzr2bNnhfdnZGRoALTXXnutWuu50VdffaUB0AIDA7XQ0NAK5eWPOaryIz2yMrV/lTl48KA2ceJEbfPmzdqePXu0Dz74QGvZsqXm6uqqbdu2rdJlpMcrmqZpx44d02677Taj7Z4wYYLhnlcpv4YXLFhg9Hrbtm01ANqmTZsMrxUXF2uhoaHaQw89ZHjN3OusfL3Tp0+vUNahQwctLCzM6DOlpKREi4uL0xo0aGDYr/I63njjDaPlx40bp3l5eVVp/4ksxZ6OG3h4eGDOnDlISUnBJ598YvJ9HTp0gLu7O/z8/NCrVy8EBgbi888/r/AtNDg4GN99951h8F///v3x66+/YurUqWjVqpVFURwffvghfvjhB6Ofyh7jmOvKlSvYvXs3HnnkEYSGhlqtXgAYPny40d+PPPII3NzcsGvXLly7dg3ffPMNBg4cCB8fH6Nv0X369MG1a9eMIhcqc9ddd2HNmjWYM2cO9u/fb1YPU1X23xrbeqNp06Zh7dq1WLRoEeLj4yuUqx4hVPXxws1yc3PxxBNPoH///hg3bhzOnTuHjIwMo/fEx8dXuNZM/agGM0v7d7N27dph8eLFGDBgALp06YLRo0fj+++/R7169fD888+btb+nTp1Cv379EBwcjE8//RS7d+/GG2+8gTVr1uDJJ5+scj0PPPCA0d8tW7aETqdD7969Da+5ubkhJiam0qi4ctW9zwYNGmT0d0FBAQ4cOICHH34Yvr6+htddXV0xYsQInD59Gr/88ovRMg8++KDR361bt8a1a9eQnZ0trp/IUoxeucnQoUMxf/58vPzyy3jooYcqfc+HH36Ili1b4tKlS9iwYQPee+89PProo/jyyy8rfX9CQoJh3EdxcTFeeOEFLFq0CG+88UaFZ75V1bJlS0OdNytvfJSWllZaXlJSYni2bcrFixdRWlqKBg0amLV9KhEREUZ/u7m5ITg4GOfPn8f58+dRUlKCt99+G2+//Xaly0uNtQ0bNmDOnDl4//33MW3aNPj6+mLgwIF44403KqzblKrsvzW2tdysWbMwZ84cvPrqqxg/fnyF8vLjc7MLFy4AAIKCgqq0npslJiaiuLgYK1euxJ49ewCUPWK5sfHg6+uLtm3bVqk+Uw1faf+qKiAgAA888ACWL1+Oq1evwtvbu1rLv/jii8jPz0daWprhUco999yDkJAQPPHEExg5ciS6du0q1nPz8fbw8ICPjw+8vLwqvJ6fn2+ynureZzdHvVy8eBGaplUaDVN+Dm++boKDg43+9vT0BABcvXq1SttAZAk2Om6i0+kwb9483H///VixYkWl77nxH3737t1RWlqK999/H59++ikefvhhZf3u7u6YMWMGFi1aZPFzeFPCw8MBAGfOnDH8Xk7TNGRmZppssJQLCgqCq6srTp8+rXyfl5eX0aDKcqp/tllZWahfv77h75KSEpw/fx7BwcEIDAw0fEtLTEysdPno6GjlNoWEhGDx4sVYvHgx0tPTsWXLFrz44ovIzs7Gtm3blMuWq8r+W2NbgbJ/yDNnzsTMmTPx0ksvVfqeVq1aYd26dSgpKTH6x/7jjz8CMA4hrapPP/0UH3/8Mb744guEhobijjvuAFDW6OjTp4/hfbt3764QsmrKyZMnKwwsrsr+VYemaQDM691JS0tDbGxshbEbd955J4CysTFVaXRYS1Xvs3I373NgYCBcXFyQmZlZ4b3lPVYhISGWbyiRlbDRUYm//e1vuP/++zF79mw0bNhQfP8bb7yBjRs3Yvr06XjooYcMkSmZmZmVfgM5duwYANvl1bj33nuh0+mwYcMGwz+Sctu2bUN+fj7+9re/Kevw9vZG165d8e9//xuvvvqqyQ+uxo0bIzs7G2fPnjU0cIqKivDVV1+ZrHvt2rVG3euffPIJSkpK0K1bN/j4+KB79+44dOgQWrduLUYTSKKiojB+/Hh888032LdvX5WXq8r+W2Nb//GPf2DmzJl45ZVXMGPGDJPvGzhwIFauXImNGzdiyJAhhtf/+c9/IjIyEu3bt6/WerOysjB27Fg89dRThkcF0dHRCAgIwKFDh4zeW/54pSpuvqarun9VdfHiRfzf//0f2rZtW6FXoarbd/ToUVy+fNnocURycjIA2KRnT6Wq95kpderUQfv27bFp0ybMnz/f0POj1+vx0UcfoUGDBlbLu0JkDWx0mDBv3jzEx8cjOztbHGkfGBiIqVOn4vnnn8fHH3+Mxx57DADQs2dPNGjQAP369cNtt90GvV6PtLQ0LFiwAL6+vpg4caJNtr1p06YYP3483nzzTeTm5qJPnz7w9vY2jC1JSEioUjKphQsX4u6770b79u3x4osvIiYmBmfPnsWWLVvw3nvvwc/PD0OGDMH06dMxdOhQPPfcc7h27RqWLFli8tEOAGzatAlubm64//778b///Q/Tpk1DmzZt8MgjjwAA3nrrLdx9993o0qULnn76aTRu3BiXLl3Cb7/9hi+++EIZXZSXl4fu3btj2LBhuO222+Dn54cffvgB27ZtM/m4zJL9t2RbFyxYgOnTp6NXr17o27dvhfEfHTp0MPzeu3dv3H///Xj66aeRn5+PmJgYrFu3Dtu2bcNHH30EV1dXo2V1Oh26du2Kb7/9ttJ1jxkzBoGBgVi4cKHR6+3atasQweLn5yf2jFm6f7t378Z9992H6dOnY/r06QCAYcOGISoqCgkJCQgJCcHx48exYMECnD17FmvWrDGq68svv0RBQQEuXboEoCyq6NNPPwUA9OnTBz4+PgCASZMmYcCAAbj//vsxefJkhISEYP/+/Zg7dy5iY2ONxmTYS1WuM5W5c+fi/vvvR/fu3TFlyhR4eHjg3XffxdGjR7Fu3Tqzx/sQ2UQND2StcapkW8OGDdMAiMnBNE3Trl69qkVFRWnNmjUzJGrasGGDNmzYMK1Zs2aar6+v5u7urkVFRWkjRozQfvrpJ7O2rSrJwTStLGHSsmXLtISEBM3Hx0fz8PDQmjVrpr3wwgtVipwp99NPP2mDBw/WgoODNQ8PDy0qKkp7/PHHtWvXrhnes3XrVq1t27aat7e31qRJE+2dd95RRq+kpqZq/fr103x9fTU/Pz/t0UcfrZBI6uTJk9oTTzyh1a9fX3N3d9dCQ0O1Tp06aXPmzFFu77Vr17SxY8dqrVu31urWrat5e3trLVq00GbMmGFIplTV6JWq7r+529q1a1dlJMjNLl26pE2YMEGLiIjQPDw8tNatW2vr1q2r9H0AtKFDh1a63pUrV2qurq5acnJyhbKkpCRNp9Np+fn5ym2viursX3mEzI1RJ3PnztXatm2r+fv7a66urlpoaKg2cOBA7b///W+FdTVq1KjKETU7d+7UevTooUVERGje3t5a8+bNtWeffVbLyckR98lUZNioUaO0OnXqVHoMyj8/NM3860yKSPvuu++0e++9V6tTp47m7e2tdejQQfviiy+qtO2OkiyPbg06TfvrASkROYWtW7figQcewOHDh9GqVaua3hwiIgOGzBI5mV27dmHo0KFscBCRw2FPBxEREdkFezqIiIjILtjoICIissC7776L6OhoeHl5IT4+Ht99953J92ZmZmLYsGFo0aIFXFxcMGnSJPttqANgo4OIiMhMGzZswKRJk/Dyyy/j0KFD6NKlC3r37o309PRK319YWIjQ0FC8/PLLaNOmjZ23tuY57JgOvV5vyGrp4+PDWHMiolpG0zRcuXIFQFlm1PLEibZah7VU539O+/btcccdd2DZsmWG11q2bIkBAwZg7ty5ymW7deuGtm3bYvHixZZsbq3isMnBcnJyKqTwJiKi2uns2bMICwuzer1Xrlwxyi5rDRkZGUap8j09PQ1z1NyoqKgIqampePHFF41e79GjB77//nurbpOz4OMVIiKiG0RGRsLf39/wY6rHIicnB6WlpRW+IIeHhyMrK8sem1rrOGxPR3naYgDoAsDVxPtMJ9sGYoR1SP0o6nlYgYuKMtOTWZcxtT/lzJsztMwFoVw1gbU0m4U0E430faNIUaY6l1WhykohTRoufTz8JpSr9kudxBooEMp/V5RJ3xqkdUud0rmKMukaviSUFyvKpP2SngmrPth8FGUAoLdw3ZcVZdK9KX0mqc5nxVmejEmfZ9J1qDpupmasKQawpHx5H+nIW+7s2bMVJvOrqoKCAkPjobKeDpWbH8VomsYhASY4bKPjxhPmCvkDrjLSTSZNzyUtryqXDqy0P9K6VSxZt6XbJZWrPrAt7XZTfSxIjSn1R4pl+yVdZ6oGC6A+n9Ixk64FW14rUrmqkWlpo8OS7Zb+VUjrVm27tF+WHFPpXErXsCXXSlWmOrTHP+E63p6o4y3dzSboSwy/1q1bt0qNl5CQELi6ulbo1cjOzubwABP4eIWIiJyDvsSyn7/Ex8cjNjYWS5cuVa7Ow8MD8fHx2LFjh9HrO3bsQKdOnWyyi7Wdw/Z0EBERVctNjYdqL/uX1NTUKj+mSUpKwogRI5CQkICOHTtixYoVSE9Px9ixYwEAU6dOxZkzZ/Dhhx8alimfyfny5cs4d+4c0tLS4OHhgdjYWPO2vRZho4OIiMhMQ4YMwfnz5zF79mxkZmYiLi4OW7duRaNGjQCUJQO7OWdHu3btDL+npqbi448/RqNGjXDq1Cl7bnqNYKODiIicg5V6OuLj4+Hi4oLExEQkJiaKi44bNw7jxo2rtGzNmjUVXnPQ9Fh2wUYHERE5hxp4vELVUysaHSWQR4xXRgpPk8IJpbBVVYiadGBPCeWqsNaqjBRXUY2pliIpvIVyaXy6JaHAZ4RyVQizFDLrL5RLy+cpyqTrTDpmqm2Toh2ka0X6eFbVL41CDxDKVdea6ngCcmi2KlJDOmaq6wiQI6FUxyVCWDZAKA9UlEnXmXSuVSHMgDoU+LyZddKtp1Y0OoiIiET6Ugt6OizNFERVwZBZIiJyDnYOmaXqY08HERHRDTimw3bY6CAiIudgpYGkZDtsdBARkXNgo8PhcUwHERHRDTimw3ZqRU/HBZhuHQUolrsm1KsKAQOAukK5KmxPCo2TwvZOKsqk0Dhpv1Thx1JIrLTd0gWlCuGUQkelY6p6Aittl7TfqlBFQB3+Kc22KoVoqr4ZSCGx0n5JIdKqYyqdL0vWLW2XdD5Ux0UVjl4V0iy0loTrSvETqplgpTlcpXVL30BVx9RUWLd0Hq2OeTocXq1odBAREYn4eMXhsdFBRETOQbMgT4fGPB32wDEdREREZBdsdBARkXNgcjCHx8crRETkHDiQ1OGxp4OIiIjsgj0dRETkHBi94vBqRaPDD3KMeWWkGPFcoVxap6eiTJp6u55QrnJBKFfl4ZBI47dNTWFdTsohosrdIB1vKQ+BKn9ClrCspflH6ivKpDwc+UK56pxcEZaVOoil61R1rQULy0pTwJ9RlEnnQ8qVoSqX/rVI65aOqepzQbo3pdwnqmMu5bEpFMqDhHLVtWIqN0ltzdNBtsPHK0RERGQXbHQQEZFzYPSKw6sVj1eIiIhEeguSg+mvP8hk9IrtsKeDiIiI7II9HURE5Bw4kNThsdFBRETOgY0Oh8dGBxEROQc2OhxerWh0BMD0hqpi4qWcE1JMfJ5QflZRpspHAcix/qq8D5KrQrnqpEvHTIq7l27bXEWZqVj/ctIxCVGUSeda2m9pv1R5IVTXKADUFcpVeTykfBXFQrk0qEuVa0M6X9IxVW2bdL6k/CSqfBhSLhnpQ1HKsaPKZyFtt5Q3JUBRZsnxBuQ8OIGKMlP5e+yep4McHgeSEhGRc2DIrMOzaU/H3LlzsWnTJvz888/w9vZGp06dMG/ePLRo0cKWqyUiolsRJ3xzeDbt6di9ezcSExOxf/9+7NixAyUlJejRowcKCgpsuVoiIiJyQDbt6di2bZvR36tXr0ZYWBhSU1Nxzz332HLVRER0q9EsSA6mSaNiyBrsOpA0L69saGZQUOVTCxUWFqKwsGxaIvaGEBFRtehLAL0504OC0St2YreBpJqmISkpCXfffTfi4uIqfc/cuXPh7+8Pf39/REZG2mvTiIiIyA7s1ugYP348jhw5gnXr1pl8z9SpU5GXl4e8vDxkZGTYa9OIiMgZWCl6hWzHLo9XnnnmGWzZsgV79uxBgwYNTL7P09MTnp5lWQ1cXa93kRXDdD4CVT4MaeyxFLfuIZSbik0H5NwLUt2qeH0pX4XUksxWlFn6VPOiUK7ab+mYhArlqk5VqcPVzA7ZKi0v5SqwZW4GKd9FoVAu5QGxpO4IRZmUz0Kiugek81H5w9/rgoVyVY4QibRu1bVyTlhWujelz0OVyzao0yx8vOLwbNro0DQNzzzzDDZv3oxvv/0W0dHRtlwdEREROTCbNjoSExPx8ccf4/PPP4efnx+ysrIAAP7+/vD2lnJ2EhERVQN7OhyeTRsdy5YtAwB069bN6PXVq1fj8ccft+WqiYjoVqO3IGRWz5BZe7D54xUiIiK70JcAejPjI25Kg+7i4oLExEQkJiZaaeMIqCUTvhEREdkL06DbDhsdRETkHKzU00G2UysaHXqYDv2TwsBUpHBCaYp41XClPyxYFgAaKsr8hWWlk6q6JdOFZaXbUpruXAoJVFFNsw6oj6l0rqVhzdLHmOpaka4jaUrxRoqyXGFZS0NPVdOZS+cyVyhX5RyWHsxK30FV16EUyiuFxErXoWq/TIWWlpPCxi0J7ZbCn68J5arjlmvi9ZoJmWWjw5FxansiIiKyi1rR00FERCRiT4fDY6ODiIicAxsdDo+PV4iIiMgu2NNBRETOQbMgOZjG5GD2wEYHERE5B30JoJdi1RTLks3x8QoRERHZRa3v6bikKJM6y6Q8A1JMvKrFJuUCOGPBuqW5ekMsqFua9jtTKPcUylU5RqQWsJRHQJWzpa6wrHQtSDdKlqLsvLCsn1Cuyjkh5fiQrkOpXFV/vrCsdL4uKMqke0/KZ6E6ZtI1KuWzkHKIqL4vS9eRdJ2q9kva7iihXOofUK3b1OewJXlFzMKeDodX6xsdREREANjoqAXY6CAiIufARofD45gOIiIiC7z77ruIjo6Gl5cX4uPj8d133ynfv3v3bsTHx8PLywtNmjTB8uXLjcr/97//YdCgQWjcuDF0Oh0WL15sw623LzY6iIjIOehLLPsxw4YNGzBp0iS8/PLLOHToELp06YLevXsjPb3ymaxOnjyJPn36oEuXLjh06BBeeuklTJgwARs3bjS858qVK2jSpAlef/11REREmLVdjoqPV4iIyDnoSy14vGJeno6FCxfi73//O5588kkAwOLFi/HVV19h2bJlmDt3boX3L1++HFFRUYbei5YtWyIlJQXz58/HoEGDAAB33nkn7rzzTgDAiy++aNZ2OSr2dBAREZmhqKgIqamp6NGjh9HrPXr0wPfff1/pMsnJyRXe37NnT6SkpKC42O7z8todezqIiMg56Evk2GHVsn/Jz89Haen1ng9PT094elYMts7JyUFpaSnCw8ONXg8PD0dWVuWB9FlZWZW+v6SkBDk5OahXr56ZO1A71IpGhxtMb6hqB6QndFL+BFVOCQBQPWkz97ovV/nTwDJSTglpu30VZVIsv5Sj4KpQror1l/JVeAnlqhwiUj4Kad0SVT4CKS9EHaFclZNClesCAHKFcuk6tSSni9TJrbp/pHtX6qJVbVuosGykUC4dM9V+NxCWlfKuqNYtXWeq+x6Qj2mGoizPxOt2jwexUqMjMtL4KpgxYwZmzpxpclGdzvisa5pW4TXp/ZW97oxqRaODiIjIXjIyMlCnzvWvA5X1cgBASEgIXF1dK/RqZGdnV+jNKBcREVHp+93c3BAcHGzhljs+jukgIiLnYKXolbp16xr9mGp0eHh4ID4+Hjt27DB6fceOHejUqVOly3Ts2LHC+7dv346EhAS4u6v6gp0DGx1EROQcrNToiI+PR2xsLJYuXSquMikpCe+//z4++OADHDt2DJMnT0Z6ejrGjh0LAJg6dSpGjhxpeP/YsWPxxx9/ICkpCceOHcMHH3yAVatWYcqUKYb3FBUVIS0tDWlpaSgqKsKZM2eQlpaG3377zYoHq2bw8QoREdENUlNTjR6vqAwZMgTnz5/H7NmzkZmZibi4OGzduhWNGjUCAGRmZhrl7IiOjsbWrVsxefJkLF26FJGRkViyZIkhXBYoe7zTrl07w9/z58/H/Pnz0bVrV3z77bfW2ckawkYHERE5BysNJK2ucePGYdy4cZWWrVmzpsJrXbt2xcGDB03W17hxY8PgUmfDxytEROQctFLzH61o10Nkq/N4haqnVvR0FMJ06JVq2I0UIiY1iKVy1dTdUnintwXr/lNYVpr2WxVSKyXclVqp0raplpeGUFkSEijlGpSmFPeWYmpNze0N+SaTvs+oguikADtLA/CuKMqk+yNQKFd1XkvXgiWh2dI1Lq27QChX3dtSbIIUhpyjKMsVlrXk8wxQh8xWnpFCvu+sTl8C6M3sIbghI2l1Hq9Q9bCng4iIiOyCjQ4iInIONRC9QtVTKx6vEBERifh4xeGxp4OIiIjsgj0dRETkHKzU00G2w54OIiJyDhzT4fDY00FERHQDjumwnVrR6LgE01Nsq6aplvJVSNNIq6b1lkj5EaSp1lXbpkgJUaVyVR4BKX+INC24dMxUsf5S56aUPyFIUSblILgolBcIB1WVy1A6ZtK1oMrNIOWrkDqapXtEVS7lb1TlgwHUeTyKhWUlAYoyabsuW1iuyrUhHW+JqSnkAeC8sKx0rUg5QtIVZabuXfvn6Si14PGKualMqTpqRaODiIhIpC8B9GaOGmCjwy44poOIiIjsgo0OIiJyDhxI6vD4eIWIiJyDlR6vcCCp7bDRQUREzoFjOhweH68QERGRXbCng4iInINWan6PhWZmqC1VS61odDSE6Q1VXV5SHg4pf4KUsyJXUSbFp0t5PFTx+JLfhXJVPH4DYdkQodxPKFfl8bA0P4KnokyVwwOQz5eUS0N1PqXrSMo/osoxIn28SvslPbVWHTcpr4N0zH0VZVeEZaX7o0BRJp0P6V+PtN+W5KaQcmmoyqXjLeXQkXKIqO4vU3XbP09HCaCXPl1NLctGhz3w8QoREdENGL1iO7Wip4OIiEhkpZ4ORq/YDhsdRETkHPh4xSJHjhyp8ntbt25t1jrY6CAiIiK0bdsWOp0OmolBteVlOp0OpaXmjdhho4OIiJwDezoscvLkSZuvg40OIiJyDmx0WKRRo0Y2XwejV4iIiKiCf/3rX+jcuTMiIyPxxx9/AAAWL16Mzz//3Ow6a0VPRxxMx5BfUCwntVvrCuWquHRAHa9/XljWkoS7luaUOKsok3KbqHIrAHIOhDChXEWVewFQ5ymQ9svS3A0qUn4E6VpQfTOQzkexheUlijJpWVvms/ASylX3tnRf5wrlUj4MVb4L6TqTqM63dI3nCuX5QrnqHjCVQ8f+eTr05n+4Mgu6kWXLlmH69OmYNGkSXn31VcMYjoCAACxevBj9+/c3q172dBARkXPQW/jzF+bpAN5++22sXLkSL7/8Mlxdr391SkhIwI8//mh2vbWip4OIiEh0U+Oh2sv+hXk6ygaVtmvXrsLrnp6eKCiQ+p1NY08HERERGYmOjkZaWlqF17/88kvExsaaXS97OoiIyDlYqaeDgOeeew6JiYm4du0aNE3Df//7X6xbtw5z587F+++/b3a9bHQQEZFzYKPDakaPHo2SkhI8//zzuHLlCoYNG4b69evjrbfewtChQ82ul40OIiIiqmDMmDEYM2YMcnJyoNfrERZmSfxhmVrR6AiG6TA3VdhejlDvJfM2x0AVGidN8S5NYa0ihUlKw59UA3lUIciAPP11PaFcFeoYKiwrhd+ppjuXwhyl6eXdhB3XLPiWVKqKS4X6Jr1NqFsa7mUq1LEq5VL4Z0OhXLVf0nUopTBSnS4p/FkKqZXuAVW5FEqfK5SrLjPpXJ4WyqX7S/WZZarM7p0H7OmwmmnTpmHmzJlwdXVFSEiI4fW8vDyMHTsW69atM6teDiQlIiLnoMH8cFkmJDXy4YcfonPnzjhx4oThtW+//RatWrXCqVOnzK6XjQ4iIiIycuTIETRu3Bht27bFypUr8dxzz6FHjx54/PHHsXfvXrPrrRWPV4iIiER8vGI1/v7+WL9+PV5++WX8v//3/+Dm5oYvv/wS9913n0X1sqeDiIicg5UyklKZt99+G4sWLcKjjz6KJk2aYMKECTh8+LBFdbLRQUREdAOmQQd69+6NWbNm4cMPP8TatWtx6NAh3HPPPejQoQPeeOMNs+vl4xUiInIOTINuNSUlJThy5AgiIyMBAN7e3li2bBkeeOABPPnkk3j++efNqtemPR179uxBv379EBkZCZ1Oh88++8yWqyMiolsZH69YzY4dOwwNjhv17dvXcSd8KygoQJs2bTB69GgMGjTI7HpU14OqLZor1CvlZhDSJyinkvYXlpVI04arSFOpq7ZbJywr5TaRcjcEK8qk/CNS+UVFmZSvQjpmnsLF4qXY8Xwh6YQqv4gkUCiPEMqvCeWZijJpKnXpHlDlfZCuI6lu1bcpKX+PNMW7lMdD9bmRISyruoYBdS4N6fNKunel3EKqc2LqOisB8LNQr1VxIKld3Ji3o7ps2ujo3bs3evfubctVEBERkRUEBQXh119/RUhICAIDA6HTmf4aeuGClMKvcg41pqOwsBCFhYUAYNHUuUREdAtiT4dFFi1aBD+/sj6vxYsX22QdDtXomDt3LmbNmlXTm0FERLURGx0WGTVqVKW/W5NDNTqmTp2KpKQkAGU9HZUNYiEiIiLbKy0txebNm3Hs2DHodDq0bNkS/fv3h5ub+U0Hh2p0eHp6wtOzbJiWq6s0tI+IiOgG7OmwmqNHj6J///7IyspCixYtAAC//vorQkNDsWXLFrRq1cqsepkcjIiInANDZq3mySefxO23347Tp0/j4MGDOHjwIP7880+0bt0aTz31lNn12rSn4/Lly/jtt98Mf588eRJpaWkICgpCVFSULVdNREREZjp8+DBSUlIQGHg9MD8wMBCvvvoq7rzzTrPrtWmjIyUlBd27dzf8XT5eY9SoUVizZk2V6wkE4GWirK5iOSmSWMr7IMXrq7qJPIRlpQOfqyhT5TcALMvTIcX6FwnlqrwOgDo3SqiwbD2hXLVfUo4CS3KyAIC34mIqEC6ky9JBVTB1X5SztCtTdY9IOUKka1yVn0Q6X9I1riqX7h9p3VK5aob0s8KyUs4W1fmQPnMChHJpedX9aSrWULqvrI6PV6ymRYsWOHv2LG6//Xaj17OzsxETE2N2vTZtdHTr1g2aproFiYiIrISNDqt57bXXMGHCBMycORMdOnQAAOzfvx+zZ8/GvHnzkJ9//dtU3bqqr//GHGogKREREdW8Bx54AADwyCOPGJKElXci9OvXz/C3TqdDaakqV64xDiQlIiLnoMH8QaQWdMq/++67iI6OhpeXF+Lj4/Hdd98p3797927Ex8fDy8sLTZo0wfLlyyu8Z+PGjYiNjYWnpydiY2OxefNmo3Jbz222a9cuw8/OnTuxc+fOSv/euXNnteplTwcRETmHGni8smHDBkyaNAnvvvsuOnfujPfeew+9e/fGTz/9VGnAxMmTJ9GnTx+MGTMGH330Efbt24dx48YhNDTUMEdZcnIyhgwZgn/84x8YOHAgNm/ejEceeQR79+5F+/btAVhvbjNTunbtavU6AUCnOeigi4KCAvj6lg2bWgnTA+ZUA5Wk8Xm36kBS1cRO0sAvKTm9dDGpYpakgaTSZHSWDEyUJhhrIMxyHajY+Bxhlq8/hAtV9VkoTegmdWVK51N13CwdSKoaVHleWDZAKFfdA9KMEWeEconqHjglLGvLgaSFQrm0vOoyVQ0kLf8efPnyZZtMF3/j/4rLg4A6Zn6VLigBfDeW/V6dbW3fvj3uuOMOLFu2zPBay5YtMWDAAMydO7fC+1944QVs2bIFx44dM7w2duxYHD58GMnJyQCAIUOGID8/H19++aXhPb169UJgYCDWrVtXoU6dTofNmzdjwIABVdrmmsTHK0RERDfIz883+imfE+xmRUVFSE1NRY8ePYxe79GjB77//vtKl0lOTq7w/p49eyIlJQXFxcXK95iqszap9Y9XVNNMhwnLSr1p0vTyqm9TUiijNOxG9a1e6oGp+pCeiqQLQuqNkL5Nqb7dSt8rIoPU5b6KHb8kfIW8LKy7VLgYvBXdR4Hh6mWL/lSXuytOip9wTEqFrqti8yaKrBKpx03VIyD1allSt9RT6C+US58bqnVLdUs9harL0NJwd+lzRdUTYmpZSz6LzGKlxys3T8ExY8YMzJw5s8IiOTk5KC0tRXi48U0eHh6OrKysSleTlZVV6ftLSkqQk5ODevXqmXyPqTprk1rf6CAiIgIATV/2Y+6y5TIyMower5RPz2HKzVPAl0d1VOf9N79e3TqtSdM0pKenIywsDN7e0sPn6uHjFSIiohvUrVvX6MdUoyMkJASurq4VeiCys7Mr9FSUi4iIqPT9bm5uCA4OVr7HVJ3WpmkamjVrhtOnT1u9bjY6iIjIKej1lv2Ui4+PR2xsLJYuXapcn4eHB+Lj47Fjxw6j13fs2IFOnTpVukzHjh0rvH/79u1ISEiAu7u78j2m6rQ2FxcXNGvWDOfPS0O6q4+PV4iIyClY6/FKampqlaNXkpKSMGLECCQkJKBjx45YsWIF0tPTMXbsWADA1KlTcebMGXz44YcAyiJV3nnnHSQlJWHMmDFITk7GqlWrjKJSJk6ciHvuuQfz5s1D//798fnnn+Prr7/G3r17De+x9dxmb7zxBp577jksW7YMcXFxFtdXjo0OIiIiMw0ZMgTnz5/H7NmzkZmZibi4OGzduhWNGjUCAGRmZiI9Pd3w/ujoaGzduhWTJ0/G0qVLERkZiSVLlhjl2ujUqRPWr1+PV155BdOmTUPTpk2xYcMGQ44OwHpzm5ny2GOP4cqVK2jTpg08PDwqjO24cMG8Uei1Pk+HahS7paPQpY4lW0av5CjKpPHLUt2qifCkVqh0TKXolWuKsibCslL0iioTr6XRK8FCEoOGLRR156qXPVuD0SvnhM8N1T0g5elQD71T58O4KCwr7LYyCkRImyLeP5ZEr6QrygB1fh5AHb2imkwRkKN2VPcmoI5eMXW+SgH88Nfv9sjTcaG3ZXk6gv5Ki9GiRQu4uLggMTERiYmJVtrS2uWf//ynsnzUqFFm1cueDiIicgo18XjFWZnbqJDUikbHHzDdylb1KEgtfykPh/TNQPVNT8ojIB141RhlqW7pVlHF60vBUVIWV+mY5yrKpKyh166oy1W9DUHCyTyfqS73FA6Mm2LHvYQTEiQklFF9iHr5qJeV5Ak9HQGKsrpC788VITGEqidEyrQq7bbq3paucen+Ei5DZS+l1BMorVvVGyHl2ZC2W+rhUZ1uUz2FnLi1djtx4gRWr16NEydO4K233kJYWBi2bduGhg0bVpjyvqoYvUJERE7BWtErVDYpXatWrXDgwAFs2rQJly+XNS2PHDmCGTNmmF0vGx1EROQU7B0y68xefPFFzJkzBzt27ICHx/V+ru7duxvmiDFHrXi8QkREJNE0C8Z03DACmGM6gB9//BEff/xxhddDQ0Mtyt/Bng4iIiIyEhAQgMzMioPdDh06hPr165tdLxsdRETkFMqjV8z9oeuGDRuGF154AVlZWdDpdNDr9di3bx+mTJmCkSNHml0vGx1EROQUOKbDel599VVERUWhfv36uHz5MmJjY3HPPfegU6dOeOWVV8yul2M6iIiIbsAxHYC7uzvWrl2L2bNn49ChQ9Dr9WjXrh2aNWtmUb21otFxFqbzP6jSL0gZ9qRLSsq+qYprF1IUKPMfAOr8I8HCsp5COtRLigMjbbeUp6OOn7pcUyTjyBbqLhROaJGUWEUhSJi8sUhIsHDhrOkyF6E/MaKxutxVkQPkWoF62YuK7QKAugHq8kDFJ4RO2K8i4YSq8j4IKUDEDLKqBLTSZSLlmhFSuiizuAqnS7z/VNsWICxbVyiXcoio8hKZypIs5f6wNk0PaGbO/s7HK5Vr2rQpmjQpyxet05l5cG/AxytEROQUmKfDulatWoW4uDh4eXnBy8sLcXFxeP/99y2qs1b0dBAREZH9TJs2DYsWLcIzzzyDjh07AgCSk5MxefJknDp1CnPmzDGrXjY6iIjIKVjr8Up8fPwtP+HbsmXLsHLlSjz66KOG1x588EG0bt0azzzzDBsdRER0a9PrAb2ZjY4bH69wIClQWlqKhISECq/Hx8ejpESYvlqBYzqIiMgpcEyH9Tz22GNYtmxZhddXrFiB4cOHm10vezqIiIioglWrVmH79u3o0KEDAGD//v34888/MXLkSCQlJRnet3DhwirXWSsaHa4wPeVzkGI5KfxMiCwVl1f14mmKMkCe8lm1bfWEXj9X4awWKUJPpQvCXapbiEdUhUJKIcpSh94VRThu/Rj1sqqp6QGgUNivM8cVdavmcAfgJZxPH0WsY6lqDncAxUIcpLcQ4qw6Lrnn1MuqwlYB9fmU7j1TIZrlVGGrwiETWRKKL3XYS9vmrSjLFZaVnjpIx1y17WEmXje/E95MFozpED+UbzFHjx7FHXfcAaBsinugbN6V0NBQHD161PC+6obR1opGBxERkcRaYzoI2LVrl03q5ZgOIiKiGzANuu2wp4OIiJyCtUJmGb1iO2x0EBGRU2AadMfHxytERERkF+zpICIip6DXzB8QqpdCDm8xBQUFNnnExJ4OIiJyCpresh+6Ljw8HE888QT27t1r1XprRU/HVZiO91bFrUstKml6bCkeP1QoV5HyDKjWXSzMjy1Nza26t6Sp66XpzN2FnBShwabLrgonxFvYODdFEhApD4eUr6KOkEREKlfxiRTeoLhLXYTp469dUZdLOUJUH8TStSCxJM+Nqbw9VaHKFQPI94D0uaDKTSGtWzhdyv2W/mdK08xbkhvFVN12z9NBVrNu3TqsWbMG9913Hxo1aoQnnngCI0eORGSk9IGlxp4OIiJyCkyDbj39+vXDxo0bkZGRgaeffhrr1q1Do0aN8MADD2DTpk1mz7/CRgcRETkFazU6mKfjuuDgYEyePBmHDx/GwoUL8fXXX+Phhx9GZGQkpk+fjitXpP45Y7Xi8QoREZFE08uP5lTLlmOejuuysrLw4YcfYvXq1UhPT8fDDz+Mv//978jIyMDrr7+O/fv3Y/v27VWuj40OIiIiMrJp0yasXr0aX331FWJjY5GYmIjHHnsMAQEBhve0bdsW7dq1q1a9bHQQEZFTsFZPBwGjR4/G0KFDsW/fPtx5552VvqdJkyZ4+eWXq1UvGx1EROQU9HrzJ4vlQFJjmZmZ8PHxUb7H29sbM2bMqFa9HEhKRERERvz8/JCdXTEm//z583B1NT9ovVb0dFyE6Q31UiwntagumLc5Bqoxu+r2oZxnQJXDIFdYVsoFoEiVIR6zUiFKykOVOAWAu+KE1QlQLyvllMg/b7os95x6WVWODwDwaSGUq/KASIlTvISjfs30V7CgCPWiIUJIfVGhutxTcT6layHIgq80lnZ1X1SU5UvrtrBclQ9DcYkCKMtJpKL6XJDybEh5OqS8RarTbep0Seu0Nj5esR5Nq/xIFhYWwsND+i9jWq1odBAREUn4eMVyS5YsAQDodDq8//778PW9ni6vtLQUe/bswW233WZ2/Wx0EBEREQBg0aJFAMp6OpYvX270KMXDwwONGzfG8uXLza6fjQ4iInIKfLxiuZMnTwIAunfvjk2bNiEwMNCq9XMgKREROQVmJLWeXbt2Wb3BAbCng4iIyMitmpE0KSkJ//jHP1CnTh0kJSUp37tw4UKz1sFGBxEROQVNM/8xiYlgjVvKoUOHUFxcbPjdFJ1OFUelxkYHERE5BU0P6M38f8hGR9kjlcp+t6Za0ejIgem8FkJaCKUCoVwa8CLFtatIeTz8FWXSPhcL5ap4fimu3lO4YvLy1OWXFOUNmqiXLRF2rPia6bIrl9TL+gaoy3081eVo3cd0WYt+wsKC41+aLHIt3aJctHmouuris+ryK4qkFlKOjxIhccQ1RaKbK8J1JP1fkU6XSq5QbkkeD0vzVqhyA0m5f1SfKYA65xGg/lwx9Vlp3uTn5tP0gMZGh03k5+dj586duO222ywKmeVAUiIiIjLyyCOP4J133gEAXL16FQkJCXjkkUfQqlUrbNy40ex62eggIiKnYK3oFQL27NmDLl26AAA2b94MTdOQm5uLJUuWYM6cOWbXy0YHERE5BU1v2Q9dl5eXh6CgIADAtm3bMGjQIPj4+KBv3744fvy42fWy0UFERERGGjZsiOTkZBQUFGDbtm3o0aMHAODixYvw8pJGAJlWKwaSEhERSfQWRK/oOZDUyKRJkzB8+HD4+vqiUaNG6NatG4Cyxy6tWrUyu142OoiIyCkwesV6xo0bh7vuugt//vkn7r//fri4lD0YadKkiUVjOmpFoyMXpp8DWRJaKkT8icurwvKkcFzpwKuiDaVQXalcNX22FOLmKs2fLVA9NtX/rl42vJ66PDDcdJkUbusfoi5HHeFJZKtHTZe1fkyoXOCiuFp+/1q9bIkqyBJwz1UvXqw4365CjGapcMjOKcJipctMuk5VE29L4bTSo33pc0EVFiuFrSqivgEAmYoy6ZhJz9Kl46I6pqbOh/kppMgRJCQkICEhwei1vn37WlRnrWh0EBERSfh4xXpKS0uxZs0afPPNN8jOzob+pvCenTt3mlUvGx1EROQU2OiwnokTJ2LNmjXo27cv4uLiLEp9fiM2OoiIiMjI+vXr8cknn6BPH0W2ZTPYJWT23XffRXR0NLy8vBAfH4/vvvvOHqslIqJbSE3l6aju/7jdu3cjPj4eXl5eaNKkCZYvX17hPRs3bkRsbCw8PT0RGxuLzZs3V3u9mzZtQs+ePRESEgKdToe0tLQq75OHhwdiYmKq/P6qsnmjY8OGDZg0aRJefvllHDp0CF26dEHv3r2Rnp5u61UTEdEtpCYaHdX9H3fy5En06dMHXbp0waFDh/DSSy9hwoQJRqnFk5OTMWTIEIwYMQKHDx/GiBEj8Mgjj+DAgQPVWm9BQQE6d+6M119/vdr79eyzz+Ktt96CZuWwHp1m7Rpv0r59e9xxxx1YtmyZ4bWWLVtiwIABmDt3rsnlCgoK4OvrCwCIhenWUWPFuiOEbcsVyqVR6mGKMmkkuTSKXVVeR1hWil5RTTYnBXFIk0pJVPd1kLCsLaNXGjZXl7u2FdrnA/5puszS6JWD75su+3KietkL6ugVnFYX55xRVK0KpYA8IVy6om7p/hH2CucUZTnCsqroLsCy6BXpw9aW0SvSZ470P1e17aaiV0oA/Puv3y9fvow6daRPr+q78X/FlwC8zRx6cFUDev/1e3W2tbr/41544QVs2bIFx44dM7w2duxYHD58GMnJyQCAIUOGID8/H19+eX2yx169eiEwMBDr1q2r9npPnTqF6OhoHDp0CG3btq3Sfg0cOBC7du1CUFAQbr/9dri7uxuVb9q0qUr13MymPR1FRUVITU01ZDIr16NHD3z//fcV3l9YWIj8/HzDDxERkb3d+H8oPz8fhYWVt6Cr+z8OKOvFuPn9PXv2REpKCoqLi5XvKa/TnPVWV0BAAAYOHIiuXbsiJCQE/v7+Rj/msulA0pycHJSWliI83PgraHh4OLKysiq8f+7cuZg1a1aF1yNg3jds6QGO1OKSpqFWfeORppdXxbwD6ph56duQtG5Vc07x5RMA0FQobyCUq76hSsfbVbha6yjugzp1hbqDhZUXCN8DC1TfrS109aLpMqknQ5i6PuM3dbmqt8JD+Mp/RfjeoOoIkb61S70NqtPpKywr9RRK+SxU17FUt3TvqnphLgjLKq4iAHJPo4qpddt9anuYn+TrxsUiIyONymbMmIGZM2dWWKa6/+MAICsrq9L3l5SUICcnB/Xq1TP5nvI6zVlvda1evdoq9dzMLtErN4faaJpWafjN1KlTkZSUBKCsy+zmE09ERGSKHvJjItWy5TIyMower3h6qpuaVf0fp3r/za9Xpc7qrre6SkpK8O233+LEiRMYNmwY/Pz8kJGRgbp16xoeaVWXTRsdISEhcHV1rdDyys7OrtBCA8pObPnJdZXSHRIREdlA3bp1qzSmo7r/4wAgIiKi0ve7ubkhODhY+Z7yOs1Zb3X98ccf6NWrF9LT01FYWIj7778ffn5+eOONN3Dt2rVKI26qwqZjOjw8PBAfH48dO3YYvb5jxw506tTJlqsmIqJbjN7Cn3Lx8fGIjY3F0qVLlesz539cx44dK7x/+/btSEhIMAzWNPWe8jrt8b914sSJSEhIwMWLF+Htff2B5sCBA/HNN9+YXa/NH68kJSVhxIgRSEhIQMeOHbFixQqkp6dj7Nixtl41ERHdQjTIEUKqZculpqZWOXpF+h83depUnDlzBh9++CGAskiVd955B0lJSRgzZgySk5OxatUqQ1QKUPYP/5577sG8efPQv39/fP755/j666+xd+/eKq8XAC5cuID09HRkZGQAAH755RcAZT0pERHq+M69e/di37598PAwHoHYqFEjnDkjjf4zzeaNjiFDhuD8+fOYPXs2MjMzERcXh61bt6JRo0a2XjUREZFNSf/jMjMzjXJnREdHY+vWrZg8eTKWLl2KyMhILFmyBIMGDTK8p1OnTli/fj1eeeUVTJs2DU2bNsWGDRvQvn37Kq8XALZs2YLRo0cb/h46dCgA0wNjb6TX61FaWnFY9OnTp+Hn51e9g3QDm+fpMNeNsdf3wnT0imp4j6WzLkoj5EMVZdIodGmkeIBQriKtWzXQSoogsWX0iip/CAA0bKgur9fEdJkUveIhPQaVTtiAhabLOk4WFhbse9N02abn1ctK0SvH1OVSrg2VXGHdJxUJM6R7V4r+Us3yLO2SI0evqNKqSNEr0jG1JHol28TrJQC++ut3e+Tp+ByAl5n1XAPQ/6/fW7RoARcXFyQmJiIxMdEKW1n7DBkyBP7+/lixYgX8/Pxw5MgRhIaGon///oiKijI7uoVzrxARkVOwVvRKdR6vOKtFixahe/fuiI2NxbVr1zBs2DAcP34cISEhRo+CqqtWNDpKYPo5nSo4SPpWIbXspfQnqph5KfZG2jbVtxLpW57UQ2NJ3VIwlpTNUfWBIPU8FQh5H65eEipQ8JC+vkpfQZMXmL9yT6Gr8thG02XC1/ZzvwvlQkZSb0VUnCoDLABcFS7yMMX5vCp8LZfuL9XplHo6pLotyaUh1S39w1Ttl3T/CJ19Fm2bqTILpjOhGhYZGYm0tDSsX78eqamp0Ov1+Pvf/47hw4cbDSytrlrR6CAiIpJYq6cjPj7+ln+8smfPHnTq1AmjR482GhdSUlKCPXv24J577jGrXjY6iIjIKdRE9Iqz6t69OzIzMxEWZjzLWF5eHrp3717pINOqYKODiIicgrUaHWQ6u+n58+ctapCx0UFEREQAgIceeghAWYr1xx9/3CgFfGlpKY4cOWJRAjI2OoiIyClYa0zHrax8BllN0+Dn52c0aNTDwwMdOnTAmDFjzK6fjQ4iInIK1nq8cisPJC3Pv9G4cWNMmTLF6mNbakWj4ypMh3O5K5aTwj8tndpeFeho6aQ2uYqyEGFZad2qcikQKk8ol6jOiTQNdoGwchchPFTFU7ivmrRSl/ucUqQFzk9SL+winLEixXcwYVEPIVOSnxA3XnzNdJkUwuwixFfXb2a67LyQZTkvV12uIs2NKYXESpHZuYoyKTJbuv9Uyfekb+pScjDpn4EqpNbUZ5J0LB0VB5KWZS21BZtO+EZERGQv1prwjYCzZ89ixIgRiIyMhJubG1xdXY1+zFUrejqIiIgkjF6xnscffxzp6emYNm0a6tWrV2kkiznY6CAiIrrBrTymo9zevXvx3XffoW3btlatl40OIiJyCpx7xXoaNmwIW8wHyzEdRETkFDQLf+i6xYsX48UXX8SpU6esWi97OoiIiMjIkCFDcOXKFTRt2hQ+Pj5wdzeOFb1w4YJZ9bLRQUREToHJwaxn8eLFNqm31jc6VM+HpKnrpSd2qhwgAOCvKJNyfEgx86p4fumkWTDDu3hMFGkbAMjP61S5TaRjIuUwUE2lrhM2LD9XXR5cT13uoypXzXUOoPic+uNOlcbDVcrDEaguLxJOaNYp02XXrqiXLRUSr6jOiV7o65buL1VAn5Q7IlQoV13DACCkL1GSpp9X7bcqhwcACKdLukyVn4em7l175+lg9Ir1jBo1yib11vpGBxEREcCp7S2Vn5+PunXrGn5XKX9fdbHRQUREdINbNXolMDDQMJ19QEBApbk5ymef5dT2RER0S+OYDsvs3LkTQUFlAxN27dplk3Ww0UFERE6DYzPM17Vr10p/tybm6SAiIiK7YE8HERE5BUavOD42OoiIyClwTIfjc+pGh9RyLRDKpbHLirQQyjJAzkmhWrel26066dKNFyyUS7lRVHykdQs5KeoqNq5ISELgLiQU0IQDU5RjuqxASNxwybzEfgAAP+GAe6gSvkCd2wRQ5+LIzVMvK90DAYr9lu4PKeeE6hqXPhc8hetMSlbTQlEmpC4Rc3yolg+xsG4pB0+uUE5UFRzTQURETsFac6/Ex8cjNjYWS5cutd/G3yKcuqeDiIhuHZxl1jLt2rWrNDdHZQ4ePGjWOtjoICIiIgwYMMDm62Cjg4iInAKjVywzY8YMm6+DjQ4iInIKjF5xfGx0EBGRU2BPh/WUlpZi0aJF+OSTT5Ceno6iIuOYsgsXzAu7Y/QKERERGZk1axYWLlyIRx55BHl5eUhKSsJDDz0EFxcXzJw50+x6a0VPxyWYbh2p4tYDhXqlnS8UylW5BKT596RyV0WZlCsjTChXtTSFtA7wFcqlXBuq/fIVmsABoeryYsUJOys0yv0svBMu55ouyz2nXrZQyCFiiRIL8494K074ZSFPhyJ1CQD1/ae6TgD5/lHVLRwSXBESVkjfiIMUQQ9XhSQ70qWgWrd070qXuLdQ/rui7KKJ182bh9R8fLxiPWvXrsXKlSvRt29fzJo1C48++iiaNm2K1q1bY//+/ZgwYYJZ9bKng4iInIK18nQQkJWVhVatWgEAfH19kZdX9i3jgQcewH/+8x+z62Wjg4iIiIw0aNAAmZmZAICYmBhs374dAPDDDz/A01PqVzONjQ4iInIKegt/6LqBAwfim2++AQBMnDgR06ZNQ7NmzTBy5Eg88cQTZtdbK8Z0EBERSawVvRIfHw8XFxckJiYiMTHRCltW+7z++uuG3x9++GE0aNAA33//PWJiYvDggw+aXS8bHURERDe4VdOgq3To0AEdOnSwuB42OoiIyCkwesW6fv31V3z77bfIzs6GXm98hKZPn25WnbWi0eEO0yF0/orlpOmxpWltpHAvU2FiAFBXWFaiWrc0hEeamVsVjiiFKkrll4RyD0VZkXDXu55VlxcqTrg0rXexNOf4EXWxj+KEewixiFIosJcqBFOYP/58prrc3V1druKhOpkAgoQbUBXFLA02k+4BVRe7dI1K5dI/p1BFWKwUhi+FzKoOufRYQfVZWZVy1WeSqTB96TPY2pgczHpWrlyJp59+GiEhIYiIiDCaCE6n0zl3o4OIiIjsZ86cOXj11VfxwgsvWLVeNjqIiMgp8PGK9Vy8eBGDBw+2er0MmSUiIqegwfxwWT5eMTZ48GBDbg5rYk8HERE5BY7psJ6YmBhMmzYN+/fvR6tWreB+0wAwc9Ogs9FBRERERlasWAFfX1/s3r0bu3fvNirT6XRsdBAR0a2NPR3Wc/LkSZvUy0YHERE5BQ4ktQ1NK2uS3Rg2a65a0ejwgun8EKqccdIU1hJplK3qIhVSL4jbpsqHIW2XFBtvyXGxNLeJKg+BlNvERdgxVbF0TKT8I9ekPB6KRCBSLo26QeryUsUJy81WL5snTD8vLK68jiPqqZd1FZKjZCnyWZxXLypew6rTJX3onRPKhdnplfendJ1J94+qXLrG/YRy6d6ONGPd14Q6ybF9+OGHePPNN3H8+HEAQPPmzfHcc89hxIgRZtdZKxodREREEvZ0WM/ChQsxbdo0jB8/Hp07d4amadi3bx/Gjh2LnJwcTJ482ax62eggIiKnwDEd1vP2229j2bJlGDlypOG1/v374/bbb8fMmTPNbnQwTwcREZEZ9uzZg379+iEyMhI6nQ6fffaZVerdvXs34uPj4eXlhSZNmmD58uVG5WvWrIFOp6vwc+2a9R5oZWZmolOnThVe79SpEzIzhfkVFNjoICIip6BZ+FNdBQUFaNOmDd555x3LN/4vJ0+eRJ8+fdClSxccOnQIL730EiZMmICNGzcava9u3brIzMw0+vHykmbeqrqYmBh88sknFV7fsGEDmjVrZna9fLxCREROwd5jOnr37o3evXubLC8qKsIrr7yCtWvXIjc3F3FxcZg3bx66detmcpnly5cjKioKixcvBgC0bNkSKSkpmD9/PgYNGmR4n06nQ0REhBlbXTWzZs3CkCFDsGfPHnTu3Bk6nQ579+7FN998U2ljpKrY00FERGQDo0ePxr59+7B+/XocOXIEgwcPRq9evQzRIJVJTk5Gjx49jF7r2bMnUlJSUFx8PW7r8uXLaNSoERo0aIAHHngAhw4dsuq2Dxo0CAcOHEBISAg+++wzbNq0CSEhIfjvf/+LgQMHml0vezqIiMgpWGsgaX5+PkpLrwcoe3p6wtPTs1r1nThxAuvWrcPp06cRGVkWcDxlyhRs27YNq1evxmuvvVbpcllZWQgPDzd6LTw8HCUlJcjJyUG9evVw2223Yc2aNWjVqhXy8/Px1ltvoXPnzjh8+LBFjz5uFh8fj48++shq9QG1pNGhytOhGjZjaZ4Od6FclV5BiseXqE6MdNKkHAdXFGUthGX9hXIhLYQyf4K3sGxNXqyqYwYAxYodk/I6uJ9SlzdUnBQ3D/WyUt4H6ZiqUoxcuijULdxAQYp+VhehrztHXazcbumYSClZpKfmqvQk0rLSZ5bqsEjXqCU5QAAgwIxlLf0Mri5rPV4pbySUmzFjBmbOnFmt+g4ePAhN09C8eXOj1wsLCxEcHAwA8PX1Nbz+2GOPGQaM3pyE6+bkXB06dECHDh0M5Z07d8Ydd9yBt99+G0uWLKnWdt4oPz8fdevWNfyuUv6+6qoVjQ4iIiKJtXo6MjIyUKfO9dST1e3lAAC9Xg9XV1ekpqbC1dX4a2h5YyMtLc3wWvk/8YiICGRlZRm9Pzs7G25ubobGys1cXFxw5513Kh/bVEVgYCAyMzMRFhaGgICASjOQapoGnU5n1BNUHWx0EBER3aB79+5wcXFBYmIiEhMTzaqjXbt2KC0tRXZ2Nrp06VLpe2JiYiq81rFjR3zxxRdGr23fvh0JCQkVZnotp2ka0tLS0KpVK7O2tdzOnTsRFFTWh79r1y6L6jKFjQ4iInIK1nq8kpqaatTTYcrly5fx22+/Gf4+efIk0tLSEBQUhObNm2P48OEYOXIkFixYgHbt2iEnJwc7d+5Eq1at0KdPn0rrHDt2LN555x0kJSVhzJgxSE5OxqpVq7Bu3TrDe2bNmoUOHTqgWbNmyM/Px5IlS5CWloalS5eaufdlunbtavg9OjoaDRs2rPRRz59//mn2OtjoICIip2DvjKQpKSno3r274e+kpCQAwKhRo7BmzRqsXr0ac+bMwbPPPoszZ84gODgYHTt2NNngAMr+2W/duhWTJ0/G0qVLERkZiSVLlhiFy+bm5uKpp55CVlYW/P390a5dO+zZswd33XWXGXthejvKH7Xc6MKFC4iOjjb78YpOKx+h4mAKCgoMz726wPTATNU4upocSBoiLGvJAL9AYVlpIJwlA0mFucnEgaS5irL6wrIBQrlqv6XJsCTCeE3ltSINJI0SZuJSDSTNFWYnS/9DXS6dL9V+hQujIqWBpNcUBybXwoGkqmtBuveE+fnEidFU14qlA0lVdUsDSaWB2rYYSHoNwJi/fr98+XKVeg+q68b/FVMg36umFAGY/9fvLVq0sPjxSm3n4uKCs2fPIjTUeNrHP/74A7GxsSgokD7ZKseeDiIicgr2frzijMp7a3Q6HaZNmwYfHx9DWWlpKQ4cOIC2bduaXb9NGx2vvvoq/vOf/yAtLQ0eHh7Izc01q56rMN3ToWrd+yjKAMu/dagIs3rDVyhvqCiTvrFcEMpVPSWqqcwBQBrDbUm4oVR3rlCu+qYmhTBL5VeFctW3X7HuS0K54qt3kbBh0rc+6TpUVZ8pTPPgI5TXUYTMSveuFKx3RlFmSe8OIH9oqtYtfaYECOWqkHWpN0/6TJKOS+VxE2VM3buFQp3WxgnfLFeeZEzTNPz444/w8Lj+KeLh4YE2bdpgypQpZtdv00ZHUVERBg8ejI4dO2LVqlW2XBURERFZqDxqZfTo0XjrrbfMzsdhik3ToM+aNQuTJ0+2OIyHiIhIorfwp1x8fDxiY2MtjgapzRYvXoySkop91xcuXBATh6lw7hUiInIKGsxvcNz4eCU1NRU//fTTLTuIFACGDh2K9evXV3j9k08+wdChQ82u16EaHYWFhcjPzzf8EBERkf0dOHDAKBy4XLdu3XDgwAGz6612o2PmzJnQ6XTKn5SUFLM2Zu7cufD394e/v3+F3PdEREQqmoU/dF1hYWGlj1eKi4tx9ao0tN60ajc6xo8fj2PHjil/4uLizNqYqVOnIi8vD3l5ecjIyDCrDiIiujVZq9HBMR3AnXfeiRUrVlR4ffny5YiPjze73mpHr4SEhCAkREp9ZZ4bpw++eYIcIiIiFebpsJ5XX30Vf/vb33D48GHcd999AIBvvvkGP/zwA7Zv3252vTYNmU1PT8eFCxeQnp6O0tJSw4x6MTExRlP6SnJguktG1VUj5X2QtkDKvmlJboZ6QrkqX4Z0U0mZBVUJMKW4eilbo7R89edqrPq6zyrKpI8P6VqQslCqBAjl0vk6/avpsnzhYpCuFel8qLqcVccbkPPghCs2TsrTIZ1PVS4bVR4NaVmgLMumiuo6lZaV7h9Vng/pXFqy3YA6f4mpHB58ZFF7de7cGcnJyXjzzTfxySefwNvbG61bt8aqVavQrFkzs+u1aaNj+vTp+Oc//2n4u127dgDK4oC7detmy1UTEdEthsnBrKtt27ZYu3atVeu0aaNjzZo1WLNmjS1XQUREBMB6j1fI2NWrV1FcbNzPZm7SMIcKmSUiIqppHEgKXLlyBePHj0dYWBh8fX0RGBho9GMuTvhGREROwVqPVziQFHjuueewa9cuvPvuuxg5ciSWLl2KM2fO4L333sPrr79udr1sdBARkVPg4xXr+eKLL/Dhhx+iW7dueOKJJ9ClSxfExMSgUaNGWLt2LYYPH25WvXy8QkREREYuXLiA6OhoAGXjNy5cKJu//O6778aePXvMrpeNDiIicgrMSGo9TZo0walTpwAAsbGx+OSTTwCU9YAEBASYXW+teLziAtOtI3/FctITuQALy1W5OKQLWMrjcVFRpoqXB+R4fNVJr5j01liOUF4glKuGH0mJdaV1q46ZlHtB6lqVzqcqf4LUspfWrcrFIZ2vIqFcyu2guk6ldUuzJ6mOWYSwrIdQbipvBKC+TgBAyoVsSZ4cKXeJlLPlvKJM+jCXyqVYhCuKMlPbLV1/1sbHK9YzevRoHD58GF27dsXUqVPRt29fvP322ygpKcHChQvNrrdWNDqIiIjsJT4+Hi4uLkhMTLxlZ5qdPHmy4ffu3bvj559/RkpKCpo2bYo2bdqYXS8bHURE5BQYvWIdxcXF6NGjB9577z00b94cABAVFYWoqCiL62ajg4iInAIfr1iHu7s7jh49Cp3OkgkgKseBpERE5BQ4kNR6Ro4ciVWrVlm9XvZ0EBERkZGioiK8//772LFjBxISEio8bjJ3MCkbHURE5BT4eMV6jh49ijvuuAMA8OuvxlNdW/LYhY0OIiJyCpxl1nK///47oqOjsWvXLpvUXysaHaEwnS9ANb5Yys0g5RmQ8kaoYu6ldmChBeWXhWV9hXLVMZPi6qU8HNIxM3+aIPlDQZVz4pywrJTbRMpnodpvKTeDKtcMoB54ZekHpXQdqnJpSB8e0rap8q6ockIAQAOhXHXvRwvLSn4TylXHVPpMko6p6jNLyosikbbtkqLM1P0jfcaS42nWrBkyMzMRFhYGABgyZAiWLFmC8PBwq9TPgaREROQU9Bb+lLuVZ5nVNOOvC1u3bkVBgfR1s+pqRU8HERGRRIP5YzOYp8M+2NNBREREAMoGid48UNSa+TrY00FERE6BA0ktp2kaHn/8cXh6lo1ku3btGsaOHVuh52fTpk1m1c9GBxEROQU2Oiw3atQoo78fe+wxq9bPRgcREREBAFavXm3T+mtFo8MbpkNmVTsghSpaGtaqCl1tLCwbIJSrQjClcEIp9M1PUSaFzEqDgKRy1TG7ICwrnS/VfuUKy0rHVKJatxSOa8nTUh+hXAqfPiuUq8IkpenlpdBuS65xab9U2ybdHxFCubRtpxRlUki5dD5VYbF5wrLSut2FctVnraleAnuHzDI5mOOrFY0OIiIiCR+vOD42OoiIyCmwp8PxMWSWiIiI7II9HURE5BT4eMXxsaeDiIicAtOgOz72dBAREd2AadBth40OIiJyCny84vhqRaOjBKYviLqK5SyZ/hqwbKpoabpyKc9AqaJMtc+AnD/BnCmqy6m2qyrlqmMq5UWRngWq8hBI+yXldDGVJ6ac6pxIuU+kfBaqCaX9hQ2/Juy4lLvhnGrdwrLStaD6kJeWvSiUeyrKpPtHKo8SylW5KU4Ly0q5alQ5QqRzKZVL17jqfJmK/LB3RAijVxwfx3QQERGRXdSKng4iIiIJH684PjY6iIjIKfDxiuPj4xUiIiKyC/Z0EBGRU+DjFcfHRgcRETkFPl5xfGx0EBGRU2BPh+OrFY0Ob5je0EDFcgFCve5CuSreHlDndpBazeeFclVMfbGwrERVd66wrJTvQto26ZiqSB8KoYoy6XxIeR+kdatyVqhyKwBy/gRV/hJP4YToLKgbsOxak/JdqPLFSMdbyn2SLZSrSDl0pPwk9RRlqn0G5P1W5dKQjomUh0MqV91DpvKqsPeAblYrGh1EREQSDeY3dNjTYR9sdBARkVPQQ+7dUy1LtseQWSIiIjPs2bMH/fr1Q2RkJHQ6HT777DOr1Lt7927Ex8fDy8sLTZo0wfLlyyu8Jzc3F4mJiahXrx68vLzQsmVLbN261SrrtyX2dBARkVOw90DSgoICtGnTBqNHj8agQYPMXLOxkydPok+fPhgzZgw++ugj7Nu3D+PGjUNoaKhhHUVFRbj//vsRFhaGTz/9FA0aNMCff/4JPz8/q2yDLbHRQURETsHejY7evXujd+/eJsuLiorwyiuvYO3atcjNzUVcXBzmzZuHbt26mVxm+fLliIqKwuLFiwEALVu2REpKCubPn29odHzwwQe4cOECvv/+e7i7l4VENGrUyIw9sD8+XiEiIrpBfn6+0U9hoRTnVbnRo0dj3759WL9+PY4cOYLBgwejV69eOH78uMllkpOT0aNHD6PXevbsiZSUFBQXl8WTbdmyBR07dkRiYiLCw8MRFxeH1157DaWl0vzMNa9W9HQEwfSGqsLXpPA0iRQeak4IWTkpRFM13bkUiiiFYKrqlqZZl8I/pQvKvFu3TIBQrvqmIoV+Sq1vD6FcVb8Umm1JGLI0aK6OsGNBwug51fm0NKxVdX9K14l0f6mOaYGwbH2hXAqZVXVwRwjLSsdMdVykUF9pvy2J3nCkqe2tMZA0MjLSqGzGjBmYOXNmteo7ceIE1q1bh9OnTxvqmzJlCrZt24bVq1fjtddeq3S5rKwshIeHG70WHh6OkpIS5OTkoF69evj999+xc+dODB8+HFu3bsXx48eRmJiIkpISTJ8+vVrbaW+1otFBREQksVajIyMjA3Xq1DH87enpWe36Dh48CE3T0Lx5c6PXCwsLERwcDADw9b3eVHzssccMA0Z1OuO90DTN6HW9Xo+wsDCsWLECrq6uiI+PR0ZGBt588002OoiIiGqTunXrGjU6zKHX6+Hq6orU1FS4uhqnXitvbKSlpRmtEwAiIiKQlZVl9P7s7Gy4ubkZGiv16tWDu7u7Ub0tW7ZEVlYWioqK4OEh9c3WHDY6iIjIKVhrIGl8fDxcXFyQmJiIxMREs+pr164dSktLkZ2djS5dulT6npiYmAqvdezYEV988YXRa9u3b0dCQoJh0Gjnzp3x8ccfQ6/Xw8Wl7Bnqr7/+inr16jl0gwNgo4OIiJyEtRodqampVerpuHz5Mn777TfD3ydPnkRaWhqCgoLQvHlzDB8+HCNHjsSCBQvQrl075OTkYOfOnWjVqhX69OlTaZ1jx47FO++8g6SkJIwZMwbJyclYtWoV1q1bZ3jP008/jbfffhsTJ07EM888g+PHj+O1117DhAkTzNx7+2Gjg4iInIK9M5KmpKSge/fuhr+TkpIAAKNGjcKaNWuwevVqzJkzB88++yzOnDmD4OBgdOzY0WSDAwCio6OxdetWTJ48GUuXLkVkZCSWLFlilAekYcOG2L59OyZPnozWrVujfv36mDhxIl544QUz9sK+dFr5CBUHU1BQYHjuNQSmW0dhijqqP/THmBRVoJrwTVq3LaNXpImbnDV6RXVcpO2WhAjlAYoyH2HZfKFcdT7bCctK0SunhU9a1bZJHxxShMnvijLpOpG+g6oihqRoIkujV1TRY2eEZdOFctVxyROWtWX0iqmJHEsB7Pnr98uXL1s8TqIyN/6vuAvy558ppQD++9fvLVq0sPjxClWOPR1EROQUrNXTUdXHK1R9taLRoZoOXdXyl6ZRl74ZSAdH6nFQkbryVF9QpW+QlkxnLh0Tad3StyXVfkvHM0AoV33Dkb7lSXk6pDwfuYoy6Zh5C+WqYWHnhWUvCReadFxUuR+knBJnhXJVng6pbqk8UFGWIywr9Yo1EcpVxyxIWFY6H1mKMikBtvS5IB1Tc/Lg2Lsb3d4ZSan6mJGUiIiI7IKNDiIicgqahT/l4uPjERsbi6VLl9pv428RteLxChERkYRjOhwfezqIiIjILtjTQURETsHeeTqo+tjTQUREToFjOhwfezqIiIhuwDEdtlMrGh3+MJ1FMFyxnNRdJsXES1lFVXHxUm4GKSuiKi+EVLe0X6osk1KeDmkqISk3iqprTZXJEVDndQDU+S6kYybl4bAkj4e0X6qcEoA6f4mUQVbabmnbLMldYEn2WWm7pXtbdR1L16gqF0ZVNFKUqbIYA+oMy4D6mJ4WlrV0GjBzcgdJx9ramKfD8dWKRgcREZFEg/ljM9josA82OoiIyCmwp8PxcSApERHRDTiQ1HZs1ug4deoU/v73vyM6Ohre3t5o2rQpZsyYgaIiKcM/ERFR9ekt/CmXmpqKn376iTPM2oDNHq/8/PPP0Ov1eO+99xATE4OjR49izJgxKCgowPz58221WiIiukXx8Yrjs1mjo1evXujVq5fh7yZNmuCXX37BsmXL2OggIiK6Bdl1IGleXh6CgqTJnSsKhOlwr2DFcqrQ0PJ6VaSwVlUI5zVhWSmsVVW3FIYmTc2tCgVWTQ9vDapvE1KIsjQFvCXhedJDPylcV7Xt0rLnhHLVMZNuYGkkvyUfANL5kq4lVci5FLotfStV3X/SuZYyWqYL5Sr1hXIprFWVIkC67/8UyqVjqjpupo63FKpubZZkFWVGUvuwW6PjxIkTePvtt7FgwQKT7yksLERhYVkkekGB9G+ZiIjoOj5ecXzVHkg6c+ZM6HQ65U9KSorRMhkZGejVqxcGDx6MJ5980mTdc+fOhb+/P/z9/REZGVn9vSEiIrIQo1dsR6dpWrUaeDk5OcjJyVG+p3HjxvDyKsu9l5GRge7du6N9+/ZYs2YNXFxMt3Nu7ukob3g8C9PdjqruSunxynmhXHq8ouoKlR6v5ArlNfV45YKwrNRlLm2bqutaevAWYcG6pf2SzpfUBJYeNahI+x2iKLO0q9KSzLnSPkuPIS4qyix9vKLaL0sfr0jd8FGKMksfr6gy0J4Ulq2pxyvf//X75cuXbZJavKCgAL6+vgCAJjA/JFMP4Pe/frfVtpIZn1khISEICVF9DF535swZdO/eHfHx8Vi9erWywQEAnp6e8PQs+yhzdbX16AIiInImfLzi+Gw2piMjIwPdunVDVFQU5s+fj3Pnrg+Xi4iQvrMSERGRs7FZo2P79u347bff8Ntvv6FBgwZGZdV8okNERCRi9Irjs1lG0scffxyaplX6Q0REZG2ahT9ke7ViwjcdTA/uylYsJ+VHsGTqbUA9UFUapCoFBKvKpWWlk6oakSPlLvERyqVjLg0QVJEGyKpII4Sk1rd66LR6MKg0IFka2Kj6BmZpHgRLcp+o7j1AHpBZV1EmXUdSuYq0XRJpsOcpC+qWHjyr7u1QYVnpGs4SylVMXYfM00E344RvREREN2DIrO3Uip4OIiIiibV6OlJTUxkyayNsdBARkVNgyKzj4+MVIiIisgv2dBARkVNgT4fjY6ODiIicAqNXHB8frxAREZFd1IqejqswnS9AlbtBiktX5QkA5EnAVLkCzgjLShOjqeLbVZM+AXLuBdV2hwvLSq1UaRIw1bZJdUuTtqkuZim3grTdUp6Pq4oy6XxIdavOt1R3Q6HcSyhXTcomXeOWfHNUHU9AvgdUpNwRUh4P6Xz5KsqkSdek3EGWTP4nTSwo3V+qY+4oeTo0mH/d8fGKfdSKRgcREZHEkoYDGx32wccrREREZBfs6SAiIqfAng7Hx54OIiJyCnoLf8oxDbrtsKeDiIicgh7mT+h3Y08H06DbDns6iIiIyC7Y00FERE6BYzocX61odJTCdJeMKp7/nFCv1HnmI5QXKMqkLqR8oVy1vBSHXiyUq3KEqPKeVKVud6FcdcFJF2ORBeVS3gcpb4q03/6KskBhWemYqc63lFtBys0g5ZzIU5RJ95fq/gCAUEWZtN1Sng7V/aPKPQLI14KUG8XPgrp/F8pVy0t5h1T5QwCgvlCuyjFi6v4y91GHudjocHx8vEJERER2USt6OoiIiCTWGkhKtsNGBxEROQU+XnF8fLxCREREdsGeDiIicgp8vOL42OggIiKnwMcrjq9WNDrCYXpqclVYnhRCJoWvqaaAB9Rhr1J4pxSuG6AoyxSWlUIwVdNNS9sthSpK3zJU4YbS+cgRylX7JYVgmrq+yknhoarQ02vCspaET0vHRArXlY6LJaHAnkK5ihTK6yWUq0KcpXMthVdLIbOqdUv3jzQVfJaiTNruEKFcda4B9babCrXnP3K6Wa1odBAREUn4eMXxsdFBREROgY9XHB+jV4iIyClYa5bZqtqzZw/69euHyMhI6HQ6fPbZZxbvAwDs3r0b8fHx8PLyQpMmTbB8+XKj8m7dukGn01X46du3r1XWb0tsdBAREZmhoKAAbdq0wTvvvGO1Ok+ePIk+ffqgS5cuOHToEF566SVMmDABGzduNLxn06ZNyMzMNPwcPXoUrq6uGDx4sNW2w1b4eIWIiJyCvR+v9O7dG7179zZZXlRUhFdeeQVr165Fbm4u4uLiMG/ePHTr1s3kMsuXL0dUVBQWL14MAGjZsiVSUlIwf/58DBo0CAAQFGQ8BHz9+vXw8fFho8MSmnb9ElBFVKhGikvREFK5dBGqlpdGoVtSLnUDWlK3dEykuqVBXKr6LV23JfsldflZsm5psjipXHVMpWULhXIpska1vBTpJG2b6jqWjrcl97YtrzNL1y3d26rlpeNt6fky5/Puxtdv/Ey3FQ2ONTZj9OjROHXqFNavX4/IyEhs3rwZvXr1wo8//ohmzZpVukxycjJ69Ohh9FrPnj2xatUqFBcXw929YnziqlWrMHToUNSpI8VF1jyHbXRcuXI9CGtuDW4HERFZ7sqVK/D1lea6dQz5+fkoLb3eZPL09ISnZ/WCwE+cOIF169bh9OnTiIyMBABMmTIF27Ztw+rVq/Haa69VulxWVhbCw8ONXgsPD0dJSQlycnJQr149o7L//ve/OHr0KFatWlWt7aspHNNBRER0g8jISPj7+xt+5s6t/lffgwcPQtM0NG/eHL6+voaf3bt348SJEwBg9PrYsWMNy+p0xv2b5b1EN78OlPVyxMXF4a677qr2NtYEh+3pCAkJwdmzZwEAPj4+lR7smpCfn4/IyEhkZGSgbl0p/diti8epanicqobHqWoc7ThpmmbotQ4JkdKTmcfHxweXL0tp16qusLAQrq6uRv9zqtvLAQB6vR6urq5ITU2Fq6txurvyHp+0tDTDa+XnKyIiAllZxmngsrOz4ebmhuDgYKPXr1y5gvXr12P27NnV3r6a4rCNDhcXF4SFhdX0ZlRQ3uVWp06dWvH8rKbwOFUNj1PV8DhVjSMeJ1s/UtHpdFbdV2vV1a5dO5SWliI7OxtdunSp9D0xMTEVXuvYsSO++OILo9e2b9+OhISECuM5PvnkExQWFuKxxx6zyjbbAx+vEBERmeHy5ctIS0sz9FicPHkSaWlpSE9PR/PmzTF8+HCMHDkSmzZtwsmTJ/HDDz9g3rx52Lp1q8k6x44diz/++ANJSUk4duwYPvjgA6xatQpTpkyp8N5Vq1ZhwIABFXpAHJnD9nQQERE5spSUFHTv3t3wd1JSEgBg1KhRWLNmDVavXo05c+bg2WefxZkzZxAcHIyOHTuiT58+JuuMjo7G1q1bMXnyZCxduhSRkZFYsmSJIVy23K+//oq9e/di+/btttk5G2Gjo5o8PT0xY8YMs57x3Up4nKqGx6lqeJyqhsfJvrp166YMBXZ3d8esWbMwa9asatXbtWtXHDx4UPme5s2b2yUM2dp0Wm3caiIiIqp1OKaDiIiI7IKNDiIiIrILNjqIiIjILtjoICIiIrtgo8MCp06dwt///ndER0fD29sbTZs2xYwZM1BUJE2tdGt59dVX0alTJ/j4+CAgIKCmN8dhvPvuu4iOjoaXlxfi4+Px3Xff1fQmOZw9e/agX79+iIyMhE6nw2effVbTm+Rw5s6dizvvvBN+fn4ICwvDgAED8Msvv9T0ZhFVio0OC/z888/Q6/V477338L///Q+LFi3C8uXL8dJLL9X0pjmUoqIiDB48GE8//XRNb4rD2LBhAyZNmoSXX34Zhw4dQpcuXdC7d2+kp6fX9KY5lIKCArRp0wbvvPNOTW+Kw9q9ezcSExOxf/9+7NixAyUlJejRowcKCgpqetOIKmDIrJW9+eabWLZsGX7//fea3hSHs2bNGkyaNAm5ubk1vSk1rn379rjjjjuwbNkyw2stW7bEgAEDzJpc6lag0+mwefNmDBgwoKY3xaGdO3cOYWFh2L17N+65556a3hwiI+zpsLK8vDwEBQXV9GaQAysqKkJqaip69Ohh9HqPHj3w/fff19BWkbPIy8sDAH4OkUNio8OKTpw4gbfffttoimKim+Xk5KC0tBTh4eFGr4eHh1eYXZKoOjRNQ1JSEu6++27ExcXV9OYQVcBGRyVmzpwJnU6n/ElJSTFaJiMjA7169cLgwYPx5JNP1tCW2485x4iM3Th1NlD2D+Pm14iqY/z48Thy5AjWrVtX05tCVCnOvVKJ8ePHY+jQocr3NG7c2PB7RkYGunfvjo4dO2LFihU23jrHUN1jRNeFhITA1dW1Qq9GdnZ2hd4Poqp65plnsGXLFuzZswcNGjSo6c0hqhQbHZUICQlBSEhIld575swZdO/eHfHx8Vi9ejVcXG6NzqPqHCMy5uHhgfj4eOzYsQMDBw40vL5jxw7079+/BreMaiNN0/DMM89g8+bN+PbbbxEdHV3Tm0RkEhsdFsjIyEC3bt0QFRWF+fPn49y5c4ayiIiIGtwyx5Keno4LFy4gPT0dpaWlSEtLAwDExMTA19e3ZjeuhiQlJWHEiBFISEgw9JClp6dzPNBNLl++jN9++83w98mTJ5GWloagoCBERUXV4JY5jsTERHz88cf4/PPP4efnZ+hB8/f3h7e3dw1vHdFNNDLb6tWrNQCV/tB1o0aNqvQY7dq1q6Y3rUYtXbpUa9Sokebh4aHdcccd2u7du2t6kxzOrl27Kr12Ro0aVdOb5jBMfQatXr26pjeNqALm6SAiIiK7uDUGIBAREVGNY6ODiIiI7IKNDiIiIrILNjqIiIjILtjoICIiIrtgo4OIiIjsgo0OIiIisgs2OoiIiMgu2OggIiIiu2Cjg8iGiouLa3oTiIgcBhsdRNWwbds23H333QgICEBwcDAeeOABnDhxAgBw6tQp6HQ6fPLJJ+jWrRu8vLzw0UcfAQA++OAD3H777fD09ES9evUwfvx4Q50zZ85EVFQUPD09ERkZiQkTJhjKioqK8Pzzz6N+/fqoU6cO2rdvj2+//dZom/bt24euXbvCx8cHgYGB6NmzJy5evGj7g0FEVE1sdBBVQ0FBAZKSkvDDDz/gm2++gYuLCwYOHAi9Xm94zwsvvIAJEybg2LFj6NmzJ5YtW4bExEQ89dRT+PHHH7FlyxbExMQAAD799FMsWrQI7733Ho4fP47PPvsMrVq1MtQ1evRo7Nu3D+vXr8eRI0cwePBg9OrVC8ePHwcApKWl4b777sPtt9+O5ORk7N27F/369UNpaal9DwwRURVwwjciC5w7dw5hYWH48ccf4evri+joaCxevBgTJ040vKd+/foYPXo05syZU2H5hQsX4r333sPRo0fh7u5uVHbixAk0a9YMp0+fRmRkpOH1v/3tb7jrrrvw2muvYdiwYUhPT8fevXttt5NERFbCng6iajhx4gSGDRuGJk2aoG7duoiOjgYApKenG96TkJBg+D07OxsZGRm47777Kq1v8ODBuHr1Kpo0aYIxY8Zg8+bNKCkpAQAcPHgQmqahefPm8PX1Nfzs3r3b8EinvKeDiKg2cKvpDSCqTfr164eGDRti5cqViIyMhF6vR1xcHIqKigzvqVOnjuF3b29vZX0NGzbEL7/8gh07duDrr7/GuHHj8Oabb2L37t3Q6/VwdXVFamoqXF1djZbz9fWtUv1ERI6EPR1EVXT+/HkcO3YMr7zyCu677z60bNlSHLDp5+eHxo0b45tvvjH5Hm9vbzz44INYsmQJvv32WyQnJ+PHH39Eu3btUFpaiuzsbMTExBj9REREAABat26trJuIyJGwp4OoigIDAxEcHIwVK1agXr16SE9Px4svviguN3PmTIwdOxZhYWHo3bs3Ll26hH379uGZZ57BmjVrUFpaivbt28PHxwf/+te/4O3tjUaNGiE4OBjDhw/HyJEjsWDBArRr1w45OTnYuXMnWrVqhT59+mDq1Klo1aoVxo0bh7Fjx8LDwwO7du3C4MGDERISYoejQkRUdezpIKoiFxcXrF+/HqmpqYiLi8PkyZPx5ptvisuNGjUKixcvxrvvvovbb78dDzzwgCH6JCAgACtXrkTnzp0NvRZffPEFgoODAQCrV6/GyJEj8eyzz6JFixZ48MEHceDAATRs2BAA0Lx5c2zfvh2HDx/GXXfdhY4dO+Lzzz+Hmxu/TxCR42H0ChEREdkFezqIiIjILtjoICIiIrtgo4OIiIjsgo0OIiIisgs2OoiIiMgu2OggIiIiu2Cjg4iIiOyCjQ4iIiKyCzY6iIiIyC7Y6CAiIiK7YKODiIiI7IKNDiIiIrKL/w8+s0fXFR1EkwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "index = 20\n", + "webbpsf.display_psf(cube, ext=3, cube_slice=index, \n", + " # Note that currently the default plot title isn't very informative for datacube modes\n", + " # so we can specify a better title directly here:\n", + " title=f'NRS IFU cube slice {index}, $\\\\lambda$={cube[0].header[\"WAVELN\"+str(index)]*1e6:.4} micron') " + ] + }, + { + "cell_type": "markdown", + "id": "03938bc2-f831-4dbd-9188-f9c4ef3dc33c", + "metadata": {}, + "source": [ + "## MIRI MRS example\n", + "\n", + "This is mostly the same, except the selection of the IFU bands is done via setting `band` directly. \n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "b9196827-456e-4f09-b40a-831550af07d7", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Band is 2A\n" + ] + } + ], + "source": [ + "miri = webbpsf.MIRI()\n", + "miri.mode = 'IFU'\n", + "miri.band= '2A' \n", + "print(\"Band is\", miri.band)" + ] + }, + { + "cell_type": "markdown", + "id": "b7612fcf-5793-497a-a4f5-38843792d31b", + "metadata": {}, + "source": [ + "Note that selecting an IFU band automatically configures the simulated pixelscale to match the default scale used in pipeline output datacubes: " + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "2bc6face-54c4-42ad-b4a5-d80de8d120fb", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Pixelscale for band 2A is 0.17 arcsec/pix\n", + "Pixelscale for band 3B is 0.2 arcsec/pix\n" + ] + } + ], + "source": [ + "print(f\"Pixelscale for band {miri.band} is {miri.pixelscale} arcsec/pix\")\n", + "\n", + "miri.band = '3B'\n", + "print(f\"Pixelscale for band {miri.band} is {miri.pixelscale} arcsec/pix\")" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "4e131678-a7b0-476f-b458-5ad3e731ecf5", + "metadata": {}, + "outputs": [], + "source": [ + "# let's get a subset for a faster PSF sim runtime\n", + "waves = miri.get_IFU_wavelengths(nlambda=50)\n", + "\n", + "# convert waves from microns to meters\n", + "# (this is a work around for a current issue with the poppy library upstream)\n", + "waves = waves.to_value(u.meter)" + ] + }, + { + "cell_type": "markdown", + "id": "4d747fa7-1851-4176-b1f7-8222f550a866", + "metadata": {}, + "source": [ + "Compute a datacube:\n", + "\n", + "(Note, for MIRI MRS you may see a warning message about the valid region of the field dependence model; this is a benign warning and can be mostly ignored.)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "e18451d5-5144-4e4f-b643-c811840362e5", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/mperrin/Dropbox (Personal)/Documents/software/webbpsf/webbpsf/opds.py:1759: UserWarning: For (V2,V3) = [-8.40143167 -5.31995333] arcmin, Field point -8.254199999999999 arcmin, -2.4800466666666674 arcmin not within valid region for field dependence model of OTE WFE for MIRI: -8.254199999999999 arcmin--6.21738 arcmin, -2.557224 arcmin--0.5632056 arcmin. Clipping to closest available valid location, 0.14723166666666643 arcmin away from the requested coordinates.\n", + " warnings.warn(warning_message)\n" + ] + } + ], + "source": [ + "cube = miri.calc_datacube(waves)" + ] + }, + { + "cell_type": "markdown", + "id": "695ff9ec-44f3-481c-8752-f6a7220cf77d", + "metadata": {}, + "source": [ + "The output FITS file has the same extensions as a regular PSF calculation, but each extension contains a 3D datacube rather than a 2D image: " + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "5986841b-97c2-443d-9677-ac48adc3fdc8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Filename: (No file associated with this HDUList)\n", + "No. Name Ver Type Cards Dimensions Format\n", + " 0 OVERSAMP 1 PrimaryHDU 990 (240, 240, 50) float64 \n", + " 1 DET_SAMP 1 ImageHDU 992 (60, 60, 50) float64 \n", + " 2 OVERDIST 1 ImageHDU 1097 (240, 240, 50) float64 \n", + " 3 DET_DIST 1 ImageHDU 1098 (60, 60, 50) float64 \n" + ] + } + ], + "source": [ + "cube.info()" + ] + }, + { + "cell_type": "markdown", + "id": "e6436fd6-e24e-4f01-af0a-c31f74c215fa", + "metadata": {}, + "source": [ + "The display_psf function works with datacubes, but you have to specify which slice of the cube to display. " + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "c73ffde0-7f85-49fc-9552-d0daabd124f2", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh0AAAHHCAYAAAAbLeozAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACDmElEQVR4nO2deXRUVdbFd2UGEggkTAECQUCJgGCCTCLQtkxKIyIyTyr90QYZ4oA4ATaCAygikyBC0yqgDSi2QEMLgiLSECZR26lBlABhCoEACUne9wempFK1X/IqVZWk2L+1aq3k1Kv77r3vvlen7j37HpthGAaEEEIIIbxMQElXQAghhBDXBnI6hBBCCOET5HQIIYQQwifI6RBCCCGET5DTIYQQQgifIKdDCCGEED5BTocQQgghfIKcDiGEEEL4BDkdQgghhPAJcjqEEEII4RPkdAghhBDCJ3jM6ViyZAlsNhtsNhs+/fRTp/cNw0CDBg1gs9nQsWNHp8/t2rXLZVk2mw1BQUGoWbMm+vXrhx9++MHlea/+vLfqWJT65LNjxw706tULsbGxCA0NRfXq1dGmTRs88sgjpnWcNGkSbDYbTp48WWh7fEl+vQpj7969uPPOOxEbG4ty5cqhSpUqaNOmDd5++22H4wr2qc1mQ9WqVdGxY0f885//9FYzLFOarkd+nx06dMjl/75i06ZNuP/++3HDDTegQoUKqFWrFnr27ImUlBSnY8+fP4+xY8ciJiYGYWFhaN68OZYvX17sOly6dAlBQUGIiIjAww8/XOzy8jl37hwef/xxdO7cGVWrVoXNZsOkSZOK9Nk333wTNpsN4eHhhR5rpQ8//fRTp3sl//Xll19abaJlSmqcCf/E4zMdERERWLRokZN9y5Yt+OmnnxAREVHkshYvXozt27fj3//+N0aNGoU1a9bg1ltvxZkzZ0qkjkWtz8cff4y2bdsiIyMDL730EjZs2IDXXnsN7dq1w4oVK4pV99JOeno66tSpg6lTp2Lt2rVYunQp6tWrh8GDB2PKlClOx+f36RdffIEFCxYgMDAQPXr0wEcffVQCtS9b3Hnnndi+fTtq1qzp0/POmzcPhw4dwpgxY7B27Vq89tprSEtLQ+vWrbFp0yaHY++55x787W9/w8SJE7Fu3Tq0bNkS/fv3x7vvvlusOthsNmzatAm33HILZs+eje+//75Y5eVz6tQpLFiwAFlZWbj77ruL/LkjR47g0UcfRUxMTJGOt9KH+UydOhXbt293eDVp0qTIdXSXkhpnwk8xPMTixYsNAMaDDz5olCtXzjh79qzD+4MGDTLatGlj3HjjjUaHDh2cPrdz505Tm2EYxuTJkw0AxltvvVXosd6uI6uPYRjGbbfdZlx33XXG5cuXneqQm5trWseJEycaAIwTJ04U2h5fkl8vd2nVqpVRp04d+/+sTy9cuGCEhoYa/fv3d/tcnqQ0XY/8Pjt48GCJ1uP48eNOtnPnzhnVq1c3br/9drvt448/NgAY7777rsOxd9xxhxETE2Pk5OQUuy7r1693eQ53ycvLM/Ly8gzDMIwTJ04YAIyJEycW+rm77rrL6NGjhzF06FCjQoUKhR5f1D40DMPYvHmzAcB4//33i9aIUkxmZmZJV0GUMB6f6ejfvz8AYNmyZXbb2bNnsXLlStx///3FKjsxMREAcPz48WKV46k6svqcOnUK0dHRCAoKcvpMQEDRuvyXX37BPffcg4oVK6JSpUoYNGgQTpw44XDMjz/+iOHDh6Nhw4YoX748atWqhR49euCrr75yOC5/ieDrr79G//79UalSJVSvXh33338/zp4963Tujz/+GM2bN0doaCji4uIwffr0ItXZDNYfBQkLC0NISAiCg4PdPtd///tf9O/fH9WrV0doaChiY2MxZMgQZGVlAQCGDRuGevXqOX3ObAmpKNcDAH744QcMGDAA1apVQ2hoKBo3bow5c+YUWucTJ07gz3/+M+rUqYPQ0FBUrVoV7dq1w7///W/6GTbtXVj7i1NPAKhWrZqTLTw8HPHx8fjll1/sttWrVyM8PBx9+vRxOHb48OFITU3Fjh07inQ+M2644QYAV5b1PEH+soUV3n77bWzZsgVz584t8meK2oeeIH9c79+/H3369EGlSpVQpUoVJCcnIycnB9999x26du2KiIgI1KtXDy+99JLD512Ns6KMsfzz7t69G/feey8qV66M6667DgDw+eef4/bbb0dERATKly+Ptm3b4uOPP6Z1L+qzS5R+PO50VKxYEffeey/eeustu23ZsmUICAhA3759i1X2wYMHAQCNGjUqFXVk9WnTpg127NiB0aNHY8eOHbh8+bLlOvbq1QsNGjTAP/7xD0yaNAkffPABunTp4lBWamoqoqKi8MILL2D9+vWYM2cOgoKC0KpVK3z33XdOZfbu3RuNGjXCypUr8cQTT+Ddd9/FuHHjHI755JNP0LNnT0RERGD58uV4+eWX8d5772Hx4sWW6p+Xl4ecnBycOHECc+fOxb/+9S+MHz/e6bjc3Fzk5OTg8uXL+PXXXzF27FhkZmZiwIABls6Xz759+9CyZUt8+eWXeO6557Bu3TpMmzYNWVlZyM7OdqtMoGjX45tvvkHLli1x4MABzJgxA//85z9x5513YvTo0Zg8ebJp+YMHD8YHH3yAZ599Fhs2bMCbb76JP/7xjzh16pSlehal/cWpJ+Ps2bPYvXs3brzxRrvtwIEDaNy4sZOz2axZM/v7xeXFF18E4Ox0GIaBnJycIr2KQ1paGsaOHYsXXngBtWvXLlZZrvrwapKSkhAUFISKFSuiS5cu+Pzzzy2Vf9999+Gmm27CypUrMWLECLz66qsYN24c7r77btx5551YvXo1/vCHP2D8+PFYtWoVLcfqPXbPPfegQYMGeP/99zF//nxs2bIFf/jDH3D27FksWrQIy5YtQ0REBHr06EGXn4vy7BJlBE9NmVw9XZ4/HXjgwAHDMAyjZcuWxrBhwwzDMCwtr3z55ZfG5cuXjXPnzhnr1683atSoYdx2220OyxbuLK+4W8ei1McwDOPkyZPGrbfeagAwABjBwcFG27ZtjWnTphnnzp0zrWP+dP64ceMc7O+8844BwHj77bfpZ3Nycozs7GyjYcOGDp/PL/Oll15yOP6hhx4ywsLC7NPJhnFlGSQmJsa4ePGi3ZaRkWFUqVLF0vLK//3f/9nbHxISYsydO9fh/fw+LfgKDQ11OtYKf/jDH4zIyEgjLS2NHjN06FCjbt26TnZXS0hWrkeXLl2M2rVrOy3bjRo1yggLCzNOnz5N6xQeHm6MHTvWrGlOyyuulluK0v7i1JMxcOBAIygoyNi1a5fd1rBhQ6NLly5Ox6amphoAjKlTp1o+z9X861//MgAYlStXNqpWrerwXv79XZQXW64qyvJK7969jbZt29rvoaIur7jCVR8ahmHs3r3bGDNmjLF69Wpj69atxltvvWU0btzYCAwMNNavX19oufljeMaMGQ725s2bGwCMVatW2W2XL182qlatatxzzz12W8FxVpQxdvV5n332WQd769atjWrVqjk8C3NycowmTZoYtWvXdngeWXl2ibKBVySzHTp0wHXXXYe33noLX331FXbu3OnW0krr1q0RHByMiIgIdO3aFZUrV8aHH35YpGl6b9SxqPWJiorCZ599hp07d+KFF15Az5498f3332PChAlo2rRpkZQQAwcOdPj/vvvuQ1BQEDZv3my35eTkYOrUqYiPj0dISAiCgoIQEhKCH374Ad9++61TmX/6058c/m/WrBkuXbqEtLQ0AEBmZiZ27tyJe+65B2FhYfbj8n+FWOHJJ5/Ezp078fHHH+P+++/HqFGjXC7TLF26FDt37sTOnTuxbt06DB06FElJSZg9e7al8wHAhQsXsGXLFtx3332oWrWq5c+bUdj1uHTpEj755BP06tUL5cuXd/gl3b17d1y6dMlUaXDLLbdgyZIlmDJlCr788ku3ZseK0v7i1tMVzzzzDN555x28+uqrSEhIcHjPbKnC6jLG1aSnp+P+++9Hz5498dBDD+HEiRNITU21v5+QkGAfV4W9ihr8WZCVK1fio48+wsKFC4vVFsC8D1u0aIGZM2fi7rvvRvv27TF8+HB88cUXqFmzJh5//PEin+Ouu+5y+L9x48aw2Wzo1q2b3RYUFIQGDRrg559/dlmGO/dY79697X9nZmZix44duPfeex1UPoGBgRg8eDB+/fVXl7O0hT27RNmh+N/eLrDZbBg+fDhmzZqFS5cuoVGjRmjfvr3lcpYuXYrGjRvj3LlzWLFiBd544w30798f69atK5E6Wq1PYmKiPe7j8uXLGD9+PF599VW89NJLTuumBalRo4bD/0FBQYiKinKYbk9OTsacOXMwfvx4dOjQAZUrV0ZAQAAefPBBXLx40anMqKgoh/9DQ0MBwH7smTNnkJeX53RuV/UpjNjYWMTGxgIAunfvDgCYMGEChg4d6vCwaty4sb2PAKBr1674+eef8fjjj2PQoEGIjIws8jnPnDmD3NzcYk9zu6Kw63Hq1Cnk5OTg9ddfx+uvv+6yDDNnc8WKFZgyZQrefPNNPPPMMwgPD0evXr3w0ksvFbnvi9L+4tazIJMnT8aUKVPw/PPPY9SoUQ7vFRyv+Zw+fRoAUKVKlSKfpyBJSUm4fPkyFi5ciK1btwK4ssSS70CEh4ejefPmRSrLnR8x58+fR1JSEh5++GHExMQgPT0dAOzLC+np6QgODkaFChUKLcusDxmRkZG46667MH/+fFy8eBHlypUr9DMF+zskJATly5d3+IGRb8/IyHBZhjv32NWqlzNnzsAwDJdKmPxr52rMFPbsEmUHr20ONmzYMJw8eRLz58/H8OHD3Soj/wupU6dOmD9/Ph588EGsX78e//jHP0qkjsWpT3BwMCZOnAigaGvZx44dc/g/JycHp06dcrj53n77bQwZMgRTp05Fly5dcMsttyAxMdHtPSUqV64Mm83mdG5X9bHKLbfcgpycHPzvf/8r9NhmzZrh4sWLlmWQVapUQWBgIH799VfT48LCwhwC3vIx67fCrkflypURGBiIYcOG0V/U+c6XK6KjozFz5kwcOnQIP//8M6ZNm4ZVq1Zh2LBhpm25mqK0v7j1vJrJkydj0qRJmDRpEp588kmn95s2bYpvv/3WKW4iP9DZXbnnP/7xD7z77rtYtGgRqlatiptvvhmAY1zHli1bEBwcXKSXO/tPnDx5EsePH8eMGTNQuXJl+2vZsmXIzMxE5cqVnWbHXFFYH5phGAaA4s0YWaWo99jVXF2//B9GR48edTouf6YqOjq6+BUVpRavzHQAQK1atfDYY4/hv//9L4YOHeqRMl966SWsXLkSzz77LO65554iK0G8VUdWn6NHj7r05POXPIoynfvOO+84TLO+9957yMnJcdi0zGaz2T3+fD7++GMcOXIEDRo0sNyeChUq4JZbbsGqVavw8ssv238BnTt3rtj7ZmzevBkBAQGoX79+ocfmf3lYXSIpV64cOnTogPfffx/PP/88fXjVq1cPaWlpOH78OKpXrw7gyi/Uf/3rX7Tswq5H+fLl0alTJ+zZswfNmjVDSEiIpbpfTWxsLEaNGoVPPvkE27ZtK/LnitJ+T9Xzr3/9KyZNmoSnn37a7kwXpFevXli4cCFWrlzpEKD9t7/9DTExMWjVqpXl8x47dgwjR47En//8Z/tyQVxcHCIjI7Fnzx77cfnLK0XBneWVGjVqOCx15vPCCy9gy5YtWLduXaFfnkXpQ8aZM2fwz3/+E82bN3eaqfAmRb3HGBUqVECrVq2watUqTJ8+3T5Dk5eXh7fffhu1a9cutlBAlG685nQAV25AT1K5cmVMmDABjz/+ON59910MGjSo2GUWp46sPl26dEHt2rXRo0cP3HDDDcjLy8PevXsxY8YMhIeHY8yYMYWWvWrVKgQFBeGOO+7A119/jWeeeQY33XQT7rvvPvsxd911F5YsWYIbbrgBzZo1Q0pKCl5++eViLS/89a9/RdeuXXHHHXfgkUceQW5uLl588UVUqFDBPi1uxp///GdUrFgRt9xyC6pXr46TJ0/i/fffx4oVK/DYY485ORIHDhyw/xI+deoUVq1ahY0bN6JXr16Ii4tzONZms6FDhw4ud5PN55VXXsGtt96KVq1a4YknnkCDBg1w/PhxrFmzBm+88QYiIiLQt29fPPvss+jXrx8ee+wxXLp0CbNmzUJubi4ttyjX47XXXsOtt96K9u3b4y9/+Qvq1auHc+fO4ccff8RHH31EN306e/YsOnXqhAEDBuCGG25AREQEdu7cifXr1+Oee+4prMstt9/deuYzY8YMPPvss+jatSvuvPNOpxiQ1q1bAwC6deuGO+64A3/5y1+QkZGBBg0aYNmyZVi/fj3efvttBAYGOnyuKNd3xIgRqFy5Ml555RUHe4sWLRxmOiIiIhyW7aywbt06ZGZm4ty5cwCuqH3yZzO7d+9uX5K4+gdAPkuWLEFgYKDTe1u2bMHtt9+OZ599Fs8++2yR+xAABgwYgNjYWCQmJiI6Oho//PADZsyYgePHj2PJkiVutbE4FGWMmTFt2jTccccd6NSpEx599FGEhIRg7ty5OHDgAJYtW+bTmRtRAngqIrWoKpLibA5mGIZx8eJFIzY21mjYsKGRk5PjtnqluHVk9TEMw1ixYoUxYMAAo2HDhkZ4eLgRHBxsxMbGGoMHDza++eYb03PnR2unpKQYPXr0MMLDw42IiAijf//+ThsKnTlzxnjggQeMatWqGeXLlzduvfVW47PPPjM6dOjgUH+2wRXbbGrNmjVGs2bNjJCQECM2NtZ44YUXirw52FtvvWW0b9/eiI6ONoKCgozIyEijQ4cOxt///neX5776ValSJaN58+bGK6+8Yly6dMnh+HPnzhkAjH79+hVah2+++cbo06ePERUVZW/DsGHDHMpcu3at0bx5c6NcuXJG/fr1jdmzZ5uqV4pyPQzDMA4ePGjcf//9Rq1atYzg4GCjatWqRtu2bY0pU6bQ+l66dMkYOXKk0axZM6NixYpGuXLljOuvv96YOHGiw2ZKRVGvFLX97tQznw4dOpiqQa7m3LlzxujRo40aNWoYISEhRrNmzYxly5Y5lVmU67tw4UIjMDDQ2L59u9N7ycnJhs1mMzIyMgqtf2HUrVvXstIlH6ZeyVfT5CthrPThtGnTjObNmxuVKlUyAgMDjapVqxq9evUy/vOf/xSpPez+Z3Xt0KGDceONN9r/dzXOijLGzDbW++yzz4w//OEPRoUKFYxy5coZrVu3Nj766KMi1720bJQnrGMzjN8WBoUoxaxduxZ33XUX9u3bh6ZNm5Z0dYSH0fUV4tpAWWZFmWDz5s3o16+fvpD8FF1fIa4NNNMhhBBCCJ+gmQ4hhBBC+ASvOx1HjhzBoEGDEBUVhfLly6N58+ZISUnx9mmFEEIInzB37lzExcUhLCwMCQkJ+Oyzz+ixR48exYABA3D99dcjICAAY8eO9V1FSwFedTrOnDmDdu3aITg4GOvWrcM333yDGTNmWNplUgghhCitrFixAmPHjsVTTz2FPXv2oH379ujWrRsOHz7s8visrCxUrVoVTz31FG666SYf17bk8WpMxxNPPIFt27aZen2MvLw8+w6R5cuXl3ZbCCHKGIZh4MKFCwCu7DRa3A0dCzuHp7DyndOqVSvcfPPNmDdvnt3WuHFj3H333Zg2bZrpZzt27IjmzZtj5syZxalumcKrm4OtWbMGXbp0QZ8+fbBlyxbUqlULDz30EEaMGOHy+KysLPv21CdOnHBrV00hhBClj+PHj6NatWoeL/fChQsOyeM8QWpqqkPenNDQUKfdn4ErOxmnpKTgiSeecLB37twZX3zxhUfr5C94dXnlf//7H+bNm4eGDRviX//6F0aOHInRo0dj6dKlLo+fNm0aKlWqhEqVKsnhEEIIUSLExMTYv4sqVapEZyxOnjyJ3NxcezqFfKpXr17sfFX+ildnOvLy8pCYmIipU6cCuLJV8ddff4158+ZhyJAhTsdPmDABycnJAK5kcaxVqxYAoA2AwALH5pFzMj+6OrEDQOH5GR1hScfZBJ/rfI2AWVLmE8R+ntgvEXsmsbuTm5H9lmB9znKjsmwN5Yk9h9gBIJ3Ys4mdDfiKxM7GjVnuUHbusxaPZ9eIjQ1WPsD7lrWPXVPW7oL3Zz5m1461zzkVn3lZbPN6dq+YfRWwNGbs/mJZT5idTdiz/gOAYGJn7WZ1Zf3N+ok95wD+LCg4zgz8Pl7Ll2ej0HMcP368SJl9XZGZmWl3HlzNdJhRcCnGMAyFBBC86nTUrFkT8fHxDrbGjRtj5cqVLo+/egrr6rwMgXC+KdnlZA1iNy4AuJ+aq2jnYHUye9CwKShmZ/3hyWHPymJ1Yu1z5xoxWFnsgcyOt1onszHDgqRYWcyBtjpuzKYtPXUtrI5xs/HHnAjWHwyrzwKz+87qfWT1nrB6vNl7bJz54hlR1HZcfS198SVcoVwoKpQzdxAoeb+PyIoVKxbJeYmOjkZgYKDTrEZaWprT7Ie4gleXV9q1a4fvvvvOwfb999+jbt263jytEEKIa5G8nOK9fiMhIQHx8fGYM2eO6elCQkKQkJCAjRs3Otg3btyItm3beqWJZR2vznSMGzcObdu2xdSpU3HffffhP//5DxYsWIAFCxZ487RCCCGuRQo4D5Y/+xspKSlFXqZJTk7G4MGDkZiYiDZt2mDBggU4fPgwRo4cCeBK2MCRI0ccYhnzMyKfP38eJ06cwN69exESEuK0MuCPeNXpaNmyJVavXo0JEybgueeeQ1xcHGbOnImBAwd687RCCCGET+jbty9OnTqF5557DkePHkWTJk2wdu1a+4z+0aNHnfbsaNGihf3vlJQUvPvuu6hbty4OHTrky6qXCKU290pmZqZdBvVHOHtHbE29FrFHmZyLrVOzwCsWnMkCspidBcwBPLiQBXexsk4TO6tTJKsQgKrEzoIwIywezwJ63QkkZe1mg52t81cmdrPfQFZjFdhYZuWwa222VsrKYqvs7FqwurLjWUAlwPuDjU0WKMuOTyf2VFYhmAd3u8JqjBCzm8UIsQgFdr2tBrqzwFN35gtcBfwf/e3v8+fPux3kacbV3xXnT/5arEDS8OjaAGDfLTQpKQlJSUkeq6vw8kyHEEII4TNKYHlFWEMJ34QQQgjhEzTTIYQQwj/Iyy3GTAcT2gtPopkOIYQQ/oGPJbPCOprpEEIIIa5CMR3eo0w4Hdlw3mGSKSNYtLzZLoQs2tvqrpZWd4o0236dbRjM1AasDUyZ48528ZWI3VP7DLozKcquN+s/tiX9OWJPJ3Z3drhl278zdUIksbPrYKaAYEqsdGJn14LZWflm9x2DjWV2DqbKYCozM9UYGzdWdy1mdWVtM1P5sD5kSiyr186q0gbgfX6mwP8+l0Z6KJBUeI8y4XQIIYQQhSKno9SjmA4hhBDiKhTT4T000yGEEMI/0D4dpR45HUIIIfwDLa+UeuR0CCGE8A+MYuzTYWifDl+gmA4hhBBC+IQyMdNRHs4VDSfHMumZmXfFkqsxqZqZXNIVTIbH5HMAl/tZ9RLZqiSTAFYxKYu1g/Ufkw0y2O8MM4kjg40DJuFjSa+OEbvZGKhL7EwizfqVSWOZ3axOrD9Y0jWWMM/qtWYSZTOYfJn9fmX3KRv7cSbnZvcFazdr3yliZ2PZTDLL3mNjltlZXVn/mUnhmWTWnXvVo3hoeSUhIUEJ37xEmXA6hBBCiEJRIGmpR8srQgghhPAJmukQQgjhH0i9UuqR0yGEEMI/kNNR6tHyihBCCCF8QpmY6YiGc2S+1YRvZg21qmhg56hs8dwsAhzgSg5WJ6ZCYNH4VhNYAbxO7PcB6yem4mBqHqbkAXi0PDsHU3iwurLjmYIE4AogNmZZwjd2TZlSxOx3Grt2TLnA+o+NWaY4KZgA7GrYWGb3C1OseSphHsDvCzbOThC7VZUP6wuA/zJk52B1ZfeRO7tTsHFT8NoZbpbvNlKvlHrKhNMhhBBCFEpeMTYHy/vdPZJ6xXtoeUUIIYQQPkEzHUIIIfwDBZKWeuR0CCGE8A/kdJR65HQIIYTwD+R0lHrKhNNRDs4qAhaxzmBKADNY1DU7N4uuZ4EzZrlXWIQ9O/dZYmfR+Oz2MsuXwqLiWTtYfzCs5tAwOzfLM2FVFWRVgWMGUxswWH+w8WQ2xpnqhPUfezCwPDGsridpjbhqrBqxsz5n9wSrq9mzg73H+oP1ORtnzM5UMABXALExy5REbPyxZ0QNWiP+XsGxmQtgp0k54tpDgaRCCCH8g/yZDndfv5GQkID4+HjMmTOnBBvjn5SJmQ4hhBCiUJTwrdSjmQ4hhBBC+ATNdAghhPAPjGJsDmb4dO/UaxY5HUIIIfyDvBwgj4U1F+GzwuuUCacjBM7qFTasmMLCbB2JRfDbiJ1F3bPocBbhblanKGJn7WO3C1NxsDabKWqs9gdTDzA7W0FlbQZ4bhQWwc8i/q0+bsyUIukWz12R2FmuFqY2MHvUMlXSeWJn/cHawOxsbAC8fZHEzsYNw51nAbuubOyzaxFN7OyBy8YxABwidva7nOX+sZqfKZZVCEBNYi/YvhxIvSIcKRNOhxBCCFEomuko9cjpEEII4R/I6Sj1SL0ihBBCCJ+gmQ4hhBD+gWY6Sj1yOoQQQvgHecWQzOZJMusLyoTTEQjnyHxP5UIArK8xWVVGsOhws9wkLOcHaweLrk8n9kxiN8spwtrB6srUBqzdTL1iNkiZOqI6sTOF0VFiP2ZybgbrJ6vHM0WDO2OcqZLYOGD5PtjxLMeKWV+w68ruI9Zu9hXDxjJTCwG8z1n/WVVusWeN2dckaze7X9i52dhnCiZ3nk8Fn9M+X7/PywHy3DxrgW3QAwICkJSUhKSkJA9VTgBlxOkQQgghfIW2QfcecjqEEEL4Bx6a6RDeQ06HEEII/0BOR6lHklkhhBBC+ATNdAghhPAPNNNR6ikTTkcezHNdXI07w41FrJtFb7uCRbizurN8GADPG8HUA2eJ/RSxXyB2szazsCqmUGC3sFm7XWGWv6MysYdaPDfrD6ZUYrlxAJ47hI1NVleWv4NdI7O8Oew9q4okVg5ThJhdO3ZfsM8wdQ7rD5bbhbXN7D127Zjyg+0UwZ41ZioflkuF2Zlqh6lXThO72bOUlVUQn3+Ny+ko9Wh5RQghhBA+oUzMdAghhBCFYhRjczBDm4P5AjkdQggh/IO8HCCPbZVYhM8Kr6PlFSGEEEL4BM10CCGE8A8001HqkdMhhBDCP5DTUeops06HVbmdmXSPSRaZ9Ix1GpORFVVeVpTPeEoaa9UOANWIvSaxMymj1UFnlqiavceuN3usMCkokzhGsgqBS11ZXdm52aOT2c36iUky2TWtSuxMQpxK7Gby6ErEblWqbjWxoFm4IJOkMwk7u0+tjj8zrErVmZ3Vicmdzb662bOx4Bg0k3F7BTkdpR6fxXRMmzYNNpsNY8eO9dUphRBCCK8zd+5cxMXFISwsDAkJCfjss89Mj9+yZQsSEhIQFhaG+vXrY/78+Q7vf/311+jduzfq1asHm82GmTNnerH2vsUnTsfOnTuxYMECNGvWzBenE0IIcS2Sl1O8lxusWLECY8eOxVNPPYU9e/agffv26NatGw4fPuzy+IMHD6J79+5o37499uzZgyeffBKjR4/GypUr7cdcuHAB9evXxwsvvIAaNWq4Va/SitedjvPnz2PgwIFYuHAhKldm+0cKIYQQxSQvtxhOh3v7dLzyyit44IEH8OCDD6Jx48aYOXMm6tSpg3nz5rk8fv78+YiNjcXMmTPRuHFjPPjgg7j//vsxffp0+zEtW7bEyy+/jH79+iE0lAUAlE287nQkJSXhzjvvxB//+Edvn0oIIYTwGdnZ2UhJSUHnzp0d7J07d8YXX3zh8jPbt293Or5Lly7YtWsXLl/2eRSMz/FqIOny5cuxe/du7Ny5s0jHZ2VlISvrSshWZiYL6RJCCCFckJdT9ERdrj77GxkZGcjN/X3mIzQ01OWMw8mTJ5Gbm4vq1as72KtXr45jx465PM2xY8dcHp+Tk4OTJ0+iZk0Wmu8feM3p+OWXXzBmzBhs2LABYWFmKZZ+Z9q0aZg8ebKT3YBz5LWZGsUqrCw2dq1OwrHWm3U+i/q3GrHOosxZoiozrJ6DJdyKJHbW30w5AAAZxM5+L7DjjxM7q5OZIomND3a9Wb+yMcBWns3GJftMpMlnXMF0Ae4kDYshdqZqsXrfMfWK2XcS6yfXXx9c7cLuCTYuzZ4FVifX2fOMjTO26M2UWwC/FgXtJaNeKcZnfyMmxnF0Tpw4EZMmTaIftdkc7wzDMJxshR3vyu6PeM3pSElJQVpaGhISEuy23NxcbN26FbNnz0ZWVhYCAx0FVhMmTEBycjKAKzMdBS+8EEII4W1SU1NRocLvYmUWVxEdHY3AwECnWY20tDSn2Yx8atSo4fL4oKAgREWZ5a/2D7zmdNx+++346quvHGzDhw/HDTfcgPHjxzs5HIDjFJar94UQQgiKh2Y6Klas6OB0MEJCQpCQkICNGzeiV69edvvGjRvRs2dPl59p06YNPvroIwfbhg0bkJiYiOBgs/kl/8BrTkdERASaNGniYKtQoQKioqKc7EIIIUSx8ZDTkZCQgICAACQlJSEpKcn0Y8nJyRg8eDASExPRpk0bLFiwAIcPH8bIkSMBXJnBP3LkCJYuXQoAGDlyJGbPno3k5GSMGDEC27dvx6JFi7Bs2TJ7mdnZ2fjmm2/sfx85cgR79+5FeHg4GjRo4GYDSwdldkdSIYQQwhukpKQUaaYDAPr27YtTp07hueeew9GjR9GkSROsXbsWdevWBQAcPXrUYc+OuLg4rF27FuPGjcOcOXMQExODWbNmoXfv3vZjUlNT0aJFC/v/06dPx/Tp09GhQwd8+umnnmlkCWEz8iNYShmZmZkID78SCpYE56AmFqjFMGskC7Bi2yuzhR8WXMWC6dwJJGUBbT8TOwuQPGFybkYssV9P7HWJPZLY3QkkZYFqVgNJfyV2VifWNgBwvZJrPSCQwSZgzQItWeArG/uMdGI/ZPG8AB9PJRlIyu4Ldn95KpDU7PnEnjfMXrSw/d9hX62eCiSd89vf58+fL/IXuRWu/q44v/B2VAhz77d05qUchI/45Eo5XqqrKCMzHefhfAOwGF+reQoAfvNanaVjzgU7t1nUCnsIMUUIy6HBHg4sSt/sQWM1lwVzFi5ZLIcdD/A+t3rtWK4R9iVhlqPmDLGz/mNj1qpjbfalzO6Xon555MOuBXOSzcY4G2us3eyasnO7k/OIwfqP3UesTuy+NnsQsz60Op5YG9gzxcx5YeP/ZIH/3dtuqxgYue4vrxi/19bK8oqwRplwOoQQQohCycsB8tycvL9qR1IryyvCGj5L+CaEEEKIaxs5HUIIIfwDDyV8S0hIQHx8PObMmWNyMuEOWl4RQgjhH2h5pdSjmQ4hhBBC+IQyMdORC2fvyGr+CbOGsvfYOVg0OZP6sShzs8h+Fv3OYEoRVg6TA5spJpgKh3mup4j9rMk5XMHqCnD5LetzlkqJKSlY/xWM0r8aNgaZmoKNA3Yt2Hg1y3PBxgf7DDueSUqZdNRMJnzO4mesjj+rcmrAuuKFlcXKYcebPQuYsqoKsVckdvY8Y5jViY3BgmPf5+oVD810CO+hmQ4hhBD+gWI6Sj1lYqZDCCGE8BWK6fAecjqEEEL4B3m5xVhecXdXMWEFOR1CCCH8g7wcIM/NqAE5HT5BMR1CCCGE8AllYqajBpwj11nFWdS4O0nDGCy6nqkNmJKC5UIAeNS3VaVDDWJnbWB1BXifM+UCS07HjmfKEpZADeCR+qw/mJ1532Z5XxhMtcPGGVNxWM3TYaZ4Yr/hmGqHnSOK2FnuH6v3FsDbwcpiK+9M9WSm4rCqWLOqGmNtMHsWRBO71USSrA0sjwoblwC/Lwra3bn+xcJDMx3KveI9yoTTIYQQQhSKh5wOBZJ6DzkdQggh/APFdJR6FNMhhBBCCJ+gmQ4hhBD+gZHr/oyF4abUVlhCTocQQgj/IC8HyDMLyzX7rJwOX1AmnI4oOKst2PA4Q+xm6hUGUzpYza3B7GY7/Z+3WBarK8sHw/JYmKlXmNKB1dVqn9cjdha9D/ABbPXxwVQL7DeTWV4KRgaxM6XDaYvlm/2+s5oXiKmCIok9jNjNcpmwe9Wq8oOdg41xMzx1b1vN4WL2NcnUXmw8sbHJzsHqaqaGKuo18rl6xUNIveI9yoTTIYQQQhSKh2Y6pF7xHnI6hBBC+AdaXikW+/fvL/KxzZo1c+sccjqEEEIIgebNm8Nms8EgQbX579lsNuTmmgUIcOR0CCGE8A8001EsDh486PVzyOkQQgjhH8jpKBZ169b1+jm0OZgQQgghnPj73/+Odu3aISYmBj///DMAYObMmfjwww/dLtPvZjqYrNNM4shkXcwjY8enEzuTpJklnmKSNNaOKsTOkl5ZTShn9h7rD5ZUzqpcl8kxAS6NZX3LrimTm7Lr4E5yNTYOWGI8loiLXQeWjA3gfc7GB7MzKS27RmZjnLWDSYXTid2qFNlMys2S9bGyIoid9R9rs9mvP6ahYH1udYxbPS/AJfoF73nfJ3zLs95Q+2c9WpMyz7x58/Dss89i7NixeP755+0xHJGRkZg5cyZ69uzpVrma6RBCCOEf5BXz9RsJCQmIj4/HnDlzfFj50sXrr7+OhQsX4qmnnkJg4O/ucmJiIr766iu3y/W7mQ4hhBDXKAWcB8uf/Q3t03ElqLRFixZO9tDQUGRmsi3rCkczHUIIIYRwIC4uDnv37nWyr1u3DvHx8W6Xq5kOIYQQ/oGHZjoE8NhjjyEpKQmXLl2CYRj4z3/+g2XLlmHatGl488033S5XTocQQgj/QE6Hxxg+fDhycnLw+OOP48KFCxgwYABq1aqF1157Df369XO73DLhdJyFc/Q9U2KzKH131Cts1YqNTXa8O6tfrB2RxF6V2FmUOasTU0wAXInAVBls7Y61zapaCODXgq3GlicnzyKNYKoFs4RerG/ZZ9jYZNeO4Y7yiMGu9SViZ/ej2fqtVYUHe1gdJfZfiJ0lKAS4+o2pwxhs3DC1ldmDmI5lYmcKI3Zuq8npAD4+CqritH5fthkxYgRGjBiBkydPIi8vD9WqVSt2mRoTQggh/AMPqVcE8Mwzz9hlstHR0XaH4+zZs+jfv7/b5crpEEII4R8YcN/h0IakDixduhTt2rXDTz/9ZLd9+umnaNq0KQ4dOuR2uXI6hBBCCOHA/v37Ua9ePTRv3hwLFy7EY489hs6dO2PYsGH4/PPP3S63TMR0CCGEEIWiQFKPUalSJSxfvhxPPfUU/u///g9BQUFYt24dbr/99mKVq5kOIYQQ/oFiOjzK66+/jldffRX9+/dH/fr1MXr0aOzbt69YZZaJmY6f4RxVHkmOrUzsZt4Vi2Y/Q+wsup7lDmGR7Gadz3I6WI3sZ6oClnuFqRPMPsNyWbA+ZwoBT9aJ5X3JJWH3rF/ZdWBqIYBH9rNrxMYNewYyu5nKh10jplZieVxYvg/W32a5VxhMlcHGMmvDcWI3i7+vTuys3ezcLH8MO95MqWQ1J5HV/C6sX81ULUXNeeTO9S8NJCQkICAgAElJSUhKSirp6pQI3bp1w86dO7F06VLce++9uHjxIpKTk9G6dWtMnjwZjz/+uFvllgmnQwghhCgUbYPuMXJycrB//37ExMQAAMqVK4d58+bhrrvuwoMPPiinQwghxDWOYjo8xsaNG13a77zzTiV8E0IIIeR0+Ibo6Gi3PyunQwghhBCoUqUKvv/+e0RHR6Ny5cqw2VjED3D6NItcMkdOhxBCCP9AMx3F4tVXX0VExJXw+ZkzZ3rlHGXC6TgO54oy1QKL4GfHAzzSnKlXWKQ5UwiwTmYR/wCPQD9H7CyfCVOKMEWNWbQ561tmZ22wqqgxexawPmfXqBLpqAgie4ogJwg1cfKZ4oXVlRXF+pVdI/6bhMMUEJWInSkpWE6bbBMJxFliZ8oIpr5g9xfLQVKD1oirkthYzjApyxWsTu4o2SKJnV1TNp6YosbsWcDGWsHnjc83+ZTTUSyGDh3q8m9PUiacDiGEEEL4ltzcXKxevRrffvstbDYbGjdujJ49eyIoyH3XQU6HEEII/0AzHR7jwIED6NmzJ44dO4brr78eAPD999+jatWqWLNmDZo2bepWudqRVAghhH+gHUk9xoMPPogbb7wRv/76K3bv3o3du3fjl19+QbNmzfDnP//Z7XI10yGEEEIIB/bt24ddu3ahcuXfA94qV66M559/Hi1btnS7XM10CCGE8A800+Exrr/+ehw/7pxMIC0tDQ0aNHC73DIx05ENZ2UDy5fCguWZcgDgEegs5wdTfjBlCYv4j6U1AtgGvEylwtrtydwHVvN3WM3pkGmtOqawgV2eXIxQInexEbe8SiQ/d1Vy8gtk0IaSgcbUHUyZw1QOZjAFFVOjsP7IIRfVbPyx99iznx3PxhPrJ3Y/AjzvC4Pdj5HEztpmpvJg14i1j7XBav+Zwe6vgnViqh+voZgOjzF16lSMHj0akyZNQuvWrQEAX375JZ577jm8+OKLyMj4/ZuzYkWmy3KmTDgdQgghhPAdd911FwDgvvvus28SZhhX3OMePXrY/7fZbMjNNduUwhGvLq9MmzYNLVu2REREBKpVq4a7774b3333nTdPKYQQ4lrFgPtLK8XYVGTu3LmIi4tDWFgYEhIS8Nlnn5kev2XLFiQkJCAsLAz169fH/PnznY5ZuXIl4uPjERoaivj4eKxevdrh/a1bt6JHjx6IiYmBzWbDBx984H4DXLB582b7a9OmTdi0aZPL/zdt2mSpXK/OdGzZsgVJSUlo2bIlcnJy8NRTT6Fz58745ptvrvkMfkIIITxMCSyvrFixAmPHjsXcuXPRrl07vPHGG+jWrRu++eYbxMY6L6IfPHgQ3bt3x4gRI/D2229j27ZteOihh1C1alX07t0bALB9+3b07dsXf/3rX9GrVy+sXr0a9913Hz7//HO0atUKAJCZmYmbbroJw4cPt3/Ok3To0MHjZQKAzcifL/EBJ06cQLVq1bBlyxbcdtttpsdmZmYiPPzK6uutcF4bjCSfY2usnozpYNNDbPdPd2I6Yog9itjZuq/VqSyz9V0WR8PsLKbDnbV2q9RmdpKnqDxZkswmAyfLJAAl0GJMx1E/junIMnmQs/ZZjelIJXYWY1WL1ojfXwzWBhbn5MmYjirEbjWmg+2Im84qBL57cMHbJRvAX3/7+/z58175sXn1d8X53kAFN39KZ+YA4Suv/G2lrq1atcLNN9+MefPm2W2NGzfG3XffjWnTpjkdP378eKxZswbffvut3TZy5Ejs27cP27dvBwD07dsXGRkZWLdunf2Yrl27onLlyli2bJlTmTabDatXr8bdd99dpDqXJD5Vr5w9e+UWrVLF9e2SlZWFjIwM+0sIIYTwNVd/D2VkZCAry/Wvj+zsbKSkpKBz584O9s6dO+OLL75w+Znt27c7Hd+lSxfs2rULly9fNj2GlVmW8FkgqWEYSE5Oxq233oomTZq4PGbatGmYPHmyk92VV82iolmDzH5JsNkRNtPBftmzmQ72y8Msspt9htlZvgX2i4TNaLgz08F+9DOP1qoqiLUN4DlC2PVmsw2MIPIzs7JJAo9g8pnQdNf2sz+7trNrUYspamryOhnkJ/ZlMpNTgcz85JGOPX3Utf2syYBil8KqQotda3ceblbVXmxsst/IbIyb/cSyeq+yOllVr7D7FCj6bJQn1XNFwkPLKzExjnPNEydOxKRJk5w+cvLkSeTm5qJ69eoO9urVq+PYsWMuT3Ps2DGXx+fk5ODkyZOoWbMmPYaVWZbwmdMxatQo7N+/H59//jk9ZsKECUhOTgZwZcqs4IUXQgghGEYed7CL8tl8UlNTHZZXQkPN0nPCKQV8vqrDyvEF7VbL9CSGYeDw4cOoVq0aypVjC7ru4ROn4+GHH8aaNWuwdetW1K7NVtuvXNj8ixsY6HOFtxBCCIGKFSsWKaYjOjoagYGBTjMQaWlpTjMV+dSoUcPl8UFBQYiKijI9hpXpaQzDQMOGDfH111+jYcOGHi3bqzEdhmFg1KhRWLVqFTZt2oS4uDhvnk4IIcQ1TF5e8V75JCQkID4+HnPmzDE9X0hICBISErBx40YH+8aNG9G2bVuXn2nTpo3T8Rs2bEBiYiKCg4NNj2FlepqAgAA0bNgQp06d8njZXp3pSEpKwrvvvosPP/wQERERds+tUqVKHp+yEUIIcW3jqeWVlJSUIqtXkpOTMXjwYCQmJqJNmzZYsGABDh8+jJEjRwK4EjZw5MgRLF26FMAVpcrs2bORnJyMESNGYPv27Vi0aJGDKmXMmDG47bbb8OKLL6Jnz5748MMP8e9//9shPOH8+fP48ccf7f8fPHgQe/fuRZUqVVxKda3y0ksv4bHHHsO8efNoHKY7eNXpyJcQdezY0cG+ePFiDBs2zJunFkIIIbxO3759cerUKTz33HM4evQomjRpgrVr16Ju3boAgKNHj+Lw4cP24+Pi4rB27VqMGzcOc+bMQUxMDGbNmuWw10bbtm2xfPlyPP3003jmmWdw3XXXYcWKFfY9OgBg165d6NSpk/3//HjIoUOHYsmSJcVu16BBg3DhwgXcdNNNCAkJcZooOH2aia3N8ek+HVa4WnvdBs5Kj2rkcyzc57LJuVjkOLMzDT5Tr7C6XkdrBLCFKLaPgKfUK2biDjbErKpXmJ21wR31ChNyRJHCWE4Wpl4pb7IpBlOvnE93bT9E1Cvk8DKlXjnqA/XKCWJn6guzVXE2nhjs/mIRae6oVxhsnw6WBYP130liN/tKKap6JRtA/k4Vvtin43S34u3TUeW3bTGuv/56BAQEICkpCUlJSR6qadnib3/7m+n7Q4cOdavcMpF7JQfOsjgmxbKaHAngnWBVlsvs7EvTTJLGHsisHZ6Sw7kjkzNzClzBHEC2mZNZ+cwJY85FRfIBtqFXMPEkg5iHCf4FzzbQYv3BTlGRfNtUIhufAQALes9Md223kErBFDNnn401dt+x8cfi+dn9aCYLZw40uxbMSWH3KXMuzJIdss8w593qs4A9a8yGAPvhVfCaml1/b1ASyyv+irtORWEotb0QQgghnPjpp5/w9NNPo3///khLSwMArF+/Hl9//bXbZcrpEEII4Rd4Sr0iruROa9q0KXbs2IFVq1bh/Pkrc2L79+/HxIkT3S5XTocQQgi/wNeSWX/miSeewJQpU7Bx40aEhPweadWpUyd7jhh3KBMxHUIIIURhGEYxYjquChxUTAfw1Vdf4d1333WyV61atVj7d2imQwghhBAOREZG4uhRZ1nanj17UKuWWa5mc8rETEc5OEe0swh35kWZbarOHGMmPWOSOxax7k4SOiZjY5HmVjeNZ7JEswwDzO9n/cei35nagPWHWdtYfSMqu7aHR7q255HKXiYSi4tMagMgm8hQmayUXVOqviDlM/krAIQQNU8ouagXiGQik9gvEmmEmQKCqSyYrNSqoIYdb6bQsjo2rf6otvrcArgChD0j2NC0mnytVO6lUAhGHmC4mZ7E3RkSf2XAgAEYP3483n//fdhsNuTl5WHbtm149NFHMWTIELfL1UyHEEIIv0AxHZ7j+eefR2xsLGrVqoXz588jPj4et912G9q2bYunn37a7XLLxEyHEEII4SsU0wEEBwfjnXfewXPPPYc9e/YgLy8PLVq0KHYCODkdQggh/AItr3ie6667DvXr1wcA2NgugxbQ8ooQQgi/QPt0eJZFixahSZMmCAsLQ1hYGJo0aYI333yzWGVqpkMIIYQQDjzzzDN49dVX8fDDD6NNmzYAgO3bt2PcuHE4dOgQpkyZ4la5ZcLpaATn3AcsFwKL0DaLxDZTbLiCCCNoAiamajHLS5BO7CwJEzsHy03CckaY5TlhE2ssDwNTZbByWGS/2YQe68NLpFIsxwojkFSKqT7MzsESskUSmQXL1WIj85OXTBJ45JEbgyW0YyoflnMml5QfZlInlgsknX/EJeyedyd9DFN1McVLOrEzBQ5rMzve7D3Wtex+Yf3E2szsAH/eFBxOZnluvIGnllcSEhKu+YRv8+bNw8KFC9G/f3+77U9/+hOaNWuGhx9+2L+dDiGEEKIw8vKAPDedjquXVxRICuTm5iIxMdHJnpCQgJwcqwLs31FMhxBCCL9AMR2eY9CgQZg3b56TfcGCBRg4cKDb5WqmQwghhBBOLFq0CBs2bEDr1q0BAF9++SV++eUXDBkyBMnJyfbjXnnllSKXKadDCCGEf1CMmA7L28v6OQcOHMDNN98M4EqKe+BK3pWqVaviwIED9uOsymjldAghhPALPBXTIYDNmzd7pdwy4XTUhHMkNQtGYZHeZsErTL1iNU9CuEW7WXQ9U9sw1Q7LBxNN7KzNZnWyGpHPBhdTzrCIeLPcK0y9cuKsa/t5Yi9HYsaiieKEHQ8A5cjFYIqQANLwy0T+k57m2p5lIoHIIR2VTT4TQurElDnBZGCaPf+Z6ukEsbMxy+4JpsQyU2ixscbOwe4Xptxi94qZko1hVQnD6sSeZ9VMzs2Gf8Gyyur3uNQr3qNMOB1CCCFEYXhKMiv1iveQ0yGEEMIv0DbopR9JZoUQQgjhEzTTIYQQwi/IM9wPCM0z27b6GiQzM9MrS0ya6RBCCOEXGHnFe4nfqV69Ou6//358/vnnHi23TMx0RMA5cp05pWxzVrOGViV2pjph57Cq7jCD5T1gEfws7wtTirAo/XO0RrzdzF6R2Fl/R5FGMyUFAFwgapQMcjyN1K/j2l491rW9AktqA/BOZx1CLmowSa5Rtbxr+9ljvEqZpJ8uEknDWSIhYflgTqe7tp/iVaJjzaqChN0T7uQUsQp7rrA2MHuWG+e2qpw5Q+xWlXJA0ZVEvs69IjzHsmXLsGTJEtx+++2oW7cu7r//fgwZMgQxMTHFKlczHUIIIfwCbYPuOXr06IGVK1ciNTUVf/nLX7Bs2TLUrVsXd911F1atWuV2/hU5HUIIIfwCTzkdCQkJiI+Px5w5c0quMaWEqKgojBs3Dvv27cMrr7yCf//737j33nsRExODZ599FhcusB13XFMmlleEEEKIwjDy+HJRUT6bj/bp+J1jx45h6dKlWLx4MQ4fPox7770XDzzwAFJTU/HCCy/gyy+/xIYNG4pcnpwOIYQQQjiwatUqLF68GP/6178QHx+PpKQkDBo0CJGRkfZjmjdvjhYtWlgqV06HEEIIv8BTMx0CGD58OPr164dt27ahZcuWLo+pX78+nnrqKUvllgmnIxDOUd9sYLHN6Mzyd7BJtAhiZ2OTRWqzqHGzHBBMOWM1H4zV+4hFpQOWxRc0l0o4qWwVlueEdQaAC6RSgUdd22vWd22v34ScoBaxm3UUu0jsbgsiHVKOXD2SFKOSyd0cSOp0gch8Tvzi2p5OqsTUQmbjj+ULYkoOptZgKgumUjELf7OqImFjnA0PVj7rPwA4bdHOFENEwETrms4qhKI/V9wLNXSfvDz3870okNSRo0ePonx5IpX7jXLlymHixImWylUgqRBCCCEciIiIQFqac3bJU6dOIZD9iikCZWKmQwghhCgMLa94DsNw3ZNZWVkICXF/txs5HUIIIfwCLa8Un1mzZgEAbDYb3nzzTYSH/76+nZubi61bt+KGG25wu3w5HUIIIYQAALz66qsArsx0zJ8/32EpJSQkBPXq1cP8+fPdLl9OhxBCCL9AyyvF5+DBgwCATp06YdWqVahcubJHy1cgqRBCCL9AO5J6js2bN3vc4QDKyExHNpylsMxb8mSSJ+b4XvaQ3Wy/OyZjY178JWJnEmKmQo1kycpMyCKJyUJJA1ngcwhpNLMDQDbRS0aQDHiVWLY5ljGP6YTNAqkq1nZtj6zn2h5AbsP0Q67tJ793bTfpp3DyXlWiQz1JJMdB5FqTfHkINHnCZBA9Jck1RyWfTHrOuoNJbwGeLC2d2Jlsld137D5lsnqAy2lZXdnx7DnE+s+dX6QFtw3wtWTWU1yrO5ImJyfjr3/9KypUqIDk5GTTY1955RW3zlEmnA4hhBCiMAzD/WUSIta4ptizZw8uX75s/5thszG3unDkdAghhPALjDwgz83vQzkdV5ZUXP3tSRTTIYQQwi8w8or3EpyMjAx88MEH+O9//1uscuR0CCGEEMKB++67D7NnzwYAXLx4EYmJibjvvvvQtGlTrFy50u1y5XQIIYTwCzylXhHA1q1b0b59ewDA6tWrYRgG0tPTMWvWLEyZMsXtcstETMcpOKtPWGQ6W84z865YVDxL4MbGJtuNnuTnMoVFfbNzMC2F1QRxAW6sh4aQrFdhJFeQjZyc2c0IYLnSiDwnmHWU1XNHmezI17iXa3vVG62d48TXru3f/9O1Pe0rXlae69EcHun68IpEzRNsUS3Erg8AhBKZSja5Ia1eIvaMYAniAICJt5iIiYh8qJ3d15GsQgCYaJEpYVgb0omdfQmYZdcord/PRh5gKKbDI5w9exZVqly5sdevX4/evXujfPnyuPPOO/HYY4+5Xa5mOoQQQgjhQJ06dbB9+3ZkZmZi/fr16Ny5MwDgzJkzCAszy5FuTpmY6RBCCCEKI68Y6pU8zXQ4MHbsWAwcOBDh4eGoW7cuOnbsCODKskvTpk3dLldOhxBCCL9Ayyue46GHHsItt9yCX375BXfccQcCflsrrV+/vv/HdAghhBDCtyQmJiIxMdHBdueddxarTDkdQggh/AItr3iO3NxcLFmyBJ988gnS0tKQV0Des2nTJrfKLRNOx0E4V5Ttis8i1s3CXlhZRHxBI9mZnXVyFq0Rj0xnOVOYnbWb5WHISGc14jkrmBqAqVGYqiWYXDx3VC0s5wc7h2WZT9XG/OTVb3Jtr9mCf8YVwUQWdJJsznPiW14WUa8wQsmp2bUoF2GpeABACEkSYvU+umDxvGYPPXa/MNETuyfY/cWGUySrELgij7Xbai4Vkk7HVKHC8r4UfG6Z5bnxBnI6PMeYMWOwZMkS3HnnnWjSpEmxtj6/mjLhdAghhBDCdyxfvhzvvfceunfv7tFyfSKZnTt3LuLi4hAWFoaEhAR89tlnvjitEEKIa4iS2gbd6nfcli1bkJCQgLCwMNSvXx/z5893OmblypWIj49HaGgo4uPjsXr1asvnXbVqFbp06YLo6GjYbDbs3bu3yG0KCQlBgwYNinx8UfG607FixQqMHTsWTz31FPbs2YP27dujW7duOHz4sLdPLYQQ4hqiJJwOq99xBw8eRPfu3dG+fXvs2bMHTz75JEaPHu2wtfj27dvRt29fDB48GPv27cPgwYNx3333YceOHZbOm5mZiXbt2uGFF16w3K5HHnkEr732GgwPy3pshqdLLECrVq1w8803Y968eXZb48aNcffdd2PatGn0c5mZmQgPvxKp0BvXXkwHWxNmsRtkQ0jLMR1muxBajekoT7ZGZDEdkWTr1iCTLSQz013bWexBzfqu7eWvIyeIJvYmfXmlbiTvWY3pSDvg2r7/bdf27z7iZZ12HQGQe9D14YfIZqjZZNCynU3NOE227TyV5trOYhisxnSQcBUA/H5hMVakqjhC7L6I6WC7oaYSO4vpMIM9C1zFdGz97e/z58+jQgX2tHWfq78r1gEo52bowUUD6Pbb31bqavU7bvz48VizZg2+/fb3GKyRI0di37592L59OwCgb9++yMjIwLp16+zHdO3aFZUrV8ayZcssn/fQoUOIi4vDnj170Lx58yK1q1evXti8eTOqVKmCG2+8EcEFtiNetWpVkcopiFdnOrKzs5GSkmLfySyfzp0744svvnA6PisrCxkZGfaXEEII4Wuu/h7KyMhAVpZrb9vqdxxwZRaj4PFdunTBrl27cPnyZdNj8st057xWiYyMRK9evdChQwdER0ejUqVKDi938Wog6cmTJ5Gbm4vq1as72KtXr45jx445HT9t2jRMnjzZyX4ezhVl3hKbIWO/VMzes/rLns0SsOPN6sS0BiyKnsF+8QeSjjpnUhaLWGd1iiA5NKqSjmLB0WbTnpfJlA3LEcJmWejP32D2AR+Qy+ajCAbL7AE6mHPJR5jKh10Ldu3y3JAvRJDnWQj5SV7BpNmuMLvv2FBjP57ZA5SNGvaMMBtl7DNsipoJidjXhDtfAqyfCl4ii5em2Bhwf5Ovqz8WExPj8N7EiRMxadIkp89Y/Y4DgGPHjrk8PicnBydPnkTNmjXpMfllunNeqyxevNgj5RTEJ+qVglIbwzBcym8mTJiA5ORkAFemzApeeCGEEIKRB/eT0V39udTUVIflldBQtnB/haJ+x5kdX9BelDKtntcqOTk5+PTTT/HTTz9hwIABiIiIQGpqKipWrGhf0rKKV52O6OhoBAYGOnleaWlpTh4acOXC5l/cwECz6AIhhBDCO1SsWLFIMR1Wv+MAoEaNGi6PDwoKQlRUlOkx+WW6c16r/Pzzz+jatSsOHz6MrKws3HHHHYiIiMBLL72ES5cuuVTcFAWvxnSEhIQgISEBGzdudLBv3LgRbdu29eaphRBCXGPkFfOVT0JCAuLj4zFnzhzT87nzHdemTRun4zds2IDExER7sCY7Jr9MX3y3jhkzBomJiThz5gzKlft9/blXr1745JNP3C7X68srycnJGDx4MBITE9GmTRssWLAAhw8fxsiRI719aiGEENcQBnisS1E+m09KSkqR1SuFfcdNmDABR44cwdKlSwFcUarMnj0bycnJGDFiBLZv345FixbZVSnAlS/82267DS+++CJ69uyJDz/8EP/+97/x+eefF/m8AHD69GkcPnwYqalXtEvfffcdgCszKTVq1DBt1+eff45t27YhJMQxaq9u3bo4coTpswrH605H3759cerUKTz33HM4evQomjRpgrVr16Ju3brePrUQQgjhVQr7jjt69KjD3hlxcXFYu3Ytxo0bhzlz5iAmJgazZs1C79697ce0bdsWy5cvx9NPP41nnnkG1113HVasWIFWrVoV+bwAsGbNGgwfPtz+f79+/QDwwNirycvLQ26ucyT4r7/+iogIN/Ie/IbX9+lwl6u113fC2TtiYgN3IkHYGhOLKGfhM0zFwcJ6zDqeKV7Yfhxkiwsasc7KN9s7hImYmc6f9VMs6ahoN+KGz55yba8U5dp+/S2kILZPRwUyCsxyrzS6y7Wd5WTJIzH+J8hmGd8670wIADhFcrIAwGnXeqgLJF3LsUOu7Uy9Uo5cbLavBwCcPcHfc0XWRdd2to9LALmxz5zk50gndqYmY0oYdjzrDrPnFnuPBUwyxQhTpqUTu5nwiNWpYLsvA8hfoPDFPh0fwnxPJjMuAej529/XX389AgICkJSUhKSkJA/UsuzRt29fVKpUCQsWLEBERAT279+PqlWromfPnoiNjXVb3aLcK0IIIfwCT6lXrCyv+CuvvvoqOnXqhPj4eFy6dAkDBgzADz/8gOjoaIelIKvI6RBCCCGEAzExMdi7dy+WL1+OlJQU5OXl4YEHHsDAgQMdAkutIqdDCCGEX+CpmY6EhIRrfnll69ataNu2LYYPH+4QF5KTk4OtW7fitttuc6tcOR1CCCH8gpJQr/grnTp1wtGjR1GtmmPE4NmzZ9GpUyeXQaZFQU6HEEIIv8BTTofgu5ueOnWqWA5ZmXA6LsK5oszHYpHbZlks2HSc1SyzbJWLqThYAlOzslg7WLZLli/FakZcgPcTi8hnN3EaCe2/eMi13cyfJoIGmnOGZaWtwKQ5ASSX57E9vFLZRCfw63bXdqZeOf2Ta3vGL6Qck4ll0udMEXKBNKE8UcoFupFKmeV9ySGDPJvUlV1rprQxywdiVaXClG/s/mXln6E14vcqU7KxS8HqxPrDLEdNUTfAZu0VpZd77rkHwJUt1ocNG+awBXxubi72799frA3IyoTTIYQQQhSGp2I6rmXyM8gahoGIiAiHoNGQkBC0bt0aI0aMcLt8OR1CCCH8Ak8tr1zLgaT5+2/Uq1cPjz76qMdjW+R0CCGEEFehQNIru5Z6A68mfBNCCCF8hacSvgng+PHjGDx4MGJiYhAUFITAwECHl7topkMIIYRfIPWK5xg2bBgOHz6MZ555BjVr1nSpZHEHOR1CCCHEVVzLMR35fP755/jss8/QvHlzj5ZbJpyOiih6RZny8azJZ1giJHbOWsRej9hZAuHqrELgSYuYjO00sR8l9nRiZzJhgP8SYLI4Ju9ldjZhx5LTAVy6d5pc8F+/d22vTcqpUJO8EWkyGZtJTvILsbMOYR3LBqaJFjSPDIRTJEM1S8YWTC5GdqhrexbLBgjgHBm0RKRMpdPB5FqztWMz+TxbxWdjk0m22aVg95DZGGfnZsOAlWV1QtysTqyfCva5WQJJb6DcK56jTp068EY+WMV0CCGE8AuMYr7E78ycORNPPPEEDh065NFyy8RMhxBCCCF8R9++fXHhwgVcd911KF++PIILTHOePs3m182R0yGEEMIv0OZgnmPmzJleKVdOhxBCCL9A6hXPMXToUK+UK6dDCCGEX6DU9sUjIyMDFStWtP9tRv5xVikTTkddACEFbCw6/BSxmyUoYiIBliyNRYEzxQlThLAETACPDiciAQpr20liNxsQrL4syp3d/KxOLAEeU5YAQJVKru0BRFJ+lgwQ24+u7ZWJtKkqkzABvCHs4jGFR7pr8wUyMFmSNgA48atr+6GfXdup6iDVtTks3bX9IruJwMcgO3fBZ0A+TI3CxqWZQovBnjdsLDM7u4dIHj0AXLXDVD4MdinY7gtmzyczZcvVuL+FVMlyrapXKleubE9nHxkZ6XJvjvzss0ptL4QQ4ppGMR3FY9OmTahS5Ur+4s2bN3vlHHI6hBBC+A2KzXCfDh06uPzbk2ifDiGEEEL4BM10CCGE8AukXin9yOkQQgjhFyimo/RTJpyOCDhHrluNmzUT97BIc6t5FdigZeIEk8B+uu5lVTnD2s3OzdoMWK8Ti3BndY0kdqZQAYAqJIGNjVT2ErkYOUQCwXKQXDKRDlROd20PJcHwF4nqhCltWD6TDJMNAtPJe1bXV3+1KNdghwN8DDKRD1OvsIcYG5dmuVesKtmYqJD9amZtM1ODsLLYudnzhuVtYpjVifVhwWejr3OviNKPYjqEEEL4BZ7KvZKQkID4+HjMmTPHd5W/RigTMx1CCCFEYSjLbPFo0aKFy705XLF79263ziGnQwghhBC4++67vX4OOR1CCCH8AqlXisfEiRO9fg45HUIIIfwCqVdKP2XC6XAVUc4iYJlaw8yLZVHdkRbLYoOW2d3JS2A11wM7B8tFY7aaZzXqmJ2bnYMpkrJNJDXZJDw+lCSOCCHyAVbO+XTX9tPHeZ2OHeLvuSKXyKFYXUNI24JNEvPUqmepSsg449p++Cyxk3LM1CtsPLFxwI5n52D3HVOfAfxZwNQrZgo0V7A6MUUXwPuDqVdYThsmuGL3qVmOGjbUCt5GZkohb6CZDs+Rm5uLV199Fe+99x4OHz6M7GzHO+30aRO5nAlSrwghhBDCgcmTJ+OVV17Bfffdh7NnzyI5ORn33HMPAgICMGnSJLfLldMhhBDCL8gr5kv8zjvvvIOFCxfi0UcfRVBQEPr3748333wTzz77LL788ku3y5XTIYQQwi/w1D4dAjh27BiaNm0KAAgPD8fZs1fWVu+66y58/PHHbpcrp0MIIYQQDtSuXRtHjx4FADRo0AAbNmwAAOzcuROhoSYBZIUgp0MIIYRfoOUVz9GrVy988sknAIAxY8bgmWeeQcOGDTFkyBDcf//9bpdbJtQr5+CcB4CpLyKJ3UyVwcpiuR5YRDaLomepQ6JpjXjkOBNykDQdtN1srz0zRY3VyHvm0bI2MCXABRMJREi6a3sO+QxTtZQnsqcgkoAi81depyPk3EyFUIPYy5GLFEzqFGwigQgh7+WRSl0mg7w8Ua+we8Uslw8bTyQVjWlZrmB5P8zUK+z3G7sfmVqOKUvSiN2sTmzcMKUN6z/23GLXwewXKXtOFKyTWbu8gafUKwkJCQgICEBSUhKSkpI8ULOyxwsvvGD/+95770Xt2rXxxRdfoEGDBvjTn/7kdrllwukQQgghfMW1ug26Ga1bt0br1q2LXY6cDiGEEH6BNgfzLN9//z0+/fRTpKWlIS/PsYeeffZZt8qU0yGEEMIv0OZgnmPhwoX4y1/+gujoaNSoUcMhEZzNZpPTIYQQQgjPMGXKFDz//PMYP368R8uV0yGEEMIv0PKK5zhz5gz69Onj8XLLhNMRAOdIaqtaXxZdD1hXr7AoehY17o4umUWHswtm1j5XsKlEs1wJrE5Ww61YND4RZaC8SeMqMPkAgeVYCSLnYHazPCe5RL3Cxll0Tdf2yGqu7TYiScpmcgYAF4icgvXH2ROu7ewUbGyYjUs2DpjqhNnZOdh4MhsyTGnGhEFMWMVUYyz/CREFAQBIGhyq9rL6jGB2sy8H9vwoqFYpCfWKu86Dllcc6dOnDzZs2ICRI0d6tNwy4XQIIYQQhaGYDs/RoEEDPPPMM/jyyy/RtGlTBBfQ6o8ePdqtcuV0CCGEEMKBBQsWIDw8HFu2bMGWLVsc3rPZbHI6hBBCXNtopsNzHDx40CvlyukQQgjhFyiQ1DsYxhWXzMYCyiyg3CtCCCGEcGLp0qVo2rQpypUrh3LlyqFZs2b4+9//XqwyNdMhhBDCL9BMh+d45ZVX8Mwzz2DUqFFo164dDMPAtm3bMHLkSJw8eRLjxo1zq9wy4XREwFnWxSR6TPIZYVI+k8kxVSST4jFJH5PYnqY14nJCZmcwORyTPpqoLmnSK9ZP7BoxO8nFBsODT4OLRGeYTS5SILlDzjN9NPg1Yu22kflG1m4mc72YyevEysokUtrTRAtKlLSW5awAl1Oyspg8lT3EmDSWSZcBLv9m52CTzUxiy54dZkOcXdbjxM6mr1m7qxC7WfJHNsYLPn/LasI3Abz++uuYN28ehgwZYrf17NkTN954IyZNmuS20+G15ZVDhw7hgQceQFxcHMqVK4frrrsOEydORHa2ScpQIYQQooywdetW9OjRAzExMbDZbPjggw88Uu6WLVuQkJCAsLAw1K9fH/Pnz3d4f8mSJbDZbE6vS5fMfjZa4+jRo2jbtq2TvW3btjh69Kjb5XrN6fjvf/+LvLw8vPHGG/j666/x6quvYv78+XjyySe9dUohhBDXMEYxX1bJzMzETTfdhNmzZxe/8r9x8OBBdO/eHe3bt8eePXvw5JNPYvTo0Vi5cqXDcRUrVsTRo0cdXmFhbI7NOg0aNMB7773nZF+xYgUaNmzodrleW17p2rUrunbtav+/fv36+O677zBv3jxMnz7dW6cVQghxjeLrmI5u3bqhW7du9P3s7Gw8/fTTeOedd5Ceno4mTZrgxRdfRMeOHeln5s+fj9jYWMycORMA0LhxY+zatQvTp09H79697cfZbDbUqFHDjVoXjcmTJ6Nv377YunUr2rVrB5vNhs8//xyffPKJS2ekqPhUvXL27FlUqcJWEIUQQgj/Yfjw4di2bRuWL1+O/fv3o0+fPujatSt++OEH+pnt27ejc+fODrYuXbpg165duHz596iZ8+fPo27duqhduzbuuusu7Nmzx6N17927N3bs2IHo6Gh88MEHWLVqFaKjo/Gf//wHvXr1crtcnwWS/vTTT3j99dcxY8YMekxWVhaysq6EkWVmmkTFCSGEEAXwVCBpRkYGcnN/D5cNDQ1FaKhJwiUX/PTTT1i2bBl+/fVXxMTEAAAeffRRrF+/HosXL8bUqVNdfu7YsWOoXr26g6169erIycnByZMnUbNmTdxwww1YsmQJmjZtioyMDLz22mto164d9u3bV6ylj4IkJCTg7bff9lh5gBtOx6RJkzB58mTTY3bu3InExET7/6mpqejatSv69OmDBx98kH5u2rRpLssORdHVK2yKzGxKx6rKgg09ZmfuUyqtEU/yxFbsmDqHtcFqsi0AiCR2qxH/TKXCVERmnCbxTLkkbJ5dC9YfTJ1ARB8AuFqJJftKJzKELFLZXFLZ8iaZzEJJp+cyuRc5t9UwNTP1CrsnrSo8zFQWrjBTgLH2sbpaVYGxL0SzLZdY+5iajJ3Daj+ZhfwzVUpJK0A8tbyS7yTkM3HiREyaNMlSebt374ZhGGjUqJGDPSsrC1FRUQCA8PDfNUWDBg2yB4wW3ISr4OZcrVu3RuvWre3vt2vXDjfffDNef/11zJo1y1I9ryYjIwMVK1a0/21G/nFWsex0jBo1Cv369TM9pl69eva/U1NT0alTJ7Rp0wYLFiww/dyECROQnJwM4MpMR8ELL4QQQjA8NdORmpqKChV+/zlldZYDAPLy8hAYGIiUlBQEBjq6fPnOxt69e+22/C/xGjVq4NixYw7Hp6WlISgoyO6sFCQgIAAtW7Y0XbYpCpUrV8bRo0dRrVo1REZGutyB1DAM2Gw2h5kgK1h2OqKjoxEdHV2kY48cOYJOnTohISEBixcvRkCAeQjJ1VNYBS+SEEII4Qs6deqEgIAAJCUlISkpya0yWrRogdzcXKSlpaF9+/Yuj2nQoIGTrU2bNvjoo48cbBs2bEBiYqJTptd8DMPA3r170bRpU7fqms+mTZvscZebN28uVlkMr8V0pKamomPHjoiNjcX06dNx4sTv2wp5M+JWCCHEtYmnlldSUlIcZjoY58+fx48//mj//+DBg9i7dy+qVKmCRo0aYeDAgRgyZAhmzJiBFi1a4OTJk9i0aROaNm2K7t27uyxz5MiRmD17NpKTkzFixAhs374dixYtwrJly+zHTJ48Ga1bt0bDhg2RkZGBWbNmYe/evZgzZ46brb9Chw4d7H/HxcWhTp06Lpd6fvnlF7fP4TWnY8OGDfjxxx/x448/onbt2g7v5a9PCSGEEJ7C1zuS7tq1C506dbL/nx8eMHToUCxZsgSLFy/GlClT8Mgjj+DIkSOIiopCmzZtqMMBXPmyX7t2LcaNG4c5c+YgJiYGs2bNcpDLpqen489//jOOHTuGSpUqoUWLFti6dStuueUWN1rB65G/1HI1p0+fRlxcnNvLKzajlHoAmZmZ9nWviSh+IKnrlbArVCN2torHgi3ZFtFsu3OzS8aCLT0VSMrqdIzYAaA6sbN5K+bRsnazQNIwN1zjkgwkZeeoTOw1SPvKkX2r3QkkDSGD+QJpyC8nXdtZsmv2EDFbCWfBiGyHeXbfsXOwe4LdW2ZlWQ0kPUvs7Blhdt+x91iwObsW7NnB7juz245du4JB1LkANv329/nz54s0e2CVq78rHoV58LIZ2QDyd5C6/vrri728UtYJCAjA8ePHUbVqVQf7zz//jPj4eLcVpmUi90oOnG969mXAbjiScgMAH6QsVwH7gmLnZgIBs9wrDBZNw/qDOR2sP8wi1tlgYQ9q1q+sP9zZNCaTPP3Yg9pqn7M6md1u7FrQnCKkDSHpru1MtWAWaZVNJCEsFw0by+wL22q+GcB6LhWrqjGmgjFzhNhXImsfu4+Y08EcKrMcJew5xHxM1j527cycMEZRnSpf517x9fKKP5I/W2Oz2fDMM8+gfPnfnzi5ubnYsWMHmjdv7nb5ZcLpEEIIIQpDCd+KT/4mY4Zh4KuvvkJIyO8/H0NCQnDTTTfh0Ucfdbt8OR1CCCGEAPC7amX48OF47bXX3N6Pg+HTbdCFEEIIb5FXzFc+CQkJiI+PL7YapCwzc+ZM5OQ4L5CdPn260I3DzJDTIYQQwi8w4L7DcfXySkpKCr755ptrNogUAPr164fly5c72d97771CNwg1Q06HEEIIIRzYsWOHgxw4n44dO2LHjh1ul1smYjqy4BzkYzU/wymT8pkSgcUuM3UCk8+lE3saqxB4tLzFVBkUNjlmJuNl6gsWoc7UK1YVIRdMQuBZPhOmKmB2q+W4IxZjKgQmpWUqBKZSyWbSCHClQxC5GEzRwNrAlChm6gU2ntgYZ6odVif2jGBSWrOyWPvY+GBBiaxOZpJ+1ofsOcSksWw8sX51R3lS8DPsWnoLBZJ6jqysLJfLK5cvX8bFiyzDVOFopkMIIYRfYBTzlY9iOoCWLVu6zJc2f/58JCQkuF1umZjpEEIIIQpD+3R4jueffx5//OMfsW/fPtx+++0AgE8++QQ7d+7Ehg0b3C5XMx1CCCGEcKBdu3bYvn076tSpg/feew8fffQRGjRogP3799MEdkVBMx1CCCH8AsV0eJbmzZvjnXfe8WiZcjqEEEL4BZ5aXhGOXLx4EZcvO4YFu7tpWJlwOgLgvA7EKs7Wi8xibZnqhCk5WNQ9S8CUTuxm6hVWX6bCYWoDFuHO+olF7wO8vuzcVYidRdGzuppF0bNrZDWfDjsHU6mYJehi44kpFNgvLHZLs3wmZjk02HuB5EnLHsBW83cwVZDZOdi1YKoMZmeY/aJl52b3PLvWbIyz8WeWpIydw6oyxJNfqux61yzwv1kup9JMQkLCNZ/w7cKFC3j88cfx3nvv4dQp528ed7PMlgmnQwghhCgMTy2vKJAUeOyxx7B582bMnTsXQ4YMwZw5c3DkyBG88cYbeOGFF9wuV06HEEIIv0DLK57jo48+wtKlS9GxY0fcf//9aN++PRo0aIC6devinXfewcCBA90qV+oVIYQQQjhw+vRpxMXFAbgSv3H69GkAwK233oqtW7e6Xa6cDiGEEH6BpzYHE0D9+vVx6NAhAEB8fDzee+89AFdmQCIjI90uV06HEEIIv8BTWWbFldT2+/btAwBMmDABc+fORWhoKMaNG4fHHnvM7XLLREyHq9wrLJqcRfabNZRFgbPIaxZNfpbYjxP7aVojfgMwb5xFk1cn9khiZ/0H8Mh+1j7Wf0zVwnJimEX2s/hp1g7WBqsKJjNFDRtPTJ3Drh1TrzC72UOT5Qhh9xGDXQvWNjPVGOtb1g52Dnat2TU1qxMri41l1n/smlrNswNwBRBTVrG6stQ8rD9YfwM8X0vB/pN6pewybtw4+9+dOnXCf//7X+zatQvXXXcdbrrpJrfLLRNOhxBCCFEYUq94hsuXL6Nz585444030KhRIwBAbGwsYmNji122nA4hhBB+gdQrniE4OBgHDhyAzWY23+UeiukQQgjhFyiQ1HMMGTIEixYt8ni5mukQQgghhAPZ2dl48803sXHjRiQmJjotN73yyitulSunQwghhF+g5RXPceDAAdx8880AgO+//97hveIsu8jpEEII4Rcoy2zx+d///oe4uDhs3rzZK+WXCacjHUBwARvzSpn0zKyhTAZoVdLHjmeSNHeke8xeidirEXu0xfIBLrlj8tEMYrcqRTZL6MU+w+TIrE5M5srGDes/gI8PJjNksIArqxJbwPn+yYc9aJmdjXF2TZmMHOCJEBlMOsrGrNXEjGawMWj1fmRycXeC69jzg41lq3JxNmYA3u6C40Zf5GWPhg0b4ujRo6hW7cq3R9++fTFr1ixUr842YLCGAkmFEEL4BZ7aHCwhIQHx8fGYM2eO7ypfSjAMR1dx7dq1yMxku8JYp0zMdAghhBCFYcD92Azt0+EbNNMhhBBCCABXgkQLBop6cr8OzXQIIYTwCxRIWnwMw8CwYcMQGnolQvLSpUsYOXKk08zPqlWr3CpfTocQQgi/QE5H8Rk6dKjD/4MGDfJo+WXC6SgH50jqcHIsWy8yW+ezqnhhZVmNcGfJnwCuUDBLDOUJzAYEW+G0mkyMKRpOELtZwjemmmAR/Kws1jZ2jczqxNrN1AasLKuJCM3GOLsvmKKBhY2xa82u6TFaIyCV2Jlqgk3wskSBrD+Ygsns3FaT71l9pphNXrNrx5REzM6eKax8M3UdUxIVHE9lNeHbtczixYu9Wn6ZcDqEEEKIwtDmYKUfOR1CCCH8Ai2vlH7kdAghhPALNNNR+pFkVgghhBA+QTMdQggh/AItr5R+yoTTEYKiq1eYEsBsE1cWwW81LwXDaoQ7wPN0sClAFiWeTuxMUcMi3AHryg8Wkc+uxS/EztQgABBJ7Cw3CutX1gY2zszyUrBcFkzhwdrHxuVJYmeqBYCPNXZuplJhbWPHMwUOALBMDmwspxH7EWK3qu4AuBKGwR6gTMnmztQyu1+s5p9i9zyrk9l9x8ZywWtndv29gaeWVxISEhAQEICkpCQkJSV5oGYinzLhdAghhBC+Qtugew85HUIIIfwCLa+UfuR0CCGE8AukXin9SL0ihBBCCJ+gmQ4hhBB+gZZXSj9lwunIhfOUDItwZwoLsykdq5H3LIKfEUnsZvk72HtWc40wlYXVvDIAbzeL+GdhWKwcdo3M+ompVJidqXOY0oY9iMwUEAymaGB9zsblKWI3y3NiNTE1U/kwZUQlYo8xOQerE1NGHCX20xbtZsojdr3Z2LT6jDBThzHYWGPjhtWV9Tezm/UTO3dJf3FreaX0o+UVIYQQQviEMjHTIYQQQhSGlldKP3I6hBBC+AVaXin9yOkQQgjhF2imo/SjmA4hhBBC+IQyO9NR1L3/8zHLAWA1J4ZZTgJXsChwM1WG1TwJVlUq7HgzL5TlgGC/ENg52LVg05tmygurkfpWc0Gw62AGy5nC2seUIoxzxP6ryWdYf9QgdjY2mfqCKZXM8guxccPuLzaWmbqDKUjY9QGADJP3XGE1H9FFYmfPLcD6PW/1/mKYqVfYtSh4v5i1yxsYcH+ZRDMdvqHMOh1CCCHE1eTBujz86s8K7+OT5ZWsrCw0b94cNpsNe/fu9cUphRBCCK+ydetW9OjRAzExMbDZbPjggw88Uu6WLVuQkJCAsLAw1K9fH/Pnz3c6Jj09HUlJSahZsybCwsLQuHFjrF271iPn9yY+cToef/xxxMSYbRMkhBBCFA+jmC+rZGZm4qabbsLs2bOLX/nfOHjwILp374727dtjz549ePLJJzF69GisXLnSfkx2djbuuOMOHDp0CP/4xz/w3XffYeHChahVq5bH6uEtvL68sm7dOmzYsAErV67EunXrvH06IYQQ1yi+Vq9069YN3bp1o+9nZ2fj6aefxjvvvIP09HQ0adIEL774Ijp27Eg/M3/+fMTGxmLmzJkAgMaNG2PXrl2YPn06evfuDQB46623cPr0aXzxxRcIDr4SfVO3bl03WuB7vDrTcfz4cYwYMQJ///vfUb584eFyWVlZyMjIsL+EEEIIX3P191BGRgaysqzKB64wfPhwbNu2DcuXL8f+/fvRp08fdO3aFT/88AP9zPbt29G5c2cHW5cuXbBr1y5cvnwlTHjNmjVo06YNkpKSUL16dTRp0gRTp05Fbq47CRp8i9dmOgzDwLBhwzBy5EgkJibi0KFDhX5m2rRpmDx5spM9A86R1J4M+rE6nFhUPIvUZu6WWR4G5g0yOyuLReqnm5ybwXKBsGvBBhdzJ5ndTOVzhtjZrcci8qsRe7jFcgDeT2zcWM33wRQ1ZjczU1OwHDWRFs/N+pv1BWB9bDL1FFNrsLqa3e/s3Kwsds+za+2OGoqNf1YW61emnGHlmI1xRsGx7GtFiKcCSQuGA0ycOBGTJk2yVN5PP/2EZcuW4ddff7WX9+ijj2L9+vVYvHgxpk6d6vJzx44dQ/Xq1R1s1atXR05ODk6ePImaNWvif//7HzZt2oSBAwdi7dq1+OGHH5CUlIScnBw8++yzlurpayw7HZMmTXLpGFzNzp078cUXXyAjIwMTJkwoctkTJkxAcnIygCtrZYoDEUIIUVQ85XSkpqaiQoXfXfXQUDPxt2t2794NwzDQqFEjB3tWVhaioqIAAOHhv/+sGTRokD1g1GZzbIVhGA72vLw8VKtWDQsWLEBgYCASEhKQmpqKl19+2f+cjlGjRqFfv36mx9SrVw9TpkzBl19+6XSxEhMTMXDgQPztb39z+lxoaKj9+MBAd34TCCGEEMWjYsWKDk6HO+Tl5SEwMBApKSlO32f5zsbVas6KFSsCAGrUqIFjxxxzRqelpSEoKMjurNSsWRPBwcEO5TZu3BjHjh1DdnY2QkLM5odLFstOR3R0NKKj2aTs78yaNQtTpkyx/5+amoouXbpgxYoVaNWqldXTCiGEEKZ4KpA0ISEBAQEBSEpKQlJSklvltWjRArm5uUhLS0P79u1dHtOgQQMnW5s2bfDRRx852DZs2IDExER70Gi7du3w7rvvIi8vDwEBVxZiv//+e9SsWbNUOxyAF2M6YmNjHf7P9+yuu+461K5d21unFUIIcY3iKacjJSWlSDMd58+fx48//mj//+DBg9i7dy+qVKmCRo0aYeDAgRgyZAhmzJiBFi1a4OTJk9i0aROaNm2K7t27uyxz5MiRmD17NpKTkzFixAhs374dixYtwrJly+zH/OUvf8Hrr7+OMWPG4OGHH8YPP/yAqVOnYvTo0W623ndoR1IhhBB+ga93JN21axc6depk/z8/JnHo0KFYsmQJFi9ejClTpuCRRx7BkSNHEBUVhTZt2lCHAwDi4uKwdu1ajBs3DnPmzEFMTAxmzZpll8sCQJ06dbBhwwaMGzcOzZo1Q61atTBmzBiMHz/ejVb4FpuRH6FSysjMzLTPjgyGf6pXKpucw2rkPQtzCrN4vBneVq+cIHazyUKmOqlI7FbVK1EWywF4PzGlDVOpMFXGMWI/QmvE1SuxxB5J7GxcsoeI2X1qVb3C2s1yzljNzwTw68rGAduKqQqxuxOpxpRp7P5izwhfqFcKniMbwGu//X3+/Plix0m44urvilvgXh8DVxRY//nt7+uvv77YyyvCNWVipuM0nCvKJHrsgWK2IQn7AmaDl31pphN7TWJ3x+mwKqNkx7M2sC9Gs8+wLxb2RcTKYdeUyVbNPsMcmKrEzpwUNp7MPHXmLFwgdquJydzZwcZqYjyr95c7uwOwL8FTxM7aza4ds5vJeJksl41BZre6qu6Oc8b63GpSQ4bZjAEb/wXr6utdIzw101HU5RVhnTLhdAghhBCF4esdSYV1fJJ7RQghhBBCTocQQgi/wFMJ3xISEhAfH485c+b4rvLXCFpeEUII4RcopqP0o5kOIYQQQviEMjHTEYyiV5TJxcxksUwaxiKv04mdqRAiiJ1F7wPWvXUWdc/UBlYTVQFcIcMi9VmdmIrDqqQP4HJCq7I5q4oTM9KJ/TixM2knG7Nm6gurnCZ2di3YrxSmFDGTZrMxaFV9we4vdxLjsfoypVklYmfjkrXN7FlgVaLP1C7s2jFZfenPV+qMr/fpENbRTIcQQgi/QDEdpZ8yMdMhhBBC+ArFdHgPOR1CCCH8Au3TUfqR0yGEEMIvMOB+bIacDt8gp0MIIYRfoJmO0k+ZcDouwbmiLLKfKQHSTcpnEeVWo8ZZBD87N4saB3jCLatqDatev1mdGOxmNUus5QoWdc6SZwFAHWJnSfaYUom1mylzzCL7WX+wz7BrxFQqLK+MmeqDjRt2H7HVbKbisJp8DLCukGHnYOOGlWOWF4WND3Z/Wc1/wsYGe3aYwZ51TJHE2sb6yZ18MAXb7ak8ML4mISFBCd+8RJlwOoQQQojCKI7sVZuD+QY5HUIIIfwCLa+UfrRPhxBCCCF8gmY6hBBC+AWeWl4R3kNOhxBCCL9AyyulnzLhdByC8zoQixpnEd0sDwjA82uwPAxMTcEUAixK/xdaI/6ZasTO6uSO+oLB1BRM5cOuBSuH1dVMbcBUKiwnBlMhsAcOuw7u9F9VYreat4SN5WMm52ZqA6vjw6pSKZzWiKsp2DVlbWAw5YSZooI9ENln2Fhmqh12rdmzA+B9zs5tVf1jlh+HwdpX8Fkg9YooSJlwOoQQQojCkHql9COnQwghhF+g5ZXSj9QrQgghhPAJmukQQgjhF0i9UvqR0yGEEMIv0PJK6adMOB0/wzm/gpmiwRVWI98BnmeiOrFbVa8wdQfAlQhWvXEWCsXKOWdSllUlBztHFLEzZYnZtWPKI9Zudo2Yna0/mkXlszwuVvPjsHOzh6PZQzOC2JnKh91fTLXA8p+Y3afsM6w/WFmsHDY2zMY4u0asn9g9wdQ8VnO7AHz8M7vVc7BxZqbQYu0r2B/uPHeLg2Y6Sj+K6RBCCCGuIiEhAfHx8ZgzZ05JV8XvKBMzHUIIIURhSDJb+pHTIYQQwi9QTEfpR8srQgghhPAJmukQQgjhF2imo/RTJpyOADhHqLPVNha5zaLuAZ7DgEV7W81hwMoxU0BkEDtTtbA2sDqxG4wpAQCugGC5MlieGKuKGrNrx2CqAtYfrE5WVVJm52Y5U9KJ/RSxnyF2T+a5YPeR1bwerC8AroBgeCpvzgmTc7DrzdRNrN1MOcPsZrlX2JhlKimr146pVFhOJYAr7wrapV4RBdHyihBCCCF8QpmY6RBCCCEKw4D7MxZaXvENcjqEEEL4BcVxHOR0+AYtrwghhBDCJ2imQwghhF+gmY7Sj2Y6hBBC+AV5xXzlo23QvUeZmOnIg7OUk0k7mbzMrKFWPS8mxbOaUIlJ28zeY/JRJsVjskQmc61Ba8T7liXGY4ndKhI7CwBj8mHAerI5JvVjxzNZopmckCUUY+OGXTv2yyuW2M0kqFaThrGywi3azRKZsWvHpL+sTqz/mDSWSZEBfl+w9rH+Y/JoNjYiWYXApedszLJnB5PluiPpZ+04WeB/s6Rx3sDVd0VRufp+0zbo3kMzHUIIIYTwCWVipkMIIYQoDMV0lH7kdAghhPAL5HSUfrS8IoQQQgifoJkOIYQQfoGnAkmF9ygTTkc5FH0gsYhrsyh6FrFuNZEUU5awTjZLJsai5RnpxM76rSqxVzc5B0v4xlQtVqfRWD+ZxZCza8GUDuwaFYy6z4e1wSyRFesPq6qCehbLYSoOADhN7KyfmOqAqTLYtWPqH4BfC3YOluiOtY0plcyuHevDw8TO7mHWr+xam90rVpVH7HnGrhG71maJFtl9V/D562v1ipZXSj9aXhFCCCGETygTMx1CCCFEYWh5pfQjp0MIIYRfoOWV0o+WV4QQQgjhEzTTIYQQwi/Q8krpp0w4HYFwHkhsgLCIazP1ChukLPKaRaYzOzt3FVoj61HurK6sbawcls8B4NHy6cTOlAAsup612UzlY1VdwpQOLE/HeWJnbQCAaGJniiSmCmKqHaZeMbt2rJ9YDhk2lplqgSlLzNQr7BwsN8pxYmd1Yudm6iKA14mNG9ZuNpYbmZzbKlafgewZwcYAuyfM3iuYJ8lsTHoDLa+Ufry+vPLxxx+jVatWKFeuHKKjo3HPPfd4+5RCCCGuQTyVZbaobN26FT169EBMTAxsNhs++OCDYrcBALZs2YKEhASEhYWhfv36mD9/vsP7HTt2hM1mc3rdeeedHjm/N/Gq07Fy5UoMHjwYw4cPx759+7Bt2zYMGDDAm6cUQgghfEJmZiZuuukmzJ4922NlHjx4EN27d0f79u2xZ88ePPnkkxg9ejRWrlxpP2bVqlU4evSo/XXgwAEEBgaiT58+HquHt/Da8kpOTg7GjBmDl19+GQ888IDdfv3113vrlEIIIa5hfL280q1bN3Tr1o2+n52djaeffhrvvPMO0tPT0aRJE7z44ovo2LEj/cz8+fMRGxuLmTNnAgAaN26MXbt2Yfr06ejduzcAoEoVx8X55cuXo3z58te207F7924cOXIEAQEBaNGiBY4dO4bmzZtj+vTpuPHGGwv9vGH8PgRcDQarU2FmwUXsPXYONjit2s3awN6zWicGW981263R7D1XsH5lu8a6MxhZu9k5WBtYf1i1m52D1YnZWXyB1TYD1ne19FR/mI0Zq2VZHftW7yGzsjxld6efrI4bFkPG2m31XgGKfo2u/v/qZ7q3MFC6YjOGDx+OQ4cOYfny5YiJicHq1avRtWtXfPXVV2jYsKHLz2zfvh2dO3d2sHXp0gWLFi3C5cuXERzsHNW1aNEi9OvXDxUqmO3fXEowvMSyZcsMAEZsbKzxj3/8w9i1a5fRv39/Iyoqyjh16pTLz1y6dMk4e/ascfbsWePHH3/MHz966aWXXnqV8dfx48e98l1z/vx5j9c1NTXV/l109uxZ49KlS4XWA4CxevVq+/8//vijYbPZjCNHjjgcd/vttxsTJkyg5TRs2NB4/vnnHWzbtm2z16sgO3bsMAAYO3bsKLSOpQHLMR2TJk1yGcBy9WvXrl3Iy7vi4z711FPo3bs3EhISsHjxYthsNrz//vsuy542bRoqVaqESpUqoUGDBlarJoQQopSSlWWWzaV0ERMTY/8uqlSpEqZNm2a5jN27d8MwDDRq1Ajh4eH215YtW/DTTz8BgIN95MiR9s/abI7zxMZvs0QF7cCVWY4mTZrglltusVzHksDyjPaoUaPQr18/02Pq1auHc+fOAQDi4+Pt9tDQUNSvXx+HD7tOnzRhwgQkJycDAPLy8nDo0CG0aNECR44cQaVKlaxWtcySkZGBmJgYpKamomLFiiVdHZ9wLbYZULuvpXZfi202DAPHjx9HgwYNEBHBxOHFo3z58jh/nonbrZOVlYXAwECHL/jQUDOhtWvy8vIQGBiIlJQUBAY6LnqFh18R0O/du9duyx8TNWrUwLFjxxyOT0tLQ1BQEKKiohzsFy5cwPLly/Hcc89Zrl9JYdnpiI6ORnQ024ngdxISEhAaGorvvvsOt956KwDg8uXLOHToEOrWrevyM6GhoQ4Xt379+gCuXKAysVblIXJzr6yYVqhQ4Zpp97XYZkDtvpbafS22GYB91jsgwDtiSZvN5tH+9FRZLVq0QG5uLtLS0tC+fXuXx7ia0W/Tpg0++ugjB9uGDRuQmJjoFM/x3nvvISsrC4MGDfJInX2B1wJJK1asiJEjR2LixImoU6cO6tati5dffhkAykSErRBCCGHG+fPn8eOPP9r/P3jwIPbu3YsqVaqgUaNGGDhwIIYMGYIZM2agRYsWOHnyJDZt2oSmTZuie/fuLsscOXIkZs+ejeTkZIwYMQLbt2/HokWLsGzZMqdjFy1ahLvvvttpBqRU482AkezsbOORRx4xqlWrZkRERBh//OMfjQMHDhT582fPnjUAGGfPnvViLUsf12K7r8U2G4bafS21+1pss2H4d7s3b97sMhB16NChhmFc+Q589tlnjXr16hnBwcFGjRo1jF69ehn79+83LffTTz81WrRoYYSEhBj16tUz5s2b53TMd999ZwAwNmzY4I2meQ2vboMeHByM6dOnY/r06W59PjQ0FBMnTnRrPa0scy22+1psM6B2X0vtvhbbDPh3uzt27GgqBQ4ODsbkyZMxefJkS+V26NABu3fvNj2mUaNGPpEhexqbURZrLYQQQogyh1LbCyGEEMInyOkQQgghhE+Q0yGEEEIInyCnQwghhBA+oUw5HR9//DFatWqFcuXKITo6Gvfcc09JV8lnZGVloXnz5rDZbA672Pkjhw4dwgMPPIC4uDiUK1cO1113HSZOnIjsbJYCrWwyd+5cxMXFISwsDAkJCfjss89KukpeZdq0aWjZsiUiIiJQrVo13H333fjuu+9Kulo+Zdq0abDZbBg7dmxJV8XrHDlyBIMGDUJUVBTKly+P5s2bIyUlpaSrJUqYMuN0rFy5EoMHD8bw4cOxb98+bNu2DQMGDCjpavmMxx9/HDExMSVdDZ/w3//+F3l5eXjjjTfw9ddf49VXX8X8+fPx5JNPlnTVPMaKFSswduxYPPXUU9izZw/at2+Pbt260RQB/sCWLVuQlJSEL7/8Ehs3bkROTg46d+6MzMzMkq6aT9i5cycWLFiAZs2alXRVvM6ZM2fQrl07BAcHY926dfjmm28wY8YMREZGlnTVRElTstuEFI3Lly8btWrVMt58882SrkqJsHbtWuOGG24wvv76awOAsWfPnpKuks956aWXjLi4uJKuhse45ZZbjJEjRzrYbrjhBuOJJ54ooRr5nrS0NAOAsWXLlpKuitc5d+6c0bBhQ2Pjxo1Ghw4djDFjxpR0lbzK+PHjjVtvvbWkqyFKIWVipmP37t04cuQIAgIC0KJFC9SsWRPdunXD119/XdJV8zrHjx/HiBEj8Pe//x3ly5cv6eqUGGfPnkWVKlVKuhoeITs7GykpKejcubODvXPnzvjiiy9KqFa+5+zZswDgN9fVjKSkJNx555344x//WNJV8Qlr1qxBYmIi+vTpg2rVqqFFixZYuHBhSVdLlALKhNPxv//9DwAwadIkPP300/jnP/+JypUro0OHDjh9+nQJ1857GIaBYcOGYeTIkUhMTCzp6pQYP/30E15//XWH1M9lmZMnTyI3NxfVq1d3sFevXt0pu6S/YhgGkpOTceutt6JJkyYlXR2vsnz5cuzevdut9Ohllf/973+YN28eGjZsiH/9618YOXIkRo8ejaVLl5Z01UQJU6JOx6RJk2Cz2Uxfu3btsmcpfOqpp9C7d28kJCRg8eLFsNlseP/990uyCW5R1Ha//vrryMjIwIQJE0q6yh6hqO2+mtTUVHTt2hV9+vTBgw8+WEI19w5Xp84GrnwRF7T5K6NGjcL+/ftdJrHyJ3755ReMGTMGb7/9NsLCwkq6Oj4jLy8PN998M6ZOnYoWLVrg//7v/zBixAjMmzevpKsmShiv5l4pjFGjRqFfv36mx9SrVw/nzp0DAMTHx9vtoaGhqF+/fpkMvCtqu6dMmYIvv/zSKWdBYmIiBg4ciL/97W/erKbHKWq780lNTUWnTp3Qpk0bLFiwwMu18x3R0dEIDAx0mtVIS0tzmv3wRx5++GGsWbMGW7duRe3atUu6Ol4lJSUFaWlpSEhIsNtyc3OxdetWzJ49G1lZWQgMDCzBGnqHmjVrOjyvAaBx48ZYuXJlCdVIlBZK1OmIjo5GdHR0occlJCQgNDQU3333HW699VYAwOXLl3Ho0CHUrVvX29X0OEVt96xZszBlyhT7/6mpqejSpQtWrFiBVq1aebOKXqGo7QauyO06depkn9UKCCgTK4FFIiQkBAkJCdi4cSN69eplt2/cuBE9e/YswZp5F8Mw8PDDD2P16tX49NNPERcXV9JV8jq33347vvrqKwfb8OHDccMNN2D8+PF+6XAAQLt27Zzk0N9//32ZfF4Lz1KiTkdRqVixIkaOHImJEyeiTp06qFu3Ll5++WUAQJ8+fUq4dt4jNjbW4f/w8HAAwHXXXefXvxBTU1PRsWNHxMbGYvr06Thx4oT9vRo1apRgzTxHcnIyBg8ejMTERPtMzuHDh/0mbsUVSUlJePfdd/Hhhx8iIiLCPtNTqVIllCtXroRr5x0iIiKcYlYqVKiAqKgov45lGTduHNq2bYupU6fivvvuw3/+8x8sWLDAr2YshXuUCacDAF5++WUEBQVh8ODBuHjxIlq1aoVNmzahcuXKJV014WE2bNiAH3/8ET/++KOTc2X4SVLkvn374tSpU3juuedw9OhRNGnSBGvXrvXrX4L56/kdO3Z0sC9evBjDhg3zfYWE12jZsiVWr16NCRMm4LnnnkNcXBxmzpyJgQMHlnTVRAmj1PZCCCGE8An+s1AuhBBCiFKNnA4hhBBC+AQ5HUIIIYTwCXI6hBBCCOET5HQIIYQQwifI6RBCCCGET5DTIYQQQgifIKdDCCGEED5BTocQQgghfIKcDiG8yOXLl0u6CkIIUWqQ0yGEBdavX49bb70VkZGRiIqKwl133YWffvoJAHDo0CHYbDa899576NixI8LCwvD2228DAN566y3ceOONCA0NRc2aNTFq1Ch7mZMmTUJsbCxCQ0MRExOD0aNH29/Lzs7G448/jlq1aqFChQpo1aoVPv30U4c6bdu2DR06dED58uVRuXJldOnSBWfOnPF+ZwghhEXkdAhhgczMTCQnJ2Pnzp345JNPEBAQgF69eiEvL89+zPjx4zF69Gh8++236NKlC+bNm4ekpCT8+c9/xldffYU1a9agQYMGAIB//OMfePXVV/HGG2/ghx9+wAcffICmTZvayxo+fDi2bduG5cuXY//+/ejTpw+6du2KH374AQCwd+9e3H777bjxxhuxfft2fP755+jRowdyc3N92zFCCFEElPBNiGJw4sQJVKtWDV999RXCw8Pt2TTHjBljP6ZWrVoYPnw4pkyZ4vT5V155BW+88QYOHDiA4OBgh/d++uknNGzYEL/++itiYmLs9j/+8Y+45ZZbMHXqVAwYMACHDx/G559/7r1GCiGEh9BMhxAW+OmnnzBgwADUr18fFStWRFxcHADg8OHD9mMSExPtf6elpSE1NRW33367y/L69OmDixcvon79+hgxYgRWr16NnJwcAMDu3bthGAYaNWqE8PBw+2vLli32JZ38mQ4hhCgLBJV0BYQoS/To0QN16tTBwoULERMTg7y8PDRp0gTZ2dn2YypUqGD/u1y5cqbl1alTB9999x02btyIf//733jooYfw8ssvY8uWLcjLy0NgYCBSUlIQGBjo8Lnw8PAilS+EEKUJzXQIUUROnTqFb7/9Fk8//TRuv/12NG7cuNCAzYiICNSrVw+ffPIJPaZcuXL405/+hFmzZuHTTz/F9u3b8dVXX6FFixbIzc1FWloaGjRo4PCqUaMGAKBZs2amZQshRGlCMx1CFJHKlSsjKioKCxYsQM2aNXH48GE88cQThX5u0qRJGDlyJKpVq4Zu3brh3Llz2LZtGx5++GEsWbIEubm5aNWqFcqXL4+///3vKFeuHOrWrYuoqCgMHDgQQ4YMwYwZM9CiRQucPHkSmzZtQtOmTdG9e3dMmDABTZs2xUMPPYSRI0ciJCQEmzdvRp8+fRAdHe2DXhFCiKKjmQ4hikhAQACWL1+OlJQUNGnSBOPGjcPLL79c6OeGDh2KmTNnYu7cubjxxhtx11132dUnkZGRWLhwIdq1a2eftfjoo48QFRUFAFi8eDGGDBmCRx55BNdffz3+9Kc/YceOHahTpw4AoFGjRtiwYQP27duHW265BW3atMGHH36IoCD9nhBClD6kXhFCCCGET9BMhxBCCCF8gpwOIYQQQvgEOR1CCCGE8AlyOoQQQgjhE+R0CCGEEMInyOkQQgghhE+Q0yGEEEIInyCnQwghhBA+QU6HEEIIIXyCnA4hhBBC+AQ5HUIIIYTwCXI6hBBCCOET/h+EVPRtz+FDDAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "index = 20\n", + "webbpsf.display_psf(cube, ext=3, cube_slice=index, \n", + " # Note that currently the default plot title isn't very informative for datacube modes\n", + " # so we can specify a better title directly here:\n", + " title=f'MIRI MRS band {miri.band}, cube slice {index}, $\\\\lambda$={cube[0].header[\"WAVELN\"+str(index)]*1e6:.4} micron') " + ] + } + ], + "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.12.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/jwst_matching_psfs_to_data.ipynb b/docs/jwst_matching_psfs_to_data.ipynb index a610f1a0..6bf5fc9e 100644 --- a/docs/jwst_matching_psfs_to_data.ipynb +++ b/docs/jwst_matching_psfs_to_data.ipynb @@ -27,7 +27,9 @@ "id": "42ce87d3", "metadata": {}, "source": [ - "Often one wants to generate PSFs matched to some particular science dataset or file. The convenience function `webbpsf.setup_sim_to_match_data` helps with this, using the file's FITS header to set up a simulated instrument matched to the appropriate instrument setup and date of observation. " + "Often one wants to generate PSFs matched to some particular science dataset or file. The convenience function `webbpsf.setup_sim_to_match_file` helps with this, using the file's FITS header to set up a simulated instrument matched to the appropriate instrument setup and date of observation. \n", + "\n", + "Let's call that function, providing a filename of a data file already downloaded in the current directory:" ] }, { @@ -91,6 +93,7 @@ "metadata": {}, "source": [ "This function will:\n", + "\n", " * Create a webbpsf instrument object for the relevant instrument\n", " * Configure it to have the correct filter, detector, and other relevant instrument parameters for that science data file (e.g. coronagraph masks and so on). \n", " * Load the measured telescope mirror alignment data from the closest-in-time wavefront sensing visit to that science data. " @@ -160,6 +163,7 @@ "metadata": {}, "source": [ "The difference between the oversampled and detector-sampled output products is readily apparent. The distortion effects are generally more subtle: \n", + "\n", " * In this example case, note the slightly blurred softer look of the DET_DIST output compared to DET_SAMP, or of OVERDIST compared to OVERSAMP. This aspect of the simulation is a model for charge transfer physics and inter-pixel capacitance within the detector which result in crosstalk between adjacent pixels.\n", " * Also included as part of the distortion is a model for optical geometric distortions (including for instance slight differences between X and Y pixel scales, small rotations and skews of the detector pixel axes, the very-slightly-different position angles for each NIRCam detector, etc.). This attempts to forward-model the distortions which the \"drizzle\" pipeline algorithm corrects for, using the same astrometric calibration information for the instruments recorded in the [science instrument aperture file](https://pysiaf.readthedocs.io/en/latest/). \n", "\n", @@ -488,7 +492,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.7" + "version": "3.12.5" } }, "nbformat": 4, diff --git a/docs/relnotes.rst b/docs/relnotes.rst index 64433389..7f183333 100644 --- a/docs/relnotes.rst +++ b/docs/relnotes.rst @@ -24,6 +24,50 @@ Version History and Change Log ------------------------------- +Version 1.5.0 +============= +*2024 December* +WebbPSF IS BEING MIGRATED TO A NEW REPOSITORY: STPSF (Space Telescope PSF) +To reflect its broader support for Roman as well as James Webb, WebbPSF is being migrated to a new repository: STPSF (Space Telescope PSF). +This transition is being done in such a way as to maintain back-compatibility for existing code, and existing installations will continue to run as-is. +This release will be among the final releases to WebbPSF and the initial release to STPSF. + +This release requires new data webbpsf-data-1.5.0 https://stsci.box.com/s/oyoih9ibuaowksxzaia02ua4t6g1bbp1 +For more information on this new data please view https://github.com/spacetelescope/webbpsf/pull/937 + +**What's Changed** +* Add IFU+datacubes docs page. by @mperrin in https://github.com/spacetelescope/webbpsf/pull/917 + +* Update for Photutils 2.0 by @mperrin in https://github.com/spacetelescope/webbpsf/pull/925 + +* MNT: Use hash for Action workflow versions and update if needed by @pllim in https://github.com/spacetelescope/webbpsf/pull/916 + +* Dependabot updates PLUS pyproject upper limit for photutils. by @dependabot in https://github.com/spacetelescope/webbpsf/pull/926 + +* update ifu docs to describe about coord_system and pixel scale by @mperrin in https://github.com/spacetelescope/webbpsf/pull/923 + +* update photutils upper bound by @BradleySappington in https://github.com/spacetelescope/webbpsf/pull/935 + +* Add automatic data download by @WilliamJamieson in https://github.com/spacetelescope/webbpsf/pull/932 + +* switch a test case float comparison to np.allclose for robustness by @mperrin in https://github.com/spacetelescope/webbpsf/pull/934 + +* Add support for new sensing point and target phase map by @obi-wan76 in https://github.com/spacetelescope/webbpsf/pull/937 + +* Minor update to docs about the WFI optical model by @Skyhawk172 in https://github.com/spacetelescope/webbpsf/pull/936 + +* update monthly trending plot to allow shifted date ranges by @mperrin in https://github.com/spacetelescope/webbpsf/pull/943 + +* Readme advise by @BradleySappington in https://github.com/spacetelescope/webbpsf/pull/944 + +* update MAST query for OPDs to also retrieve OPDs from FP6 by @mperrin in https://github.com/spacetelescope/webbpsf/pull/945 + +## New Contributors +* @WilliamJamieson made their first contribution in https://github.com/spacetelescope/webbpsf/pull/932 + +**Full Changelog**: https://github.com/spacetelescope/webbpsf/compare/v1.4.0...1.5.0 + + Version 1.4.0 ============= diff --git a/docs/roman.rst b/docs/roman.rst index 44ae5af0..731b1917 100644 --- a/docs/roman.rst +++ b/docs/roman.rst @@ -22,6 +22,11 @@ Wide Field Instrument (WFI) The WFI model is based on the `Cycle 9 instrument reference information `_ from the Roman team at Goddard Space Flight Center (GSFC). The reported jitter for the Roman observatory is 0.012 arcsec per axis, per `GSFC `_. +.. note:: + + The current Roman WFI optical model was provided by Goddard Space Flight Center circa 2021 (the Cycle 9 reference data); a new optical model is currently being implemented in WebbPSF. + + To work with the WFI model, import and instantiate it just like any of the JWST instruments:: >>> from webbpsf import roman diff --git a/pyproject.toml b/pyproject.toml index 9c20212b..31257c27 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,7 +24,7 @@ dependencies = [ "astropy>=5.1.0", "photutils>=1.10.0", "poppy>=1.0.0", - "pysiaf>=0.19.1", + "pysiaf>=0.23.3", "synphot>=1.0.0", "astroquery>=0.4.6", ] diff --git a/requirements.txt b/requirements.txt index d437a7ec..5a16f3d2 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,10 +1,10 @@ -astropy>=6.1.2,<6.2.0 -ipython>=8.27.0,<8.28.0 +astropy>=6.1.2,<7.1.0 +ipython>=8.27.0,<8.31.0 matplotlib>=3.9.1,<3.10.0 numpy>=2.1.0,<2.2.0 -photutils>=1.13.0,<1.14.0 +photutils>=1.13.0,<2.1.0 poppy>=1.0.0 -pysiaf>=0.22.0,<0.23.0 +pysiaf>=0.23.3,<0.25.0 scipy>=1.14.0,<1.15.0 -synphot>=1.4.0,<1.5.0 +synphot>=1.4.0,<1.6.0 astroquery>=0.4.7,<0.5.0 diff --git a/webbpsf/__init__.py b/webbpsf/__init__.py index 76358dfe..6a99574a 100644 --- a/webbpsf/__init__.py +++ b/webbpsf/__init__.py @@ -40,7 +40,7 @@ class UnsupportedPythonError(Exception): # required. If changes to the code and data mean WebbPSF won't work # properly with an old data package, increment this version number. # (It's checked against $WEBBPSF_DATA/version.txt) -DATA_VERSION_MIN = (1, 4, 0) +DATA_VERSION_MIN = (1, 5, 0) class Conf(_config.ConfigNamespace): diff --git a/webbpsf/mast_wss.py b/webbpsf/mast_wss.py index 7ac1e992..def23240 100644 --- a/webbpsf/mast_wss.py +++ b/webbpsf/mast_wss.py @@ -357,7 +357,7 @@ def infer_pre_or_post_correction(row): return 'UNKNOWN' -def retrieve_mast_opd_table(aperture_list=['NRCA3_FP1'], verbose=False): +def retrieve_mast_opd_table(aperture_list=['NRCA3_FP1', 'NRCA1_FP6'], verbose=False): """Retrieve table of OPDs from MAST. Returns : Astropy table listing available OPDs and metadata such as dates and sensing type. diff --git a/webbpsf/tests/test_detectors.py b/webbpsf/tests/test_detectors.py index bff7e9fd..72c7315e 100644 --- a/webbpsf/tests/test_detectors.py +++ b/webbpsf/tests/test_detectors.py @@ -195,7 +195,7 @@ def test_ipc_basic_effect_on_psf_fwhm(): for extname in ['OVERSAMP', 'DET_SAMP']: fwhm_ipc = poppy.measure_fwhm(psf_withipc, ext=extname) fwhm_noipc = poppy.measure_fwhm(psf_noipc, ext=extname) - assert fwhm_ipc == fwhm_noipc, f'Adding IPC should not have any effect on the {extname} data.' + assert np.isclose(fwhm_ipc, fwhm_noipc), f'Adding IPC should not have any effect on the {extname} data.' print(f'test ok for {extname}') for extname in ['OVERDIST', 'DET_DIST']: diff --git a/webbpsf/tests/test_errorhandling.py b/webbpsf/tests/test_errorhandling.py index badd592d..3a3e8a74 100644 --- a/webbpsf/tests/test_errorhandling.py +++ b/webbpsf/tests/test_errorhandling.py @@ -4,6 +4,7 @@ import logging import os import os.path +from pathlib import Path import pytest @@ -48,7 +49,7 @@ def test_invalid_masks(): assert _exception_message_starts_with(excinfo, "Instrument NIRCam doesn't have a pupil mask called 'JUNK'.") -def test_get_webbpsf_data_path_invalid(monkeypatch): +def test_get_webbpsf_data_path_invalid(monkeypatch, tmp_path): real_env_webbpsf_path = os.getenv('WEBBPSF_PATH') real_conf_webbpsf_path = conf.WEBBPSF_PATH real_webbpsf_path = real_env_webbpsf_path or real_conf_webbpsf_path @@ -57,9 +58,19 @@ def test_get_webbpsf_data_path_invalid(monkeypatch): # config says to (and env var has been unset) monkeypatch.delenv('WEBBPSF_PATH') monkeypatch.setattr(conf, 'WEBBPSF_PATH', 'from_environment_variable') - with pytest.raises(EnvironmentError) as excinfo: + + # Patch the function that gets the home directory so we don't overwrite the + # what is on the system + def mockreturn(): + return Path(tmp_path) + + monkeypatch.setattr(Path, "home", mockreturn) + + with pytest.warns(UserWarning, match=r"Environment variable \$WEBBPSF_PATH is not set!\n.*") as excinfo: _ = utils.get_webbpsf_data_path() - assert 'Environment variable $WEBBPSF_PATH is not set!' in str(excinfo) + + # Check that the data was downloaded + assert any((tmp_path / "data" / "webbpsf-data").iterdir()) # Test that we can override the WEBBPSF_PATH setting here through # the config object even though the environment var is deleted diff --git a/webbpsf/tests/test_psfgrid.py b/webbpsf/tests/test_psfgrid.py index c0944582..301649dd 100644 --- a/webbpsf/tests/test_psfgrid.py +++ b/webbpsf/tests/test_psfgrid.py @@ -107,7 +107,7 @@ def test_all_detectors(): # Case 1: Shortwave -> check that only the SW detectors are applied for the SW filter nir.filter = shortfilt - grid1 = nir.psf_grid(all_detectors=True, num_psfs=1, add_distortion=False, fov_pixels=1, oversample=2, verbose=False) + grid1 = nir.psf_grid(all_detectors=True, num_psfs=1, add_distortion=False, fov_pixels=4, oversample=2, verbose=False) det_list = [] for hdu in grid1: det_list.append(hdu.meta['detector'][0]) @@ -117,7 +117,7 @@ def test_all_detectors(): # Case 2: Longwave -> check that only the LW detectors are applied for the LW filter nir.filter = longfilt - grid2 = nir.psf_grid(all_detectors=True, num_psfs=1, add_distortion=False, fov_pixels=1, oversample=2, verbose=False) + grid2 = nir.psf_grid(all_detectors=True, num_psfs=1, add_distortion=False, fov_pixels=4, oversample=2, verbose=False) det_list = [] for hdu in grid2: det_list.append(hdu.meta['detector'][0]) @@ -189,29 +189,29 @@ def test_nircam_errors(): # Shouldn't error - applying SW to SW and LW to LW nir.filter = longfilt nir.detector = longdet - nir.psf_grid(all_detectors=False, num_psfs=1, fov_pixels=1, detector_oversample=2, fft_oversample=2, verbose=False) + nir.psf_grid(all_detectors=False, num_psfs=1, fov_pixels=4, detector_oversample=2, fft_oversample=2, verbose=False) nir.filter = shortfilt nir.detector = shortdet - nir.psf_grid(all_detectors=False, num_psfs=1, fov_pixels=1, detector_oversample=2, fft_oversample=2, verbose=False) + nir.psf_grid(all_detectors=False, num_psfs=1, fov_pixels=4, detector_oversample=2, fft_oversample=2, verbose=False) # Should error - Bad filter/detector combination (LW filt to SW det) with pytest.raises(RuntimeError) as excinfo: # Errors inside calc_psf() call nir.filter = longfilt nir.detector = shortdet - nir.psf_grid(all_detectors=False, num_psfs=1, fov_pixels=1, verbose=False) # error + nir.psf_grid(all_detectors=False, num_psfs=1, fov_pixels=4, verbose=False) # error assert 'RuntimeError' in str(excinfo) # Should error - Bad filter/detector combination (SW filt to LW det) with pytest.raises(RuntimeError) as excinfo: # Errors inside calc_psf() call nir.filter = shortfilt nir.detector = longdet - nir.psf_grid(all_detectors=False, num_psfs=1, fov_pixels=1, verbose=False) # error + nir.psf_grid(all_detectors=False, num_psfs=1, fov_pixels=4, verbose=False) # error assert 'RuntimeError' in str(excinfo) # Should error - Bad num_psfs entry (must be a square number) with pytest.raises(ValueError) as excinfo: - nir.psf_grid(all_detectors=False, num_psfs=2, fov_pixels=1, verbose=False) # error + nir.psf_grid(all_detectors=False, num_psfs=2, fov_pixels=4, verbose=False) # error assert 'ValueError' in str(excinfo) diff --git a/webbpsf/trending.py b/webbpsf/trending.py index ae530817..d1d912d0 100644 --- a/webbpsf/trending.py +++ b/webbpsf/trending.py @@ -313,11 +313,13 @@ def wfe_histogram_plot( elif ote_only is True: opd_data = fits.getdata(full_file_path, ext=1) mask = opd_data != 0 + sensing_apername = fits.getheader(full_file_path, ext=0)['APERNAME'] + + # Get WSS Target Phase Map for the sensing aperture + # Note that the sensing maintenance program changed field point from NRC A3 to A1 around Dec 2024. + was_targ_file = webbpsf.utils.get_target_phase_map_filename(sensing_apername) + - # Get WSS Target Phase Map - was_targ_file = os.path.join( - webbpsf.utils.get_webbpsf_data_path(), 'NIRCam', 'OPD', 'wss_target_phase_fp1.fits' - ) target_1024 = astropy.io.fits.getdata(was_targ_file) target_256 = poppy.utils.krebin(target_1024, (256, 256)) / 16 wf_si = target_256 * mask # Nircam target phase map at FP1 @@ -607,6 +609,7 @@ def single_measurement_trending_plot( opd, opdhdu = _read_opd(filename) mask = opd != 0 opd[~mask] = np.nan + sensing_apername = opdhdu[0].header['APERNAME'] # hdr_rmswfe = opdhdu[1].header['RMS_WFE'] visit = opdhdu[0].header['OBS_ID'][0:12] @@ -628,18 +631,31 @@ def single_measurement_trending_plot( print(' Previous OPD is:', prev_filename) prev_opd, prev_opd_hdu = _read_opd(prev_filename) + prev_sensing_apername = prev_opd_hdu[0].header['APERNAME'] if subtract_target: if verbose: print(' Subtracting NIRCam SI WFE target phase map') - # Get WSS Target Phase Map - was_targ_file = os.path.join(webbpsf.utils.get_webbpsf_data_path(), 'NIRCam', 'OPD', 'wss_target_phase_fp1.fits') + # Get WSS Target Phase Map for the sensing aperture + # Note that the sensing maintenance program changed field point from NRC A3 to A1 around Dec 2024. + was_targ_file = webbpsf.utils.get_target_phase_map_filename(sensing_apername) + prev_was_targ_file = webbpsf.utils.get_target_phase_map_filename(prev_sensing_apername) + + target_1024 = astropy.io.fits.getdata(was_targ_file) target_256 = poppy.utils.krebin(target_1024, (256, 256)) / 16 # scale factor for rebinning w/out increasing values + if prev_was_targ_file != was_targ_file: + prev_target_1024 = astropy.io.fits.getdata(prev_was_targ_file) + prev_target_256 = poppy.utils.krebin(prev_target_1024, + (256, 256)) / 16 # scale factor for rebinning w/out increasing values + else: + prev_target_256 = target_256 + + opd -= target_256 - prev_opd -= target_256 + prev_opd -= prev_target_256 # Compute deltas and decompose deltatime = get_datetime_utc(opdhdu, return_as='astropy') - get_datetime_utc(prev_opd_hdu, return_as='astropy') @@ -961,9 +977,13 @@ def vprint(*text): opd, opdhdu = _read_opd(opdtable[which_opds_mask][0]['fileName']) mask = opd != 0 + sensing_apername = opdhdu[0].header['APERNAME'] + + # Get WSS Target Phase Map for the sensing aperture + # Note that the sensing maintenance program changed field point from NRC A3 to A1 around Dec 2024. + + was_targ_file = webbpsf.utils.get_target_phase_map_filename(sensing_apername) - # Get WSS Target Phase Map - was_targ_file = os.path.join(webbpsf.utils.get_webbpsf_data_path(), 'NIRCam', 'OPD', 'wss_target_phase_fp1.fits') target_1024 = astropy.io.fits.getdata(was_targ_file) target_256 = poppy.utils.krebin(target_1024, (256, 256)) @@ -1106,12 +1126,16 @@ def filter_opdtable_for_daterange(start_date, end_date, opdtable): return opdtable -def filter_opdtable_for_month(year, mon, opdtable): +def filter_opdtable_for_month(year, mon, opdtable, shift_dates_offset=None): """Filter existing opdtable for a given month This includes the last measurement in the prior month too (if applicable), so we can compute a delta to the first one """ start_date, end_date = get_month_start_end(year, mon) + if shift_dates_offset: + start_date += shift_dates_offset * u.day + end_date += shift_dates_offset * u.day + return filter_opdtable_for_daterange(start_date, end_date, opdtable) @@ -1129,7 +1153,7 @@ def get_opdtable_for_daterange(start_date, end_date): return opdtable -def get_opdtable_for_month(year, mon): +def get_opdtable_for_month(year, mon, **kwargs): """Return table of OPD measurements for a given month. retrieve the opd table and filter according to month/year designation """ @@ -1137,7 +1161,7 @@ def get_opdtable_for_month(year, mon): opdtable0 = webbpsf.mast_wss.retrieve_mast_opd_table() opdtable0 = webbpsf.mast_wss.deduplicate_opd_table(opdtable0) - opdtable = filter_opdtable_for_month(year, mon, opdtable0) + opdtable = filter_opdtable_for_month(year, mon, opdtable0, **kwargs) return opdtable @@ -1196,7 +1220,8 @@ def get_dates_for_pid(pid, project='jwst'): return -def monthly_trending_plot(year, month, verbose=True, instrument='NIRCam', filter='F200W', vmax=200, pid=None, opdtable=None): +def monthly_trending_plot(year, month, verbose=True, instrument='NIRCam', filter='F200W', vmax=200, pid=None, opdtable=None, + shift_dates_offset=None): """Make monthly trending plot showing OPDs, mirror moves, RMS WFE, and the resulting PSF EEs year, month : integers @@ -1209,6 +1234,8 @@ def monthly_trending_plot(year, month, verbose=True, instrument='NIRCam', filter Image display vmax for OPDs, given here in units of nanometers. opdtable : astropy.table.Table Table of available OPDs, Default None: as returned by retrieve_mast_opd_table() + shift_dates_offset: int + number of dates to shift later the time period, for showing trends around the month divisions """ def vprint(*text): @@ -1217,6 +1244,10 @@ def vprint(*text): start_date, end_date = get_month_start_end(year, month) + if shift_dates_offset: + start_date += shift_dates_offset * u.day + end_date += shift_dates_offset * u.day + # Look up wavefront sensing and mirror move corrections for that month if pid: pid_dates = get_dates_for_pid(pid) @@ -1228,11 +1259,11 @@ def vprint(*text): try: if opdtable is None: vprint('obtaining opdtable for month') - opdtable = get_opdtable_for_month(year, month) + opdtable = get_opdtable_for_month(year, month, shift_dates_offset=shift_dates_offset) else: vprint('filtering opdtable for month') opdtable = check_colnames(opdtable) - opdtable = filter_opdtable_for_month(year, month, opdtable) + opdtable = filter_opdtable_for_month(year, month, opdtable, shift_dates_offset=shift_dates_offset) except ValueError as e: print(e) raise @@ -1342,7 +1373,8 @@ def basic_show_image(image, ax, vmax=0.3, nanmask=1): fs = 14 # Font size for axes labels - fig.suptitle(f'WF Trending for {year}-{month:02d}', fontsize=fs * 1.5, fontweight='bold') + title_extra = f", plus {shift_dates_offset} days" if shift_dates_offset else "" + fig.suptitle(f'WF Trending for {year}-{month:02d}{title_extra}', fontsize=fs * 1.5, fontweight='bold') # Plot 1: Wavefront Error @@ -1518,7 +1550,10 @@ def basic_show_image(image, ax, vmax=0.3, nanmask=1): for i in range(3): im_axes[i, j].set_visible(False) - outname = f'wf_trending_{year}-{month:02d}.pdf' + if shift_dates_offset: + outname = f'wf_trending_{year}-{month:02d}+{shift_dates_offset}days.pdf' + else: + outname = f'wf_trending_{year}-{month:02d}.pdf' plt.savefig(outname, dpi=200, bbox_inches='tight') vprint(f'Saved to {outname}') diff --git a/webbpsf/utils.py b/webbpsf/utils.py index 29c6ddc7..f4c0ea47 100644 --- a/webbpsf/utils.py +++ b/webbpsf/utils.py @@ -180,6 +180,39 @@ def setup_logging(level='INFO', filename=None): """ +def auto_download_webbpsf_data(): + import os + import tarfile + from pathlib import Path + from tempfile import TemporaryDirectory + from urllib.request import urlretrieve + + # Create a default directory for the data files + default_path = Path.home() / "data" / "webbpsf-data" + default_path.mkdir(parents=True, exist_ok=True) + + os.environ["WEBBPSF_PATH"] = str(default_path) + + # Download the data files if the directory is empty + if not any(default_path.iterdir()): + warnings.warn("WebbPSF data files not found in default location, attempting to download them now...") + + with TemporaryDirectory() as tmpdir: + # Download the data files to a temporary directory + url = "https://stsci.box.com/shared/static/qxpiaxsjwo15ml6m4pkhtk36c9jgj70k.gz" + filename = Path(tmpdir) / "webbpsf-data-LATEST.tar.gz" + urlretrieve(url, filename) + + # Extract the tarball + with tarfile.open(filename, "r:gz") as tar: + tar.extractall(default_path.parent, filter="fully_trusted") + + if not any(default_path.iterdir()): + raise IOError(f"Failed to get and extract WebbPSF data files to {default_path}") + + return default_path + + def get_webbpsf_data_path(data_version_min=None, return_version=False): """Get the WebbPSF data path @@ -201,7 +234,12 @@ def get_webbpsf_data_path(data_version_min=None, return_version=False): if path_from_config == 'from_environment_variable': path = os.getenv('WEBBPSF_PATH') if path is None: - raise EnvironmentError(f'Environment variable $WEBBPSF_PATH is not set!\n{MISSING_WEBBPSF_DATA_MESSAGE}') + message = ( + 'Environment variable $WEBBPSF_PATH is not set!\n' + f'{MISSING_WEBBPSF_DATA_MESSAGE} searching default location..s' + ) + warnings.warn(message) + path = auto_download_webbpsf_data() else: path = path_from_config @@ -1010,3 +1048,25 @@ def label_wavelength(nwavelengths, wavelength_slices): else: raise ValueError('Maximum number of wavelengths exceeded. ' 'Cannot be more than 10,000.') return label + + +def get_target_phase_map_filename(apername): + """Get WSS Target Phase Map for the specified aperture + Note that the sensing maintenance program changed field point from NRC A3 to A1 around Dec 2024. + """ + path = os.getenv('WEBBPSF_PATH') + if apername == 'NRCA3_FP1': + fn = 'wss_target_phase_fp1.fits' + elif apername == 'NRCA1_FP6': + fn = 'wss_target_phase_fp6.fits' + else: + raise ValueError(f"Target phase map not available for aperture = {apername}") + + was_targ_file = os.path.join( + get_webbpsf_data_path(), 'NIRCam', 'OPD', fn) + + if not os.path.exists(was_targ_file): + raise ValueError("File wss_target_phase_{}.fits, \ + not found under {}.".format(apername.split('_')[1].lower(),path)) + + return was_targ_file \ No newline at end of file diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index d043bd7c..40efa4b5 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -1053,7 +1053,7 @@ def aperturename(self, value): if not has_custom_pixelscale: self.pixelscale = self._get_pixelscale_from_apername(self._aperturename) _log.debug( - f'Pixelscale updated to {self.pixelscale}', + f'Pixelscale updated to {self.pixelscale}' + f'based on average X+Y SciScale at SIAF aperture {self._aperturename}' ) @@ -1697,12 +1697,12 @@ def load_wss_opd(self, filename, output_path=None, backout_si_wfe=True, verbose= axes[0, 0].set_xlabel(f'RMS: {utils.rms(opdhdu[0].data*1e9, ote_pupil_mask):.2f} nm rms') if backout_si_wfe: - if verbose: - print('Backing out SI WFE and OTE field dependence at the WF sensing field point') - # Check which field point was used for sensing sensing_apername = opdhdu[0].header['APERNAME'] + if verbose: + print(f'Backing out SI WFE and OTE field dependence at the WF sensing field point ({sensing_apername})') + # Create a temporary instance of an instrument, for the sensng instrument and field point, # in order to model and extract the SI WFE and OTE field dep WFE at the sensing field point. @@ -1720,10 +1720,13 @@ def load_wss_opd(self, filename, output_path=None, backout_si_wfe=True, verbose= # note that there is a slight focus offset between the two wavelengths, due to NIRCam's refractive design # Set to the sensing aperture, and retrieve the OPD there sensing_inst.set_position_from_aperture_name(sensing_apername) - # special case: for the main sensing point FP1, we use the official WAS target phase map, rather than the - # WebbPSF-internal SI WFE model. - was_targ_file = os.path.join(utils.get_webbpsf_data_path(), 'NIRCam', 'OPD', 'wss_target_phase_fp1.fits') - if sensing_apername == 'NRCA3_FP1' and os.path.exists(was_targ_file): + # special case: for the main sensing points FP1 or FP6, we use the official WAS target phase map, + # rather than the WebbPSF-internal SI WFE model. + + # Select correct target phase map based on sensing field point. + # Note that the sensing maintenance program changed field point from NRC A3 to A1 around Dec 2024. + if sensing_apername in ['NRCA3_FP1', 'NRCA1_FP6']: + was_targ_file = utils.get_target_phase_map_filename(sensing_apername) sensing_fp_si_wfe = poppy.FITSOpticalElement(opd=was_targ_file).opd else: sensing_fp_si_wfe = sensing_inst.get_wfe('si')