diff --git a/docs/index.rst b/docs/index.rst index 1ca2437e..b8bddd11 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -58,6 +58,7 @@ Changelog: - discuss MCMC autocorrelation in MCMC tutorial (@michaelkmpoon) - add time warning if OFTI doesn't accept an orbit in first 60 s (@michaelkmpoon) +- implementation of Hipparcos-Gaia catalog of accelerations fitting! (@semaphoreP) **2.2.2 (2023-06-30)** diff --git a/docs/tutorials.rst b/docs/tutorials.rst index 4a1acd01..a001d4b6 100644 --- a/docs/tutorials.rst +++ b/docs/tutorials.rst @@ -55,6 +55,7 @@ us if you are still confused). tutorials/Using_nonOrbitize_Posteriors_as_Priors.ipynb tutorials/Changing_bases_tutorial.ipynb tutorials/Hipparcos_IAD.ipynb + tutorials/HGCA_tutorial.ipynb diff --git a/docs/tutorials/HGCA_tutorial.ipynb b/docs/tutorials/HGCA_tutorial.ipynb new file mode 100644 index 00000000..cf52098e --- /dev/null +++ b/docs/tutorials/HGCA_tutorial.ipynb @@ -0,0 +1,239 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Fitting the Hipparcos-Gaia Catalog of Acclerations (HGCA)\n", + "\n", + "Jason Wang (2023)\n", + "\n", + "We will demonstrate how to fit the stellar absolute astrometry from the [Hipparcos-Gaia Catalog of Accelerations (HGCA; Brandt 2021)](https://ui.adsabs.harvard.edu/abs/2021ApJS..254...42B/abstract) to help constrain the mass and orbital parameters of the companion. The implementation of fitting the HGCA here is based on [Brandt et al. 2019](https://ui.adsabs.harvard.edu/abs/2019AJ....158..140B/abstract), but utilizes the DR3 version of the HGCA. \n", + "\n", + "## Difference to fitting IAD directly\n", + "\n", + "This method is an alternative to fitting the Hipparcos IAD and Gaia astrometry. Instead of fitting the individual epochs of Hipparcos data, which also includes needing to fit for the position, proper motion, and parallax of the star, we are only fitting for two differential proper motions: the proper motion difference between Hipparcos and the change in position between Hipparcos and Gaia; the proper motion difference between Gaia and the change in position between Hipparcos and Gaia. This simplifies the fit, but also ignores any detectable curvature in the stellar astrometry seen by Hipparcos. For planet companion cases, the acceleration should be well within the noise of the Hipparcos measurements. \n", + "\n", + "The benefit of using this technique is that the errors should be more robust. The HGCA catalog inflates the error bars from the Hipparcos and Gaia measurements on a global scale to match the true observed scatter in the data. Additionally, there are bad epochs in the Hipparcos IAD that may not be removed. \n", + "\n", + "We encourage users to try both to see how similar or different the results are, as the two methods utilize the same underlying data. \n", + "\n", + "\n", + "## Obtain the necessary data\n", + "\n", + "While `orbitize!` will automatically download the HGCA catalog, you need to obtain two other data files to reconstruct the Hipparcos and Gaia observations for your star. These files tell us when and in what orientation did Hipparcos and Gaia take data of the star for the forward modeling that happens in `orbitize!`. \n", + "\n", + " 1. The Hipparcos IAD file of the star. Follow the section in the [Hipparocs IAD tutorial](https://orbitize.readthedocs.io/en/latest/tutorials/Hipparcos_IAD.html#Part-1:-Obtaining-the-IAD) on how to obtain this file for your star of interest.\n", + " 2. The anticipated Gaia epochs and scan directions in a CSV file that is obtained from the [Gaia Observation Forecast Tool (GOST)](https://gaia.esac.esa.int/gost/). For GOST, after entering the target name and resolving its coordinates, use 2014-07-26T00:00:00 as the start date. For the end date, use 2016-05-23T00:00:00 for DR2 and 2017-05-28T00:00:00 for EDR3. You probably will be fitting the EDR3. The output CSV file is all you need. \n", + "\n", + "## Setup the `HGCALogProb` Object\n", + "\n", + "The user just needs to setup the `orbitize.gaia.HGCALogProb` object, which will handle the rest of the HGCA modeling under the hood. To setup the `HGCALogProb` object, we also need to create an `orbitize.hipparcos.HipparcosLogProb` instance, but it will not actually be used in the fitting. It is simply used to handle reading in the Hipparcos IAD to maximize code reuse. \n", + "\n", + "In the example below, we setup the HGCA fitting for beta Pictoris (HIP 27321).\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using HGCA catalog stored in /home/sblunt/Projects/orbitize/orbitize/example_data/HGCA_vEDR3.fits\n" + ] + } + ], + "source": [ + "import os\n", + "from orbitize import DATADIR, hipparcos, gaia\n", + "\n", + "# the necessary input data for beta Pic is part of the orbitize! example data!\n", + "iad_filepath = os.path.join(DATADIR, \"HIP027321.d\")\n", + "gost_filepath = os.path.join(DATADIR, \"gaia_edr3_betpic_epochs.csv\")\n", + "\n", + "# Create the HGCA and helper Hipparcos object\n", + "# we're just going to assume one planet for simplicity\n", + "hipparcos_lnprob = hipparcos.HipparcosLogProb(iad_filepath, 27321, 1)\n", + "hgca_lnprob = gaia.HGCALogProb(27321, hipparcos_lnprob, gost_filepath)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup orbit fit\n", + "\n", + "After you do this step, the rest of the orbit fit setup is basically the same as a standard `orbitize!` orbit fit. The only differences is that you need to pass the `HGCALogProb` object into the system class (but not the `HipparcosLogProb` because the Hipparcos data is already being handled in HGCA). We also recommend you explicitly set up the system to fit for the dynamical mass of the companions in the system (this should happen automatically with the inclusion of HGCA data, but we recommend being explicit so it is clear what you are fitting for).\n", + "\n", + "We will also note that that the default prior on the companion mass is a log-uniform prior between 1e-6 and 2 solar masses. If you want a more restricted mass, now is also the time to adjust that. We present one example here, but refer to the Modifying Priors tutorial for further details. " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "from orbitize import read_input, system, priors\n", + "import astropy.units as u\n", + "\n", + "# read in relative astrometry\n", + "astrometry_filepath = os.path.join(DATADIR, \"betaPic.csv\")\n", + "data_table = read_input.read_file(astrometry_filepath)\n", + "\n", + "# set up the system, passing in hgca_lnprob and setting it fit dynamical mass\n", + "stellar_mass = 1.75\n", + "stellar_mass_err = 0.05\n", + "plx = 51.44\n", + "plx_err = 0.12\n", + "\n", + "this_system = system.System(\n", + " 1,\n", + " data_table,\n", + " stellar_mass,\n", + " plx,\n", + " mass_err=stellar_mass_err,\n", + " plx_err=plx_err,\n", + " fit_secondary_mass=True,\n", + " gaia=hgca_lnprob,\n", + ")\n", + "\n", + "# adjust the prior on mass to be log uniform between 1 and 50 Jupiter masses\n", + "mjup = u.Mjup.to(u.Msun)\n", + "this_system.sys_priors[this_system.param_idx[\"m1\"]] = priors.LogUniformPrior(\n", + " mjup, 50 * mjup\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run orbitize! \n", + "\n", + "From here onwards, everything is the same as an regular `orbitize` fit, so we refer to reader to the other tutorials. We recommend using the MCMC sampler, since the orbits are generally too constrained for OFTI. To pick the right number of temperatures, walkers, steps for MCMC, check out papers that performed similar fits (e.g., Section 3.1 of [GRAVITY Collaboration et al. (2020)](https://ui.adsabs.harvard.edu/abs/2020A%26A...633A.110G/abstract) and Section 3.1 of [Hinkley et al. (2023)](https://ui.adsabs.harvard.edu/abs/2023A%26A...671L...5H/abstract))." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Starting Burn in\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/sblunt/Projects/orbitize/orbitize/priors.py:354: RuntimeWarning: invalid value encountered in log\n", + " lnprob = -np.log((element_array*normalizer))\n", + "/home/sblunt/Projects/orbitize/orbitize/priors.py:463: RuntimeWarning: invalid value encountered in log\n", + " lnprob = np.log(np.sin(element_array)/normalization)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Burn in complete. Sampling posterior now.\n", + "100/100 steps completed\n", + "Run complete\n" + ] + } + ], + "source": [ + "from orbitize import sampler\n", + "\n", + "# MCMC parameters\n", + "# for demonstration purposes only. You will need to increase these likely\n", + "n_temps = 2\n", + "n_walkers = 50\n", + "n_threads = 1\n", + "burn_steps = 1\n", + "total_orbits = 100 * n_walkers\n", + "\n", + "# create the sampler, run it, and save posteriors\n", + "this_sampler = sampler.MCMC(this_system, n_temps, n_walkers, n_threads)\n", + "\n", + "this_sampler.run_sampler(total_orbits, burn_steps=burn_steps, thin=10)\n", + "\n", + "this_sampler.results.save_results(\"demo_hgca.hdf5\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot proper motions from the orbit fit over the HGCA observations\n", + "\n", + "William Balmer (2023)\n", + "\n", + "Now, you can plot the result of your fit using the `orbitize.plot.plot_propermotion` function directly, or it's wrapper within the `sampler.results` object. The function is similar to `orbitize.plot.plot_orbits` that you may already be familiar with. These plots show the H-G proper motion (with a large x-axis error bar) in addition to the two \"real\" proper motion measurements from Hipparchos and Gaia.\n", + "\n", + "Note that the HGCALogprob fits for the time-averaged proper motions, and not the instantaneous proper motions that are shown in this visualization. This is a subtle point, but important, because the data points show over the orbits here are not exactly what is driving the MCMC solver towards the most likely orbits." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Important Note of Caution: the orbitize! implementation of the HGCA \n", + " fits for the time-averaged proper motions, and not the instantaneous proper \n", + " motions that are being plotted here. This plot is provided only for the \n", + " purpose of an approximate check on the fit.\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA44AAAGPCAYAAAATNI4YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAD4X0lEQVR4nOzdeYBcZZX4/e/da+s1naSTkJBAAmFHQEJQR5FIWGWTbUTWHzgossQBDS8GBJyAjoCyioOKC4PiCIPoBGMURQlbMIIKUTB70un0Vvvdn/ePW13QpBNISJNOcj5ODelbt566t6qXOvc8zzmaUkohhBBCCCGEEEJshL6tD0AIIYQQQgghxPAmgaMQQgghhBBCiE2SwFEIIYQQQgghxCZJ4CiEEEIIIYQQYpMkcBRCCCGEEEIIsUkSOAohhBBCCCGE2CQJHIUQQgghhBBCbJIEjkIIIYQQQgghNsnc1gcghBBCCCGEEP1c18X3/SEb37ZtUqnUkI2/o5LAUQghhBBCCDEsuK7L+PG70tXVOWTP0d7eztKlSyV43EwSOAohhBBCCCGGBd/36erq5P/+73my2YatPn65XOSYY96P7/vvOHC86667+NrXvkZHRwcHHHAAd9xxB4ceeuhG93/44Yf50pe+xLJly5gyZQq33HILxx577KD7/tu//Rvf+ta3uO2227jiiivq2ydOnMjy5csH7Dt37ly++MUvvqNjHgoSOAohhBBCCCGGlUwmRzab2+rjKqU2a/8f//jHzJo1i3vvvZdp06Zx++23M3PmTJYsWcKoUaM22P/pp5/mrLPOYu7cuRx//PE8+OCDnHTSSbz44ovsu+++A/Z95JFHeOaZZxg7duygz33DDTdw0UUX1b9uaNj6gfTmkOI4QgghhBBCCDGIW2+9lYsuuojzzz+fvffem3vvvZdMJsN3vvOdQff/xje+wdFHH81VV13FXnvtxY033shBBx3EnXfeOWC/1atX87nPfY4f/ehHWJY16FgNDQ20t7fXb9lsdquf3+aQwFEIIYQQQgixUykUCgNunudtsI/v+yxatIgZM2bUt+m6zowZM1i4cOGg4y5cuHDA/gAzZ84csH8cx3zqU5/iqquuYp999tnoMd58882MGDGC973vfXzta18jDMPNPc2tSqaqCiGEEEIIIYYVpZLbUIwLMH78+AHbr7vuOq6//voB27q6uoiiiNGjRw/YPnr0aF599dVBx+/o6Bh0/46OjvrXt9xyC6Zpctlll230OC+77DIOOuggWltbefrpp5k9ezZr167l1ltvfbtTHDISOAohhBBCCCGGFaXUZq9HfKfjAqxcuZLGxsb6dsdxtvpzDWbRokV84xvf4MUXX0TTtI3uN2vWrPq/999/f2zb5tOf/jRz5859z471rWSqqhBCCCGEEGJYUUN4A2hsbBxwGywYa2trwzAM1q1bN2D7unXraG9vH/S429vbN7n/U089RWdnJxMmTMA0TUzTZPny5Xz+859n4sSJG309pk2bRhiGLFu2bKP7DDUJHIUQQgghhBDDy1BHju+AbdscfPDBLFiwoL4tjmMWLFjA9OnTB33M9OnTB+wPMH/+/Pr+n/rUp3jppZdYvHhx/TZ27FiuuuoqnnjiiY0ey+LFi9F1fdBKru8VmaoqhBBCCCGEEIOYNWsW5557LocccgiHHnoot99+O+VymfPPPx+Ac845h3HjxjF37lwALr/8cj784Q/z9a9/neOOO46HHnqIF154gfvuuw+AESNGMGLEiAHPYVkW7e3t7LnnnkBSYOfZZ5/liCOOoKGhgYULF3LllVdy9tln09LS8h6e/UASOAohhBBCCCGGlaFe4/hOnXHGGaxfv545c+bQ0dHBgQceyLx58+oFcFasWIGuvzGJ8/DDD+fBBx/k2muv5ZprrmHKlCk8+uijG/Rw3BTHcXjooYe4/vrr8TyPSZMmceWVVw5Y97gtaGoo3hEhhBBCCCGE2EyFQoGmpiYWLPgr2ezWb3hfLhc58sh9yOfzA4rjiLcnaxyFEEIIIYQQQmySTFUVQgghhBBCDCuq9r+hGFdsGck4CiGEEEIIIYTYJMk4CiGEEEIIIYYVpZLbUIwrtoxkHIUQQgghhBBCbJJkHIUQQgghhBDDimQchx8JHIUQQgghhBDDynDp4yjeIFNVhRBCCCGEEEJskmQchRBCCCGEEMOKTFUdfiTjKIQQQgghhBBikyTjKIQQQgghhBhWZI3j8COBoxBCCCGEEGJYkamqw49MVRVCCCGEEEIIsUkSOAohhBBCCCGE2CSZqiqEEEIIIYQYVmSN4/AjgeNWFMcxa9asoaGhAU3TtvXhCCF2MEopisUiY8eORdd3zAkj8ntUCDHUdobfpTsCWeM4/EjguBWtWbOG8ePHb+vDEELs4FauXMkuu+yyrQ9jSMjvUSHEe2VH/l0qxFCQwHEramhoAJJfRI2Njdv4aIQQO5pCocD48ePrv2t2RPJ7VAgx1HaG36VCDAUJHLei/mlVjY2N8oFHCDFkduQpnPJ7VAjxXtmRf5fuCGSN4/AjgaMQQgghhBBieFG121CMK7aIBI5CCCGEEEKIYUXV/jcU44otI6WkhBBiB/P73/+eE044gbFjx6JpGo8++uiA+5VSzJkzhzFjxpBOp5kxYwb/+Mc/BuzT09PDJz/5SRobG2lububCCy+kVCpt8nld1+Wzn/0sI0aMIJfLceqpp7Ju3bqtfXpCCCGE2AYk4yiEEDuYcrnMAQccwAUXXMApp5yywf1f/epX+eY3v8kDDzzApEmT+NKXvsTMmTP529/+RiqVAuCTn/wka9euZf78+QRBwPnnn8/FF1/Mgw8+uNHnvfLKK/nFL37Bww8/TFNTE5deeimnnHIKf/zjH4fsXIXYUUVRRBAE2/owtkuGYWCapqxhFGIrk8BRCCF2MMcccwzHHHPMoPcppbj99tu59tprOfHEEwH4/ve/z+jRo3n00Uc588wzeeWVV5g3bx7PP/88hxxyCAB33HEHxx57LP/5n//J2LFjNxg3n89z//338+CDD/LRj34UgO9+97vstddePPPMMxx22GGDHo/neXieV/+6UCi8q3MXYkdQKpVYtWqVFPF4FzKZDGPGjMG27W19KGILSR/H4UcCRyGE2IksXbqUjo4OZsyYUd/W1NTEtGnTWLhwIWeeeSYLFy6kubm5HjQCzJgxA13XefbZZzn55JM3GHfRokUEQTBg3KlTpzJhwgQWLly40cBx7ty5fPnLX96KZyjE9i2KIlatWkUmk2HkyJGSNdtMSil832f9+vUsXbqUKVOmoOuyMmt7pGKFiodgjeMQjLmzkJ8kIYTYiXR0dAAwevToAdtHjx5dv6+jo4NRo0YNuN80TVpbW+v7DDaubds0NzdvdNzBzJ49m3w+X7+tXLlyc09JiB1KEAQopRg5ciTpdJpUKiW3zbil02mampoYO3YscRzj+/62fkvFDuCuu+5i4sSJpFIppk2bxnPPPbfJ/R9++GGmTp1KKpViv/3245e//OVG9/23f/s3NE3j9ttvH7B9S2oNDDUJHIUQQmwzjuPUezZK70Yh3iCZxndHsozbv/6pqkNx2xw//vGPmTVrFtdddx0vvvgiBxxwADNnzqSzs3PQ/Z9++mnOOussLrzwQv70pz9x0kkncdJJJ/GXv/xlg30feeQRnnnmmUGXgHzyk5/kr3/9K/Pnz+fxxx/n97//PRdffPHmHfxWJj9VQgixE2lvbwfYoNrpunXr6ve1t7dv8AcxDEN6enrq+ww2ru/79PX1bXRcIcR7xwsUF93dw0V39+AFMjVPiC116623ctFFF3H++eez9957c++995LJZPjOd74z6P7f+MY3OProo7nqqqvYa6+9uPHGGznooIO48847B+y3evVqPve5z/GjH/0Iy7IG3Ndfa+C//uu/mDZtGh/84Ae54447eOihh1izZs2QnevbkcBRCCF2IpMmTaK9vZ0FCxbUtxUKBZ599lmmT58OwPTp0+nr62PRokX1fX7zm98QxzHTpk0bdNyDDz4Yy7IGjLtkyRJWrFhRH1cIIYR4p5LsoBqCWzJ+oVAYcHtzobZ+vu+zaNGiAev3dV1nxowZLFy4cNDjXrhw4YD9AWbOnDlg/ziO+dSnPsVVV13FPvvsM+gYm6o1sK1IcZwhEEXxtj4EIcROrFQq8dprr9W/Xrp0KYsXL6a1tZUJEyZwxRVXcNNNNzFlypR6O46xY8dy0kknAbDXXntx9NFHc9FFF3HvvfcSBAGXXnopZ555Zn06zerVqznyyCP5/ve/z6GHHkpTUxMXXnghs2bNorW1lcbGRj73uc8xffr0jRbG2RTPC7fKayHEzmbFihV0dXXhh4quFUmV4j/9qRHbTKa+trW1MWHChG12fN/73ve44oorNpidIMR7bfz48QO+vu6667j++usHbOvq6iKKokHrArz66quDjtvR0bHJOgIAt9xyC6Zpctlll210jM2tNfBekMBxCIShBI5CiG3nhRde4Igjjqh/PWvWLADOPfdcvve973H11VdTLpe5+OKL6evr44Mf/CDz5s2r93AE+NGPfsSll17KkUceia7rnHrqqXzzm9+s3x8EAUuWLKFSqdS33XbbbfV9Pc9j5syZ3H333Vt0DtKGQIjNt2LFCqZO3YtqtTJg+6O3vPHvdDrDq6++MiTB43nnnUdfXx+PPvrogO1PPvkkRxxxBL29vZxxxhkce+yxW/25xY6nP0M4FOMCrFy5csC6esdxtvpzDWbRokV84xvf4MUXX9zu1jJL4CiEEDuYj3zkI5v8Y6tpGjfccAM33HDDRvdpbW3lwQcf3Oj9EydO3OA5UqkUd911F3fdddfmH7QQ4l3r6uqiWq3wkXPvpbl9jw3u7+v4O08+8G90dXVts6xjOp0mnU5vk+cW25eh7uP4TgqytbW1YRjGJusCvFV7e/sm93/qqafo7Owc8DMYRRGf//znuf3221m2bNkW1Rp4L8gaRyGEEEKIHUhz+x60TThgg9tgweR77Xvf+96Atj3XX389Bx54IN/61rcYP348mUyG008/nXw+v+0OUoga27Y5+OCDB6zfj+OYBQsWbHT9/vTp0wfsDzB//vz6/p/61Kd46aWXWLx4cf02duxYrrrqKp544on6GJtba+C9IBlHIYQQQgixzbz22mv85Cc/4ec//zmFQoELL7yQz3zmM/zoRz/a1ocmtqmhmaoKmzfmrFmzOPfccznkkEM49NBDuf322ymXy5x//vkAnHPOOYwbN465c+cCcPnll/PhD3+Yr3/96xx33HE89NBDvPDCC9x3330AjBgxghEjRgx4DsuyaG9vZ8899wTeWa2BbUECRyGEEEIIsVU8/vjj5HK5AduiKNrkY1zX5fvf/z7jxo0D4I477uC4447j61//urTzEdvcGWecwfr165kzZw4dHR0ceOCBzJs3r14AZ8WKFQP6hh5++OE8+OCDXHvttVxzzTVMmTKFRx99lH333Xeznvftag1sCxI4CiGEEEKIreKII47gnnvuGbDt2Wef5eyzz97oYyZMmFAPGiGZphfHMUuWLJHAcSc21MVxNsell17KpZdeOuh9Tz755AbbTjvtNE477bR3PP6yZcs22PZ2tQa2BQkchRBCCCHEVpHNZpk8efKAbatWrdpGRyO2Z0NdHEdsPgkchRBCCCF2IH0df9+s7dvaihUrWLNmTX3t1jPPPIOu6/X1XkKI4UECRyGEEEKIHUBbWxvpdIYnH/i3je6TTmdoa2t7D4/q7aVSKc4991z+8z//k0KhwGWXXcbpp59en6b6yCOPMHv27I02XBdCvDckcBRCCCGE2AFMmDCBV199ha6uLvxQccvPCgB84ZRGbDNpNN7W1rbNejhuzOTJkznllFM49thj6enp4fjjj+fuu++u35/P51myZMk2PEKxLQynNY4iIYGjEEIIIcQOYsKECUyYMAEvULS90AvA+97XgmNpQ/7c3/ve9wbd/pGPfKT+Yf28887jvPPO22CfSy65hEsuuWTQx2/sMWLHJoHj8COBoxBCCCHEDsaxNL79mdZtfRhCbDnF5rZcfOfjii2iv/0uQgghhBBCCCF2ZhI4CiGEEEKIbeL6669n8eLF2/owxDCkhvAmtoxMVRVCCCGEEEIMKypWqHgI1jgOwZg7CwkchRBCCCGEEMOKUsltKMYVW0amqgohhBBCCCGE2CQJHIUQQgghhBBCbJJMVRVCCCGE2MF4fsxn5i4B4O7Ze+LYkisQ2xfp4zj8SOAohBBCCCGEGF5kkeOwI5efhBBCCCGEEEJskgSOQgghhBA7oNWvP82j95zMggW/fk+e77zzzuOkk07aYPuTTz6Jpmn09fVt8vF/+tOfOOOMMxgzZgyO47Drrrty/PHH8/Of/1ymF+6E+hOOQ3ETW0YCRyGEEEKIHYxSihd/ewc965Zw3XVzhn3g9b//+78cdthhlEolHnjgAV555RXmzZvHySefzLXXXks+n9/WhyjETm+7DxxvvvlmNE3jiiuuqG9zXZfPfvazjBgxglwux6mnnsq6des2OY5Sijlz5jBmzBjS6TQzZszgH//4xxAfvRBCCCHEu+P58Qa3X/7fE6xf9Wf2OexcXnj+WX7xi3kb7DNclMtlLrzwQo477jh+8YtfcNRRR7Hbbrux1157ceGFF/LnP/+ZpqambX2Y4j3WXxxnKG5iy2zXxXGef/55vvWtb7H//vsP2H7llVfyi1/8gocffpimpiYuvfRSTjnlFP74xz9udKyvfvWrfPOb3+SBBx5g0qRJfOlLX2LmzJn87W9/I5VKDfWpCCGEEEJskf7qqf2UUjx+/zWM3OUADp15NZ2rFnPBJbM5/sJd0TStvt/91+31Xh/qoH71q1/R3d3N1VdfvdF93nzcYucgtXGGn+0241gqlfjkJz/Jt7/9bVpaWurb8/k8999/P7feeisf/ehHOfjgg/nud7/L008/zTPPPDPoWEopbr/9dq699lpOPPFE9t9/f77//e+zZs0aHn300Y0eg+d5FAqFATchhBBCiG1p9et/ZP2qP/O+j3wWTdN430c+y/pVf2b16xu/gL61PP744+RyuQG3Y445ZpOP+fvf/w7AnnvuWd/2/PPPDxjj8ccfH9LjFsORGsKb2BLbbcbxs5/9LMcddxwzZszgpptuqm9ftGgRQRAwY8aM+rapU6cyYcIEFi5cyGGHHbbBWEuXLqWjo2PAY5qampg2bRoLFy7kzDPPHPQY5s6dy5e//OWteFZCCCGEEJvn7tlvBFxKKT70ofMZPf4Axu3+AQDG7f4BRo8/gJ4l9/O//3XBkGbvjjjiCO65554B25599lnOPvtsAHK5XH372Wefzb333jvoOPvvvz+LFy8GYMqUKYRhODQHLIYtyTgOP9tl4PjQQw/x4osv8vzzz29wX0dHB7Zt09zcPGD76NGj6ejoGHS8/u2jR49+x48BmD17NrNmzap/XSgUGD9+/Ds9DSGEEEKId82x35hA9sQTT/DC889y1Nn31QNETdM44MOf5Vc/vJjfPflrZs6cOWTHks1mmTx58oBtq1atqv+7PxgEaGxsBJLAEGDJkiX1C/yO42wwjhBi29ruAseVK1dy+eWXM3/+/G2+9tBxHBzH2abHIIQQQggB/YX+rqOxdQKpTAtda/5Wvy+VaaGxdQJz5lzHUUcdtc3WDA4WDB511FG0trZyyy238Mgjj2yDoxLD0VAVspHiOFtuuwscFy1aRGdnJwcddFB9WxRF/P73v+fOO+/kiSeewPd9+vr6BmQd161bR3t7+6Bj9m9ft24dY8aMGfCYAw88cEjOQwghhBBia/J9n5WrVlHoWc1j95026D6rVsX4vj+sLnzncjn+67/+izPOOIPjjjuOyy67jClTplAqlZg3bx4AhmFs46MUQmx3geORRx7Jyy+/PGDb+eefz9SpU/nCF77A+PHjsSyLBQsWcOqppwLJ1IcVK1Ywffr0QcecNGkS7e3tLFiwoB4oFgoFnn32WS655JIhPR8hhBBCiK3BcRyefWYh69evxw9i5n5nOQCzL9gV20qms44aNWpYBY39Tj75ZJ5++mluueUWzjnnHHp6emhqauKQQw7hoYce4vjjj9/WhyjeY7LGcfjZ7qqqNjQ0sO+++w64ZbNZRowYwb777ktTUxMXXnghs2bN4re//S2LFi3i/PPPZ/r06QMK40ydOrU+HaK/D+RNN93EY489xssvv8w555zD2LFjOemkk7bRmQohxNApFotcccUV7LrrrqTTaQ4//PAB68bPO+88NE0bcDv66KPfdty77rqLiRMnkkqlmDZtGs8999xQnoYQ4i3Gjx/PQQcdxPvedxBtY/embezevO99B3HQQcltl112GbLn/t73vjdoNfqPfOQjKKU2qD/xVocccggPP/ww69atIwgCurq6mDdvHmeccYa049gJqVgN2W1zbe7ftocffpipU6eSSqXYb7/9+OUvfzng/uuvv56pU6eSzWZpaWlhxowZPPvsswP2mThx4gZ/h2+++ebNPvatabvLOL4Tt912G7quc+qpp+J5HjNnzuTuu+8esM+SJUvI5/P1r6+++mrK5TIXX3wxfX19fPCDH2TevHnbfB2lEEIMhf/3//4ff/nLX/jBD37A2LFj+eEPf8iMGTP429/+xrhx4wA4+uij+e53v1t/zNtlKX784x8za9Ys7r33XqZNm8btt9/OzJkzWbJkCaNGjRrS8xFCDOTY+rDp0yjE9mxz/7Y9/fTTnHXWWcydO5fjjz+eBx98kJNOOokXX3yRfffdF4A99tiDO++8k912241qtcptt93GUUcdxWuvvcbIkSPrY91www1cdNFF9a8bGhqG/oQ3QVOyQnSrKRQKNDU10dnZzciRrdv6cIQQO5j+3zH5fL5ejXBLVKtVGhoa+N///V+OO+64+vaDDz6YY445hptuuonzzjuPvr6+Tfayfatp06bx/ve/nzvvvBOAOI4ZP348n/vc5/jiF7846GM8z8PzvPrX/dWp163rYtSoEVt2gkJsx1zXZenSpUyaNEkuXr8Lm3odt9bvUjE0+t+f++/7PZlM7u0fsJkqlRIXXvwvrFy5csD7v7Gil5v7t+2MM86gXC4P6D162GGHceCBB260/Uz/Of/617/myCOPBJKM4xVXXMEVV1zxbk53q9rupqoKIYR4d8IwJIqiDT5MpdNp/vCHP9S/fvLJJxk1ahR77rknl1xyCd3d3Rsd0/d9Fi1aNKAfrq7rzJgxg4ULF270cXPnzqWpqal+k5ZGQgghABSqXll1q95Icmbjx48f8Pdn7ty5GxzDlvxtW7hw4YD9AWbOnLnR/X3f57777qOpqYkDDjhgwH0333wzI0aM4H3vex9f+9rXtnk/0x1yqqoQQoiNa2hoYPr06dx4443stddejB49mv/+7/9m4cKF9VL5Rx99NKeccgqTJk3i9ddf55prruGYY45h4cKFg1Y37OrqIoqiQfvhvvrqqxs9FumHK4QQYjBDXRxnsIzjW23J37aOjo531Bv+8ccf58wzz6RSqTBmzBjmz59PW1tb/f7LLruMgw46iNbWVp5++mlmz57N2rVrufXWWzfrfLcmCRyFEGIn9IMf/IALLriAcePGYRgGBx10EGeddRaLFi0C4Mwzz6zvu99++7H//vuz++678+STT9an0WwN0g9XCCHEYIa6j2NjY+M2nap8xBFHsHjxYrq6uvj2t7/N6aefzrPPPltfN/nmi6r7778/tm3z6U9/mrlz526zv5syVVUIIXZCu+++O7/73e8olUqsXLmS5557jiAI2G233Qbdf7fddqOtrY3XXntt0Pvb2towDIN169YN2L6pHrpCCCHEcLYlf9va29vf0f7ZbJbJkydz2GGHcf/992OaJvfff/9Gj2XatGmEYciyZcu27GS2AgkchRBiJ5bNZhkzZgy9vb088cQTnHjiiYPut2rVKrq7uxkzZsyg99u2zcEHH8yCBQvq2+I4ZsGCBRvtoSuEEEIMZ1vyt2369OkD9geYP3/+2/4tjON4QLG4t1q8eDG6rm/TKuUyVVUIIXZCTzzxBEop9txzT1577TWuuuoqpk6dyvnnn0+pVOLLX/4yp556Ku3t7bz++utcffXVTJ48mZkzZ9bHOPLIIzn55JO59NJLgWRazbnnnsshhxzCoYceyu233065XOb888/fVqcpxE7L80I+fckTAHzrnpk4jnzkE9uZIZqqurkLJ9/ub9s555zDuHHj6sV1Lr/8cj784Q/z9a9/neOOO46HHnqIF154gfvuuw+AcrnMV77yFT7+8Y8zZswYurq6uOuuu1i9ejWnnXYakBTYefbZZzniiCNoaGhg4cKFXHnllZx99tm0tLRsxRdj88hvESGE2Anl83lmz57NqlWraG1t5dRTT+UrX/kKlmURhiEvvfQSDzzwAH19fYwdO5ajjjqKG2+8ccC6itdff52urq7612eccQbr169nzpw5dHR0cOCBBzJv3rwNigQIIYQQb2eoi+O8U2/3t23FihXo+huTOA8//HAefPBBrr32Wq655hqmTJnCo48+Wu/haBgGr776Kg888ABdXV2MGDGC97///Tz11FPss88+QLL+/6GHHuL666/H8zwmTZrElVdeOWDd47YgfRy3IunjKIQYSjtD77H+c5Q+jmJntbX6OFYqHqeceiuVai/XXjOTI4/8yKAVkYej733ve1xxxRX09fVt8RjSx3H71f/+fOvu35JOb/0+jtVqiU9/5gh5/7eArHEUQgghhNiB/OxnP2Pq1D14Yt41PPW7rzFz5gx22213fvaznw35c3d0dHD55ZczefJkUqkUo0eP5gMf+AD33HMPlUrlHY1xxhln8Pe//32Ij1QIsblkqqoQQgghxA7iZz/7GZ/4xCfYZfz7Ofa4y2hu2ZW+3uW8/PJP+MQnPsFPf/pTTjnllCF57n/+85984AMfoLm5mf/4j/9gv/32w3EcXn75Ze677z7GjRvHxz/+8bcdJ51Ok06nh+QYxfZjqNtxiM0nGUchhBBCiB1AFEVceeUsdhn/fo746LWMHDUVy0ozctRUjvjotewy/v1ceeXniaJoSJ7/M5/5DKZp8sILL3D66aez1157sdtuu3HiiSfyi1/8ghNOOAGAW2+9lf32249sNsv48eP5zGc+Q6lUqo/zve99j+bm5vrXr7/+OieeeCKjR48ml8vx/ve/n1//+tdDcg5i+Ohf4zgUN7FlJHAUQgghhNgBPPXUU6xYsZz99jsdTRv4EU/TdPbd7zRWrFjGU089tdWfu7u7m1/96ld89rOfJZvNDrqPpmkA6LrON7/5Tf7617/ywAMP8Jvf/Iarr756o2OXSiWOPfZYFixYwJ/+9CeOPvpoTjjhBFasWLHVz0MMIxI5DjsSOAohhBBC7ADWrl0LQHPLroPe31Lb3r/f1vTaa6/VW/y8WVtbG7lcjlwuxxe+8AUArrjiCo444ggmTpzIRz/6UW666SZ+8pOfbHTsAw44gE9/+tPsu+++TJkyhRtvvJHdd9+dxx57bKufhxg+1BDexJaRwFEIIYQQYgcwZswYAPp6lw96f29te/9+74XnnnuOxYsXs88++9Sbm//617/myCOPZNy4cTQ0NPCpT32K7u7ujRbPKZVK/Pu//zt77bUXzc3N5HI5XnnlFck4CvEek8BRCCGEEGIH8KEPfYgJE3bl5Zd/glLxgPuUivnLyw8zYcJEPvShD2315548eTKaprFkyZIB23fbbTcmT55cL3azbNkyjj/+ePbff3/+53/+h0WLFnHXXXcB4Pv+oGP/+7//O4888gj/8R//wVNPPcXixYvZb7/9Nrq/2DEks0rVENy29ZltvyRwFEIIIYTYARiGwW233cqqlc/z29/cRGfnKwRBhc7OV/jtb25i1crnue22rw9JP8cRI0bwsY99jDvvvJNyubzR/RYtWkQcx3z961/nsMMOY4899mDNmjWbHPuPf/wj5513HieffDL77bcf7e3tLFu2bCufgRh2ZK7qsCOBoxBCCCHEDuKUU07hpz/9KagO/u8XV/HgD0/n/35xFRqdQ9qKA+Duu+8mDEMOOeQQfvzjH/PKK6+wZMkSfvjDH/Lqq69iGAaTJ08mCALuuOMO/vnPf/KDH/yAe++9d5PjTpkyhZ/97GcsXryYP//5z/zrv/4rcRxv8jFCiK1P+jgKIYQQQuxATjnlFI4++jhOOfVWKtVerr1mJkce+ZEhyTS+2e67786f/vQn/uM//oPZs2ezatUqHMdh77335t///d/5zGc+QyaT4dZbb+WWW25h9uzZ/Mu//Atz587lnHPO2ei4t956KxdccAGHH344bW1tfOELX6BQKAzpuYhtT/o4Dj+akldvqykUCjQ1NdHZ2c3Ika3b+nCEEDuY/t8x+XyexsbGbX04Q6L/HNet62LUqBHb+nCEeM+5rsvSpUuZNGkSqVRqi8fxvJBPX/IEAN+6ZyaOs3PlCjb1Ou4Mv0u3Z/3vz523/5p0evDWLu9GtVrm0itmyPu/BXau3yJCCCGEEDsBxzH53neO29aHIcQWG6qWi5Iy23KyxlEIIYQQQgghxCZJxlEIIYQQQggxrMgax+FHAkchhBBCCCHEsCJTVYcfmaoqhBBCCDHMSFbk3ZHXT4itTwJHIYQQQohhor9lhu/72/hItm+VSgUAy7K28ZGILaWG8Ca2jExVFUIIIYQYJkzTJJPJsH79eizLQtflGv/mUEpRqVTo7Oykubl5yHtXiiE0VFGeRI5bTAJHIYQQQohhQtM0xowZw9KlS1m+fPm2PpztVnNzM+3t7dv6MMS7IMVxhh8JHIUQQgghhhHbtpkyZYpMV91ClmVJpnEHIMVxhh8JHIUQQgghhhld10mlUtv6MIQQok4mzgshhBBCCCGE2CTJOAohhBBCCCGGFVnjOPxI4CiEEEIIIYQYVhRDtMZx6w+505CpqkIIIYQQQgixEXfddRcTJ04klUoxbdo0nnvuuU3u//DDDzN16lRSqRT77bcfv/zlLwfcf/311zN16lSy2SwtLS3MmDGDZ599dsA+PT09fPKTn6SxsZHm5mYuvPBCSqXSVj+3zSGBoxBCCCGEEEIM4sc//jGzZs3iuuuu48UXX+SAAw5g5syZdHZ2Drr/008/zVlnncWFF17In/70J0466SROOukk/vKXv9T32WOPPbjzzjt5+eWX+cMf/sDEiRM56qijWL9+fX2fT37yk/z1r39l/vz5PP744/z+97/n4osvHvLz3RRNyUTfraZQKNDU1ERnZzcjR7Zu68MRQuxg+n/H5PN5Ghsbt/XhDIn+c1y3rotRo0Zs68MRQuyAdobfpduz/vfnazfPI53KbvXxq26Zq7549Dt+/6dNm8b73/9+7rzzTgDiOGb8+PF87nOf44tf/OIG+59xxhmUy2Uef/zx+rbDDjuMAw88kHvvvXfQ5+g/51//+tcceeSRvPLKK+y99948//zzHHLIIQDMmzePY489llWrVjF27NgtOfV3TTKOQgghhBBCiGFFDeENkmDtzTfP8zY4Bt/3WbRoETNmzKhv03WdGTNmsHDhwkGPe+HChQP2B5g5c+ZG9/d9n/vuu4+mpiYOOOCA+hjNzc31oBFgxowZ6Lq+wZTW95IEjkIIIYQQQojhZYgjx/Hjx9PU1FS/zZ07d4ND6OrqIooiRo8ePWD76NGj6ejoGPSwOzo63tH+jz/+OLlcjlQqxW233cb8+fNpa2urjzFq1KgB+5umSWtr60af970ggaMQQuyEisUiV1xxBbvuuivpdJrDDz+c559/vn6/Uoo5c+YwZswY0uk0M2bM4B//+Mfbjru5BQSEEEKIwQx1xnHlypXk8/n6bfbs2e/RmSWOOOIIFi9ezNNPP83RRx/N6aefvtF1k8OFBI5CCLET+n//7/8xf/58fvCDH/Dyyy9z1FFHMWPGDFavXg3AV7/6Vb75zW9y77338uyzz5LNZpk5cyau6250zM0tICCEEEJsK42NjQNujuNssE9bWxuGYbBu3boB29etW0d7e/ug47a3t7+j/bPZLJMnT+awww7j/vvvxzRN7r///voYb/3bGYYhPT09G33e94IEjkIIsZOpVqv8z//8D1/96lf5l3/5FyZPnsz111/P5MmTueeee1BKcfvtt3Pttddy4oknsv/++/P973+fNWvW8Oijj2503FtvvZWLLrqI888/n7333pt7772XTCbDd77znY0+xvO8DdaZCCGEECiFGoLb5jSHtG2bgw8+mAULFtS3xXHMggULmD59+qCPmT59+oD9AebPn7/R/d88bv86y+nTp9PX18eiRYvq9//mN78hjmOmTZv2jo9/a5PAUQghdjJhGBJFEalUasD2dDrNH/7wB5YuXUpHR8eAxf1NTU1MmzZtk4v7N7eAAMDcuXMHrDEZP378uzw7IYQQO4L+GG8obptj1qxZfPvb3+aBBx7glVde4ZJLLqFcLnP++ecDcM455wyY5nr55Zczb948vv71r/Pqq69y/fXX88ILL3DppZcCUC6Xueaaa3jmmWdYvnw5ixYt4oILLmD16tWcdtppAOy1114cffTRXHTRRTz33HP88Y9/5NJLL+XMM8/cZhVVQQJHIYTY6TQ0NDB9+nRuvPFG1qxZQxRF/PCHP2ThwoWsXbu2vvB+c4oBbEkBAYDZs2cPWGOycuXKd3l2QgghxNZzxhln8J//+Z/MmTOHAw88kMWLFzNv3rz637sVK1awdu3a+v6HH344Dz74IPfddx8HHHAAP/3pT3n00UfZd999ATAMg1dffZVTTz2VPfbYgxNOOIHu7m6eeuop9tlnn/o4P/rRj5g6dSpHHnkkxx57LB/84Ae577773tuTfwtzmz77Frjnnnu45557WLZsGQD77LMPc+bM4ZhjjgHAdV0+//nP89BDD+F5HjNnzuTuu+/e4MPMmymluO666/j2t79NX18fH/jAB7jnnnuYMmXKe3FKQgjxnvvBD37ABRdcwLhx4zAMg4MOOoizzjprwLSY94LjOIOuKxFCCLFzq08tHYJxN9ell15azxi+1ZNPPrnBttNOO62ePXyrVCrFz372s7d9ztbWVh588MHNOs6htt1lHHfZZRduvvlmFi1axAsvvMBHP/pRTjzxRP76178CcOWVV/Lzn/+chx9+mN/97nesWbOGU045ZZNjbkkRCCGE2J7tvvvu/O53v6NUKrFy5Uqee+45giBgt912qy+835xiAFtSQEAIIYQQ24/tLnA84YQTOPbYY5kyZQp77LEHX/nKV8jlcjzzzDPk83nuv/9+br31Vj760Y9y8MEH893vfpenn36aZ555ZtDxtrQIBEhRByHE9i+bzTJmzBh6e3t54oknOPHEE5k0aRLt7e0DFvcXCgWeffbZjS7u35ICAkIIIcTGDJc1juIN213g+GZRFPHQQw9RLpeZPn06ixYtIgiCAcUZpk6dyoQJEzZanGFLikD0k6IOQojt1RNPPMG8efNYunQp8+fP54gjjmDq1Kmcf/75aJrGFVdcwU033cRjjz3Gyy+/zDnnnMPYsWM56aST6mMceeSR3HnnnfWv366AgBBCCCG2X9vdGkeAl19+menTp+O6LrlcjkceeYS9996bxYsXY9s2zc3NA/bfVHGGLSkC0W/27NnMmjWr/nWhUJDgUQixXehvdrxq1SpaW1s59dRT+cpXvoJlWQBcffXVlMtlLr74Yvr6+vjgBz/IvHnzBlRiff311+nq6qp/fcYZZ7B+/XrmzJlDR0cHBx544IACAkIIIcQ7NlTpQUk5brHtMnDcc889Wbx4Mfl8np/+9Kece+65/O53v3vPj0OKOgghtlenn346p59++kbv1zSNG264gRtuuGGj+/QXKXuzTRUQEEIIId4piRuHn+1yqqpt20yePJmDDz6YuXPncsABB/CNb3yD9vZ2fN+nr69vwP6bKs6wJUUghBBCCCGEEGJnsl0Gjm8VxzGe53HwwQdjWdaA4gxLlixhxYoVGy3OsCVFIIQQQgghhBBDR4rjDD/b3VTV2bNnc8wxxzBhwgSKxSIPPvggTz75JE888QRNTU1ceOGFzJo1i9bWVhobG/nc5z7H9OnTOeyww+pjTJ06lblz53LyyScPKAIxZcoUJk2axJe+9KUNikAIIYQQQggh3iuqdhuKccWW2O4Cx87OTs455xzWrl1LU1MT+++/P0888QQf+9jHALjtttvQdZ1TTz0Vz/OYOXMmd99994AxlixZQj6fr3/9TopACCGEEEIIId4bssZx+NGUkpdvaykUCjQ1NdHZ2c3Ika3b+nCEEDuY/t8x+XyexsbGbX04Q6L/HNet62LUqBHb+nCEEDugneF36fas//25Yc7jpFLZrT6+65aZc8Px8v5vgR1ijaMQQgghhBBCiKEjgaMQQgghhBBCiE3a7tY4CiGE2HlEUYymaei6ttmPVUoRRwrD3LJrpFEYoxsamrb5zy2EEOLdUbX/DcW4YstI4CiEEGJYCoOIcslHNzQaGjevWJlSilLRI44U2QYHczODR88NcasBtmOQztib9VghhBDvnhTHGX5kqqoQQohhyfciAOJIEceb95c+CmPiKHlMGESb/dxuNRhwDEIIId5jaghvYotI4CiEEGJYit90WTiK4s167JuvKPcHkFt8HJsZtAohhBA7IpmqKoQQYlh6c7eoOFJgvfPHvjnYi+LNCzrfGijGsdqiNZZCCCG2nKxxHH4kcBRCCDEsqTcFcJub9Xtz0Nk/jueFhEGEaRrohoZSoGlgWcZGH/vW4xDbP88NCaMYXQdN0zBNHdM03v6BQoj3lKxxHH4kcBRCCDHsKKUG/HHf3ODtzYGmUlCtBoRBRBBGVCoBUS2LaFk62ayDbb8ROLz1uaIoxkICix1BtRpQrfpEoSKqvc+6rtGQs6UI0k4ijhWFvLutD0O8A8nfgSHIOErkuMUkcBRCCDHsvPXv+uZON33zBwPfj4hRBH6MYeok3TUUGuB5EUp5NDY69azTW59bPmTsGKrVgELeRVHLMmsxKq4FEgUPTdNIpTdjPrTY7sRxTG9vlaobbOtDEWK7JIGjEEKIYefdThftf7jvRwRBRBQrbMfAcUwsy8AwNMIgplTyCPyIajUgk9EwDH2D55biONs/3wsplzwUkEqZZLM2pqmjFBQLLlU3pFBwsWwdw5Ds8o6qkHcJghhD1iwLsUUkcBRCCDHsVCs+ugpwUiaxUnjVENsxSaUtlFJUSj62Y2LZg3/IV0oRBhFRFBOFMbqhYVkGqZSJYSQFxQ1DJ4pjKmUfz43QtIBMxqqXTdB1jTje/FYgw1UcvzHtS9O0nabgj+9HlCs+YaRwHIPGRgddT74HNA0aGh2CMCIMFfm8R3Nz+h29Nv3TqftfU13X0LSd4zXdHpVKHq4XocFm94UV24ascRx+JHAUQggx7JQKPrFXRR+RIQhCKgUfTYP2cU141RDfS24tbdlBH69ihe9H6KaGQmEYOrZt1IPGfo5jEoUxrhsShjHVakh/zKAbSeC4IxTHCWpB9Jvpuo5l6Tt0sBPHCtcN8L0I2zLIZu160NhP13Wam9P0dFeTILPkvW1gEUUxwVv6g2pacnFiZwnItydhGFEu+0ByoSAMt/+f6Z2BrHEcfqSPoxBCiOGn9oddARrJB/GwNqU0DKNNPDDh+REKiCIwTQMNBhTA6WcYepLJTJn4fkgUxYRhEmD1BxhKbb/TVZVKAuj+oFHT3siKxXGM520YUO5IfD/Ed0M0XSOVNrHt5Hr5Wz84mqZBQ6ODBlSrIb4fbnTMgUFj/+up1V7rcIOAUmx7hYKLUuDYBhkpgiTEFtuijONjjz222Y/52Mc+Rjqd3pKnE0IIsZMxzSSgsy0DZehUKwGxSoKdIIhRSqFpWv2/bxZFcTJNNYxJZUwCX2FZxkYza7ZtEIYxpmkkmblQYRoaoOrTVVWsYDvMJAVBTFwrLJSs7UyC4ThWBEGEUqoe6Lw1G7u9SzLIAWGsSKeTta3VakAUxcSxwjR10mmr/n2RSpm4VQPPT7JTg33PhGFcv3BhGHq9lYtSijCMk6nR0Ruvt3jvKJVMK3/rNOxy2cf3Y5miuh1SDNFU1a0/5E5jiwLHk046abP21zSNf/zjH+y2225b8nRCCLFD2dKLbzsV7U3/0TRMQ0+KmcSAUgR+hGUbqFihGQM/3LtukKxNVDGWbYJKMk4bfSpNw3FMwjCip6dKFCbFM+I4CQY0IJ2xMMy3D6x8L6RU9GrTFpNjtlPmBtMj3wtJgLRh0AjJejzbNuqBZRBEO9y6R88NqFQCQBGFBn19VYJgYHa1c22RdMYilbawTINU2iIIkkys54WkUm9UWe0PDmFg0AgDp6m+eVqwBI9vL661xnk7/T9baGBbxoCfrSCIcKvBgCBD0zQMEyrlpIJqNpcURBJiS9x111187Wtfo6OjgwMOOIA77riDQw89dKP7P/zww3zpS19i2bJlTJkyhVtuuYVjjz0WgCAIuPbaa/nlL3/JP//5T5qampgxYwY333wzY8eOrY8xceJEli9fPmDcuXPn8sUvfnFoTvId2OI1jh0dHYwaNeod7dvQ0LClTyOEEDucLb341tbWNjQHNBzVPgAqxYBFFZap43vJ9sBPqqXqb/ps3p8x7M8gGoYOVhJgDpad7BeGEYW8i+9GVCs+pmlgO2aS1YwUpZJHk6ljGHq9yEySfYqJoxhN04jjiELe2+BytmXqjBjd8J5+aA3DiHIpQNPBcTZc2wnUg9sgSDK5vh9i2+YOETz6fsT6rjKVsodtmwRelFyAMHVSKRNNg86OMq7r09sLmbSFricVdS3HwElZlEoeTu17AKgFjaoeJPbrz96iwDCTyrz9mUet9pxicP3TfnVdH3QqeT+34rN+fRm3GuC5ARoahm2QTps0t2TRtKRlTxwpwiAmipPX3nMD0CCXdchmZYrq9ma4FMf58Y9/zKxZs7j33nuZNm0at99+OzNnzmTJkiWDxkJPP/00Z511FnPnzuX444/nwQcf5KSTTuLFF19k3333pVKp8OKLL/KlL32JAw44gN7eXi6//HI+/vGP88ILLwwY64YbbuCiiy6qf72tY6otChzPPffczZp2evbZZ9PY2LglTyWEEDskufi2aar+X4VB/1pDBVqSxfHdsB48vnm6oO+FxKp/jaJGOm1SqU01jaIYwxhYDCaOYwoFj0rZB03DdkyiOEbXNFSsMCwN14uwAoO+3iq6nvR+jGOFoeu1aq06biUk31dF05LsZK7RIQoVvhcRhDHd64q0jsphmkNbjCaKYnwvwnWDenZMxQrfi7Hs5MP5m7Of/cGj76v6ekjH2fi03uHA95NCRqEfJ4FgZmBGNwgiOtYUKBU9dENL0tZardBGrPDcELfqE4ch1WpIOm3geSGmZaBrGm4lpFz2yeUcbMugsSlNXPv+AeqBYBwnmW/fD9/4IOollVo1TUM3NMIwQtN2vGnAW0M94H4bxWKVzrUlCvkqYRTXZx6EZZ+uzojOjiItI7M0Nabr7XaUUlQqAaWij0Lh2CZhGEsGeLujGJqJpZs35q233spFF13E+eefD8C9997LL37xC77zne8Mmv37xje+wdFHH81VV10FwI033sj8+fO58847uffee2lqamL+/PkDHnPnnXdy6KGHsmLFCiZMmFDf3tDQQHt7++ae4JDZosDxu9/97mbtf88992zJ0wghxA5JLr69PfVG5FifZpqsNQRd0zAtnTBI1psl2UUd34uSYNJLvjZMHdM00LUQP4zo667Wp7dlsjZxrOjpruLVCqHkMjZ2ymDd2phCLQi0LR3XjaiUPXINaXwvxLSN2gfQJPtULPpUK0nVVydl0dCUwrIM0s3Jc6xfV8T1Qpb+vYvWtiyZnE1mCLIfUZS0FumfYmvWMqSaVsuo9lZA08g1ODS1vPH9p2nJtFXPi4Dkw3x/EZmh0N8ipb/VSf8aUqVULdNrYL7lA36l4lPscwnDCPWWWj6GqZPNOeQabVw3pHNdkd6eCr4XkcvZYOp4XoCTsojCmJ58Bd9LMlKjR+eSYMTSSWccVBhRqQT45ZBC3kVBkoWtBX66rteDwGrFJ4qSb1TD0NB1vb5+sj8IN4xkHa7jWDtEJndriWtVjyG5wGNZgwfWPd0V1qzopVhwiRXJz07axnYMymWPaldAX4+P64Zo42DEqBzpdPJau25AKm0Rq5hYQXdXmRFt2fqFJzH8DXXGsVAoDNjuOA6O4wzY5vs+ixYtYvbs2fVtuq4zY8YMFi5cOOj4CxcuZNasWQO2zZw5k0cffXSjx5TP59E0jebm5gHbb775Zm688UYmTJjAv/7rv3LllVdimtuuKca7euYgCDj66KO59957mTJlytY6JiGE2KFt6cW3t/6R25EFfphkiEgyOJBkzmIV19trqNrF6DB8o8hLGEVoRvLBwDT0WmAZke+tYFoGtm0S1tb1haHC80N0TaOpOYVp6LjVgDhKprlaloHnh1RLyboqtxrR2JzCsnSaW1JYlkml5FEueQA0N6VpacvgeyFRpKiUfdIZi7ZROZb/owvPC1i7qo+m5gx2OmlCb9lJoLQ11kD2B85KQSZrYZoGjmMSxzE968tAMu222OcSBBEtIzL1IKg/ePT9sJ4J2trZmThWeLXWGBsTBFEyddFIsr+GoVHMe5SLXn0fTU+qoJpWcrEgCmPWryuwanlIFCtcLyLwAwzdwPeT7HBjYwrbNigWPNxKEjS0jcph2gZp06JU9OnrLpPJ2ji2iWlq9PR4lAou3XbyXqXSZj3bmFSqTS4cpGqFd/qFYYznhfXvySCMiCNFJmsP60zueynJNL4x7ffNr0sUxQS1tiiv/KUDt+JjWiYjxzQwojWD7ZhUKgGpKKaxMU3V8nCrET1dJeyUBQoUMWGkaGh0aGywKRQ9wlDRtb6EYxt0rS9uu5MX79wQJxzHjx8/YPN1113H9ddfP2BbV1cXURQxevToAdtHjx7Nq6++OujwHR0dg+7f0dEx6P6u6/KFL3yBs846a8BF4ssuu4yDDjqI1tZWnn76aWbPns3atWu59dZb38lZDol3FThalsVLL720tY5FCCF2KnLxbeMK+SpEFgpFU0uaWCl0TSOO3lhfGMcK10uChUolQNeTKqK6lnwwdT0NsxxQLiUZCQdIZ2zcSkBPTxnN1HAci+a2LIauUa0kRTT02ofZajUkjpJCJ24te1XMu2QbHOJIgaXwvBDb1tF1g0zOJvAjnJSJ5yYBWKXsYxgaqYxNueShGTq93SWaWzOYtfWSnhti2QaWrZPvdeldXyKOwUmb2I5JS0uGdG7wDGUUxrXMV0y5nKxpzGRtjFoxIYAwiDENA7PBwNA18n0ubjkgr7tkG2xs26gXxrEso17cRde1rTbF0vdCXDes98Q0LR1d19F1rX4DcN0Qt+ITVCJ6uksUe5Og3XFMGppTNDSmBgRgcRzT3VmmXPYolTzKZQ/fj2lpSZNOW6AllTRtJ8ky9fVWCfyIbM6iUvIxbZ3AjSiXXEoFHwzIpm3QNAxDwy0n7/GYsY1oegrHMetTogFsx9ggwDZNHdO06wVbXDekUgmSQKbB2Skzj74X4tWyfMmFiYBKOcQt1wpJOcn3ei7roGoZ8r//bR2lggtoNLfY6LqGUuAHIVGYFC9KZ0x0TWHoAZVKyKql3TSPyCTZyaxN06gc6ayDnbLoXFukZ32ZwI/QzR23BY1451auXDkgUHtrtvG9EAQBp59+OkqpDWZovjlruf/++2PbNp/+9KeZO3fuNjlWeJeBIyRTqO6//35uvvnmrXE8Qgix05CLbxvXP92yVPDx3Aiv6gM6saawTbO+3sxOmbX1fBFhpEg5BkqBZRoYuo5hJNVR4yjGMHVSWYtqNaBY9NB1jewYB7fi16fNGbqGZSdr3qIgyW7GCnINdtIHzjFZ31Gk2pjC0AEFKcdmxKgsQS275bkhqbSF5yZr8dZ3FNE1jbbROTwvxnMDojBCIwkywjCma12B3vUVDFOrT4uslCIqJY/erhINTWlGjMpi2xamlayT9KoBbjUJdqvVABXFqEjDrQSoWNWnm/YHOamURSpjoRk6xT4Xt+JjmBphEJFKWxjGG8V/+ouWvNtKq1GtJYbvhoRBhGkZ5BodbGfgx48wSAKBKEymjVYqPr1dSYVb3dDJ5HRs2yQKY4oFD9NMqpr6frI2MZNzCKOYSiVI1py6Ac1tGZqa0vhuiO8l6x59PyKbtTH05PXrWF1FEWPqBoaZrF/tLJawbR2NJOvp+SGGnkx77n/e/rYPmyroYtUqfxqmQbnk4XshhViRzdk7zVq7OFa4lYDgTX0xOzuK9HSVURGksyaWnVys6OuqsCaIyDQ6FApVCnkXQzcYM6GJdCr5vvf9pLJqtRpiGKDipPhTyrbwg5hYKXq6y6TSyfpUatnzKEp+B0RhnPRllYYM2welNui5urXGBWhsbHzbZSBtbW0YhsG6desGbF+3bt1G1x62t7e/o/37g8bly5fzm9/85m2PZdq0aYRhyLJly9hzzz03ue9QedeBYxiGfOc73+HXv/41Bx98MNlsdsD92zKdKoQQw51cfNsITSOMFH4QgwalSoBt6lTKBqEZoxkajmMkAZimo+uQtkyclEnaNglqAYhp6gRuSBQq3IrPqqXd9KyvUip56DpEboBpm2imTiZtkck5lIpJsRzT0okihWXrGIZBW1uGfJ9LEEV0rysQBjHZhhQTd2utF+XoX/fmVgM0DQo9ZUp5F93Q2WViMy0jTLq6ypQLHtVqQK4pRaXisX5tERWDUgbt7VlyTSn8akix4LF+bYFy0aens0T7uGYyOXvAhynd0HFSZn3dThBERIGiVHDJ5GxCv9YawkmClUzWJvRDSgWfSsEjnbOJIoXtGAOKi/RPWe3PSG4utxbYetUQ3wvQDSMJ7koelm1ip0w0koxeGCaFi3RDx3cDPDekpTWDbmikUkmvxTCIk+O0dQI/ptDnEkVR0mJDTz5gNjWn6e2pEEewZlkfvaky2UaHjlV9VIoBVlrHcQxs26JrXZEwSCqlOq0muYYMrhviVUKCIMC0koBPKQ1Nh8ZGO2nVYhmk0xZO6u1fF03TyGQsDEOjVPIIgqherTWV2rHXPcaxolz0iGsXeZy0RW9XmUJfNbmw4Rg0t2Rw0gZuJSQKPUpFl/VdJQp5j+aWNOMmNjJydEOyzlnB2lV5uteXqJQ9TNMgmzEp9imiIMZzQ6quj26ZZLIBpq4olVNUyj7FgottmbS2ZUDTKJVL2/jVEdsL27Y5+OCDWbBgQb0iehzHLFiwgEsvvXTQx0yfPp0FCxZwxRVX1LfNnz+f6dOn17/uDxr/8Y9/8Nvf/pYRI0a87bEsXrwYXdffcWG9ofCuA8e//OUvHHTQQQD8/e9/H3CfzOUXQohNk4tvg9N0HdsxSacMNF0nm0mCG8fSSWcsKmWf1csKydRQx0TTNRqaU8SNKXQ0NC3JRnSuzdOxso9SsUoq4xDGMUEYoymwDKM2PTHAtg3iMMbzw2QKo2Um/eUMDd+LyeaSFg3tY226u0qsXtGH60aYlkm5EpBKWdgpk3TGploJcN2ANUt7CYOQYt4lnbMp9FXJ5lLYhk5V18n3VfEqPtWqV6v6Ck4KujtLqBhyjSkam3XKfVVWLOvFNjUKvVXGTmzCsW1My6ClLZmWp0iKjCRFg5J1fypWFHqqxLHCcsz6tFNdT9Y6Rn5EpNWmjtYyZ2EYk0olwaPn9a8djTeZWXsr3w3oWJ2nkHep1or1pFM2pq2jAbqpk0rFuBUfTaeeYdVMqJb9pMItSYDbNiqHUSt81J85jVVM6CuiOMZ1Q4Iwolz2Cf0Yt+rh2Mn+KKiUPfJ9FeIoxkol00rzPRVWvuaSyto4jk62MY1b8SkGLpat46QMdENLghNdw/NjSn0uHWaR0WNBRRpaKzQ0vvOpYklbjyQz7Psxupa8R07KHNJCRNuSW0nWC2u6RiZrUyq6dK0rYZgGo8Y0MmJUlkLepdTn4nsRMQona7B2ddJvJ45jnJSF74ZYjknXuiJ9PRUqpYAoiImDCBXFOI6VrIU2dLxAEZSquNWAtK3jVSPCIKahyUFv0mhvb8aydZYtrWzrl0e8A6pWOGsoxt0cs2bN4txzz+WQQw7h0EMP5fbbb6dcLterrJ5zzjmMGzeOuXPnAnD55Zfz4Q9/mK9//escd9xxPPTQQ7zwwgvcd999QBI0fuITn+DFF1/k8ccfJ4qi+vrH1tZWbNtm4cKFPPvssxxxxBE0NDSwcOFCrrzySs4++2xaWlq24quxed71b6vf/va3W+M4hBBip7QtLr5FUcT111/PD3/4Qzo6Ohg7diznnXce1157bf05zzvvPB544IEBj5s5cybz5s3b5Nib2yR5YxobbHJZBydlYtkmPV0REKOhUeytsm5Nkd6uEp4fYVkmmq7INaQZt2tTEhTFMT0dRXrXV6iWPSKg6kVkczatI7KMHJ2hu7NCuejipC3SGZvIjygXXCrlAMMMGDO+CV03cJ2AONbw/IimpmSdWzbngOYTBiF9PRVQSQGWTC4JRpYu6aNa8YhiRa4hhW7oVEsBKoaGxjSm5dPTWWFpVxnL0Mk1Opi6Tm9XlSCIWLakC8s0qJR9Aj8inbYolAOqVZ/X/tZJ+7hGJkwZiWaAZVtYtoFpJu+dk0qmnVZKXq3vXUjjm6ZG+m6YTNQzNHQtqbqqlMKrBDgZi0o5wLIMLNsgDKNaIaGkIM3b6euusPL1bvJ9LpViFd9T+F5A6Ec0taZpbs2SbbRxKyGWo5NrSGFaBk0tKUoFj8AL0XWdbM4mm7Wp1oJI00rWgAZ+hFdJMpSxSjLPfjmk2FelWgqSICVtMnpMAyoGv+rT1VXFcnRGjMzh2Ab/eKUTzUjarWSyNipSlAsuhT4Xx0kuAMQKgjDGdX3iMEJpGtVqMiV69NgmKmWfctEjnbVrbT5q7T5qmdr+7K+mgUaSsVQKDF3HthWBF2HaOm41JAiSYH1Hatnhe2F9emoma+N6AWtX5olihWmA70Ys/3sXlZJPoc+tVzYul5LqxHEQYlsGa1YUaGlNoaHR1+sSeBGOrZNpzdUq3CYzD7INKcoVj3SDTW9XhdiPWLOyF8M0aWxJo4hJpS3KJY9M1iaV2jGD9R2NGqKpqps75hlnnMH69euZM2cOHR0dHHjggcybN69eAGfFihUDCpwdfvjhPPjgg1x77bVcc801TJkyhUcffZR9990XgNWrV/PYY48BcOCBBw54rt/+9rd85CMfwXEcHnroIa6//no8z2PSpElceeWVG1Rr3VJhGLJmzZoBrT/eiXf9k3Puuedy4YUX8i//8i/vdighhNjpbIuLb7fccgv33HMPDzzwAPvssw8vvPAC559/Pk1NTVx22WX1/Y4++ugBFWDfbjH+5jZJ3pTAj7Bbk+mjaMl0xb7eKp3rSoTVgMAPsG0LJ2WDCYEXke9OpsE1NNhEKKolnziIsGwL2zHQzOQDZkPOJnCTKa62Y5FKmTQ1OZTLPpVuHz8IyVoOQSUk12Jh2Tb5vIfrBkmVz6KHaehkMxaVosfyv3fS1ZBm7IQmGhodqhUfFcdYtsWIZocwUBimhqZptTYEAbmszWt9VbrXF3EMi8Zmh0yDTVzy6Fqdr0+XrZQ8MhmHtjE50lmLatnFr0as/Gcv1bJP77hGTNtkxKgcYyc2YZpJAKLrGtkGh0rJB5V8kPe9MFmz54VJ2ffWNHGkKBWqGKGBpoOmKex0UtglDONkzaWu4fshgZ9si6KkH146Yw+YalksuKxZ0Udfd4VCvkoqbWJEEWExIp02iENFpeLjBbXANDAIw5gxYxvpWlcmjhWpjE1Ta5psrdCQ70fEYUxY6/dXqa2X8/2IMFYEbkhfdyXJsGox6XSKttFZ0hmHVMpk/dqYVNoELanCWyp6tI7M4lVC7JSOW47o6ymilCKXc5KiOKaOpWsEYQgxxJGG5/oUe6t0rSnQ01liwsQW+rrK5Boc0plk/atp6ZiWUV+DOpj+Nh2gqJYitFphILfqk05bpNL2dj99NYrieqGppB1GxOqlfRRLPtW8S2NLBqJk2m61EiaZ+qxNT1cJv/a9OXJ0I8W+Kr5XZn2HQUOjg2boOLbBiJGNNDQ65JrSyevphpRKLpqu09SUpqkpxZplvZRKHoQ+RhGyWYtSwaNS8UmlTKqeZBzF5rn00ks3OjX1ySef3GDbaaedxmmnnTbo/hMnTnzb4PWggw7imWee2ezjfKf++te/ctBBBxFFb99L9c3edeCYz+eZMWMGu+66K+effz7nnnsu48aNe7fDCiHETmFbXHx7+umnOfHEEznuuOOA5I/Yf//3f/Pcc88N2M9xnM1qPLy5TZIBPM/D895otdDfcqS7s0zsGTQ0Z8jmLOIgwvd81q0oojSFoes0NqfINtg0tWYo9FTp7iyRX1/mtRV9uJWAbINFW3sTzSMzxDHYloHjGFRKHsWCj2lppNMWgR/RsTZZ72ZaBrmcQxQovGqIYXjEelJsJ/AiVq3vRank+W3TpLtUolLxKfRWWb86abVhpgyyaYsxE5rp6axSdX0yaRt0iIKYStEjn68ShxFBNQbDp2tdmWopqbzZ2JZBt3X6eisEfkxJedi9Bk41oqkljZ1K1g3me6rke6qMHNdEtRTQ01Vm7Phmxk1sBpL1d7ZjJOvyTK0+bTSqBWGZrE3n6gJRCBgRkZeccyqIyTQ6KAVhoKi6yZTDYq9L6Ef1NWu6kRSHSWccdB1Wr+hj/doChT6X1hFpUo6JChRt7TmctE1Dc7LerFSo4rsRtpcUnHl9yXrslEU6azOqvYGGxlTy/ZfSk96LtfYM+d4KoR8SBjFuLTvc21UGXRFHGiPHNDFiRAZDSwoMdXeXCaKYbIODpkGp4FLsq5JtcBg3oYmudWUqVRfd0Mg2pEnnDIo9VXw/REUKTSl0PSYOk+IqpaJPHMbkuyt0LOtj1LhG2kblGD2uiVxjiihMXpvQ03EyyXpb3dBrWZNalkMlU4U9L8QwQdN04ijGrUS4lRDdcElnbDIZC9MytssgMrlYoZIg2tZZ8VoP69cWKearjGzPYRga2SaHfE+ZTKOJZRhUKj5+GJNusBg1qoGmlgxrV/TRVS3TtaZAx4qYEaNz7HfIOEa1N9Dcksa0kyJZvV1lyiWfaikg22CjAsg1pZKWOJUAwzAoFlyCQKGpmEopoFyRNY5CbIl3HTg++uijrF+/nh/84Ac88MADXHfddcyYMYMLL7yQE088EcuytsZxCiHEDmlbXHw7/PDDue+++/j73//OHnvswZ///Gf+8Ic/bLCe8sknn2TUqFG0tLTw0Y9+lJtuummjC/i3pEkywNy5c/nyl7+8wfZK0YXAxPciin0mhk7yod71CaMYFSl6e0qMHJXD0HXaRmdpG51j6T+6Wb++hOsGxLHCtMtYKYOW9gYgmfJYyLt4rg+aTaXo44ch1WqIaRhousJyLJqaHBqaU3iVANcPKOY9oigm1+Cg6zoNTQ7Vss+Y8c2Uyx4dawrke6vk+1zSGYuW1gzVaoSKY9C0pJ9iFFMueHSuKlKteIQxjGzL0NNbYeVrXSiSLOCoMY20jc6wflUe343IZJJsSbHPpVK0aR2dwzIMurtKaLpGx6o+mlpSZIopqgWP7s4SY3ZpZsToDCiNdMbGThkEfkQp7xIGEemMje+FmLXA0nJM0jmDUj5p1RGGMY0tKYJQEfoh3etL+OUQyzZJZy28akChxyX0Y3RTw3V9OlcVqJaTD+olw00qqpZDbMegZUQG3/Xxyl6yIDOOqZRC1neUsGutGMZObK4HeJmsXV/7aBg6fhzipCyCUKHFyTrKfG+5VrUzpLUtjV0rWpNO25QKLoGfrHPUNDBMnUrZR9OT9Y+rl/eRTls0tqTI5Ww8N2L9mkI9yxn4Yb31i2HqmIZOQ6NDMe8TuEmQHrgBlbyH5wZM2W8Mo9obCWrrQlHgVZOKr6lMMpX4zcIwxveTKrJJ4Z+ISiUgjpKCMpWyXyugk1Qd3V6CSM8Nk3WNmkY6Y7FmVR8r/tlNKe/R1p6loTlNY6PDmpV5giDCMHQC16djdRHPC2hqStPclCbwI5rasgRhhFvxiVRIteix4h89NDenickAyTrKfJ8LQLbRRsUQqoiGxhSNTSmCIKpP/873VQmqQTKdWG1elkVsI0Pcx3FH1r8EZmOq1eoWjbtVJnmPHDmSWbNmMWvWLF588UW++93v8qlPfYpcLsfZZ5/NZz7zGelRJoQQg9gWF9+++MUvUigUmDp1KoZhEEURX/nKV/jkJz9Z3+foo4/mlFNOYdKkSbz++utcc801HHPMMSxcuBDD2HCt25Y0SQaYPXv2gDUbhUKB8ePHU857aIGZNPNuTIqrlEs+SimclEmp6BO4Ed3rK8QxlAoeqYxJttFm/MRmOtaUccse+b4quqmhGUmxmTVLK4ShoqHZhiCi4saUiy6aCaER49gGZc9DhQrDMDF0RSXv0ddZolj2cVIGo8c0ksmYpDMWqazNmFQTDc0Ory/pplhyKfd5uG5IxqmgmyajxzZgGDqaAt8PyfdWKZdcbNOgO4jwA59qNcK0dNA0MjmTvzyzmmLepWlEhnG7teAHMcW8i1/18aoBucY0o8c10tVZJvRDCr0VLMeiWvIpF3361ldoakmRa0ozclyOTC6DWw2olv2kn56mYUcGtp0UZ9E0jVTGItvg0LW2iF8N6SwXUTpJwFcKkumYdtKuxNChWg2pFjziOKS7o0yl7FHp88g2OFitKfq6K4RBjOVk6OlM7tf0pNKtbmgEfoQiCRSbmlMUuqv4lYDRuzQlvRYbHCzbwK0EVKoBvhdiWRr5Pg+vEmCZJlEc0joyh2nqxJGi2OfiV0PCOCaVtqmWXdJZm/VrS1imTtX1KJaTSrGeHzFmlyb6usoEboRp6Fg5k7amZJ1cueBh6Ek/x3TWqrX60ChXdFQcE0SKfD758FXq8xmzazN7HdBOrimNVw2Sqb1BRCkfJdVhMxa2k7zWpqmj1BvBYCZr09SSplr1qZSCeg/IwI+wHSPpy2npSb9Pa8uq3A61OIrr7WFSGYtCvsKSl9biVUKaWtK0j20iimJWr8gThhGptEkqa9PdVcXOmjSNzNA2Ik2+r0qx4NHbUwZg1LgcxElQWi15LPrDCkaNaSDX5KAZGlFMEpgriOI4ubDTnKKpKUWp4GOYJp1r+igX/eTnq8EmljaO2wUVb34hm3c67o7ub3/7G2eeeSaTJk0a9P61a9duUFfhndiqq4PXrl3L/PnzmT9/PoZhcOyxx/Lyyy+z995789WvfpUrr7xyaz6dEELsEN7ri28/+clP+NGPfsSDDz7IPvvsw+LFi7niiisYO3Ys5557LgBnnnlmff/99tuP/fffn913350nn3ySI488cqsdi+M4g66d9N2QUuChSJqBVyyPaslDodOQdcjlUvh+QFBbc9fdWSRVm9ZoWhYtI9PkTaiWQgp9LmGg0PQChqEzYlQOwpjG5gyuG2ClNXrXV0lZGs1taQI/mc7Zvb4IsYahgR8ka+2qBcV6VaTU59LcnCWbszFr2bwJk1pZ/vf1eJpHpZhM18vkNErFpMm55wasXZ1HtyCdM3ELAX09FXQNdMeksTlNy6gMr/11PaViFV3XGT+ljdZROfLd5aTiaa06qm5q2LbDLhNtutaVkp6XPVXMERlyLSm8os/y7jKObdLTlWO3qSMZMTqHZZt4bohf9Ql9g6aWNKmsRbng41YDUmkL2zboXleir6uM6wWAhmXq5JpT+H5MoaePSjUgk7JpanEoFpN1gSpUpBttGprS5FpSlCsBtqOTythJ0OtH6JrCdiwCL8L3Q3QTWkZkiSNFV0eRxtY0vh/R0pYl15hCacnawRgNTVN0dpSToCyMMGyDES0pMo6VTItV4HkBhYKLW/UJQ8WIUTn6uitUqx5uOUJDI9eQIoxiHMdgyeK1OCmDdNpm9NgmWsc00LO2SEl3SadNwlBhpYykqIsdoRsOucYUfb1ViBSGZSTZaxWz7BWftct7GTuxhUmTR9I4IlULHmPiCCrFZDqlnTJJpfvbniSl/X0/CWYzGYdUysb3w6Rwjh/huiGGrmHHBmEQUyXAsg3sWiZyuKhWAqhlaKM4ZPHCVXiVkGyjw4Q9RlDJu5SKPipWtLTlaGhN07Gqj3LRI5dzGDEyi46OZUfke/uoVgJsyyDyFZOnjiKbs/nnP9aT76my/B/dOI6B4Zi0jEzT1JpBI8ncZrM2Ohr5bhffDyn3VvHcCF3T0W3QgHgoegOKrU8yjlts3333Zdq0aVxyySWD3r948WK+/e1vb/a47zpwDIKAxx57jO9+97v86le/Yv/99+eKK67gX//1X+uNLB955BEuuOACCRyFEGIT3u7i24UXXrhVnueqq67ii1/8Yj043G+//Vi+fDlz586tB45vtdtuu9HW1sZrr702aOC4JU2SN6WxJY2ubHRNh0jR3VWkt7eKk7FIN9o4pkHbqBwKnZiIwI8JvJBCl4tfDZIshKmRaksThiGVQlJUJZWzyfdWME0Npbk0tqbB1LHzfpJR0wxaxqQo9LpUSj7lqk+p6KJDUv20xSHWtCTbVumjsSmVtOJwTHw3wNI1LF0j0iGOFKEbEWVj8n1lOlcX8VyfTC6FZRsUI48ojol1yNkGhqFY83ovlYJHjMauk1sYMSpLtewR+hG2qRMZOsTJNMdKNcC2dEaNbaDQl7S+KJU8giikbWQO3TEo9Fbwl/dRKSRFSRpbUjS2pnGrSSsD349wshZ+EFIueLhlHydtkWm0k+xXNUBFIalUKmlZ4kUU+1yiWKGjkWt0SNeyr1GUprHFYcwuLYBGFIKug5O10HQdy9JBh56OCtVqgGHoSZuGvEdc6+vXuapAqtbvz3IMco1p7JSRZPFqlTqrFR/DMrFsnVzWIZdzSKctfDcgiKKkomvJJeWYrF7Wm0wZNRWGDk2tWTI5B01XvPZyJ2hgmCbjJjWTa87glQP6eisopbAMHUNXFHpd/CAk3ZCmsdkkChUYOsVeF88LsNIWVtpAxRrFgsvfX17L2qU95JrTjBidI5O1AWo9I6m9FhpO2iKVNtEMPZlKq+mk0ia6pqHpGikn+bdhasSRStp46DGmaeBWfX7zmz/S2dnBLruM44iPfnibLg3yvTApYKRp2CmDZ3+3lErJw7ZNJu7ZSqXPpVDwsAydxvYMuq6zelkP3evL2LZJtsEmDhVBHJDPV0g3JN+DjQ1pWlpTtQJCIZOmjKSnq8TqFXn6equkUBT6NCrlgOaWNON3bwWlUSq6lPIu+Z5K0oszVCgUBiZRpIhDmaoqdmwf+MAHWLJkyUbvb2ho2KLaCu86cBwzZgxxHHPWWWfx3HPPbVBWFuCII46gubn53T6VEELscDbn4tvWChwrlcqA0uEAhmEQb2L+1qpVq+ju7mbMmDGD3r8lTZI3xU4bOIZDOpW0hyjmPXw3RNc11i3tAy2mY5lB04gMI8c00Doim0xL1XVilWQUVAy5Jgtdt3BskzUr8hS6ypiaxlo3JJNxKJVcWkZkaWxJE3ohlbJHseiRzVk0Nacol5OgsVoJaWpNk2tIYZgaPd1l3FJIvq+Cn7Ixa0VQQi+kqTWLU/Hpy7toWtKOolLyCcMY2zHJpK0koPEiMg1Jf8PQi1i3soiuKyI/BhNKRZ98X5WWtjSmbRD6EUGg8KoepaJPuRKQ90OsLp1Mg00mY9PXVaFcUFQLPrmcQyZrE8WKvi6X3u4qURCTylqMHtdIOmsTBCHr1yYFiQq9VXSSi/FtYxowbZPmERncakgYJC08Ai/CtgwMWyedtam6IaWuMoEXYuqQa0ozfo+RrF7ajenoZHI2pmWSyTmEYYTnhrSOyjBqXA4VKvK9VTJpjTBWlItJZrhU8NBMyDWkqZR8GltSdK4uAooISOcctNBHUxa+GeCbyWsTa1Ct+Di2gd2SJQxDSl1Jq5RSOWbsLo20jMiCFrFiSS+ZjE2u0aZ5dJZQQdf6cjJFmmQ6dK4xhe+GOGmbYsHFdyMKPREpx8TWk8AvjCJ68xWKeQ3TMLBtHS3W6V1XxrR0bMsg05gUcMrmUvXv7zhWSasOLanAatg6dtoklbZIpwxMy6y1IDFA0wiCpC2Khsb/PfE4N/7HHNasXVkfb9zY8dx0w1xOOeVkbCcpyqMnvUDQ9SQQHarprfGbqqg6KZOXnltJsaeKbuhM2nMElXxIseBhmjqpnE3oRZSKVfp6qliWQSZnkU5ZKKUo5gM0wyCTMxk3sZHdp4wi9CO6O8t4VR+3EuBWQ0xbp3lkhihI+oE2tibB6JqleSxbp1rxqZR8XDdA06GhOYVpakm7Dx2q1Z0g5bQDULX/DcW4O7pvfOMbm7x/991336Kq7u86cLzttts47bTTSKVSG92nubmZpUuXvtunEkKIHc62uPh2wgkn8JWvfIUJEyawzz778Kc//Ylbb72VCy64AIBSqcSXv/xlTj31VNrb23n99de5+uqrmTx5MjNnzqyPc+SRR3LyySfXA8O3a5K8OVb8o4e0XiaTdShVfYqFCqEXE7gh5YKL0sCxDXq7KwR+RBDEhF6EYSTr8NxCSBhGFHoAdKqFKqEXoSno6y5hGAZ5VcZaa7Ium2f0Lk1ksg757gpRpPCrJnbKJpW2MUyDTJNF4Cn8IET3oZL3CcIIzXbw4whDT4If0zFpaE6hkcFaW6RU8ulYlqdc9ojDmLbROdat6qNc8NA00FMGbj5Zv6dQqFARa6AFGn1dVVy3A9MwaGxJ0zwiQ67BwXc14jgidAPKJRelNAp9VWzLwE6bxKGir7tKT1eFhqYUI8flaGixKRd93EqE3xtS6nOxUjojRzeRbbEhUmiahpVK+t5Vyz6VoocirlUyVZSqSTY33WCTa0phOQZ+NaSv18Ur+5imTn5dmSUvrKRjVQFd12huy+BXA8rlANPQMTSNVM4G2ySKIhzHJPB8gnJAHEMUxkRxSFBWVIs+6YzJ08/+ngcfvYd/PeWzvG//afhhRGNjMjVRR8et+KBrFPsqlEo+tmXiODrlUoim6RS7SzSOyNDbU6K7s0hvR4k4jmkd3YCTStG3vky1kvQTLPW6BGFIQ1OKfF8FW9fxgxivElAoVHCLAZ4bQaRQmiLStWTtXcUDX6HpSWCbzZioWEepmFSPybpleayUSVOtGqimIEYRRzG6BigIY4Wma0mm065laDU96VdYa+fy9KLfMPebV7PLLu/n2OMup7llV/p6l/PySz/hgos+RceKO/jYkcdi2wZOysKyk2nNmmagG0kvTt3UMPSkZYtm6MlUaUOvtwbRaus6dV2vBZyb7ilbKfv1Kaqv/W0dnauLxHHM+N2aiQJFIV8lVuBYRhLEl4LkfUol71NLS5Zy0cereESaQjc1Ro5pZMLEEbX2JgbjJtn0dZd5/ZX1uNUQ0DAApUGuOYWGorujiG5olAsurhsRRRFtYxoZMbaJhiYHyzAIo5iudSW8yN3s30liG0iuAg7NuDsR13V56aWX6OzsHHCBWNM0TjjhhM0aa4sDxzlz5nDiiSfyqU99akuHEEKInd7mXHzrb1Xxbt1xxx186Utf4jOf+QydnZ2MHTuWT3/608yZMwdIso8vvfQSDzzwAH19fYwdO5ajjjqKG2+8ccB6xNdff52urq7612/XJHlzLH91PblsFkPX8fyYwEumNqo4xi36KA3KhoGT8ajkPZb/rQs3SBrAW46VtFCIwTR1DANCLyKutdEIA0UYeahYEavkw3LnyjwNTTboOlEAuq5AB9MyaR2RQdN1ymUfpWLCMFm3pps6gRdi6Rpr+z88axphqHAcHc0w8Mo+fXkPz0vW+61bVaBa8AjjmJaWDG7JJ4wUugZ22qRUSip1anotk1PyiIHedSV6mx0aW7M4GYvAjYhiDUMzcF2PahCjwrgWaIBhaARRRKm3Qm9niREjG7DTJn6pSqHgohs6qbSNV/QxUjamqZHNWcmU0EpAb75C5MXYtk4UKgIvSDK/PVVUTNILrxJQLFTxKj5eJUrW3WWSqZzVUjL119AUStOwHCMJWgww9WQaaxzHRJEiCKMkYIwUQRgSxxAHUdJLMlb8bOFdrOt9nR/++G7sytgkADU10o6FaRqEQUilHBD6yfRDpRShH6OimCCMsdPJR50oSC48BFGMYxk4WQvD1NFIgqKgGuD6IVFMUg23tt3QNTQ0VBwTExN5Cj+OUCHEKkJP2RiGlvS4dEMMPekRaNkalqVjmEYSHKJhmAamoWPWXo/kgCFWMbECrda6I+nbkTyvYWigacRE3PXw1xi3y/s54shr0bRk1sDIUVM54shr+c2Cm/jq7TdiFHbB0JNKxJqeZEJNx0z6TNoapqZjOAaGodXeEz256RqGrqFbSbbSsCwMHXRTq9+vG8m/DUND03WiKHkddEuje12JNSvyKKVoak2T76rQ29WV9F9M1y4QeCrp8+gHWIZOnLbIr+sELckmpxyLXGMaI1L0rivSHdd+DqoB6zuKFIs+bsVHV8n3iFKKVauK+F5IFCXvnWkb2E7Su7Wp0SZl6oRuiBf4VCo+oRvgF4PN/p0k3ntJG5uhGXdnMW/ePM4555wBf6v7aZr23vVxXLVqFccccwy2bXPCCSfw8Y9/nCOPPBLbtrd0SCGE2Glsy4tvDQ0N3H777dx+++2D3p9Op3niiSfedpxly5ZtsG1TTZI3h+uGqNBHNzWiIEKFMbapgwZxrKGZYBqg/Ihy3qXHS3ro6aaG0jWIkyb1jm1h2zphpIhjBbpOXOtjGIURKIXnRRBrdKcNTFMDXUOFgKXh6Dodr3VBbd0ZGsSahq4lBWMMw6DqBknvP6WwDSM5DstAaTHKj1G6hg6EmLjVWoVSTSVZxkARRxG5BhvLtij0VihVoiSQdDx0XSOIkmm3hqnXgwrT0jFrH+KV0ojDGM8PCGqvg1LUW4FomoaKkq81XQEaCpWsHyVZi6kZSdsL0zTRDIVSGqYOhm4QhRGhH+FHMcSAlkxojeOkv2FS+jbJZtm2QRBGhGFyDsmzaWiGgaElF/pVrFAqTuKyECIAFRMlb0fSD5HkPDvLS1jX+zr7HHYuf33mAf7wuycZ27gX9fyXrtXOkXoAqFQMKrkoYOj92xTEECXxGIaZvJf942iAH4KKo9pB90/z1JPTrb2H/RO8o1gRRnHS7qPqYaZ1NE1HaeD5IQRAJcmW6LqOoSVj6loSzCWBYfLBTdd0altqZx2jlEryqRqAjgLW9C0hX+rkAx++uh409tM0nf32P43/+8VVPPnEb2lv2qP/xUfXNGJN1V4uHb3/rPXaK6bV1l4mO0NtD60WeCbvX+2u/tdaB10z6sccBCG+H6Fi0I3kwkIUKXRLx9ANTLP2XsWKMPkxBJWcu24lYzoZG9vWsQ0zCax1DV0HP4Kq69Uv8jhWUiCKKE7WvJZ9quUAP4iS49WTQHXU6Bx/XpdUUQ6CkKiWaUmKKJXf/peQEDuAz33uc5x22mnMmTNniy7ivtUWB47f+c53iOOYP/7xj/z85z/niiuuYO3atXzsYx/jxBNP5Pjjj6e1tfVdH6AQQuyI5OLbplUKVbwoqDWsT/rCGbW2BJqmofsaUfJpHMMgyU5FCqVBHCTTAC1Dw7RrLQ+Ulnwgj0mmgoYqCWKiJLOlatMENU0nVkm2yTL1JOAL+5u3KwxLSz4sx7UPz8mn6eT/TANNJVM+UYogit8ITBREKukLGEdxfU2pipPsYLUQEoVJ4GcenxQf8rfSa6k28u/4LdvD2m1zvLmmpwK82r/7wxrtTfdtbGzjLf+tj6cUf73/O4zc5QAOnXk1nasW81f3j0z65GeHZL3exnP+G9KBjZWi2dwPVhtbWayoBdU15Zfz8DI0t+w66P4tte1eXE4K5cTJ95tSMVqUfN/GqCS7SX92843Xsb//ZHIp4I1ylv2tC7SkHCmangSysUp6akZERCFQu1BhGBr52vRWpSVXDzSSgF5pOhrJhY0kJq6t8TRqQayh1aYgJ2tAMTXQTfp/zIjjpHgQgA4qfuP4NS05H0USrPd2FDFqhYfM2s+2ZiYBcqQ8hNgZrFu3jlmzZm2VoBHe5RpHXdf50Ic+xIc+9CG++tWv8sorr/Dzn/+cb33rW1x88cUceuihfPzjH+ess84a8obWQgixPZGLb5vmlXy0MPnQGiuVZD2qEbX8FYYJlmOgGUkWLopiVKwwdI0gSD6YYuvEtayQAkwj2R+SAFOFMUF/1kjTkimSYVwLHBVBLXLTFEn2R4coSj6/apqGYemEoULFEbptQJw0P0+qNsZouo6ha0RRTByoegsAXdcIDYVhAZpOlPQcIfCTrOhW7ZO1HVv9+h9Zv+rPHHX2fWiaxvs+8ll+9cOLWf36H9ll8ge39eG95zK5kQD09S5n5KipG9zf27scgH0O3J09JowjjpJpwmgaKkrWqgZ+RBQrolAlwWMUo/Tke1zVArYoqmXmoiTIVLWpu/3BomkYhDHEUUS1EuBVk/WphqFhWxoYSfY1jMJaVJxkzYmT6bi6BkZ/RtfQ0QwNPVK1qw1JgKtiRWxQD2xVGEOYPD6Ma6Ft/1WQ5BQxNA1NV2gquZAU+BEYScbVtsxk2roJhq7jR1vrsowYSqp2sW0oxt1ZfOITn+DJJ59k99133yrjbdW/T3vttRd77bUXV199NevXr+exxx7jscceA+Df//3ft+ZTCSHEdk8uvm1crjmNHidZk5Dkw2AcQ+An0yUNSwc7yWtpKsZAR09pSRP32odQxzFBT9Yc6oaWTH0zk6lsuqGhqZggJAn0VJxU5gyjpHS/SjKKSoFl61iWmeRfdA1NgW5oWLZBEMYQJlPy0EjGAFBJARtiRRwm46O0+riWaRDXpj8aVnIefpA0F1HPPZNMXzQ1lB+DDjo6mqbQbTM5D8tANxSEyQdjTA3b0UlnbJRSeEGyppMoxq3ERHEEQVzL8IBpGWhWMgXRMA0cXRGhEYdJ24dIKQwUYazVMr5gOwZWysRA4VZjAj8mDIKkx6IBtm1BrfJnttEmnbExHYtMxkLXNZRRm5qoAxrJ+rgYiMALk6IzYRgTx0klxfkLbmXULgcwbvcPADBu9w8wctz+vPx/X2PqERl0K2lf0dycJtvooNCoFjxMK1m/l83ZaCo5Nkhec8PUSaUdLDt5fTVdo1r0cb0I0zHINaZoG5OjaUQG4phyOcT3Q1Dguz49HWWqXrKWNp2zk0xaDMQK0zaIY0W54uFWQ0pFH98Pif0QN4iJ/AgtToof6Sr5/kzZJrqlESsNTana66NhmjrZBpt0g4NlG1iGhr7PQSz65VhefuknA9Y4QpLN+8tLDzOqbRwfO+ZjRBG1ixjJ66lBvcCOhp4EWLXMndIUqGS6KJpWW0MKhm1iWwZKB0PTQUG55OJ5IcU+j7Wre6gWApSCTM6isTWbBJe1wC5WSZsOzw2J4xjHschmTUzbwrB1MqlamxYz+UGrlD2KFR+vGhEHEX6k0GKFFsdoWrJuNfAiokAl64z9KPn5TOnompGEnHE9T0oQxcRh8j2laq8BybfbG+tLxfAmxXHetTvvvJPTTjuNp556iv3222+Dtj2XXXbZZo33rgPHarWKUopMJgPA8uXLeeSRR9h777258MILt1r5eCGE2NHJxbc3TNhrFEZsodDJZC0yDRZosHp5X9KuQdcwtKS/XRiE6MrAtpMP45GmEcVJ3z6lIA6TSo22ZYCuY1qgDANTT9ajBWHSoF1FimrJRUVg2MkUOts2cDIGlmGRyiQtA1wvxjI0NEMRRkmiRLP0pECHH6GbWlJkx48IgwhN02kekaGp2aHQ56JpGqmMiR/UPrzrGqW8h++FxERJkKjXWqSEyZj9mVcVJYVUMDVsw8AyDSKSTKWKFIVCMn4SvCaFXay0RlxW+JEi8nx0SyeMIlRFYZo6dsrCzFhEKkYpjVTWxLZMiBWVSoCfpGSxzGSNXhjG6KYiZWkoHKJAkUrpGI4GgQmGRipjo5samayJnU6m/EahIgg8/DBZxxjFisAL8WvrJFVcK78fw+vLX2Tl2lfr2UZIsrzvO+JSfvXDi1nT8xfGte5FseJT7Kli2xrprEOuwQZdkWvMJIF8FEMUo2vJWjtD09At0EwwbA2vEhBGEYZVWzMbRfhlHz9tYhgQBQGaiomDiEpfBdOE1lyaTKODX43w/QjPTQoBeW4StZiGhuPoGJpDGFtJD07PJ/TBdwPiICmwVPUDgihEr0ImbWNaBiqOCANFHOmEYUS+r4pjG6TSFnbG5LzTLudr93yR3y64iX33P42Wll3p7V3OX156mFWrnufaWf+JmTIxSS5cxJFK1qFqOqaZ3GxHx7SSwElp1N6bWnGe/imhClRELVMIXuTjuSG+F1EpefR0V9AMk0yTSa7Voak5je9FBEHy8+RVQ1wvmart5GyyGYtcYwonY5JK20k2VKmkiJMX4nsKn2S6d6rBwDZ1UikbJ22gQgiiEDRIOw62Y6L0mHTGorE5jWWbxHHyMxBFMVqsiAKFHwSEXkSl7FMqeHiVgErVx3dDTFvBI+/97zUh3mv//d//za9+9StSqRRPPvnkgGn+mqa994HjiSeeyCmnnMK//du/0dfXx6GHHopt23R1dXHrrbdyySWXvNunEEKIHZZcfBtc04gUhBbFvEex6qE7oBkm6UYHJ+uQzhlkMzZu3iXf7Sbrmkwj+RDqhui6jpO1kqIwTpLJ0ZQiiCLcEDQtJjS0eqXVTCYpcmM5JnEtoLNTJoapYegGtpPM5zNMnYaUiWVp+NUQtGS9opVKqmRqWY1MQxIUlgsuoa/INtmMGJkh05jCTtmkczYjxuQIwuRDrq5rdK7IJ9VF46RYj2npeF5EtRygawrLNogjhVsN8CphctymhkJHQyXnE8aga/V1oWi1NZ6GTjpjoRkRvlsrEhNGWKaJMkCzNIIwSjJdjolhaUk2TjNw9FrRkrSOZVtUKx5eWSMmaYkA4KSS/nxKAUZyPnEYEUc6xdjFquj1IjkqVoRxf2YxJtYgjkgqi2oKE4Vh6Dz5wo9pbBlPKtNC15q/1b8vUpkWGlrG84c//w+nf/j/Q4s1wjCg7ELgJZmo9vHNSXElRT17FmgKLY6wHDtZaxlDpc8n8EMiFWNaJp4fJhVg/RDVU2vL4ke4bkDsJ2toLUcnbZu4lYgwDImi5Pun6kf1qaBxnMyNjuOYIFBJtheAKOnZqYOuKUI/wvOTFjGuV03q5RhJARi9P7Ota8navTiZwtmS3odzTr6Gxxb8F//3i6vqr0trczsX/+t17DLyIFYs602mX+skxWwApWmoMCaMo3pWrnZ37f+9sQ43DlVS8TZO0vxRlLyuQZAUo4nCWg9KXSOdMwl6Y7o6S0Sxqk1npX4Otm1gWBqxAeWKT6HsEasSppZMjVUKwijJIsYqRtd1bFMn1DSqbkC54ieVYXUNyzGpeD4VN8A0NcJI4brJ+mfdqH0YVv3/ietTWY20Rc4yyLYqWmJFGERUqlIcZ3ugGKKqqlt/yGHr//v//j++/OUv88UvfnGD/s1b4l0Hji+++CK33XYbAD/96U9pb2/nT3/6E//zP//DnDlzJHAUQohNkItvg3O9ABXpuH4SAPT1ekRhFd+PkmDATtHYbNI0qpF0a5bIj7Fsk2JflUJPBc3QaB6RIZW2QFc4aQcVxVSKfhJURHGtb52O5STBYRjEeEGY9AWsVU5VtdRdgIZF0rMuilSy3tFPKoOm0xZhnLSu0DQN34swDJ1cY5p0LvkzW3VjyqUipm2CG7J+bQnDNshkbXq6OuitdGM5BgqFqwK0WMcLAsqBR1NjKy3tY1ARlEs+gRuhohCvklSKjMKQlGMl/SNrgVgchdgpOykIUit1ogO+GyYVKL2QKIxrLRoMND1Zu6abBlYt4xhGCk0li9LSOQfT0kmndSqpkDCIiEOVVKqMkoIpoRdiGmCYJlqtooph6knWNCI5Ak0RKQ0MA90w0IB0GizbwrYNdMCtuhTKPRQK63nsvtMG/wbJtqFrMYaVHHPkK6IgpNwTscrtYcToDKPGJL05lUqm6KazNnam1kNRKdINFg1mCssxqRZ9dEMj15yiZWQOXQPXDSj0uYDCSNs4GZOG1gymptWnQ0IyJTSMk3WtURjjekmw6VVDwiDAD+Kkz2gYE4URfhijQoUfhklvUj8k8hWhSqrURkqhRbXp19SmWGpJhdhYxew55XCu2m0aKzpeoer20drSxuTd98c0TLT+skT14LN2pBpopomBkaxvrFWZjWsVaOuFcWJFRH8AmFSO9f2YKEqyjfXMJEnV3VI5BC2s94DUDQPd0DA0HSdjYZoQxuCXQ7QwTqZnaxoKDSOpe4Nu6NgZE8sysCwNx7HfOCZdw7JMnNr7pmpTDG3HxDBr56pUPfgFahlT/Y2quZqGTa1KLMl03VRVpqqKnYPv+5xxxhlbJWiErRA4VioVGhoaAPjVr37FKaecgq7rHHbYYSxfvvxdH6AQQuzI5OLb4MI4yTzkmtPJ+jFLp5ivYqJhGAqv7NO5KpkGajsm2VYbK23hBCGpwAKStY35vIfp6GiGSSpj05JziKOolhWDUl+FMAbd1ElnTBqMNJgk/fhMg3LJo1IOCPwoCQBqU+xCX5GyDVINDnbaJKiE5HtdYqWwbYNU1iSVcjBraxKLJZ9iyUNFPlYqqeZhmjqGXeXzc87A9zfekNxx0nz7jv+ldeRo0qUAt+ITexFRS5Q0o4/jpFVBrV+i6q8wGyeVJsMIQjfJUlqOQcbSsXyT0E8yUFGg0BVEWjJ1VYUaSkvWehqmjuNYmGbSiN73krWalm3jBUkvTMPQUHGEbZtoVjIlNQ4i7IyD7RjoGQfNVKDX+hkaBpauSGUssjm7VmE2ye7FUUQUZrnpmm9TKPYl6zGN5BhsO6lEW8q76HGGlNGIWwrwq8mUz2oxoFryiaoh+fVVUqaBbZk0tKRpbEkxYnQjhqkn2dZ0Eni4bkDkKxpyDk0tKcZNaiXXlMJ3A1Yt7yVlG9iOSVNrhua2DIah86b+2W+yYQ4jCEKCIMYr+1SqPuW8R1/Bo1x28apJb1LPj3DdsNYbNCJwwyTgDEPCoDbVGoUOBHGMFkdEXoCmwZjWKTiWiWbqeCUf0opU2sRJJesikwxmbT1j/wq/WkuQqL89Tf+xRjFxFOP7ydToahCiPC3JTqqYwI/qU7I1gyRgx0jGM5KgNtmuY1tJD0VT1+tldS07WV9qmgaapuoXbQxDw7KSbUlrDoNUxsKxTJyUQTpjJ/1Sg6R3KkAqY2I7Vv/p1CSB6Jvpen97kdr6zVpYrWkaxeLW6YkrhpYUx3n3zj33XH784x9zzTXXbJXx3nXgOHnyZB599FFOPvlknnjiCa688koAOjs7aWxsfNcH+FZz587lZz/7Ga+++irpdJrDDz+cW265hT333LO+j+u6fP7zn+ehhx7C8zxmzpzJ3XffvclStEoprrvuOr797W/T19fHBz7wAe655x6mTJmy1c9BCCH6ycW3waXTFk3NOXSgpTWNbVusWtaNW/DRbB0/iImjCPP/b+/e46Mo7/2Bf55nbrubewgkgICoVFAUKGqMYrWaY1B7Kki9HatIrVQrKKb1QsvF2zlUT1W0ovzsUaqnUi3nKLW1xYNRezOiIAhURaUotpJwzW0vc31+fzwzkywJISRZcuH7fr22dWdnZ5+Z3Qzzne/zfJ+IP36qIALOORxLg2UZsiubJbNxTsKFY7qASEDhQlaP9OejEwxQOWBDA+OAI2TZ/kh2BMITyClQoUcsJJvk1CBayoRteeBcjvnLztbBVQUK5BhHAQY9wsPxenKcpIqcbNm11U46EI4LAQamKti1axcsK4Vzpi9FfslXWh2HupqP8eYzNyBlNSEWPQrZ2YYcFykEGutMOJ4L25JZK9sWcp5IVxYS4VzWyXRsGynLhZOywTUFRkSBEtGhGgJ2yoKregC4nGqEy2IprisQb7Sg6ypigwwYhoJk0oXnuf78kJ7MNgqGrDwDriP7RnIOpDwXQuWwHQ92o4ArbChcQDU4FE1DTq6AHovCYwJNjZbs1ikYwIRfWJOhuHgwjho2DIah+tOpAJomj6vnygI8qaQMhpnrIpVw4NoO6vcmsfPLBiSaLOzalYDleBgZGYTskQWI5cipbhjT4LoOEg0WLNuBnZLj51IJBXV7E6jbk4Rl2kgmbCiagqLiHOQNkEHjodB1/xIrPxYu8zwPluWiqclEY52JxoYEkkkbqaSDZNJCMu4ilbJgJi2kki5sy/GDOcD1bAhPhSoEPEduK266EElZrEepl1NMyMJKHLoWBGscqioDNX9eDAh/TkvPk+Mtg66pwpNjfV3Xg207EB4DUzhUTYWiKTAMBbEcA3pEhR6RWUJFVaDrTHbZVlUAHhRNhaZxGLoKPapCUxVZiViVRZZ0XYOhKRCMwUpaaGy04NgeDEPOJapHFOgRGRyGF/kCMqiMaOGyjl7/7z97iwfjkL5L0kP8aZAyst0jhOu6eOCBB/Dqq6/i5JNPblUc56GHHjqk7XU5cFywYAH+7d/+DbfeeivOPfdclJWVAZAXQBMmTOjq5lv54x//iJtuugmnnnoqHMfBj370I5x//vn44IMPkJWVBQC49dZb8corr2DFihXIy8vDrFmzcMkll+Cvf/3rAbf7wAMP4NFHH8UzzzyDkSNHYv78+aioqMAHH3yASORQZncihJCOO9w33/qKxr1JOAkGParCdWTF02SjDHyKj8pDQ10Sjg3k5msYflwRolEVibiDHdgHCIFotgGuMJhxG3pEFpHxXA/JlINkgwnHcSFcAVVTEYtpgMblmLOkDdVRYLkywNQNFbkFURgRDVbCBvNceLYNJhhisQjgT96elx/B0GG5yM43sHNHI5oaLGgqR0Tl4ApHdo4ip+awXSi6AstyoUcUNMRlMJJf8hUUDR934OPRkED9viRiWQbyC2KymmhBDIkmC6bl+OPlHIBlIdFkQuEcXOMwEzaScQvxhhSSCRWKqiASk9k2zxVwoips04YrANt0YZk2VF2B4wBQOTzGkIhb8Jgu553kDGqWBs8ViBVwKFyVBW2YHNuoaH4BGiHHSLqOA8WVGS/H8uBZFhptB8l6G0xhUBUFmh4UbdGgaQpURUDT5OdZSdlxUtNkV2LGAE1RkFNgYMRXipCbH4URlYWTrKSN+j1J1HyxD59+uBO7v2yEY7uob0iiblcTBg7Okd+j5cBKMCgag+IwmJ6AYaiIZslCNoCcziGSpaNgYAz5hVn+fJ1dxzlHJMIRiWgoKsqGEIVIJm3Emywk4xYs25XtMx2YpoNU0ka80UIi6cBM2jBTJuykA8uWRXQcx4PrF6RxbQ/ChQwALVm4B/48iuHnM78aThCLhQGY8MccemAKh8IYdMOQN0iiCrKyDGTlajCiBoyIimhMQzRLR1ZMhxbRZOVWJmBEVEQMDYrC4Lkyswm/CI7sispbVTRVVQ7NkJnFWEz1s7rNGVHOOTRdkVlMrZu6mHK7e7ZDMksgMwMSj5y4EZs2bQrjsc2bN6e91pn5cLscOH7rW9/CpEmTsGPHDowb1/yP3nnnnYepU6d2dfOtrFq1Ku35L37xCwwaNAjr1q3D1772NdTX1+Opp57C8uXLce655wIAli1bhjFjxuDtt9/G6aef3mqbQggsXrwY8+bNw8UXXwwAePbZZ1FcXIyVK1fiiiuuaLMtpmnCNJsnkW1ooK4PhJBD0/Lm23nnnZfxm299hSxKw+CZHpyIBzPlIBJVkFeUjYLCGFSFw7ZcGDENnHHkF8ZQVKzJ5UNcDBgUg6brcgyXJ6CpClSDw9AVNCUcNO5NoGFfAhCAHtHgOgLC85CIW+E4Ssfx4DmyIEpE19HgxRFv5IhGNGTlR5CVawB+BVbbkmMmNV3H0ccNRLzBkvM1Oh7MpAVF44gV6cgrisHzgHhDCq4j0Fic06HjkVuQhVhuBApXwDVAcA6uANmFUUT88YoCAslGC1lZGlRNRX5RFIwx1O1NonFvQnb38+RUI67twfE8OCkbrhuBmbTR1GQhFRewkh64yqFHVESzDEBlaKgzYZsOoMjxiAJyCgdF5dA0jixdRXZBBK4tx/lFYiqycyOyCqvlwrXlJPGcC6gKh/DkpO9qREEsokNRZSVZznhaZk8IWajFEzJTlZVtILswiqyYDiOiwohoaUHdwCG5OOrYQow8oRifbtqB7Z/uQWODhS+21yGRcDB4RC4KinP87KWc+iK3IIKc3CiGHF0AxuW0LynTgaYpyMrWuy1obAtjDLGYjpg/jYrjyO6iVsqFaduwTVlEyLJdWCkHpunCsuTUJcmkBcuU2VbHceC4cgylvAEgu3Y6tg3XFn6Xv+Yxf8H/ya7DDJwr4JyDc3kDQDdUaBENWVk6svMMGFEd2Tk6IoYKTVfBOIfK5VymnMv9kO/Tmqu3+tPYqKqca9VzPdmOoKsvk91JEwkLUTD5fRrNl6WeP58kP8RML+k/qKtq173xxhvdur1umcexpKQEe/fuxf/93//BsponVa2pqcHo0a0nqe1O9fX1ABBOlL1u3TrYto3y8vJwndGjR2P48OGorq5uM3Dctm0bampq0t6Tl5eH0tJSVFdXHzBwXLRoEe6+++7u3B1CyBHmcN986yuGHVOAnFg2XEcgJy+CXTsaoakcA4fkIBrTkZ0bQXa+AZUzeIJD1XU59122Dsf2kJUdgaqryMo1oOuKrIAKyIqmKseebAMDS/KQPyCK7DzZq8S2HezZHYeVsMEZw55dcdTvSUAIIBJRoemyOqkAw6CjcpCbE4HtCPzzs71INMSxO2EimbBhRHQUDcpGwaAsNNUnkWxSYCZsOfbPkV2HAAbDkFmYjhhYnI2hR+XBdT0oCoMR1cAVLouFRDXY/nx5WsRDstGBYsjxiJwBisKRlRNBJEuHbihwXAHX8QAmYJsuGvYlkYzbyDUtJBptpOIpJOMuFIUhElHAGQd0D6oix65Fc3RYcRtckVOl5OTGcPTxA5CTq2PvriQYYygZmovsXFl4hkFO9+A4HnL8Y+04HhIJO60Sqew+LMK74IrKoeuKfERlYGXbMijyPOEXnDGhaXIcoqxAyhDN0nHUyEIUleRg6DF7se2DWtT+sx71++JwbQ81nzdAj6qIZcmCN9lGBEcdW4hITHZlTSVtaJ4ix6AexqCFsSArp8DvQAXX9eTYR0fOyWhZrlzmeBDBeFt/bGIq6Xdr9QvxeH6JUwEB07Rl8G4LOK6MIIUnwJic01Q3AE3TYBiya2nEUMH84kUC8rtRFA7NUBHxu5pyhUEIf45Iv3KvrittFuFgDOC89W/dtuXfAmNo9bdAASMhvU+XA8e///3vmDp1KjZt2gTGWBjFByd++Q9kZniehzlz5uDMM8/E2LFjAchgVdd15Ofnp61bXFyMmpqaNrcTLN9/DGR77wGAuXPnorKyMnze0NCAYcOGdWZXCCFHsJKSEpSUlKQtO+2003qoNb1Ddm4EBXnZ0HQFtu0hmiUnlJddGRmMmIbBwwvgmi7ijRZSSQuJBhMN++KIxiJy3kFDBhRZWTo4Y0glZPEUVVdgpRwIAFqLDAfnHNGIjljMQG6eLJSyZ2cTUnELnHMMP64Iu3c2IJmw/ekjZDe7AUVyLKbtCjQ1mOCKC6ZyJJM2cgdko2gox+4djbBSsviLosoMaf6AGHbWf9ah45FfGMOIYwbI7qWW7L4ZydLlnI2enOQ9lbSx658NiEY12d1PCDi2gMo4jNwI8gujiGQbMFNyvCaYnPeSCaD2ywaYKRuRiBzr+eXWvajbmwRXGCIxDdFYAVQFiGZFYKYccAaYpoPCoiwMKM5B4aAs2JYHuByaoaB4aF56N6j9hpSpKkdOjg7HUeG6zdVmWFA0qI2gIRgvqmmKzMyZTjhlhmW5fgCphO+NRDWMGluMggFRfP73faj9Zx0a9yWgKnJeRNUPcrLyZREjAH5AKq9bdKPnK28qijwWLQ9fkD30PJnddT2EvwHbktNmhNNpCOFnEtGcaWyxLcZYOP6P+YE3h/xbAJddW4NMqOfJ6TIUP9OocB5Og6GqMtA+VJbfNVjTlYxmdknfJDw0Z6i7ebukc7ocON5yyy0YOXIkqqqqMHLkSLzzzjvYs2cPfvCDH+CnP/1pd7TxgG666SZs3rwZf/nLXzL6OQdiGAYMgwZYE0K67oMPPsD27dvTem0AwDe/+c0ealHPUnUN2XkRRGMaEnELuXkDUTQ0B54rYKYcqCpHotFCLEdHvh5F/b4EahtMWKaHrDx/jjk/AwIA0SwdEDKb5FiurOrI0Dz/G5rrJQRLGGPIK4hBaZFBGTqiEE0NKezzM5E5+QaycnUMG1WIeKOJpoYUHEeOnczJi0KLKEjGTURiBnLyFHBVbn3YyHwwwZG/PatDx4OBQdMVGIaKpnrZzdW1PWTnGmGgpHIOfXgBbFsGUa4fXAkhUDAwhkhMD2/wJhI2XNfzgwqGWLaORKMFQ5dFTIpLcrF3VxOSTTbMlI3cgigGFGeDcYbdXzaCqwy5+VFwzuXchP6k9lxhyMoxOjR2JsiwqarsFuz6XRltWwY/qiqzXPtvq2VmznH8bpmO57/PlWPmdPm6EAJ5RVkYwQE9wmElcyFcD/mFMeT607Uoql/tFTIDJgRkN1q15wPHtgTBHudKm20UQvhzL3poOZiLc3k82wrQOvJ9BRlPMyWzm8E2Y1lyTOOhcvzpSxhrUUiIENKrdfkvtbq6Gq+//jqKior8f4A4Jk2ahEWLFuHmm2/G+vXru6OdrcyaNQu/+93v8Kc//QlHHXVUuLykpASWZaGuri4t61hbW9vqjn7L9wTrDB48OO0948ePz0j7CSEE6NleG72ZAFA0JAfCE9AiGvSIiqLiHOyqaYQeUWGlHGhCjumLxjSZOcrSACFQUJSNrFwDmp+dCkSz5Vg6M2XL96gy4AnWCeaIYy0urDVdCT8PkN0nBw8vAFM4UgkbmqEiGtOhahxGTIcelYVjPFfATFhIxYFUygLjHJFs2aU0EtHgOkB2roacPFlxs67m4zaPQ8vluq6CMYas3Aia6lPwXIHGuhQiUQ2ekN1OGWMoHJgFMMCx5G11PaKmBQuMMWRl6Ugm5TQQcl5KjgGD5HyYtu0iJxpB4aBs7N0ZR/3uBCzLwd5dCcCVrwsTKCxSkJVrIJalw7JcKClHBmyHGAS0DiBlRQzHkcEPY3I8ZBDktiQzXYoMakwZOAZBJCDAuCKr3+ZFMTJqYO/ORqQSDlRDhWW6UDUFuQXyO5BZzCDb2HcDGZm1lcfM9SukyrGiAp4nu4bK48nC7r0HE5yX5DhGBYwBtu1BUeX4Xse2oOkcuq52KHPoeQJmShaoUTXKNpIDoKqqvU6Xz4yu64al5IuKivDll1/i+OOPx4gRI7Bly5YuN3B/QgjMnj0bL730Et58802MHDky7fWJEydC0zRUVVVh2rRpAIAtW7Zg+/btYdGJ/Y0cORIlJSWoqqoKA8WGhgasWbPmiJ1DjRByePRkr43eTFeZrICZcmBEOCJRDZoux6Ul4hb0iPznK8iuCAFk5xooLM5BTn4E2dkGIlGt9Xb9aQQitotEkwXHcQHI9bwWF8ctxbJ1WZzSE4jEZPCZkx+Fosg5JBWF+ZkThqLiHHAux0d6wkMq6fpj9DSoigI1i0PhCur3JZFMWCgcUIBINIY3n7nhgMciEo1h8ODisF2cM2TnRZBoMuFYHlKJ5gqRkZgWVp5sL2Mmi7JocBw5N2LQJRI6EG2xXvHQXMSydOyuaUCi0YIAwBhHVp4GM+UgJz8CrnBEohyO7cLzGFStc2PTmgNI4c8zKLtfyq6SAoAXricDHiAYHyeE391RYTBTsuqoTLZ5MBlDNKoiElMx5OhC1O2JI9FowbWF7F4LmWlMJW2ZbfSLu/R1QQCpqrJKqet6frdg4QeUzevJYwq0PJ4AwoCzZTERmWXUoSgsDLY9T8Dyi/KoGg9vBLQVlHqeQCJhwXPleFajDwfpJLOoOE7v0+W/1rFjx+L999/HyJEjUVpaigceeAC6ruPJJ5/EMccc0x1tTHPTTTdh+fLl+M1vfoOcnJxwDGJeXh6i0Sjy8vJw3XXXobKyEoWFhcjNzcXs2bNRVlaWVhhn9OjRWLRoEaZOnQrGGObMmYP77rsPo0aNCqfjGDJkCKZMmdLt+0AIIYGe6rXR2w0amgvDkFNxaP4k7IzJ7EQ0piMZt2BEVVlYRXhyLGOOnFNRVflBi84EgYHnZ2SUoNAMWnfbY4whlp0+LEHXFbAcw+8yyKFBQDdUqJoMxIyoJqdPMG0wziFc/0KFy4vyeKOJeNyFrhbgrb+sRW1NLZg//Uc0psH1/AI2AAYPLsaor6T/e8o5Q3ZuBJbptMiSKYecKZPHof0gKSc/gqxcHYlGE6bpQNPkcRdC+IGCA8ZYOEG70sWgKwh4AB7eGGienkGEyw5EN+QUH44jg08I+T0n4rbMtubIiri2n0VOJm3Y/lg7VeOIRrVOlanvzWR2UQmDSBlINh9PGSAebCusVfdh3c8y2rYL25LZXlnNVXZBVdTmsZBBZVXblhlQ7ncnp2wjIX1HlwPHefPmIR6PAwDuuecefOMb38BZZ52FAQMG4IUXXuhyA/f3xBNPAADOOeectOXLli3DtddeCwB4+OGHwTnHtGnTYJomKioq8Pjjj6etv2XLlrAiKwDcfvvtiMfjmDlzJurq6jBp0iSsWrWK5nAkhGTU4e610VfoEQ05+VFwzppHafnXl5quQAiZWbNSLhzLk11TNSWc+40d5GKUMZkZCy5yFYWndcc7GE1TYPtjJbNz04NKeYEug8tsEQkzNnLz8r+jMQ0N+1JwbA8DCopRVFAMxuX4wJaFYSLR5nGabR4nQz0s3So558jOiyK7xbJk3IKZcpBoah6Xa0TVbq1E2hxESsGxDL7n4JgGP46WU0MEwU3QddWxm8dQhisDYdCoaQoiUbXfBY0tBd1+FQXQtJYZRaT9RoPjyVjz8Wzv7yIYcyqrwHqwHVcWcrI9OGgdkSoKQzSW2alOSD/Qi+ZxXLJkCf7zP/8TNTU1GDduHH72s5+1W8RuxYoVmD9/Pj777DOMGjUK999/Py688EIAgG3bmDdvHn7/+9/j73//O/Ly8lBeXo6f/OQnGDJkSLiNvXv3Yvbs2fjtb38bxjWPPPIIsrOzD/SxGdflf20qKirC/z7uuOPw0UcfYe/evSgoKMjIybcj6eVIJIIlS5ZgyZIlHd4OYwz33HMP7rnnni63kRBCOupw99roK8Jqjy36zfEW/6bIDKMBMynnMPQcD5Fcw88qdSz4U/1J5W3bhRFRwyxWR96rqByMpWcsD7QPbW3PMDRkxXQ0NZp+dVBPdqP1K31y5cDVRXuLqF/V1UzJLp5cYW12D+5O8ljKwKejgnkFRSTorinHoHqenHszyIh1dGqU/iQIJLtLcxVYNewaG2SMGeRvhHM592d/DtBJ9+gtQxxfeOEFVFZWYunSpSgtLcXixYtRUVGBLVu2YNCgQa3Wf+utt3DllVdi0aJF+MY3voHly5djypQpeO+99zB27FgkEgm89957mD9/PsaNG4d9+/bhlltuwTe/+U2sXbs23M5VV12FHTt2YPXq1bBtGzNmzMDMmTOxfPnyrh6CTmOCOvp2m4aGBuTl5WHnzj0YOLCwp5tDCOkDXn31VcTjcVxyySX49NNP8Y1vfAMff/xx2Gvj3HPPDdcNzjH19fXIzc3twVZnTrCPH2z8DKNGH4V4oxlml/IKoqjflwzXLSjKQjJhIdFkQTMURKNy/KOi8FZZwLZ4nkBjfQoAkJMXQVNDSo6VzDE61N0y0WSFQWdXAibXlWP5+mrgIoTMLClq2xU7CeltjoRzaV8WfD83fftZGHqs27dvWgks+eU1Hf7+S0tLceqpp+Kxxx4DIKe/GTZsGGbPno0777yz1fqXX3454vE4fve734XLTj/9dIwfPx5Lly5t8zPeffddnHbaafj8888xfPhwfPjhhzjhhBPw7rvv4pRTTgEArFq1ChdeeCH+8Y9/pGUmD6du6d+SSqWwceNG7Ny5U44naOFILSVPCCEdcbh7bfQVQQAij0EwT0YbmbuIBst0IbzmLoeK2rHjxnlzd9WgMAqQPkVHezRDkRPSWzJ47Oz3JbM0nXprr8CYnCqEEEK6VYZTjg0NDWmL25pmz7IsrFu3DnPnzg2Xcc5RXl6O6urqNjdfXV2dNs87IP+tX7ly5QGbVF9fD8ZYOCNEdXU18vPzw6ARAMrLy8E5x5o1azB16tSD7mYmdDlwXLVqFa655hrs3r271WuMsSO2lDwhhHQU3XxrjSlIn6n8ADiXXQyDOfyA9quJ7k83VDi2lVYcpaMBoOpn2DxPIBG3oKoKjAhViCSEkL5g2LBhac8XLlyIu+66K23Z7t274bouiouL05YXFxfjo48+anO7NTU1ba4fFPTcXyqVwh133IErr7wyzIDW1NS06garqioKCwsPuJ3Docv/ws2ePRuXXnopFixY0OogEUIIad+qVatw9dVXY8+ePa1eO5JvvgXBW0diOCOqwnHkxO2Ms0Oq6qlpSph1DJ4fShsjUQ2JuAXHlmO5KHAkhJC+4Ysvvkjrqrp/tvFwsG0bl112GYQQYQHQ3qzLo+5ra2tRWVlJQSMhhHTC7Nmzcdlll2HHjh3+vHXNjyM1aAwwpE9OzljbSUhF4YhlG7JgTtahV2qMxnQ5TjGmHXKFUk1XEMvS5VQY1F2TEEK6jfBExh4AkJubm/ZoK3AsKiqCoiiora1NW15bW4uSkpI2211SUtKh9YOg8fPPP8fq1avTgtiSkhLs3LkzbX3HcbB3794Dfu7h0OXA8Vvf+hbefPPNbmgKIYQceejmW/dQVY5oTOvUHIKcy8xhZyci13Q5t2SmK4oSQsiRJBjimIlHR+m6jokTJ6Kqqipc5nkeqqqqUFZW1uZ7ysrK0tYHgNWrV6etHwSNn3zyCV577TUMGDCg1Tbq6uqwbt26cNnrr78Oz/NQWlra8R3oZl3uU/PYY4/h0ksvxZ///GecdNJJ0LT0fzhvvvnmrn4EIYT0W8HNt2OPPbanm9I77Z88bDE9ByGEEJJplZWVmD59Ok455RScdtppWLx4MeLxOGbMmAEAuOaaazB06FAsWrQIAHDLLbfg7LPPxoMPPoiLLroIzz//PNauXYsnn3wSgAwav/Wtb+G9997D7373O7iuG45bLCwshK7rGDNmDCZPnozrr78eS5cuhW3bmDVrFq644ooeq6gKdEPg+Ktf/Qr/93//h0gkgjfffHO/bkWMAkdCCGkH3Xw7sP3HNwbzIlLcSAgh5HC5/PLLsWvXLixYsAA1NTUYP348Vq1aFfYU2r59Ozhv7u1yxhlnYPny5Zg3bx5+9KMfYdSoUVi5ciXGjh0LAPjnP/+Jl19+GQAwfvz4tM964403cM455wAAnnvuOcyaNQvnnXceOOeYNm0aHn300czvcDu6PI9jSUkJbr75Ztx5551pB+1IRPM4EkIO1VNPPYUbbrgBkUgEAwYMaHXz7e9//3v4/EiYeyzYx48/2o5jjxuKZNyGZTkA5LyN9XsT4byOBUVZPdlUQkgfdSScS/uy4Pu58bJlGZvH8Ylfz6DvvxO6nHG0LAuXX375ER80EkJIZ/z4xz/G3XffTTff2tJGSdW0eR0JIYT0WwICXcxvHXC7pHO6fJUyffp0vPDCC93RFkIIOeLQzbcDYwwdmsuREEIIIZnX5SsV13XxwAMP4Oyzz8bs2bNRWVmZ9iCEEHJgPXHzzXVdzJ8/HyNHjkQ0GsWxxx6Le++9N+3OrhACCxYswODBgxGNRlFeXo5PPvnkoNtesmQJjj76aEQiEZSWluKdd97pUls7Mo8jIYSQ/qc3VFUl6brcVXXTpk2YMGECAGDz5s1przH6F58QQtoV3Hx79dVXcfLJJ7cqjvPQQw91+2fef//9eOKJJ/DMM8/gxBNPxNq1azFjxgzk5eWFxXgeeOABPProo3jmmWcwcuRIzJ8/HxUVFfjggw8QiUTa3O4LL7yAyspKLF26FKWlpVi8eDEqKiqwZcsWDBo06JDbyRiDoqTf31RUDtf1Dn2nCSGEENIlXQ4c33jjje5oByGEHJF64ubbW2+9hYsvvhgXXXQRAODoo4/Gr371qzA7KITA4sWLMW/ePFx88cUAgGeffRbFxcVYuXIlrrjiija3+9BDD+H6668PS5QvXboUr7zyCp5++mnceeedbb7HNE2Yphk+b2hoAABEY3IiZj2iwvMEVE2Ry7N0MCaXE0II6cc8IR+Z2C7plE79y7tx40aMHTu2w2Ny/va3v+H444+HqtI/9IQQ0lJP3Hw744wz8OSTT+Ljjz/GV77yFbz//vv4y1/+EmY3t23bhpqaGpSXl4fvycvLQ2lpKaqrq9sMHC3Lwrp16zB37txwGecc5eXlqK6uPmBbFi1ahLvvvrvVct2QgSJjDNEsvcU2GWLZxqHvNCGEkL5FZKY4DvVV7bxOjXGcMGEC9uzZ0+H1y8rKsH379s58FCGEkG5255134oorrsDo0aOhaRomTJiAOXPm4KqrrgKAcCLiYI6qQHFxcfja/nbv3g3XdQ/pPQAwd+5c1NfXh48vvviiK7tGCCGEkAzpVApQCIH58+cjFuvY3CqWZXXmYwghpF/qbK+N7vLrX/8azz33HJYvX44TTzwRGzZswJw5czBkyBBMnz692z6nIwzDgGFQBpEQQkg64QmIDHQrzcQ2jxSdChy/9rWvYcuWLR1ev6ysDNFotDMfRQgh/c6ECRNQU1ODgQMHdmj9srIybNiwAUVFRd3y+bfddluYdQSAk046CZ9//jkWLVqE6dOno6SkBABQW1uLwYMHh++rra3F+PHj29xmUVERFEVBbW1t2vLa2tpwe4QQQgjpuzoVOL755pvd3AxCCDly9HSvjUQi0SrbqSgKPE9WKx05ciRKSkpQVVUVBooNDQ1Ys2YNbrzxxja3qes6Jk6ciKqqKkyZMgUA4HkeqqqqMGvWrG5tPyGEkCOA8B+Z2C7pFKpWQwghh1lP99r413/9V/z7v/87hg8fjhNPPBHr16/HQw89hO985zsAZEGaOXPm4L777sOoUaPC6TiGDBkSBoUAcN5552Hq1KlhYFhZWYnp06fjlFNOwWmnnYbFixcjHo+HVVYJIYQQ0ndR4EgIIYdZZ3ttBFNVdNXPfvYzzJ8/H9///vexc+dODBkyBN/73vewYMGCcJ3bb78d8XgcM2fORF1dHSZNmoRVq1alzeG4detW7N69O3x++eWXY9euXViwYAFqamowfvx4rFq1qlXBHEIIIeRgRIaqqmakUusRggk6et2moaEBeXl52LlzDwYOLOzp5hBC+pngHFNfX4/c3Nyebk5GBPtYW7sbgwYN6OnmEEL6oSPhXNqXBd/Pd7/xJHSt+2ukWHYS//W7mfT9d0KnpuMghBBCCCGEEHLkoMCREEIIIYQQQki7aIwjIYQQQgghpFcRQj4ysV3SOZRxJISQXmD16tXhtBubN2/Ga6+9hpqamh5uFSGEENJDgsgxEw/SKYccOP75z38GAPz1r39td71du3Z1rkWEEHIE+uEPfwhd1/Hoo4/i0ksvxbJly3D++efjsssu67ZqqoQQQkhfITyRsQfpnEMOHP/whz+guroar7zySrvr3XjjjZgzZw7++c9/hssee+yxQ28hIYQcAXRdBwAsX74c7733Hp577jls3LgRF1xwQThPIiGEEEJITzmkwPHuu++G4zg499xz4bou7rnnngOue9FFF+Gll17Csccei3POOQeTJk3CihUrutxgQgjpjwYNGoRbb70VtbW14Lz51Dxjxgxs3ry5B1tGCCGEEHKIxXEWLlyIn//857j33nuRn5+P7373uwdcd9GiRXj33XcxaNAgbNy4EQsWLMA3vvGNLjeYEEL6o2effRYvvPACnn/+eUybNg3f+MY3cMIJJ2DDhg1QFKWnm0cIIYQcVlQcp/c55K6qjuPghz/8IVzXbXe9oUOHIplMAgBOPvlkvPTSS3jwwQc710pCCOnnBgwYgO9///soLS3FL3/5S9TX1+OJJ57Ahx9+iOeff76nm0cIIYQcXlQcp9c55Ok4brzxRgDA9773vXbXe+ihhzB58mR87Wtfw7hx4/Dll18iEol0rpWEENJPffHFFxg2bFjasvz8fNxxxx091CJCCCGEkNYyNo/jhAkTsG7dOrz66qvYuHEjsrKy8PLLL2fq4wghpE8aMWIECgsLMW7cOIwfPz58WJaFRx99FM8880xPN5EQQgghpOuBY1t3ywOxWAxTp07F1KlTu/oxhBDSL23btg3r16/Hhg0bsH79evz617/Gl19+CQDIzc3t4dYRQgghPSNTU2fQdBydd8hjHPc3YsQIFBUV4bzzzsMPfvAD/Pd//zc2bdqEdevWYfr06d3RRkII6bdGjBiBKVOm4K677sJvfvMbfPHFF/jLX/6CY489Fk888URPN48QQgjpEQIZGuLYibYsWbIERx99NCKRCEpLS/HOO++0u/6KFSswevRoRCIRnHTSSfj973+f9vqLL76I888/HwMGDABjDBs2bGi1jXPOOQeMsbTHDTfc0InWd58uB47btm3Df/3Xf+Gss87Cp59+ih/96EcYP348TjvtNOqaSgghnVBWVoZHHnkEP/3pT3u6KYQQQkjP6CXFcV544QVUVlZi4cKFeO+99zBu3DhUVFRg586dba7/1ltv4corr8R1112H9evXY8qUKZgyZUra1FrxeByTJk3C/fff3+5nX3/99dixY0f4eOCBBw6p7d2ty11VR4wYEd4xD1RXV2P69OntzvNICCEEsCwLuq63Wj5q1Cj87W9/64EWEUIIIf1fQ0ND2nPDMGAYRqv1HnroIVx//fWYMWMGAGDp0qV45ZVX8PTTT+POO+9stf4jjzyCyZMn47bbbgMA3HvvvVi9ejUee+wxLF26FABw9dVXAwA+++yzdtsYi8VQUlJyyPuWKV3OOLaF7pYTQkjHZGdnY/z48ZgxYwYeeeQR/OlPf8Knn36Kn/3sZygvL+/p5hFCCCH90rBhw5CXlxc+Fi1a1Gody7Kwbt26tH+POecoLy9HdXV1m9utrq5u9e93RUXFAddvz3PPPYeioiKMHTsWc+fORSKROORtdKcuZxzpbjkhhHTe66+/jvfffx/vv/8+nnvuOcydOxepVAoAMHnyZCxYsAAnnXQSTjrpJAwZMqSHW0sIIYQcHsKTj0xsF5AFPlsWoWsr27h79264rovi4uK05cXFxfjoo4/a3H5NTU2b69fU1BxSO//t3/4NI0aMwJAhQ7Bx40bccccd2LJlC1588cVD2k536nLgmJ2djRNOOAETJkzA+PHjMWHCBAwZMoTulhNCSAdMmjQJkyZNCp97noctW7Zgw4YN2LBhA9555x38/Oc/x86dO7Fv374ebCkhhBBy+AghIA5xPGJHtwvIyuW9uXr5zJkzw/8+6aSTMHjwYJx33nnYunUrjj322B5pU5cDx0O5Wz569OguN5gQQvozzjnGjBmDMWPG4MorrwyX19bW9mCrCCGEkCNPUVERFEVp9W9wbW3tAccelpSUHNL6HVVaWgoA+PTTT/tu4Hgod8td1+3qxxFCyBGpuLi41UB+QgghpN/qRAXUDm+3g3Rdx8SJE1FVVRUWAvU8D1VVVZg1a1ab7ykrK0NVVRXmzJkTLlu9ejXKysq60upwyo7Bgwd3aTtd0e3FcYK75VdeeSXuv/9+rFq1Cjt27AgntO4Of/rTn/Cv//qvGDJkCBhjWLlyZdrrQggsWLAAgwcPRjQaRXl5OT755JODbvdQ52ghhBBCCCGE9F+VlZX4+c9/jmeeeQYffvghbrzxRsTj8bDK6jXXXIO5c+eG699yyy1YtWoVHnzwQXz00Ue46667sHbt2rRAc+/evdiwYQM++OADAAiTbsE4yK1bt+Lee+/FunXr8Nlnn+Hll1/GNddcg6997Ws4+eSTD+Pep8tIVdW27D9ItCvi8TjGjRuHJUuWtPn6Aw88gEcffRRLly7FmjVrkJWVhYqKirALbVsOdY6W9ijKYTushBBCCCGE9Du9ZBpHXH755fjpT3+KBQsWYPz48diwYQNWrVoVxjbbt2/Hjh07wvXPOOMMLF++HE8++STGjRuH//mf/8HKlSsxduzYcJ2XX34ZEyZMwEUXXQQAuOKKKzBhwoRwug5d1/Haa6/h/PPPx+jRo/GDH/wA06ZNw29/+9suHtWuYSITo04PI8YYXnrppTB9LITAkCFD8IMf/AA//OEPAQD19fUoLi7GL37xC1xxxRVtbqe0tBSnnnoqHnvsMQAyDT1s2DDMnj27zTla2tLQ0IC8vDzs2bMPhYX5Xd43QghpKTjH1NfX9+oB/V0R7GNt7W4MGjSgp5tDCOmHjoRzaV8WfD/fPutn0NVot2/fcpL45Z9n0/ffCf0uNbZt2zbU1NSkVXTNy8tDaWnpAedP6cwcLQBgmiYaGhrSHoQQQgghhJCuCaqqZuJBOqffBY5B3+BDmT+lvTla2ptzZdGiRWkThw4bNqyLrSeEEEIIIYT0mr6qJNTvAsfDae7cuaivrw8fX3zxRU83iRBCCCGEEEK6XZen4+htgjlSamtr08rV1tbWYvz48W2+pzNztACAYRgwDKPrjSaEEEIIIYSEesFsHGQ//S7jOHLkSJSUlKCqqipc1tDQgDVr1hxw/pSWc7QEgjlaujrnCiGEEEIIIYT0dX0y49jU1IRPP/00fL5t2zZs2LABhYWFGD58OObMmYP77rsPo0aNwsiRIzF//nwMGTIkrLwKAOeddx6mTp0azqlSWVmJ6dOn45RTTsFpp52GxYsXp83RQgghhBBCCDk8hCcgvO5PD2Zim0eKPplxXLt2LSZMmIAJEyYAkEHfhAkTsGDBAgDA7bffjtmzZ2PmzJk49dRT0dTUhFWrViESiYTb2Lp1K3bv3h0+P9gcLYQQ0l8cffTRYIy1etx0000AgHPOOafVazfccEO72xRCYMGCBRg8eDCi0SjKy8vxySefHI7dIYQQ0g9RVdXep8/P49ib0DyOhJBM6q65x3bt2gXXdcPnmzdvxr/8y7/gjTfewDnnnINzzjkHX/nKV3DPPfeE68RisXY/8/7778eiRYvwzDPPhD09Nm3ahA8++CDtpt3B0DyOhJBMo3kce7fg+7mi9OGMzeP4/Jpb6fvvhD7ZVbW3cxwPnifAOevppvRp8q4QwjtDitInE+S9xv7Hk3OZSSJHnoEDB6Y9/8lPfoJjjz0WZ599drgsFou1WxysJSEEFi9ejHnz5uHiiy8GADz77LMoLi7GypUrccUVVxzwvaZpwjTN8DnNh0sIIQQAIPxHJrZLOoWuxDPAcVxYlgPLcikd3kmO48I0HViWA9t2Ydvyuet6Pd20Psl1PZim2+p4Oo5Hv9FOksfUgdfHx0pYloVf/vKX+M53vpN2I+G5555DUVERxo4di7lz5yKRSBxwG9u2bUNNTQ3Ky8vDZXl5eSgtLUV1dXW7n0/z4RJCCCF9A2UcMyDINHqeB9P0oKoKVJVi9I4QQsC23fBiPLiQDTJltu3CdQU0jVO2rIMcx4XjBAE3Q3DYhBBwHBeuy6BpCmXIO8jz5O8wCLgdx4OuKz3cqs5buXIl6urqcO2114bL/u3f/g0jRozAkCFDsHHjRtxxxx3YsmULXnzxxTa3UVNTAwCtxoQXFxeHrx3I3LlzUVlZGT5vaGig4JEQQggVx+mFKHDMAF1XoetqeHHpOHIsEQWP7RNCpGVpWwbc8jh6cF0PnufBtgU0TaHg8SBkoC2DRkXhUNXmgNtxvDDjaFkudJ2Cx4PxPAHLcsLnqsr7fBfqp556ChdccAGGDBkSLps5c2b43yeddBIGDx6M8847D1u3bsWxxx7brZ9P8+ESQgghfUPfvuLpxThnMAw1vKiUWR/3IO86stl20G2SQdfVtECbMZkV03V5r0NewFNX4Pa4rhcGjaqqtAq0VZXDMIJl8nj29W6XmRQE2ADAOYdhqFDVvn3z4vPPP8drr72G7373u+2uV1paCgBp0yC1FIyFrK2tTVteW1vb4XGShBBCCOndKHDMME1rzpoFGR7SmiwoJI9Ne5kvzpkfPLKw6yoFj60F3SmB5kxjWxhj0HUKHjtCHk/h38ToH12lly1bhkGDBuGiiy5qd70NGzYAAAYPHtzm6yNHjkRJSQmqqqrCZQ0NDVizZg3Kysq6rb2EEEKOHHKYUmYepHMocDwMWna5lGPKKHhsyfNaduc9eHdJGTwq4Xttm45nS0FADTRnatvTVvBIwXg6x2kOqPtLF2nP87Bs2TJMnz4dqto8amHr1q249957sW7dOnz22Wd4+eWXcc011+BrX/saTj755HC90aNH46WXXgIgf0Nz5szBfffdh5dffhmbNm3CNddcgyFDhmDKlCmHe9cIIYT0C5mKGukap7NojONhoqoKhJDdB23bBWOsW8eTCSHSCsr0lrFqntc80eqBpn8IghzOD5wZ2x/nMiCShXQ82DYOGiAdart6SlCt03X97zN4gQFM/g80jYfddtt6f3OX344dkyB4DILGYMxjbzouPUXe2JA3J/pTEaHXXnsN27dvx3e+85205bqu47XXXsPixYsRj8cxbNgwTJs2DfPmzUtbb8uWLaivrw+f33777YjH45g5cybq6uowadIkrFq16pDmcCSEEEJCnoBgGQjyqGdVp1HgeBhpmgwePc87pGIknieQStpwbA+cAwAD44DCORSNw7bdVl1gg8Cqpy50bduFZbXOrioKh6bxMGvTPB2E7P53KILxo0EBGMbYQQPPYJqUlt0xGWNQFOaPA+y5LoiyCm/r7zJoqfAEzJT8HQByjGI0qkJROYyIBs7ZfkHOoe1LEDyaphtmLftLdq0rWt7Y6OuFcFo6//zz28wsDxs2DH/84x8P+v7938sYwz333IN77rmn29pICCGEkN6DAsfDTNM4LEuEWZ3m4iRtc2wXibgdlg6WcZj877hlwXE8RLM0KEpzkBBkH+U8iK5fVObwBZCplB0WEQFkEBtMpxEUbDFNF7rG4QkRBnydCVAUhUOI5uJDjKHNi3shBJJJOy0oa9kux5EBl2nKbG00qh224xVU6mzZ5Tb4zhiTvSpcx0UibkHhDFA5bMuBZXlwbBe6ocK2XMSydLhecxa1M0FOc+bRCbsB9+WpJrqq5TyXh3pjIxOCccCEEEIIIYcbBY4Z0LJc//462iVQCAEz6cBMyW1xhSES1cJAwrQcOAmZObNSDrJyDBiGCsYYGPwusf70FcGk75GIekjdOTsjCFaBYFoSpcW8lnIsY5Dxa4pbAIBYTOvSVCWqysOg9EDdgJNJOdl9cPxbHnPXlUWLgvkjPQ84HEm24Pu3LTfMKgaVTgEWZkUt04aZcsEZg6oyZOdqgBBIJB24jgfLdP1996CoDEZEhaZ1/k+7q92A+4v0qXS6nnmVU8nIGyWMHXoXact00FhvdqkNhBBCSF+RqUI2VMah8yhwzIBEwkI8biEW09q8MDxYl0DP85BosuA68petRxQYET8oZDKgcJMCmq7CMh04tou6vQmomoJIi0yZwjkYANMP1BJNJhSVIxLRwotXReXdNhedzNjJQDcS0VplqoKKqLquIpmwYVkOhAAsywPndnhx3plMX3vdgC3L8bORMpO4f5CqKHL/DUNNu7jPFNlGmS0Mzl3Bx1kpB/EmM+x6m0rYcFwPwhPgnCGWrcMTMsDUNQUuZ1BUDst00NRkIuIHncFNhAMJtn+gY90yk9vRbsD9TZAB7uy+u66HVMqBmbJhmS4EmscgB39/mq4gYmjQDeWAf4NCCCTjNmzLpfH8hBBCjhwUOfY6FDhmQDJhY9/eBJoaOGJZBgxDAVeau0UKL+i2KbsoCghwJgM4x3GRjNt+8AIYUQ3MZrBtmZ3zhECiyYJpOlA1jmhMh2O7cGwXnuMhFbeg+RnNVEoGTAICriODISOiwra9MKh1HA8m/K6NKk+bPuRQJRImGutMcEUGQMIT0A0VekRBNKqH6wkhwPzKqPX7UmiqT8F1BVRNBrWKymHoKoyoAk1TOxRICiGPl2W6EJ4HM+kXhmEMyWRzZlNR2t+WDCI7tfsHJbO/zWMYhRBwbReOI28eeMGJTAikUg4STaYskCMAParCMFQkEjYScRueK8e7aoYGK+WgqSklt5NykJ2tYU9MR3ZOBLEsPfxuXc8DRFAYSH5UML5TUThULf0GQstKwM2Zt94RPHp+lli0HODOmoMyzhkY73yRKHkDobmbaqLJ8p8HrwtYKRsCAkKWLEIkqkGPyN9sfV0SDXVJOG7LYw2oCgfjDIrfvlTSQYNIQVHkTZVIVG2uLMzkd5WMW/BcAc459MiRl/klhBBCSO9AgWMG1O9NwLZUCL/7Y3aOgWhMh6oAtuVCM1QoqgIEXTdNF4rKwfwLX0CuZxgqbMWVwQ5jsFIuTMtGKunAEx5UVYNjuvA8F1bSRv2+FEzLgQCgaSpUncMTHhgAzuR4woZ9CbgewBmQna0jEtUQiepQVfkZskwNg6IAwu+y6QoB5gFclYVtGBgUf7yXEEDDvgTq9ibRUJ8CA6CqgJl0ZMBjqIhkaVA1FbEsHbEcA4o/Ri8Rt6H6gZxjy+NgJh3ohoJGIduoqhyRmIZoVIMeCcZyIsy8mkmZiXGDYAxCdv0UQDIuvw/TdOA6LqyoAa4koRkqVJWDcw6uyOCCK9z//0OvsOp5Aq4/D6Xs6tpc4VYIGbTb/nQOaeM9Hbl+EPBwAEzhMhPpeFAVDkDAMDgcVyDemIJrC7ieB9d1kIzbiMdtmEkHlmnBsgWE68GxPNm1WeeIRCPIyVMRjUUQy9URixkwYhpiMR2aroIpgOf5NxBMgHH5W1EUWQwGAFzHg+vK46ppClRNAWcA62Q1Wtd18ec//xk7duzA4MGDcdZZZ0HpQLRupmzYfvGgtICxHYxBfrfB9xx8536w3Fb7gy7ElunAc7y04LO+LoGmeguJuAUh0scbarqKRNJCImGDMQ7NYIhEDESjKlSNAeCwUg5sS4619YQ8xq7nwrE92Q3dE1BVBYqiQFEAx/LCG0jRbA1mKtGh/SaEEEL6OpGhqqodvYYgrVHgmAGW6aCoSJUZB8eTgaTpoHFfEq7rQdflRaDws0YQDLYjwBRAU1WoKkNWji67WzKGeIMJz++j1lBnwUxZ4BxwEg5M00FDfRKJuA3LcWAlXdi2A1VTYEQ1RCIq4I97FGBwXQ/JhOkHh3JKB91QEI1p0AwVwnXheQzJhAUr5fhBJ8BUDs7kexSVQ9EUuLaDZJMNKHIfFD+wbNxnyW6ynCOZTKK+LoGsnAjMhI1dNY2Ixy0wAFk5EWRla8jKjUBVZDDpOB6slA0GBtuWwRj2AICAoijQDDk+0UzYSDTafrYWUDQFmsERjepQVC6zegxoajRhJmxoqoJU3Pan75EX51xlUP0uwkFmiQHgKgPAwjFoXGHNXXpVJSw+5PpjSNMrtAKuJ+C5HlwXzd0TgeZMcxBoCk9mkQwFmsqRTLmo2x1HU2MKru1CCAaVA6Yt4LgOAAbHcpFKWmFXX8vyYFo2hCu3CyEr7jLIgFbxA2I1okFTObJzDOQVRFBYFEVWThRZOREoGodwZXEg7ndf1vx9DboxCyHCjKiqNGcmwwyf0jLTJwMizgHG/YDcD75efPFF3HprJbZv/zw8ZsOHj8DDDz+ESy65pM2/J9tykYxb4c2B4DgHAWFLws+mep4I/9t1PLgHGHbMuMy4tgwobUsWIgJj0P3KtImEiX27E/73IgvlqJoGrjBYlot4k4XGeBxNDRYcR45HNSIq7LiDJtXvgm4L2J4HO+nAclw4KQem7cK2BCA82TNAeBBuMP4Rsvu5oSKWbUDTOEwnecDzDiGEEEJIJlHgmAHbPtoJs8FDKm7LLoiuQE5BDJ7tAozBc+JQNAW248K1Hbh+RU/XERCukOOeslRomiYDEzBEohxGTAWgAFxAU+UYtHijDTNlwfEYGDw/O6nAdTwkGkwkG03ougotKscWOqYNzwVsxwbjHE2NKXge4LgemJ/lcxwhX3cAweQyhTNoigLGActyYaYceBBQuALOBVRNhaZzPwMqoGsK9IgcM8g9BtPdCytpy6I+tgdFURDJkmO7dF1FLEtDVnYUWlSDoshuf5qqwPODNK5wmMKDU5dEY10Kru1B1TiMiAx8FduDZTLE60wwhYMrgO14iDeYAGPIyYmAK0wWk7FdWIkkbNOGqivwXAbLdmAl7eZgFTILKrNbLgQ4OBNgjIOpLAxcdF2FasgxqJGoBk1XwRUFqiq7va559y+4/8F7ceftC3HaKWfCMh1YKReuJ8CYrK5bvzeFvXvjMBtN2JYHD4DiR+wMABMCri2DINuRlXRdW2Yq/TgYjAO6rslAmDM4rgXXApKmC7ge3EYTXAB7uYBgABOym2s0piErP4KcHNm1VY/ofgDF5XeuyoBKVRRwVd7Y4IqcBkTTlANm7bDfIs4ZXvn9b/Gd66/GUcNOxYUXzUZ+wQjU7fscmzb9Gt/61rfw/PO/xiVTp/rvZ/BcD8m4Jcf2wS8QFdGg+uMBD1ZcRggBzw0ywJ78b1eEwb7wH44nAMjvPJm0EG8yIQRgGCqSDKjfl4RjubLLr6YgryCCWLYOxhjqG0wkzSa4jotkwoGqyGPmWh4a9iZgJuXv3XFlYAhPdm0FRHiDgoHB8s8Truv6XVsFFM6hKrK4VTxuw9AVWC5lHAkhhBwZaIhj70OBYwZsfusLfKrvhuPILBHgZ2T8C0LH9ABFXuwrjMucFGcyi2gJmYlCMMYJcDwPDLL7INM4DF2RRVBUBY7rj4XUOQxDg6LJi3xVYfAEIDiHyhiYIlviMEBTOJSICithIWW5sJIOXNeF48qCIJ7tgDMGzS+iI9viT/GRcOB5LhgY6ht3I2k1QgiZXXNtB57twfUEFBWIaDnI1gph2U548Q/IQyJkQhCqwsBUDlVVoegMqqpCMxREoyqiMR05+VHoQTVYDsTrTdiWAzCGaJaOaJYhg2yDyyyrpkHVFTAAiaQFM+FAj6hwdAWu68jiM/UmzKQcR6ppClIpF5ZpAwBs00FTg4lkwoJry2BDFvAxkTI9uJYLV8BPH/rfLpffH2McOufguhyvpugcT/z6bmz/8iPcc/dd+PHNjyE7Pwrd0CE8D/EmU2aKk7Y/llGAK4Cqa+Aq989ssotoojGFlOnCsWQ3TU3hAOfgGve79DJ4totU0oZty2q8tunBcR24LgA/iAJk2wWTvw+VA4xxOYelzqBqGhRDQUxXoEZV6FENmsKhGfLGA2cMzOCIqCr0qCKzpYr8zoLuoIwBqqqCczkuT1UVCHj40fw7cdRRp+Lr584DYzJjOXDQaHz93Hl44/X7UHnrD1A64WtQNQ1CCKQSNhgEmKIgElURieky8DbbSB/6ASRjLf67jedgMnuswB9vHGSPbQ/xJlnUCn4X1L21TWhqMOEJD/CAWI4OITxs29KAvbuSqK9LyuJFtgMBBk1VwPwMs/DHQ9q2rHbrODKbzjmHoity3lBFgaYxQOHQOOCy4F9ImfHNzo/AEQJ20oXruGiybDTFKeNICCHkCCGQmaJwndjmkiVL8J//+Z+oqanBuHHj8LOf/QynnXbaAddfsWIF5s+fj88++wyjRo3C/fffjwsvvDB8/cUXX8TSpUuxbt067N27F+vXr8f48ePTtpFKpfCDH/wAzz//PEzTREVFBR5//HEUFxcf+g50EwocM6B2+z7oPNVcNsMPLDzXg0xmMSgAmD+siyn+5aYrAFe+SzCB/btgMzBwlYErCpjCwIP3ewDzAxnO5DrMDzyFfCPgAi4A4coxj1D8cXVc9msUQsBzZKDEGAPTZaEcRVPAITOSniv7mgsBJFJ78cuqH8FxDjw9gKoYuKTsbmRHBwCCNf+dimAkJWDZLoTngjFbdhdlAIcKrjUnrRgHWFBQUgVULgMWI6qCqxyqrkLTAZVrftAmx4jJqUwYVE1eiHueBwgBRVWhKgoURcjPh59ltRyYcQu268J1ZVDhmC5SSReWa8ssnycQXN8LwcJjB7QIUvzg6Ytdm7D9y49w4unT8be3n8HKF1/BiEEny8JHihy/yVUZFGXFdESzIuBZCjxHwLQdmEkbqYTlj+tz5Y0Af9xpwnWax5/6vyfH8WQVVLu5iAuEzDDKLtHyt+AJ2W0SArAA+Vvz5HNw7gd+siqvoilQOQNTVWhRBkVRwfxqvZrqd09W5ffFFdm9VXZzZTAMBaqmAgz48OP12LHjH7jwojlh0Bj+rhnH2JMuxR9euQ0rlr+M8SedCtdyITiT3bqz5PfFNQVcUcAVDsOQmVAjosngVJNdo8OxjIc4zYWZdJA05XjRxoYUbNORbfBkBth1XPzj831o2JuEbTkyewkBJvy/Vy6z7SoXcP2usoD8e2RMyGrHCoPnyDGqHIBwHMRTAsJzAQ/QojqiEQ5dV6BFNRQMiMCyPCQ1B8mkDTPpQThuu/tCCCGE9B+9I3J84YUXUFlZiaVLl6K0tBSLFy9GRUUFtmzZgkGDBrVa/6233sKVV16JRYsW4Rvf+AaWL1+OKVOm4L333sPYsWMBAPF4HJMmTcJll12G66+/vs3PvfXWW/HKK69gxYoVyMvLw6xZs3DJJZfgr3/966HvcjehwDEDGvcloMBrHoflX8PKzKL8by4AMC7H6HEuowHBwot7AH7RFHnRLwBwjUN4HK7nQFgCwvbCC36Fpb8PLdL74Vgv4XdHVYLsjAxwPA+y76InwMDD7CRXuHwvwkRIGJDuTeyE45g4Z/pS5Jd8pdUxqKv5GG8+cwNSJxyH7CEnhM3a3/7NDrQ3zbkDINXO66Gs/Z63jFe8/T5E9R+xjmz44IQQWPvUYgw8ahxOq7gdO/+xAe/+owrFFdeF3Sstf90mALv33wDz29JN7TkUAs1t69Sbbf/RIjm29QvZxTK/YESbbyvwl1f97j188Z6swMsgb6pw+N1xWYssov93woK0L4LiRko4ppJx2cWaBZlQRQmfMya7ADtBZtDx4Nge4DGAC7iOK7uOCw+eI/92ZIZSfjbXZddtVVPAOQfjMiAXfnZXyEo3UIIxoGow9pRBeAKu7SBhuvAcV/79ieAUIHseKLoCrioyMwrA81x4rgfTpIwjIYSQI4Pw5KVxJrZ7KB566CFcf/31mDFjBgBg6dKleOWVV/D000/jzjvvbLX+I488gsmTJ+O2224DANx7771YvXo1HnvsMSxduhQAcPXVVwMAPvvsszY/s76+Hk899RSWL1+Oc889FwCwbNkyjBkzBm+//TZOP/30Q9uJbkKBYwYwt3ncl8xMhK/IAE3I7B9jHrgCuMwNC6fI8Wr+xa0f4HFVAddkFVHh2PAsWZQFHvO73wHuAYaZCeGPU/QzfADCSp7BHRchZBEXz3PDLqTB/womZIAlgj80GUgmTZlpzC/5CoqGj+u+g9dP/HPrX7HrH+/j/G8/CcYYJpxzE/7vlzPxz61/xVHHTerp5h12seyBAIC6fZ9j4KDRrV7ft08Wy7EbNOxGI8CEHw/6v9sWAxLCn7qAn51PXyZvmngt0u0yOGQC8g6LP6TRE0J2KQ1u7Ij9bm4wFnZzBZOZeM7kNhhkQCgACFc0B7BBFprJrukMftDpeeENHAYBz2vedNDsVpjcttwfuYbtHTjDTwghhJCOa2hoSHtuGAYMw0hbZlkW1q1bh7lz54bLOOcoLy9HdXV1m9utrq5GZWVl2rKKigqsXLmyw21bt24dbNtGeXl5uGz06NEYPnw4qqurKXDsT2IFBhQeAWNyXBnzgiImcgwb40xebDJAgPvBmgtP+Bem/jQIXFH8b0hA2LLAh2t7YTdDAGF0KNDi6tnPxvhhYZglBIOco1D2qZQX18H1NRgURQtf86/K/avtIEBl4UsppXlexvbYa9bCzN4ZdqMF8xN94UV9cKUPP2UTPJXdYllQIhTwK3fKTKlMiwbHEH4lzba68bXYH39PmB9MM3+sYnN7RFqWNliJc/jjPHna8WWMhZVGw464TL53/eYHMXDoyRh67JkAgKHHnomBQ0/G+pd/gpKTfwgviOpdQPhfgvCruzI/1cz8uSk1XQc3FCh+qk0AYJ7MZrkQcC00l5ZWAdWfA1B+CMKKvHDkd8gU+Bkw5ndj9f+fyalCwhsdnmyjX9YFwvXnIA2OExNg4a3AICXdHAIJT4RHfqDwkB0pwsaNv8a55zWPcQRkkLdp4wpkR4owuPh48PA1mbVjAASXfWz9RF6YmGdgUIJuuJw134CBnO5Gtrbl99Y85Qxv/lmE4x1F+DsQ4XIIQAGT1XABcL/LNnhwM8X/voDwxpD8DXvNB1Ow5jjUY8GfuF9tNmhE86EM2sHUIHCVBYvUjHTZIYQQQnqfYPqyTGwXAIYNG5a2fOHChbjrrrvSlu3evRuu67YaV1hcXIyPPvqoze3X1NS0uX5NTU2H21hTUwNd15Gfn9+l7XQ3ChwzoPCYfGgwZEbDvxIXYNB0FXpEgRpRwMGg+GPUXNsLx6VxlUHTFcgxgQKM+0Gj48A0hezaZntyYvHgIt6WUzt4QmY+RNDlVbDm8W3+tB9BptHzL+oh/BiMcSiMyaIsXIFwXTh+IMaEvOLlnENRAQ4OHm8APjz4sSgenochxUVgngyyXD+rKbMvzfMdKgoD4xwQDJ7nwHFkUOX44w2BoEciC6u/CgByeJgAPA8Cqt99MQiU/a7AjAEKAw/mTOSQhWX8eRtdPwsk/MBVuPKYwJ9mIri4VzTFDzoZRHixL8fTMfiBMQf+/o/12NXwd5z/zSfDzDNjDBO+Pgv/98uZiKufYVjJyXCC4juuzEAJz4PnyvSuwuXUJ1zl0HQOIyLn24wY3J9zUgWDB8sf/2hbcr7PaMyAEVFlZVkOqKoC1VAhXP/3pXNEDUX21A2mC/GCMEwEYan8Xbryd+d6TFYAduTNDdeV80x6EHKKCb9bp+t6fndND26QYfN/Y8IDzppwFf5Q/Sher7oPJ518KQoKRmDfvs+xaeMK/OOLd/Evp3wf2flReVj9mwnB+1vGS6Ll/4eRfvMjuCchbzD4YWIYn7EWmT75fQU3JoLvlEP+PfD9pmkJ/n5Yi+MkX/PXCfqUQ8jjL794qEx2uVVVmfdnnEEzFEQiGhRdgWoo0FVZXIhx+bcgx5oyuI4AV/2bTYLDSsXxwrqD/90RQgghpH1ffPEFcnNzw+f7ZxtJaxQ4ZkDhwFyoqi6zN6YL13PBOIemy/FXjiMzj1xRoHgAIrLbHBOAEVEA/4JR4XJS8UhE9TNUMkPhpFw4tiMDMduFnXJgpvyCLn5AqGgMelQBZ3JOQ8t04XpyjKMIxmxBQFGaK7syRZHv01U5MbkrA1UIwHHcMMsUydLB9nVs8N2JE4fhK18ZBV1XAeHC8QQcy4VjubLipOXAdeV+uLYL03Lh2C7gCXj+JbpjObBMWclVCMB1Hbi2zDAG8wSCM8Dz4Dl+IOMCXAXA5DQVnAmkkq4fuDIoOgP3qwupisw+McjMomLIi/ZgUnqZwRJ+lkpWgWU86MnIwzkKOeQcjm//4SXkFgxDJFaA3V9+EB6LSKwAuQXD8M7HL+O0snNkUBomdxk8y4bryCDS9QSY58G1XRnkQkBVNRhZKooG5yA7x4DrCTi2TE8VFkVRMCgHOfkR5BVGoUfldC62ZSPeYKKx3oTnODCiepj99sKKn/KYe64skMT8Aj8t52RUNA4mmP+ZDhzbgWkL2KYD23RhO0ElV/k7cy0XpiuDR9cREI6HYaMqUFSSjd9W/Rf+8Mpt4XHJyxmES867FaNHnh5W8fE82Zbm8bl+xhMtMqhpAWXrO5KsObZr0eW0OQfN/Ok8ZLDm30jQGTjj4FwB48LPzHKojPvVf+XfserP7anrMuBTuQJVZ+EYYUWRlX0Vf1ylosuAH0I2RdNUaJHm069/n8NvD4OiMihclcv8sZqGocK0ksDCDv3pEUIIIX1a0BsoE9sFgNzc3LTAsS1FRUVQFAW1tbVpy2tra1FSUtLme0pKSg5p/QNtw7Is1NXVpWUdD3U73Y0CxwwYVJKF/Px8aBENdkpeuOuGClXnMJMWUqZfRVFTkJ1rwIjp0FQGPapC1RQ5PYMnkJ8fgRGVmctU0gZjsniNYzngHDBisqIkANimjaZGC2ZCzkWoawqiWQYUhcGyHKQStqwG6XmwTReqqsCIaVAUBjNlI5V0kEo5SMbl9AJcALH8KPIHRKDrGlzLxe6dcZhxE64rYBhah47FhLLhOPNrJ4bjI11PZrgcT0C4MnBJNqXQuDeFun1J2LYjK08KD4amwfVcpJIOXM8FV1VwIZBosmBaLoQ/H5/nCD+IBITHYLtynkkhBAxDScv6uZYLz//vVMqBovgTwCsKDIMjtzCGWJ6OiKEDENB0FY4f4Hqe8KedkIFlOEegkPvkCQHLthH/+T407KvFy09e2uYx0ZXB+OqZwxGNGX5hFQbFzyq5DmAlbaQaTSTjFpLJFMyEi6ZGE05K7mOiwUQqYSOSpSM3N4IRxw9EbkEUgwbnIJZjpM1taKZUMMaRnRdFLFtmIz3/uLl+pVw5L6QH13Fh+YG7bbtwbE9+b0KEwaaiMGhaFhSNQdN4WIyG+YGm53r+e124rgfHCcbtAq4rMNkZh1t+8B1Ur3kLNTW1yM8txPGjxkPTFAj/+3E9eSPB87OewhP+vItBN1k/i+h393QcB8K/syJEWMsYYEJ271U49IicSkRVAK6oUDUVRoTLiqyGAsWvJgu/eI3CAC+YG1ORvQAihpwzVNUUqLoqu6H7x0UE6fuwj7ispir87r/CcdDYYMHz5Bym0SxNZrj9pCjjDAJ+V2UmK9pyhUHTFOi6AlWVc1c2NqaPxyCEEEL6LQ/7FR/oxu12kK7rmDhxIqqqqjBlyhT5ds9DVVUVZs2a1eZ7ysrKUFVVhTlz5oTLVq9ejbKysg5/7sSJE6FpGqqqqjBt2jQAwJYtW7B9+/ZD2k53o8AxA44+tgg5OdlQFQXRLB3ZuQYcx4OZcOB6HhzTwZ59CaiKgoGDs1AyNB+qyqEbKnbVNAIACgZkoXBgFsyUg0TcQrzJgqIwxGI6OGcwstIDN8YYLFPOl2iZDhzLRSoup3IIuoNGs3VoupwywfMEUnEbruP6U3xwuJaLnTuaULcnAcaEnODeUBDN1gFwDBlRiHh9CrbtIrVxJwBZPbUtwXLGGKyk449P9DN7CgfnAh5nUDiHrsaQmxfDkJFyTKiiysxMENgIyGkzACAa1aFwBtf1YKUcpBKWnGTd8WCZNpJJG67rIZpjICtLQ6LJwt7aOFIpGwwMBQOiiGTraNyXQqLRhGW6yMrRkFuYhcFH50M3NCQa5QTwuQUR6LoKI6LKwF9TOvT9bzj9XezcubPFeDmkdcstKhqIoUOPCtfnTAYOAPxMrwfHcpBI2NhTG0e8KYV4Qwpmwka8PoVU0gGEQF5hDCVH52PwsDwUDspu1T7bcpFoksVUIjENhp/h4gqX80W2sT+eJ8fSyoDShWO6sEwXpunAddLPtLYtwBwHXJXdOhVVzqWZZciReK4foAItxmACALIx/JhvAgAMQ0UkpkFRFf/34q/i/4cQLdOKTLbNk92/w6BdBO1u/66kEDJYDqoJa7oCTVP9rCpLu8HA/LGqtu3B8QPhlttnjEHLUvzfBU8L1vf/TM8TaKxLwjA0GcTn6nKsZYvmMtY8Bljxb0yE43lbsJyO3bAhhBBCSPeorKzE9OnTccopp+C0007D4sWLEY/Hwyqr11xzDYYOHYpFixYBAG655RacffbZePDBB3HRRRfh+eefx9q1a/Hkk0+G29y7dy+2b9+OL7/8EoAMCgGZaSwpKUFeXh6uu+46VFZWorCwELm5uZg9ezbKysp6rDAOQIFjRhQPzkXUyAIYkJ1nIH9AFoQQsFIOFIVDUWSXv6a4Bc/1oOkq8gfEYCZtDCjOAWdA/gA5l0Q0SwdTGBSNQ+EM0ZgORZUXlkL42SJXThWgarL/ZHauAVVVZPbHcmXmLaKGmbKA5/ldDf1gEwwoKM6BJwTMlA0r4SDRlEIibsvxcRrHUaMKoaoKooUnwngkijefueGAxyESiWFgURFcx4PbxpztLXGFIRLVoUfUNi/C43HTb4cIuw6q2Tpi2bJIj+t54IzDiKkwTRdWyoGmcQgBjBjlwYzbaGxIhd0dFUVBdl4EisKgGyqKSnKgGyoa62XQaERURKIaNEOFqiphYNcRw4YNazXg+lAEXV9zNAVZOYbsvuq4MtNq2vBcASspM455hTHk5EXSgkAhBCzTRTIuJ9UwIiqisY4VM5JTWciiTLqhpk1p4jgeLMtp7p4adG915dyEtuUCsMP5LIP98Dw/4PSDMw4ZyGZl62HAGBSmCYvUhFVPm4vWHCwwDMjxsizM2jHm110SzQGopnEoCm93O0xhMBQeBtyOLYNo23LheR4s04FlOuBcjl/WdEWOT27BdTwkE/LvR9dV5ORFZJdVQgghhPQJl19+OXbt2oUFCxagpqYG48ePx6pVq8ICONu3bw97AALAGWecgeXLl2PevHn40Y9+hFGjRmHlypXhHI4A8PLLL4eBJwBcccUVANIL9Dz88MPgnGPatGkwTRMVFRV4/PHHD8MeHxgTmeg8fIRqaGhAXl4etn36T+Tk5MqxTaoCz/NkVzc/Wxior0+iqT7ld0eTF6eMwQ9omn+A8bgls2hRDVobWaIgk+I46VVE5Dg9Hs5r1xFCCDQ1WWE3TyvpIpWSmS6uckAIqLoCI6Lioy1bUfPlzjDj4rkeFEVmTvMKoigZXIyjjjpKZg295sAgIIvTcL+raPsX054n0ORnz2IxDRAyo+bYblomzBNAImGBc4bc3IifXfO7UTouGhtMNNan4LkejIgGMMB1RJjpUjUFRlQG8q67f996v1trMP7vECaZ78hxD7uPtgiQVFWBqvJwHduSQbFtp1eQ5X4xHcZZeLMAkN2hs3KNbm1ryza7rgfXEWGBHsfxZHAeBIJhtRr53chMdsvsHvfb3vHj2TxPIwsLHMllPO13HhxP1/VavJdDaydD2FHNQaTTKhMZ/J5dp/mzGWPIzouE32VnBeeY+vr6g47L6KuCfayt3Y1Bgwb0dHMIIf3QkXAu7cuC7+fi4xZBUyLdvn3bTeE3n86l778TKOOYAbkFURQWZh90vVhMl11GE7Y/VQdHNKa1yoQEF6YHCv4YY/4YqPTgQ2YkgwCjOeCRQU/rbnDBtnRdgWnKyqY5+RHkIIKi4mwIAViWCzNlIxm3UVIyFCUlQ6Ey2XWUc4ZYloH8ATGZrQp0w6+Mczney7ZlUZ1oVAuzbJ4/VtKxXSQTlj+FApNFdvYLsDRNQWFRcxpNVquUGTrFL3wSy9L9gK1l8CGDctcVcN3mNgUBZGcCSRFULg27WTYHIIoiC7G03CZjMjuqG6rMFluuzP753ShbTkei+OP6jANkcLtD8LtTVcDwv+TgdxeMSwyeyzGK/u/TE/D8LLkcEijguCIMumRWXlaF5X510aA7acs7evsLAtmWn9uyrZqmdPgGysGomuL//nT5PViOn4kU8gaO0xww6oYCPaJ1OWgkhBBCjiQy4ZCZ7ZLOocCxB2maAk1ToORwGIYKw2j9dQTdUIPsYXtaXninByVyigf532nvCAPIsOikP9ec43hwHC8MjlpehMturwyKJrsDRmMaGGTXxkx2wwsCR8eRwUFYzVThYZdCVVegaLKQiKErzUFEizFwqibHear7dVcMjnNLQddNVW0OhIIgsq1xdfJYph/Plieolt0xW2NQVRaOcWsP5wyGHxgKIfzslvyutTCoOfzkd8KADnx88Ntsndlt5o+QhPCY/5t2045pMI7Ua56Mc7/2yAC8uwLGtrTsoirH5QZ/O3KMcKYCd0IIIaRfo8ix16HAsYdpWpDd89DW9DHBHIaHeuErs0FBgRHujxHz4Hlokd1q3X00wBlgOx5M00YkEhTkaM5YCg8wIoCuKy1ezyxVlYGe61fubDvQlsV1IhE17P7bUe1d4DcH5YCmNRdjCTNsLSaMP5Te30G2MuhS3BmMMT8D1qm395jmoLy5q+7+2cKgy2tbQWHbMteduCOCbKnWsSGlhBBCCCF9BgWOPSwIHIOxWK27qcoU4cGyje0JMmDBvIXA/sVIgCCIDESiGkTCgRCAqvCwIE/AsmS1m/a6DmaCpikHCRy9w9Ku/ceNtj6erQPyYLqHllleykZJLQPzwMF+o8H76HgSQgghhGQeBY49LH3sntsqQAwyjorSvRfEzd0p296upilhl1BPCKj7XZBnql0Ho2kcpsn8sWRe2rixlpmqw92ugx1PcujomBJCCCFHrqBmQia2SzqHqjX0AkGlVNv20n7MQbEPoGsZx87S9eZ2tXQo4y67myxywv12pRe+aXmsKPNECCGEENJ3BUMcM/EgnUOBYy8QFO8I5mUMtBzfmMniHgdulyzssf+UBs0BWs8EZ0GgHWREW7eLftaEEEIIIYR0J7rC7iWCYMiymrNoPR0Icc7CrqAts3vN3VR7pl1BAZL06UYQBt09FdASQgghhBDSX1Hg2EsEgWPL7F5PZ/ZatqtlN9qeDmhlu+RnB4F2y6kxKONICCGEENLHUV/VXoeK4/QSMrunwHFcWJYLw2A9ntkDmrvRep6AbbvgnMPzmidr7ymyGq0bFskJgsae6tZLCCGEEEK6D03j2PtQaqYXaS5G4yKVsiGECLtl9oZ2maYL05TTcKhqzxagaVkkJ5Wyw+lBtB6a+J4QQgghhJD+jALHXkRVeVrhF6A5aOtJuq6GYwqDbqq9pV2MsbCbqgwme75dhPR2Rx99tD/dSfrjpptuAgCkUincdNNNGDBgALKzszFt2jTU1ta2u00hBBYsWIDBgwcjGo2ivLwcn3zyyeHYHUIIIf2QnMc5Mw/SORQ49jKGoYZdLRWFp81T2JMMQw0zjJqm9HgWFJDdUiOR5nbpukLdVAnpgHfffRc7duwIH6tXrwYAXHrppQCAW2+9Fb/97W+xYsUK/PGPf8SXX36JSy65pN1tPvDAA3j00UexdOlSrFmzBllZWaioqEAqlcr4/hBCCCEk82iMYy/DOUN2thHOk9hbqCpHTk7va5emKdA0pde1i5DebODAgWnPf/KTn+DYY4/F2Wefjfr6ejz11FNYvnw5zj33XADAsmXLMGbMGLz99ts4/fTTW21PCIHFixdj3rx5uPjiiwEAzz77LIqLi7Fy5UpcccUVmd8pQggh/Yt38FV61XaPAD2fNiJt6q1BELWLkP7Fsiz88pe/xHe+8x0wxrBu3TrYto3y8vJwndGjR2P48OGorq5ucxvbtm1DTU1N2nvy8vJQWlp6wPcETNNEQ0ND2oMQQgihrqq9T78OHJcsWYKjjz4akUgEpaWleOedd9pdf8WKFRg9ejQikQhOOukk/P73vz9MLSWEkJ6xcuVK1NXV4dprrwUA1NTUQNd15Ofnp61XXFyMmpqaNrcRLC8uLu7wewKLFi1CXl5e+Bg2bFjndoQQQgghGdVvA8cXXngBlZWVWLhwId577z2MGzcOFRUV2LlzZ5vrv/XWW7jyyitx3XXXYf369ZgyZQqmTJmCzZs3H+aWE0LI4fPUU0/hggsuwJAhQ3rk8+fOnYv6+vrw8cUXX/RIOwghhPQulHHsffpt4PjQQw/h+uuvx4wZM3DCCSdg6dKliMViePrpp9tc/5FHHsHkyZNx2223YcyYMbj33nvx1a9+FY899thhbjkhhBwen3/+OV577TV897vfDZeVlJTAsizU1dWlrVtbW4uSkpI2txMs37/yanvvCRiGgdzc3LQHIYQQQnqffhk4WpaFdevWpY234ZyjvLz8gONtqqur09YHgIqKinbH59DYHEJIX7Zs2TIMGjQIF110Ubhs4sSJ0DQNVVVV4bItW7Zg+/btKCsra3M7I0eORElJSdp7GhoasGbNmgO+hxBCCGmXyOCDdEq/DBx3794N13UPabxNTU3NIY/PobE5hJC+yvM8LFu2DNOnT4eqNhfYzsvLw3XXXYfKykq88cYbWLduHWbMmIGysrK0iqqjR4/GSy+9BEAWp5ozZw7uu+8+vPzyy9i0aROuueYaDBkyBFOmTDncu0YIIaQ/oMCx16HpOLpg7ty5qKysDJ83NDRQ8EgI6RNee+01bN++Hd/5zndavfbwww+Dc45p06bBNE1UVFTg8ccfT1tny5YtqK+vD5/ffvvtiMfjmDlzJurq6jBp0iSsWrUKkUgk4/tCCCGk/8nUeEQa49h5/TJwLCoqgqIohzTepqSk5JDH5xiGAcMwut5gQgg5zM4///wD/uMZiUSwZMkSLFmy5IDv3/+9jDHcc889uOeee7q1nYQQQo5QmcoOUtzYaf2yq6qu65g4cWLaeBvP81BVVXXA8TZlZWVp6wPA6tWraXwOIYQQQggh5IjXLzOOAFBZWYnp06fjlFNOwWmnnYbFixcjHo9jxowZAIBrrrkGQ4cOxaJFiwAAt9xyC84++2w8+OCDuOiii/D8889j7dq1ePLJJ3tyNwghhBBCCDniUMKx9+m3gePll1+OXbt2YcGCBaipqcH48eOxatWqsADO9u3bwXlzwvWMM87A8uXLMW/ePPzoRz/CqFGjsHLlSowdO7andoEQQgghhJAjkxDykYntkk7pl11VA7NmzcLnn38O0zSxZs0alJaWhq+9+eab+MUvfpG2/qWXXootW7bANE1s3rwZF154Yac+l36PhBBCCCGEdF4QN2bicaiWLFmCo48+GpFIBKWlpXjnnXfaXX/FihUYPXo0IpEITjrpJPz+97/fb98EFixYgMGDByMajaK8vByffPJJ2jpHH300GGNpj5/85CeH3vhu1K8DR0IIIYQQQkgf1EsixxdeeAGVlZVYuHAh3nvvPYwbNw4VFRXYuXNnm+u/9dZbuPLKK3Hddddh/fr1mDJlCqZMmYLNmzeH6zzwwAN49NFHsXTpUqxZswZZWVmoqKhAKpVK29Y999yDHTt2hI/Zs2cf+nHsRhQ4EkIIIYQQQkgbHnroIVx//fWYMWMGTjjhBCxduhSxWAxPP/10m+s/8sgjmDx5Mm677TaMGTMG9957L7761a/iscceAyCzjYsXL8a8efNw8cUX4+STT8azzz6LL7/8EitXrkzbVk5ODkpKSsJHVlZWpne3XRQ4EkIIIYQQQo4oDQ0NaQ/TNFutY1kW1q1bh/Ly8nAZ5xzl5eWorq5uc7vV1dVp6wNARUVFuP62bdtQU1OTtk5eXh5KS0tbbfMnP/kJBgwYgAkTJuA///M/4ThOp/e3O/Tb4jiEEEIIIYSQvkkIAeF1f+GQYB7iYcOGpS1fuHAh7rrrrrRlu3fvhuu6YXHNQHFxMT766KM2t19TU9Pm+jU1NeHrwbIDrQMAN998M7761a+isLAQb731FubOnYsdO3bgoYce6uCedj8KHAkhhBBCCCG9i8jQ1Bn+Rr/44gvk5uaGiw3DyMSndVplZWX43yeffDJ0Xcf3vvc9LFq0qMfaSl1VCSGEEEIIIb2LyOADQG5ubtqjrWCsqKgIiqKgtrY2bXltbS1KSkrabHZJSUm76wf/fyjbBIDS0lI4joPPPvvsgOtkGgWOhBBCCCGEELIfXdcxceJEVFVVhcs8z0NVVRXKysrafE9ZWVna+gCwevXqcP2RI0eipKQkbZ2GhgasWbPmgNsEgA0bNoBzjkGDBnVll7qEuqoSQgghhBBCepkW6cFu327HVVZWYvr06TjllFNw2mmnYfHixYjH45gxYwYA4JprrsHQoUOxaNEiAMAtt9yCs88+Gw8++CAuuugiPP/881i7di2efPJJAABjDHPmzMF9992HUaNGYeTIkZg/fz6GDBmCKVOmAJAFdtasWYOvf/3ryMnJQXV1NW699VZ8+9vfRkFBQfcdikNEgSMhhBBCCCGkVxEZGuN4iNM44vLLL8euXbuwYMEC1NTUYPz48Vi1alVY3Gb79u3gvLkT5xlnnIHly5dj3rx5+NGPfoRRo0Zh5cqVGDt2bLjO7bffjng8jpkzZ6Kurg6TJk3CqlWrEIlEAMjxls8//zzuuusumKaJkSNH4tZbb00b99gTmBCHevjIgTQ0NCAvLw+7d+/FgAE9dzeAENI/BeeY+vr6tAH9/Umwj7W1uzFo0ICebg4hpB86Es6lfVnw/ZyXMw8qi3T79h2RQlXjffT9dwJlHAkhhBBCCCG9ihACIgM5R8qZdR4VxyGEEEIIIYQQ0i7KOBJCCCGEEEJ6l0wlBinh2GkUOBJCCCGEEEJ6Feqq2vtQV1VCCCGEEEIIIe2ijCMhhBBCCCGkV+kt03GQZpRxJIQQQgghhBDSLso4EkIIIYQQQnoXIZCRnCOlHDuNAkdCCCGEEEJIr0LFcXofChwJIYQQQgghvQtNx9HrUOBICCGEEEII6VWoOE7vQ8VxCCGEEEIIIYS0izKOhBBCCCGEkN5FhP+Tge2SzqCMIyGEEEIIIYSQdlHGkRBCCCGEENKr0BjH3ocyjoQQcgT65z//iW9/+9sYMGAAotEoTjrpJKxduzZ8/dprrwVjLO0xefLkg253yZIlOProoxGJRFBaWop33nknk7tBCCGknxKeyNiDdA5lHAkh5Aizb98+nHnmmfj617+OP/zhDxg4cCA++eQTFBQUpK03efJkLFu2LHxuGEa7233hhRdQWVmJpUuXorS0FIsXL0ZFRQW2bNmCQYMGZWRfCCGEEHJ4UOBICCFHmPvvvx/Dhg1LCwpHjhzZaj3DMFBSUtLh7T700EO4/vrrMWPGDADA0qVL8corr+Dpp5/GnXfe2fWGE0IIIaTHUFdVQgg5wrz88ss45ZRTcOmll2LQoEGYMGECfv7zn7da780338SgQYNw/PHH48Ybb8SePXsOuE3LsrBu3TqUl5eHyzjnKC8vR3V19QHfZ5omGhoa0h6EEEIIRAYfpFMocCSEkCPM3//+dzzxxBMYNWoUXn31Vdx44424+eab8cwzz4TrTJ48Gc8++yyqqqpw//33449//CMuuOACuK7b5jZ3794N13VRXFyctry4uBg1NTUHbMuiRYuQl5cXPoYNG9Y9O0kIIaRPExAQIgMPihw7jbqqEkLIEcbzPJxyyin4j//4DwDAhAkTsHnzZixduhTTp08HAFxxxRXh+ieddBJOPvlkHHvssXjzzTdx3nnndVtb5s6di8rKyvB5Q0MDBY+EEEJIL0SBIyGEHGEGDx6ME044IW3ZmDFj8L//+78HfM8xxxyDoqIifPrpp20GjkVFRVAUBbW1tWnLa2tr2x0naRjGQYvuEEIIOfK4MPvUdo8EFDgSQsgR5swzz8SWLVvSln388ccYMWLEAd/zj3/8A3v27MHgwYPbfF3XdUycOBFVVVWYMmUKAJnZrKqqwqxZs7qt7YQQQvo3XddRUlKCt2seydhnlJSUQNf1jG2/v6LAMQMYYz3dBEIIOaBbb70VZ5xxBv7jP/4Dl112Gd555x08+eSTePLJJwEATU1NuPvuuzFt2jSUlJRg69atuP3223HcccehoqIi3M55552HqVOnhoFhZWUlpk+fjlNOOQWnnXYaFi9ejHg8HlZZJYQQQg4mEolg27ZtsCwrY5+h6zoikUjGtt9fUeCYAYpCgSMhpPc69dRT8dJLL2Hu3Lm45557MHLkSCxevBhXXXUVAEBRFGzcuBHPPPMM6urqMGTIEJx//vm4995707qVbt26Fbt37w6fX3755di1axcWLFiAmpoajB8/HqtWrWpVMKcjdJ3+eSKEkCNVJBKhwK4XYkIIKi3UTRoaGpCXl4f6+nrk5ub2dHMIIf3MkXCOORL2kRDSs+g8Q0jn0HQchBBCCCGEEELaRYEjIYQQQgghhJB2UeBICCGEEEIIIaRdFDgSQgghhBBCCGkXBY6EEEIIIYQQQtpFgSMhhBBCCCGEkHb1ucDxxRdfxPnnn48BAwaAMYYNGza0WieVSuGmm27CgAEDkJ2djWnTpqG2trbd7QohsGDBAgwePBjRaBTl5eX45JNPMrQXhBBCCCGEENJ39LnAMR6PY9KkSbj//vsPuM6tt96K3/72t1ixYgX++Mc/4ssvv8Qll1zS7nYfeOABPProo1i6dCnWrFmDrKwsVFRUIJVKdfcuEEIIIYQQQkifovZ0Aw7V1VdfDQD47LPP2ny9vr4eTz31FJYvX45zzz0XALBs2TKMGTMGb7/9Nk4//fRW7xFCYPHixZg3bx4uvvhiAMCzzz6L4uJirFy5EldccUVmdoYQQgghhBBC+oA+l3E8mHXr1sG2bZSXl4fLRo8ejeHDh6O6urrN92zbtg01NTVp78nLy0NpaekB3wMApmmioaEh7UEIIYQQQggh/U2/Cxxramqg6zry8/PTlhcXF6OmpuaA7wnW6eh7AGDRokXIy8sLH8OGDeta4wkhhBBCCCGkF+rVXVWfe+45fO973wuf/+EPf8BZZ53Vgy1KN3fuXFRWVobP6+vrMXz4cMo8EkIyIji3CCF6uCWZE+wbnUcJIZlyJJxLCcmEXh04fvOb30RpaWn4fOjQoQd9T0lJCSzLQl1dXVrWsba2FiUlJQd8T7DO4MGD094zfvz4A36WYRgwDCN8vnv3bgCgzCMhJKP27NmDvLy8nm5GRuzZswcAnUcJIZnXn8+lhGRCrw4cc3JykJOTc0jvmThxIjRNQ1VVFaZNmwYA2LJlC7Zv346ysrI23zNy5EiUlJSgqqoqDBQbGhqwZs0a3HjjjR3+7MLCQgDA9u3b++WJqKGhAcOGDcMXX3yB3Nzcnm5Ot+vv+wf0/33s7/sX9GoIzjX9EZ1H+7b+vn9A/9/H/r5/wJFxLiUkE3p14NiWvXv3Yvv27fjyyy8ByKAQkFnDkpIS5OXl4brrrkNlZSUKCwuRm5uL2bNno6ysLK2i6ujRo7Fo0SJMnToVjDHMmTMH9913H0aNGoWRI0di/vz5GDJkCKZMmdLhtnEuh4zm5eX125MtAOTm5tL+9XH9fR/7+/4F55r+iM6j/UN/3z+g/+9jf98/oH+fSwnJhD4XOL788suYMWNG+DyYKmPhwoW46667AAAPP/wwOOeYNm0aTNNERUUFHn/88bTtbNmyBfX19eHz22+/HfF4HDNnzkRdXR0mTZqEVatWIRKJZH6nCCGEEEIIIaQX63OB47XXXotrr7223XUikQiWLFmCJUuWHHCd/QdEM8Zwzz334J577umOZhJCCCGEEEJIv0E5+m5kGAYWLlyYVjCnP6H96/v6+z7S/vV9/X0faf/6vv6+j/19/4AjYx8JyQQmqBYxIYQQQgghhJB2UMaREEIIIYQQQki7KHAkhBBCCCGEENIuChwJIYQQQgghhLSLAkdCCCGEEEIIIe2iwLGFP/3pT/jXf/1XDBkyBIwxrFy5Mu312tpaXHvttRgyZAhisRgmT56MTz75JG2dmpoaXH311SgpKUFWVha++tWv4n//93/T1tm7dy+uuuoq5ObmIj8/H9dddx2ampoyvXsAumcft27diqlTp2LgwIHIzc3FZZddhtra2rR1emofFy1ahFNPPRU5OTkYNGgQpkyZgi1btqStk0qlcNNNN2HAgAHIzs7GtGnTWrV/+/btuOiiixCLxTBo0CDcdtttcBwnbZ0333wTX/3qV2EYBo477jj84he/yPTuddv+3XzzzZg4cSIMw8D48ePb/KyNGzfirLPOQiQSwbBhw/DAAw9kardC3bF/77//Pq688koMGzYM0WgUY8aMwSOPPNLqs3ri+wO6Zx/37NmDyZMnY8iQITAMA8OGDcOsWbPQ0NCQtp2e2sf+fi6l8yidR4Heex4F+v+59Eg4jxLSKwkS+v3vfy9+/OMfixdffFEAEC+99FL4mud54vTTTxdnnXWWeOedd8RHH30kZs6cKYYPHy6amprC9f7lX/5FnHrqqWLNmjVi69at4t577xWcc/Hee++F60yePFmMGzdOvP322+LPf/6zOO6448SVV17ZJ/axqalJHHPMMWLq1Kli48aNYuPGjeLiiy8Wp556qnBdt8f3saKiQixbtkxs3rxZbNiwQVx44YWtvqMbbrhBDBs2TFRVVYm1a9eK008/XZxxxhnh647jiLFjx4ry8nKxfv168fvf/14UFRWJuXPnhuv8/e9/F7FYTFRWVooPPvhA/OxnPxOKoohVq1b1+v0TQojZs2eLxx57TFx99dVi3LhxrT6nvr5eFBcXi6uuukps3rxZ/OpXvxLRaFT8v//3/3r9/j311FPi5ptvFm+++abYunWr+O///m8RjUbFz372s3Cdnvr+umsf9+7dKx5//HHx7rvvis8++0y89tpr4vjjj0/7G+vJfezv51I6j9J5VIjeex4Vov+fS4+E8yghvREFjgew/8XAli1bBACxefPmcJnrumLgwIHi5z//ebgsKytLPPvss2nbKiwsDNf54IMPBADx7rvvhq//4Q9/EIwx8c9//jNDe9O2zuzjq6++Kjjnor6+Plynrq5OMMbE6tWrhRC9ax937twpAIg//vGPYVs1TRMrVqwI1/nwww8FAFFdXS2EkBeFnHNRU1MTrvPEE0+I3NxcYZqmEEKI22+/XZx44olpn3X55ZeLioqKTO9Sms7sX0sLFy5s84Ln8ccfFwUFBeH+CiHEHXfcIY4//vju34l2dHX/At///vfF17/+9fB5b/n+hOi+fXzkkUfEUUcdFT7vLfvY38+ldB6V6Dw6rtXy3nIeFaL/n0v7+3mUkN6Cuqp2kGmaAIBIJBIu45zDMAz85S9/CZedccYZeOGFF7B37154nofnn38eqVQK55xzDgCguroa+fn5OOWUU8L3lJeXg3OONWvWHJ6dOYCO7KNpmmCMpU2aG4lEwDkP1+lN+1hfXw8AKCwsBACsW7cOtm2jvLw8XGf06NEYPnw4qqurAcj2n3TSSSguLg7XqaioQENDA/72t7+F67TcRrBOsI3DpTP71xHV1dX42te+Bl3Xw2UVFRXYsmUL9u3b102tP7ju2r/6+vpwG0Dv+f6A7tnHL7/8Ei+++CLOPvvscFlv2seW+vu5lM6jdB4N9JbzKND/z6VH2nmUkJ5CgWMHBSecuXPnYt++fbAsC/fffz/+8Y9/YMeOHeF6v/71r2HbNgYMGADDMPC9730PL730Eo477jgActzOoEGD0ratqioKCwtRU1NzWPdpfx3Zx9NPPx1ZWVm44447kEgkEI/H8cMf/hCu64br9JZ99DwPc+bMwZlnnomxY8eGbdN1Hfn5+WnrFhcXh22rqalJu9gJXg9ea2+dhoYGJJPJTOxOK53dv47oyDHItO7av7feegsvvPACZs6cGS7rDd8f0PV9vPLKKxGLxTB06FDk5ubiv/7rv8LXess+7q+/n0vpPErn0UBvOI8C/f9ceiSeRwnpKRQ4dpCmaXjxxRfx8ccfo7CwELFYDG+88QYuuOACcN58GOfPn4+6ujq89tprWLt2LSorK3HZZZdh06ZNPdj6junIPg4cOBArVqzAb3/7W2RnZyMvLw91dXX46le/mnYceoObbroJmzdvxvPPP9/TTckI2r+D27x5My6++GIsXLgQ559/fje2rnt0dR8ffvhhvPfee/jNb36DrVu3orKysptb2P36+7mUzqN9S3/fP6D/n0uPxPMoIT1F7ekG9CUTJ07Ehg0bUF9fD8uyMHDgQJSWloZdibZu3YrHHnsMmzdvxoknnggAGDduHP785z9jyZIlWLp0KUpKSrBz58607TqOg71796KkpOSw79P+DraPAHD++edj69at2L17N1RVRX5+PkpKSnDMMccAQK/Yx1mzZuF3v/sd/vSnP+Goo44Kl5eUlMCyLNTV1aXdiaytrQ3bVlJSgnfeeSdte0Eltpbr7F9hr7a2Frm5uYhGo5nYpTRd2b+OOND+Ba9lWnfs3wcffIDzzjsPM2fOxLx589Je6+nvD+iefSwpKUFJSQlGjx6NwsJCnHXWWZg/fz4GDx7cK/bxQPr7uZTOo3QeDbbTk+dRoP+fS4/k8yghPaF33drsI/Ly8jBw4EB88sknWLt2LS6++GIAQCKRAIBWd4wVRYHneQCAsrIy1NXVYd26deHrr7/+OjzPQ2lp6WHag4M70D62VFRUhPz8fLz++uvYuXMnvvnNbwLo2X0UQmDWrFl46aWX8Prrr2PkyJFpr0+cOBGapqGqqipctmXLFmzfvh1lZWVh+zdt2pR20bZ69Wrk5ubihBNOCNdpuY1gnWAbmdId+9cRZWVl+NOf/gTbtsNlq1evxvHHH4+CgoKu78gBdNf+/e1vf8PXv/51TJ8+Hf/+7//e6nN66vsDMvcdBueYYIxdT+5jR/X3cymdR+k82hPnUaD/n0vpPEpID+m5ujy9T2Njo1i/fr1Yv369ACAeeughsX79evH5558LIYT49a9/Ld544w2xdetWsXLlSjFixAhxySWXhO+3LEscd9xx4qyzzhJr1qwRn376qfjpT38qGGPilVdeCdebPHmymDBhglizZo34y1/+IkaNGnXYpuPo6j4KIcTTTz8tqqurxaeffir++7//WxQWForKysq0dXpqH2+88UaRl5cn3nzzTbFjx47wkUgkwnVuuOEGMXz4cPH666+LtWvXirKyMlFWVha+HpSRP//888WGDRvEqlWrxMCBA9ssI3/bbbeJDz/8UCxZsuSwlOjujv0TQohPPvlErF+/Xnzve98TX/nKV8LfRFD9r66uThQXF4urr75abN68WTz//PMiFotlvIx8d+zfpk2bxMCBA8W3v/3ttG3s3LkzXKenvr/u2sdXXnlFPP3002LTpk1i27Zt4ne/+50YM2aMOPPMM3vFPvb3cymdR+k8KkTvPY921z725nPpkXAeJaQ3osCxhTfeeEMAaPWYPn26EKK5TLOmaWL48OFi3rx5aWW2hRDi448/FpdccokYNGiQiMVi4uSTT25VUn7Pnj3iyiuvFNnZ2SI3N1fMmDFDNDY29pl9vOOOO0RxcbHQNE2MGjVKPPjgg8LzvF6xj23tGwCxbNmycJ1kMim+//3vi4KCAhGLxcTUqVPFjh070rbz2WefiQsuuEBEo1FRVFQkfvCDHwjbttPWeeONN8T48eOFruvimGOOSfuM3r5/Z599dpvb2bZtW7jO+++/LyZNmiQMwxBDhw4VP/nJT/rE/i1cuLDNbYwYMSLts3ri++uufXz99ddFWVmZyMvLE5FIRIwaNUrccccdYt++fb1iH/v7uZTOo3QeFaL3nke7ax9787n0SDiPEtIbMSGE6FhukhBCCCGEEELIkYjGOBJCCCGEEEIIaRcFjoQQQgghhBBC2kWBIyGEEEIIIYSQdlHgSAghhBBCCCGkXRQ4EkIIIYQQQghpFwWOhBBCCCGEEELaRYEjIYQQQgghhJB2UeBICCGEEEIIIaRdFDgS0gcxxrBy5cqebgYhhPRZdB4lhJBDQ4EjIYfo2muvBWOs1WPy5Mk93TRCCOkT6DxKCCF9j9rTDSCkL5o8eTKWLVuWtswwjB5qDSGE9D10HiWEkL6FMo6EdIJhGCgpKUl7FBQUAJDdn5544glccMEFiEajOOaYY/A///M/ae/ftGkTzj33XESjUQwYMAAzZ85EU1NT2jpPP/00TjzxRBiGgcGDB2PWrFlpr+/evRtTp05FLBbDqFGj8PLLL2d2pwkhpBvReZQQQvoWChwJyYD58+dj2rRpeP/993HVVVfhiiuuwIcffggAiMfjqKioQEFBAd59912sWLECr732WtoFzRNPPIGbbroJM2fOxKZNm/Dyyy/juOOOS/uMu+++G5dddhk2btyICy+8EFdddRX27t17WPeTEEIyhc6jhBDSywhCyCGZPn26UBRFZGVlpT3+/d//XQghBABxww03pL2ntLRU3HjjjUIIIZ588klRUFAgmpqawtdfeeUVwTkXNTU1QgghhgwZIn784x8fsA0AxLx588LnTU1NAoD4wx/+0G37SQghmULnUUII6XtojCMhnfD1r38dTzzxRNqywsLC8L/LysrSXisrK8OGDRsAAB9++CHGjRuHrKys8PUzzzwTnudhy5YtYIzhyy+/xHnnndduG04++eTwv7OyspCbm4udO3d2dpcIIeSwovMoIYT0LRQ4EtIJWVlZrbo8dZdoNNqh9TRNS3vOGIPneZloEiGEdDs6jxJCSN9CYxwJyYC333671fMxY8YAAMaMGYP3338f8Xg8fP2vf/0rOOc4/vjjkZOTg6OPPhpVVVWHtc2EENKb0HmUEEJ6F8o4EtIJpmmipqYmbZmqqigqKgIArFixAqeccgomTZqE5557Du+88w6eeuopAMBVV12FhQsXYvr06bjrrruwa9cuzJ49G1dffTWKi4sBAHfddRduuOEGDBo0CBdccAEaGxvx17/+FbNnzz68O0oIIRlC51FCCOlbKHAkpBNWrVqFwYMHpy07/vjj8dFHHwGQlfqef/55fP/738fgwYPxq1/9CieccAIAIBaL4dVXX8Utt9yCU089FbFYDNOmTcNDDz0Ubmv69OlIpVJ4+OGH8cMf/hBFRUX41re+dfh2kBBCMozOo4QQ0rcwIYTo6UYQ0p8wxvDSSy9hypQpPd0UQgjpk+g8SgghvQ+NcSSEEEIIIYQQ0i4KHAkhhBBCCCGEtIu6qhJCCCGEEEIIaRdlHAkhhBBCCCGEtIsCR0IIIYQQQggh7aLAkRBCCCGEEEJIuyhwJIQQQgghhBDSLgocCSGEEEIIIYS0iwJHQgghhBBCCCHtosCREEIIIYQQQki7KHAkhBBCCCGEENKu/w/U+eBp6qHOHQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "figure = this_sampler.results.plot_propermotion(\n", + " cbar_param=\"m1\", alpha=0.1, periods_to_plot=3\n", + ")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "python3.11", + "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.11.4" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/orbitize/example_data/gaia_edr3_betpic_epochs.csv b/orbitize/example_data/gaia_edr3_betpic_epochs.csv new file mode 100644 index 00000000..cd087df9 --- /dev/null +++ b/orbitize/example_data/gaia_edr3_betpic_epochs.csv @@ -0,0 +1,45 @@ +Target, ra[rad], dec[rad], ra[h:m:s], dec[d:m:s], ObservationTimeAtGaia[UTC], CcdRow[1-7], zetaFieldAngle[rad], scanAngle[rad], Fov[FovP=preceding/FovF=following], parallaxFactorAlongScan, parallaxFactorAcrossScan, ObservationTimeAtBarycentre[BarycentricJulianDateInTCB] +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2014-09-24T03:17:43.690,1,0.003366090442748648,-2.300433958482974,FoVP,-0.7177768121570808,0.716584607472628,2456924.6385198794 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2014-09-24T05:04:17.895,6,-0.0017592745186114803,-2.2989421430856525,FoVF,-0.7179415637921202,0.7163662372846105,2456924.712528852 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2014-10-01T18:42:06.590,7,-0.00650994606000498,2.541359236566953,FoVP,0.7116794191056144,0.7166374323350725,2456932.280652986 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2014-10-01T20:28:40.824,5,-0.001445323169629156,2.543117372208244,FoVF,0.7112554019851087,0.7169776440973757,2456932.3546622447 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2014-11-03T08:05:49.312,4,7.379734901303453E-4,-1.6391612276689718,FoVF,-0.6803572545190627,0.7088631454338691,2456964.839525243 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2014-12-07T01:43:30.212,7,-0.004849959376093157,-2.5453289312449967,FoVF,0.6523701313098679,0.7030559574218078,2456998.5744054615 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-01-05T23:43:47.075,2,0.0025747885979415836,-0.4942730752652883,FoVP,-0.6582104476275792,0.7000839976813767,2457028.491161975 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-01-06T01:30:21.281,6,-0.0020271807631072774,-0.49115917840781126,FoVF,-0.6594673585899907,0.6989435391249444,2457028.565168211 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-02-07T11:22:03.024,7,-0.006295401625795344,-1.4403312052347343,FoVP,0.6844806580297031,0.7039768380081998,2457060.975515688 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-02-07T13:08:37.264,5,-0.0012850130662692158,-1.438006139309293,FoVF,0.683776941405708,0.7047219521301625,2457061.0495214257 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-03-12T10:53:44.404,4,-0.0016724966096404376,0.6827836603917707,FoVP,-0.7079503473251454,0.7058933185346745,2457093.9550098255 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-06-07T09:46:48.222,2,0.0032231934804026075,0.6243368912655289,FoVP,0.6755480884910476,0.720672913448407,2457180.906924771 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-07-16T20:48:42.265,3,9.159206573612685E-4,2.8086879399064,FoVP,-0.6868433075451905,0.7281972252904217,2457220.366807382 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-07-16T22:35:16.471,7,-0.004261568960156839,2.8081985368817692,FoVF,-0.6856104978194572,0.7294175121952486,2457220.4408154115 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-08-10T02:42:59.981,1,0.0036191819871523073,1.6824763736679587,FoVP,0.7085412300520944,0.7247153497310905,2457244.6132780225 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-09-15T14:01:23.521,3,0.001089927452839157,-2.444198739257349,FoVP,-0.7196839438339331,0.7206183765561472,2457281.0852944436 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-09-15T15:47:57.737,7,-0.004085851094577186,-2.4430036981963634,FoVF,-0.7196627408057288,0.7206192439872735,2457281.1593035525 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-11-16T01:25:49.514,3,8.854859632097954E-4,-1.4264296688090004,FoVP,-0.6684715550056647,0.7086359784904921,2457342.5619674516 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-11-16T03:12:23.738,7,-0.003773896187501599,-1.4235221893706327,FoVF,-0.6695598930166724,0.7075315446893528,2457342.6359756514 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-12-20T07:03:04.009,3,9.588368448122845E-4,-2.303779872316881,FoVP,0.6488519167675069,0.7064083940105257,2457376.796352539 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2015-12-20T08:49:38.184,1,0.005480454883411049,-2.3005608291804354,FoVF,0.6475441389228678,0.7076150554246801,2457376.8703589914 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-01-19T06:48:28.296,1,0.00455958460866644,-0.2506064763631727,FoVP,-0.6711781552129196,0.6994346336977754,2457406.785930979 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-01-19T08:35:02.483,5,-1.6103865792494655E-4,-0.2476878009395681,FoVF,-0.6723150392499362,0.6983982869510355,2457406.8599365856 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-02-19T18:28:00.555,6,-0.004873371904800672,-1.2381780932616115,FoVP,0.6934653764733474,0.7070101491885915,2457438.271051817 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-02-19T20:14:34.781,4,2.6074977073553676E-4,-1.23621070452396,FoVF,0.6929968714204228,0.7075225391539418,2457438.3450571797 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-03-24T17:54:55.881,6,-0.0043337270022286226,0.8991683120726217,FoVP,-0.7074817657423257,0.7082054154742482,2457472.247171013 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-04-19T05:44:58.009,4,-0.0014381446312868695,-0.22886248328921097,FoVP,0.6963362021808351,0.712819874024422,2457497.739628123 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-04-19T07:31:32.219,2,0.003886499388137549,-0.2288081754393124,FoVF,0.6971805758323766,0.7119662135302238,2457497.813633586 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-05-27T16:50:32.073,5,-0.0029145249998517703,1.948752602913419,FoVP,-0.6680580667308067,0.725043819579978,2457536.201254772 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-06-18T22:47:24.786,7,-0.005962641445618194,0.8409130362544286,FoVP,0.6714549790241314,0.7239052814778538,2457558.4490316375 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-06-19T00:33:58.986,5,-8.288178662300212E-4,0.8400634884096743,FoVF,0.6729927004193484,0.722489770277802,2457558.5230387584 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-07-28T05:37:36.385,1,0.005979846420325473,3.021466380293874,FoVF,-0.6961663851608516,0.725506680122301,2457597.7343020677 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-07-28T09:51:16.403,7,-0.0064117397384628435,3.0209332244477936,FoVP,-0.6937088110107983,0.7279934284451889,2457597.9104628544 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-08-22T09:43:34.594,6,-0.005021369192567612,1.891477111963893,FoVP,0.7163779836170526,0.721693043409109,2457622.9056568476 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-09-26T22:53:20.282,2,0.004663702203171653,-2.2467618758998045,FoVF,-0.7127789642106954,0.7211462908866082,2457658.4550358946 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-10-28T02:48:28.158,3,6.303049850276892E-4,3.0113148092982573,FoVP,0.6912829527953595,0.7107404199561433,2457689.619063065 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-11-27T20:32:15.113,1,0.004541520064988245,-1.207663399140569,FoVP,-0.6588548698686335,0.7090410364005234,2457720.3582387106 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2016-11-27T22:18:49.298,5,-1.748119107680908E-5,-1.2045518070063896,FoVF,-0.6600772355459219,0.7078414863572481,2457720.432246054 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2017-01-01T02:09:11.754,7,-0.006902628486465427,-2.079797953973847,FoVP,0.6550633131402046,0.7064251343794189,2457754.5922265993 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2017-01-01T03:55:46.012,6,-0.002282603205303411,-2.07671056686794,FoVF,0.653834682789045,0.7075793064419984,2457754.6662335903 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2017-01-31T15:38:55.842,2,0.004129546721164321,-0.007922543428096887,FoVF,-0.6861814540849321,0.6996909360179947,2457785.154084105 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2017-04-06T20:18:02.634,2,0.003134500531205524,1.1100328332672587,FoVP,-0.7041159475800817,0.7099434017490356,2457850.346244692 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2017-04-06T22:04:36.843,6,-0.0022166543522921323,1.1103871962777563,FoVF,-0.7034680010111511,0.7105601458702862,2457850.420249916 +bet Pic,1.5153157780177544,-0.8912787619608181,05:47:17.088,-51:03:59.441,2017-05-01T09:56:47.276,6,-0.0030387516834771058,-0.018759200967237467,FoVF,0.686165682370106,0.7161488905344787,2457874.914289276 \ No newline at end of file diff --git a/orbitize/gaia.py b/orbitize/gaia.py index 39974632..ac313732 100644 --- a/orbitize/gaia.py +++ b/orbitize/gaia.py @@ -1,9 +1,19 @@ +import os import numpy as np import contextlib +import requests with contextlib.redirect_stdout(None): from astroquery.gaia import Gaia from astropy import units as u +import astropy.io.fits as fits +import astropy.time as time +from astropy.io.ascii import read +from astropy.coordinates import get_body_barycentric_posvel +import numpy.linalg + +from orbitize import DATADIR +import orbitize.lnlike class GaiaLogProb(object): """ @@ -173,3 +183,220 @@ def compute_lnlike( return lnlike + +class HGCALogProb(object): + """ + Class to compute the log probability of an orbit with respect to the proper + motion anomalies derived jointly from Gaia and Hipparcos using the HGCA catalogs + produced by Brandt (2018, 2021). With this class, both Gaia and Hipparcos + data will be considered. Do not need to pass in the Hipparcos class into System! + + Required auxiliary data: + * HGCA of acceleration (either DR2 or EDR3) + * Hipparcos IAD file (see orbitize.hipparcos for more info) + * Gaia Observation Forecast Tool (GOST) CSV output (https://gaia.esac.esa.int/gost/). + + For GOST, after entering the target name and resolving its coordinates, + use 2014-07-26T00:00:00 as the start date. For the end date, use + 2016-05-23T00:00:00 for DR2 and 2017-05-28T00:00:00 for EDR3. + + Args: + hip_id (int): the Hipparcos source ID of the object you're fitting. + hiplogprob (orbitize.hipparcos.HipLogProb): object containing + all info relevant to Hipparcos IAD fitting + gost_filepath (str): path to CSV file outputted by GOST + hgca_filepath (str): path to HGCA catalog FITS file. + If None, will download and store in orbitize.DATADIR + + Written: Jason Wang, 2022 + """ + def __init__(self, hip_id, hiplogprob, gost_filepath, hgca_filepath=None): + + # use default HGCA catalog if not supplied + if hgca_filepath is None: + # check orbitize.DATAIDR and download if needed + hgca_filepath = os.path.join(DATADIR, "HGCA_vEDR3.fits") + if not os.path.exists(hgca_filepath): + hgca_url = 'https://cdsarc.cds.unistra.fr/ftp/J/ApJS/254/42/HGCA_vEDR3.fits' + print("No HGCA catalog found. Downloading HGCA vEDR3 from {0} and storing into {1}.".format(hgca_url, hgca_filepath)) + hgca_file = requests.get(hgca_url, verify=False) + with open(hgca_filepath, 'wb') as f: + f.write(hgca_file.content) + else: + print("Using HGCA catalog stored in {0}".format(hgca_filepath)) + + # grab the entry from the HGCA + with fits.open(hgca_filepath, ignore_missing_simple=True, ignore_missing_end=True) as hdulist: + hgtable = hdulist[1].data + entry = hgtable[np.where(hgtable['hip_id'] == hip_id)] + # check we matched on a single target. mainly check if we typed hip id number incorrectly + if len(entry) != 1: + raise ValueError("HIP {0} encountered {1} matches. Expected 1 match.".format(hip_id, len(entry))) + # self.hgca_entry = entry + self.hip_id = hip_id + + # grab the relevant PM and uncertainties from HGCA + self.hip_pm = np.array([entry['pmra_hip'][0], entry['pmdec_hip'][0]]) + self.hip_pm_err = np.array([entry['pmra_hip_error'][0], entry['pmdec_hip_error'][0]]) + hip_radec_cov = entry['pmra_pmdec_hip'][0] * entry['pmra_hip_error'][0] * entry['pmdec_hip_error'][0] + + self.hg_pm = np.array([entry['pmra_hg'][0], entry['pmdec_hg'][0]]) + self.hg_pm_err = np.array([entry['pmra_hg_error'][0], entry['pmdec_hg_error'][0]]) + hg_radec_cov = entry['pmra_pmdec_hg'][0] * entry['pmra_hg_error'][0] * entry['pmdec_hg_error'][0] + + self.gaia_pm = np.array([entry['pmra_gaia'][0], entry['pmdec_gaia'][0]]) + self.gaia_pm_err = np.array([entry['pmra_gaia_error'][0], entry['pmdec_gaia_error'][0]]) + gaia_radec_cov = entry['pmra_pmdec_gaia'][0] * entry['pmra_gaia_error'][0] * entry['pmdec_gaia_error'][0] + + + # compute the differential PMs by subtracting Hip and Gaia from HG. Also propogate errors + self.hip_hg_dpm = self.hip_pm - self.hg_pm + self.hip_hg_dpm_err = np.sqrt(self.hip_pm_err**2 + self.hg_pm_err**2) + self.hip_hg_dpm_radec_corr = (hip_radec_cov + hg_radec_cov)/(self.hip_hg_dpm_err[0] * self.hip_hg_dpm_err[1]) + self.gaia_hg_dpm = self.gaia_pm - self.hg_pm + self.gaia_hg_dpm_err = np.sqrt(self.gaia_pm_err**2 + self.hg_pm_err**2) + self.gaia_hg_dpm_radec_corr = (gaia_radec_cov + hg_radec_cov)/(self.gaia_hg_dpm_err[0] * self.gaia_hg_dpm_err[1]) + + # condense together into orbitize "observations" + self.dpm_data = np.array([self.hip_hg_dpm, self.gaia_hg_dpm]) + self.dpm_err = np.array([self.hip_hg_dpm_err, self.gaia_hg_dpm_err]) + self.dpm_corr = np.array([self.hip_hg_dpm_radec_corr, self.gaia_hg_dpm_radec_corr]) + + # grab reference epochs for Gaia from HGCA so we can forward model it + self.gaia_epoch_ra = entry['epoch_ra_gaia'][0] + self.gaia_epoch_dec = entry['epoch_dec_gaia'][0] + # read in the GOST file to get the estimated Gaia epochs and scan angles + gost_dat = read(gost_filepath, converters={'*':[int, float, bytes]}) + self.gaia_epoch = time.Time(gost_dat["ObservationTimeAtGaia[UTC]"]).decimalyear # in decimal year + gaia_scan_theta = np.array(gost_dat["scanAngle[rad]"]) + gaia_scan_phi = gaia_scan_theta + np.pi/2 + self.gaia_cos_phi = np.cos(gaia_scan_phi) + self.gaia_sin_phi = np.sin(gaia_scan_phi) + self.gost_filepath = gost_filepath # save for saving + + # save the same infor from Hipparcos. we also have the errors on the Hipparcos data + self.hipparcos_epoch = hiplogprob.epochs # in decimal year + self.hipparcos_cos_phi = hiplogprob.cos_phi + self.hipparcos_sin_phi = hiplogprob.sin_phi + self.hipparcos_epoch_ra = entry['epoch_ra_hip'][0] + self.hipparcos_epoch_dec = entry['epoch_dec_hip'][0] + self.hippaarcos_errs = hiplogprob.eps + self.hiplogprob = hiplogprob # save for saving + + + def _save(self, hf): + """ + Saves the current object to an hdf5 file + + Args: + hf (h5py._hl.files.File): a currently open hdf5 file in which + to save the object. + """ + # save hipparocs IAD using exiting code + self.hiplogprob._save(hf) + + # save Gaia GOST file. avoid unicode!! + gost_dat = read(self.gost_filepath, converters={'*':[int, float, bytes]}) + hf.create_dataset("Gaia_GOST", data=gost_dat) + + + def compute_lnlike( + self, raoff_model, deoff_model, samples, param_idx + ): + """ + Computes the log likelihood of an orbit model with respect to a single + Gaia astrometric point. This is added to the likelihoods calculated with + respect to other data types in ``sampler._logl()``. + + Args: + raoff_model (np.array of float): NxM primary RA + offsets from the barycenter incurred from orbital motion of + companions (i.e. not from parallactic motion), where N is the + number of epochs of combined from Hipparcos and Gaia and M is the + number of orbits being tested, and raoff_model[0,:] are position + predictions at the Hipparcos epoch, and raoff_model[1,:] are + position predictions at the Gaia epoch + deoff_model (np.array of float): NxM primary decl + offsets from the barycenter incurred from orbital motion of + companions (i.e. not from parallactic motion), where N is the + number of epochs of combined from Hipparcos and Gaia and M is the + number of orbits being tested, and deoff_model[0,:] are position + predictions at the Hipparcos epoch, and deoff_model[1,:] are + position predictions at the Gaia epoch + samples (np.array of float): R-dimensional array of fitting + parameters, where R is the number of parameters being fit. Must + be in the same order documented in ``System``. + param_idx: a dictionary matching fitting parameter labels to their + indices in an array of fitting parameters (generally + set to System.basis.param_idx). + + Returns: + np.array of float: array of length M, where M is the number of input + orbits, representing the log likelihood of each orbit with + respect to the Gaia position measurement. + """ + # find the split between hipparcos and gaia data + gaia_index = len(self.hipparcos_epoch) + + # Begin forward modeling the data of the HGCA: linear motion during the Hip + # and Gaia missions, and the linear motion of the star between the two missions + # for now, think about only 1 model at a time to not break our brains + model_ra_hip = raoff_model[:gaia_index, 0] + model_dec_hip = deoff_model[:gaia_index, 0] + model_ra_gaia = raoff_model[gaia_index:, 0] + model_dec_gaia = deoff_model[gaia_index:, 0] + + hip_fit = self._linear_pm_fit(self.hipparcos_epoch, model_ra_hip, model_dec_hip, + self.hipparcos_epoch_ra, self.hipparcos_epoch_dec, + self.hipparcos_sin_phi, self.hipparcos_cos_phi, self.hippaarcos_errs) + model_hip_pos_offset = hip_fit[0:2] + model_hip_pm = hip_fit[2:4] + + + gaia_fit = self._linear_pm_fit(self.gaia_epoch, model_ra_gaia, model_dec_gaia, + self.gaia_epoch_ra, self.gaia_epoch_dec, + self.gaia_sin_phi, self.gaia_cos_phi, 1) + model_gaia_pos_offset = gaia_fit[0:2] + model_gaia_pm = gaia_fit[2:4] + + # compute the change in position between hipparcos and gaia + hg_dalpha_st = model_gaia_pos_offset[0] - model_hip_pos_offset[0] + model_hg_pmra = hg_dalpha_st/(self.gaia_epoch_ra - self.hipparcos_epoch_ra) + + hg_ddelta = model_gaia_pos_offset[1] - model_hip_pos_offset[1] + model_hg_pmdec = hg_ddelta/(self.gaia_epoch_dec - self.hipparcos_epoch_dec) + model_hg_pm = np.array([model_hg_pmra, model_hg_pmdec]) + + # take the difference between the linear motion measured during a mission, and the + # linear motion of the star between missions. These are our observables we compare + # to the data. Each variable contains both RA and Dec. + model_hip_hg_dpm = model_hip_pm - model_hg_pm + model_gaia_hg_dpm = model_gaia_pm - model_hg_pm + + # chi2 = (model_hip_hg_dpm - self.hip_hg_dpm)**2/(self.hip_hg_dpm_err)**2 + # chi2 += (model_gaia_hg_dpm - self.gaia_hg_dpm)**2/(self.gaia_hg_dpm_err)**2 + + lnlike = orbitize.lnlike.chi2_lnlike(self.dpm_data, self.dpm_err, self.dpm_corr, np.array([model_hip_hg_dpm, model_gaia_hg_dpm]), np.zeros(self.dpm_data.shape), []) + + return np.sum(lnlike) + + def _linear_pm_fit(self, epochs, raoff_planet, decoff_planet, epoch_ref_ra, epoch_ref_dec, sin_phi, cos_phi, errs): + """ + Performs a linear leastsq fit to determine the inferred proper motion given the stellar orbit around the barycenter + """ + # Sovle y = A * x + # construct A matrix + A_pmra = cos_phi * (epochs - epoch_ref_ra) / errs + A_raoff = cos_phi / errs + A_pmdec = sin_phi * (epochs - epoch_ref_dec) / errs + A_decoff = sin_phi / errs + A_matrix = np.vstack((A_raoff, A_decoff, A_pmra, A_pmdec)).T + + # construct y matrix + y_vec = (raoff_planet * cos_phi + decoff_planet * sin_phi)/errs + + x, res, _, _ = numpy.linalg.lstsq(A_matrix, y_vec, rcond=None) + + return x + + \ No newline at end of file diff --git a/orbitize/plot.py b/orbitize/plot.py index c3936092..1825cb56 100644 --- a/orbitize/plot.py +++ b/orbitize/plot.py @@ -11,6 +11,7 @@ import matplotlib.pyplot as plt from matplotlib.collections import LineCollection import matplotlib.colors as colors +from matplotlib.ticker import FormatStrFormatter from erfa import ErfaWarning @@ -333,7 +334,7 @@ def plot_orbits(results, object_to_plot=1, start_mjd=51544., fig = plt.figure(figsize=(14, 6)) ax = plt.subplot2grid((2, 14), (0, 0), rowspan=2, colspan=6) else: - plt.set_current_figure(fig) + plt.figure(fig.number) if rv_time_series: ax = plt.subplot2grid((3, 14), (0, 0), rowspan=2, colspan=6) else: @@ -683,4 +684,280 @@ def plot_orbits(results, object_to_plot=1, start_mjd=51544., ax2.locator_params(axis='x', nbins=6) ax2.locator_params(axis='y', nbins=6) - return fig \ No newline at end of file + return fig + +def plot_propermotion(results, system, object_to_plot=1, start_mjd=44239., + periods_to_plot=1, end_year=2030.0, alpha = 0.05, + num_orbits_to_plot=100, num_epochs_to_plot=100, + show_colorbar=True, cmap=cmap, + cbar_param=None, tight_layout=False, + # fig=None + ): + """ + Plots the proper motion of a host star as induced by a companion for + one orbital period for a select number of fitted orbits + for a given object, with line segments colored according to a given + parameter (most informative is usually mass of companion) + + Important Note: These plotted trajectories aren't what are fitting in the + likelihood evaluation for the HGCA runs. The implementation forward models + the Hip/Gaia measurements per epoch and infers the differential proper motions. + This plot is given only for the purposes of an approximate visualization. + + Args: + system (object): orbitize.system object with a HGCALogProb passed to system.gaia + object_to_plot (int): which object to plot (default: 1) + start_mjd (float): MJD in which to start plotting orbits (default: 51544, + the year 2000) + periods_to_plot (int): number of periods to plot (default: 1) + end_year (float): decimal year specifying when to stop plotting orbit + tracks in the Sep/PA panels (default: 2025.0). + alpha (float): transparency of lines (default: 0.05) + num_orbits_to_plot (int): number of orbits to plot (default: 100) + num_epochs_to_plot (int): number of points to plot per orbit (default: 100) + show_colorbar (Boolean): Displays colorbar to the right of the plot [True] + cmap (matplotlib.cm.ColorMap): color map to use for making orbit tracks + (default: modified Purples_r) + cbar_param (string): options are the following: 'sma1', 'ecc1', 'inc1', 'aop1', + 'pan1', 'tau1', 'plx', 'm0', 'm1', etc. Number can be switched out. Default is None. + tight_layout (bool): apply plt.tight_layout function? + fig (matplotlib.pyplot.Figure): optionally include a predefined Figure object to plot the orbit on. + Most users will not need this keyword. + + Return: + ``matplotlib.pyplot.Figure``: the orbit plot if input is valid, ``None`` otherwise + + + (written): William Balmer (2023), based on plot_orbits by Sarah Blunt and Henry Ngo + + """ + + if Time(start_mjd, format='mjd').decimalyear >= end_year: + raise ValueError('start_mjd keyword date must be less than end_year keyword date.') + + if object_to_plot > results.num_secondary_bodies: + raise ValueError("Only {0} secondary bodies being fit. Requested to plot body {1} which is out of range".format(results.num_secondary_bodies, object_to_plot)) + + if object_to_plot == 0: + raise ValueError("Plotting the primary's orbit is currently unsupported. Stay tuned.") + + with warnings.catch_warnings(): + warnings.simplefilter('ignore', ErfaWarning) + + data = results.data[results.data['object'] == object_to_plot] + possible_cbar_params = [ + 'sma', + 'ecc', + 'inc', + 'aop' + 'pan', + 'tau', + 'plx', + 'm0', + 'm1' + ] + + if cbar_param is None: + pass + elif cbar_param[0:3] in possible_cbar_params: + index = results.param_idx[cbar_param] + else: + raise Exception( + "Invalid input; acceptable inputs include 'Epoch [year]', 'plx', 'sma1', 'ecc1', 'inc1', 'aop1', 'pan1', 'tau1', 'sma2', 'ecc2', ...)" + ) + # Select random indices for plotted orbit + num_orbits = len(results.post[:, 0]) + if num_orbits_to_plot > num_orbits: + num_orbits_to_plot = num_orbits + choose = np.random.randint(0, high=num_orbits, size=num_orbits_to_plot) + + # Get posteriors from random indices + standard_post = [] + if results.sampler_name == 'MCMC': + # Convert the randomly chosen posteriors to standard keplerian set + for i in np.arange(num_orbits_to_plot): + orb_ind = choose[i] + param_set = np.copy(results.post[orb_ind]) + standard_post.append(results.basis.to_standard_basis(param_set)) + else: # For OFTI, posteriors are already converted + for i in np.arange(num_orbits_to_plot): + orb_ind = choose[i] + standard_post.append(results.post[orb_ind]) + + standard_post = np.array(standard_post) + + sma = standard_post[:, results.standard_param_idx['sma{}'.format(object_to_plot)]] + ecc = standard_post[:, results.standard_param_idx['ecc{}'.format(object_to_plot)]] + inc = standard_post[:, results.standard_param_idx['inc{}'.format(object_to_plot)]] + aop = standard_post[:, results.standard_param_idx['aop{}'.format(object_to_plot)]] + pan = standard_post[:, results.standard_param_idx['pan{}'.format(object_to_plot)]] + tau = standard_post[:, results.standard_param_idx['tau{}'.format(object_to_plot)]] + plx = standard_post[:, results.standard_param_idx['plx']] + + # Then, get the other parameters + if 'mtot' in results.labels: + mtot = standard_post[:, results.standard_param_idx['mtot']] + elif 'm0' in results.labels: + m0 = standard_post[:, results.standard_param_idx['m0']] + m1 = standard_post[:, results.standard_param_idx['m{}'.format(object_to_plot)]] + mtot = m0 + m1 + + raoff = np.zeros((num_orbits_to_plot, num_epochs_to_plot)) + deoff = np.zeros((num_orbits_to_plot, num_epochs_to_plot)) + vz_star = np.zeros((num_orbits_to_plot, num_epochs_to_plot)) + epochs = np.zeros((num_orbits_to_plot, num_epochs_to_plot)) + + # Loop through each orbit to plot and calcualte ra/dec offsets for all points in orbit + # Need this loops since epochs[] vary for each orbit, unless we want to just plot the same time period for all orbits + for i in np.arange(num_orbits_to_plot): + # Compute period (from Kepler's third law) + period = np.sqrt(4*np.pi**2.0*(sma*u.AU)**3/(consts.G*(mtot*u.Msun))) + period = period.to(u.day).value + + # Create an epochs array to plot num_epochs_to_plot points over one orbital period + epochs[i, :] = np.linspace(start_mjd, float( + start_mjd+(period[i]*periods_to_plot)), num_epochs_to_plot) + + # Calculate ra/dec offsets for all epochs of this orbit + raoff0, deoff0, _ = kepler.calc_orbit( + epochs[i, :], sma[i], ecc[i], inc[i], aop[i], pan[i], + tau[i], plx[i], mtot[i], tau_ref_epoch=results.tau_ref_epoch + ) + + raoff[i, :] = raoff0 + deoff[i, :] = deoff0 + + # Create a linearly increasing colormap for our range of epochs + if cbar_param is not None: + cbar_param_arr = results.post[:, index] + norm = mpl.colors.Normalize(vmin=np.min(cbar_param_arr), + vmax=np.max(cbar_param_arr)) + + elif cbar_param is None: + + norm = mpl.colors.Normalize() + + # Create figure for orbit plots + fig, axs = plt.subplots(1, 2, figsize=(8,4), facecolor='white') + + # Plot each orbit (each segment between two points coloured using colormap) + for i in np.arange(num_orbits_to_plot): + epoch_in_yr = Time(epochs[i,:], format='mjd').decimalyear + # masses (in same units, solar) + m_b = standard_post[:, results.param_idx['m1']][i] + m_a = standard_post[:, results.param_idx['m0']][i] + # dt + timestep = epoch_in_yr[1]-epoch_in_yr[0] + # dra/dt and ddec/dt + ddec_b = np.gradient(deoff[i,:], timestep) # in mas/yr + dec_b_radian = deoff[i,:]*(2.7777778e-7)*(0.017453293) # mas -> deg -> radian + ra_b = raoff[i,:] + rastar_b = ra_b*np.cos(dec_b_radian) # in mas + drastar_b = np.gradient(rastar_b, timestep) # in mas/yr + + # convert to dRA^star_star (lol) and dDec_star + mass_ratio_ = (-1*m_b/(m_a+m_b)) + ddec_a = ddec_b * mass_ratio_ + drastar_a = drastar_b * mass_ratio_ + + if cbar_param is not None: + color = cmap(norm(standard_post[:, results.param_idx[cbar_param]][i])) + else: + color = 'k' + + axs[0].plot(epoch_in_yr,drastar_a+system.gaia.hg_pm[0], + color=color, + alpha=alpha, + zorder=0 + ) + axs[1].plot(epoch_in_yr,ddec_a+system.gaia.hg_pm[1], + color=color, + alpha=alpha, + zorder=0 + ) + + + + axs[0].set_xlim(1980,2030) + axs[0].yaxis.set_major_formatter(FormatStrFormatter('%.1f')) + axs[1].set_xlabel('Epoch') + + axs[0].set_ylabel(r'$\mu_\alpha^*$ [mas/yr]') + + axs[0].errorbar(np.nanmedian(system.gaia.hipparcos_epoch), + system.gaia.hip_pm[0], + yerr=system.gaia.hip_pm_err[0], + zorder=30, + mec='k', + fmt='s', color='cornflowerblue') + + hgca_epoch = (system.gaia.gaia_epoch_ra+np.nanmedian(system.gaia.hipparcos_epoch))/2 + hgca_epoch_err = (system.gaia.gaia_epoch_ra-np.nanmedian(system.gaia.hipparcos_epoch))/2 + + axs[0].errorbar(hgca_epoch, + system.gaia.hg_pm[0], + xerr=hgca_epoch_err, + yerr=system.gaia.hg_pm_err[0], + zorder=30, + mec='k', + fmt='^', color='#6280D6') + + + axs[0].errorbar(system.gaia.gaia_epoch_ra, + system.gaia.gaia_pm[0], + yerr=system.gaia.gaia_pm_err[0], + zorder=30, + mec='k', + fmt='o', color='#5f61b4') + + + axs[1].set_xlim(1980,2030) + axs[1].yaxis.set_major_formatter(FormatStrFormatter('%.1f')) + + axs[1].errorbar(np.nanmedian(system.gaia.hipparcos_epoch), + system.gaia.hip_pm[1], + yerr=system.gaia.hip_pm_err[1], + zorder=30, + mec='k', + fmt='s', color='cornflowerblue', label='Hip.') + + axs[1].errorbar(hgca_epoch, + system.gaia.hg_pm[1], + xerr=hgca_epoch_err, + yerr=system.gaia.hg_pm_err[1], + zorder=30, + mec='k', + fmt='^', color='#6280D6', label='H-G') + + axs[1].errorbar(system.gaia.gaia_epoch_ra, + system.gaia.gaia_pm[1], + yerr=system.gaia.gaia_pm_err[1], + zorder=30, + mec='k', + fmt='o', color='#5f61b4', label='Gaia') + + axs[1].set_ylabel(r'$\mu_\delta$ [mas/yr]') + axs[1].set_xlabel('Epoch') + axs[0].set_xlabel('Epoch') + + cbar_ax = fig.add_axes([1.03, 0.15, 0.03, 0.80]) + + cbar = mpl.colorbar.ColorbarBase( + cbar_ax, cmap=cmap, norm=norm, orientation='vertical', + label=cbar_param + ) + + axs[0].set_rasterization_zorder(1) + axs[1].set_rasterization_zorder(1) + + axs[1].legend() + + print("Important Note of Caution: the orbitize! implementation of the HGCA \n", + "fits for the time-averaged proper motions, and not the instantaneous proper \n", + "motions that are being plotted here. This plot is provided only for the \n", + "purpose of an approximate check on the fit.") + + if tight_layout: + plt.tight_layout() + + return fig diff --git a/orbitize/results.py b/orbitize/results.py index 9349109c..dee48b8d 100644 --- a/orbitize/results.py +++ b/orbitize/results.py @@ -215,16 +215,32 @@ def load_results(self, filename, append=False): ) os.system('rm {}'.format(tmpfile)) + + # load Gaia data try: gaia_num = int(hf.attrs['gaia_num']) dr = str(hf.attrs['dr']) gaia = orbitize.gaia.GaiaLogProb(gaia_num, hipparcos_IAD, dr) except KeyError: gaia = None + + # alternatively load HGCA. Note requires hipparcos_IAD attribute + gaiagost_data = hf.get("Gaia_GOST") + if gaiagost_data is not None: + + tmpfile = 'thisisprettyhackysorrylmao' + tmptbl = table.Table(np.array(gaiagost_data)) + tmptbl.write(tmpfile, format="ascii", overwrite=True) + + gaia = orbitize.gaia.HGCALogProb(int(hip_num), hipparcos_IAD, tmpfile) + hipparcos_IAD = None # HGCA handles hipparocs, so don't want to pass Hipparcos also into the system + + os.system('rm {}'.format(tmpfile)) else: hipparcos_IAD = None gaia = None + try: fitting_basis = str(hf.attrs['fitting_basis']) except KeyError: @@ -338,5 +354,27 @@ def plot_orbits(self, object_to_plot=1, start_mjd=51544., plot_errorbars=plot_errorbars, fig=fig ) + def plot_propermotion(self, + # system, + object_to_plot=1, start_mjd=44239., + periods_to_plot=1, end_year=2030.0, alpha = 0.05, + num_orbits_to_plot=100, num_epochs_to_plot=100, + show_colorbar=True, + cmap=orbitize.plot.cmap, + cbar_param=None, + # fig=None + ): + """ + Wrapper for orbitize.plot.plot_propermotion + """ - + return orbitize.plot.plot_propermotion(self, self.system, + object_to_plot=object_to_plot, start_mjd=start_mjd, + periods_to_plot=periods_to_plot, end_year=end_year, + alpha = alpha, num_orbits_to_plot=num_orbits_to_plot, + num_epochs_to_plot=num_epochs_to_plot, + show_colorbar=show_colorbar, + cmap=cmap, + cbar_param=cbar_param, + # fig=fig + ) diff --git a/orbitize/sampler.py b/orbitize/sampler.py index 5251cb09..764abe50 100644 --- a/orbitize/sampler.py +++ b/orbitize/sampler.py @@ -122,24 +122,26 @@ def _logl(self, params): self.system.param_idx, ) - if self.system.gaia is not None: - gaiahip_epochs = Time( - [self.system.gaia.hipparcos_epoch, self.system.gaia.gaia_epoch], - format="decimalyear", - ).mjd - - # compute Ra/Dec predictions at the Gaia epoch - raoff_model, deoff_model, _ = self.system.compute_all_orbits( - params, epochs=gaiahip_epochs - ) + if self.system.gaia is not None: + gaiahip_epochs = Time( + np.append( + self.system.gaia.hipparcos_epoch, self.system.gaia.gaia_epoch + ), + format="decimalyear", + ).mjd + + # compute Ra/Dec predictions at the Gaia epoch + raoff_model, deoff_model, _ = self.system.compute_all_orbits( + params, epochs=gaiahip_epochs + ) - # select body 0 raoff/deoff predictions & feed into Gaia module lnlike fn - lnlikes_sum += self.system.gaia.compute_lnlike( - raoff_model[:, 0, :], - deoff_model[:, 0, :], - params, - self.system.param_idx, - ) + # select body 0 raoff/deoff predictions & feed into Gaia module lnlike fn + lnlikes_sum += self.system.gaia.compute_lnlike( + raoff_model[:, 0, :], + deoff_model[:, 0, :], + params, + self.system.param_idx, + ) return lnlikes_sum diff --git a/orbitize/system.py b/orbitize/system.py index 7c299b87..4f1dddd9 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -135,8 +135,10 @@ def __init__(self, num_secondary_bodies, data_table, stellar_or_system_mass, self.track_planet_perturbs = ( self.fit_secondary_mass and ( - (len(self.radec[0]) + len(self.seppa[0] > 0) or - (self.num_secondary_bodies > 1) + ((len(self.radec[0]) + len(self.seppa[0]) > 0) or + (self.num_secondary_bodies > 1) or + (hipparcos_IAD is not None) or + (gaia is not None) ) ) ) diff --git a/tests/test_gaia.py b/tests/test_gaia.py index 12822731..70dbce85 100644 --- a/tests/test_gaia.py +++ b/tests/test_gaia.py @@ -1,81 +1,92 @@ import numpy as np import os from orbitize import DATADIR -from orbitize import hipparcos, gaia, basis, system, read_input +from orbitize import hipparcos, gaia, basis, system, read_input, sampler, results + def test_dr2_edr3(): """ Test that both DR2 and eDR3 retrieval gives ballpark similar values for beta Pic """ - hip_num = '027321' # beta Pic + hip_num = "027321" # beta Pic edr3_num = 4792774797545800832 dr2_number = 4792774797545105664 num_secondary_bodies = 1 - path_to_iad_file = '{}HIP{}.d'.format(DATADIR, hip_num) + path_to_iad_file = "{}HIP{}.d".format(DATADIR, hip_num) myHip = hipparcos.HipparcosLogProb(path_to_iad_file, hip_num, num_secondary_bodies) - dr3Gaia = gaia.GaiaLogProb( - edr3_num, myHip, dr='edr3' - ) - dr2Gaia = gaia.GaiaLogProb( - dr2_number, myHip, dr='dr2' - ) + dr3Gaia = gaia.GaiaLogProb(edr3_num, myHip, dr="edr3") + dr2Gaia = gaia.GaiaLogProb(dr2_number, myHip, dr="dr2") - assert np.isclose(dr2Gaia.ra, dr3Gaia.ra, atol=0.1) # abs tolerance in degrees + assert np.isclose(dr2Gaia.ra, dr3Gaia.ra, atol=0.1) # abs tolerance in degrees def test_system_setup(): """ Test that a System object with Hipparcos and Gaia is initialized correctly """ - hip_num = '027321' # beta Pic + hip_num = "027321" # beta Pic edr3_num = 4792774797545800832 num_secondary_bodies = 1 - path_to_iad_file = '{}HIP{}.d'.format(DATADIR, hip_num) + path_to_iad_file = "{}HIP{}.d".format(DATADIR, hip_num) myHip = hipparcos.HipparcosLogProb(path_to_iad_file, hip_num, num_secondary_bodies) - myGaia = gaia.GaiaLogProb( - edr3_num, myHip, dr='edr3' - ) + myGaia = gaia.GaiaLogProb(edr3_num, myHip, dr="edr3") - input_file = os.path.join(DATADIR, 'betaPic.csv') + input_file = os.path.join(DATADIR, "betaPic.csv") plx = 51.5 num_secondary_bodies = 1 data_table = read_input.read_file(input_file) betaPic_system = system.System( - num_secondary_bodies, data_table, 1.75, plx, hipparcos_IAD=myHip, - gaia=myGaia, fit_secondary_mass=True, mass_err=0.01, - plx_err=0.01 + num_secondary_bodies, + data_table, + 1.75, + plx, + hipparcos_IAD=myHip, + gaia=myGaia, + fit_secondary_mass=True, + mass_err=0.01, + plx_err=0.01, ) assert betaPic_system.labels == [ - 'sma1', 'ecc1', 'inc1', 'aop1', 'pan1', 'tau1', 'plx', 'pm_ra', 'pm_dec', - 'alpha0', 'delta0', 'm1', 'm0' - ] + "sma1", + "ecc1", + "inc1", + "aop1", + "pan1", + "tau1", + "plx", + "pm_ra", + "pm_dec", + "alpha0", + "delta0", + "m1", + "m0", + ] assert betaPic_system.fit_secondary_mass assert betaPic_system.track_planet_perturbs + def test_valueerror(): """ Check that if I don't say dr2 or edr3, I get a value error """ - hip_num = '027321' # beta Pic + hip_num = "027321" # beta Pic edr3_num = 4792774797545800832 num_secondary_bodies = 1 - path_to_iad_file = '{}HIP{}.d'.format(DATADIR, hip_num) + path_to_iad_file = "{}HIP{}.d".format(DATADIR, hip_num) myHip = hipparcos.HipparcosLogProb(path_to_iad_file, hip_num, num_secondary_bodies) try: - myGaia = gaia.GaiaLogProb( - edr3_num, myHip, dr='dr3' - ) - assert False, 'Test failed!' + myGaia = gaia.GaiaLogProb(edr3_num, myHip, dr="dr3") + assert False, "Test failed!" except ValueError: pass @@ -104,70 +115,77 @@ def test_orbit_calculation(): pm_a = 0 pm_d = 0 - plx = 100 # [mas] + plx = 100 # [mas] m0 = 1 m1 = 1 a0 = 0 d0 = 0 - hip_num = '027321' # beta Pic + hip_num = "027321" # beta Pic edr3_num = 4792774797545800832 num_secondary_bodies = 1 - path_to_iad_file = '{}HIP{}.d'.format(DATADIR, hip_num) + path_to_iad_file = "{}HIP{}.d".format(DATADIR, hip_num) myHip = hipparcos.HipparcosLogProb(path_to_iad_file, hip_num, num_secondary_bodies) - myGaia = gaia.GaiaLogProb( - edr3_num, myHip, dr='edr3' - ) + myGaia = gaia.GaiaLogProb(edr3_num, myHip, dr="edr3") param_idx = { - 'sma1':0, 'ecc1':1, 'inc1':2, 'aop1':3,'pan1':4, 'tau1':5, 'plx':6, - 'm0':7, 'm1':8, 'alpha0':9, 'delta0':10, 'pm_ra':11, 'pm_dec':12 + "sma1": 0, + "ecc1": 1, + "inc1": 2, + "aop1": 3, + "pan1": 4, + "tau1": 5, + "plx": 6, + "m0": 7, + "m1": 8, + "alpha0": 9, + "delta0": 10, + "pm_ra": 11, + "pm_dec": 12, } - # Case 1: only proper motion explains Gaia-Hip offset - pm_a = 100 # [mas/yr] - pm_d = 100 # [mas/yr] - sma = 1e-17 + pm_a = 100 # [mas/yr] + pm_d = 100 # [mas/yr] + sma = 1e-17 - raoff = np.zeros((2, 1)) + raoff = np.zeros((2, 1)) deoff = np.zeros((2, 1)) myGaia.ra = myHip.alpha0 + ( - myGaia.mas2deg * pm_a * (myGaia.gaia_epoch - myGaia.hipparcos_epoch) / - np.cos(np.radians(myHip.delta0)) + myGaia.mas2deg + * pm_a + * (myGaia.gaia_epoch - myGaia.hipparcos_epoch) + / np.cos(np.radians(myHip.delta0)) ) myGaia.dec = myHip.delta0 + ( myGaia.mas2deg * pm_d * (myGaia.gaia_epoch - myGaia.hipparcos_epoch) ) test_samples = [sma, ecc, inc, aop, pan, tau, plx, m0, m1, a0, d0, pm_a, pm_d] - lnlike = myGaia.compute_lnlike( - raoff, deoff, test_samples, param_idx - ) + lnlike = myGaia.compute_lnlike(raoff, deoff, test_samples, param_idx) assert np.isclose(np.exp(lnlike), 1) # Case 2: only H0 offset explains Gaia-Hip offset - test_samples[param_idx['pm_dec']] = 0 - test_samples[param_idx['pm_ra']] = 0 - a0 = 100; d0 = 100 - test_samples[param_idx['alpha0']] = a0 # [mas] - test_samples[param_idx['delta0']] = d0 # [mas] - - myGaia.ra = myHip.alpha0 + myGaia.mas2deg * a0 / np.cos(np.radians(myHip.delta0)) + test_samples[param_idx["pm_dec"]] = 0 + test_samples[param_idx["pm_ra"]] = 0 + a0 = 100 + d0 = 100 + test_samples[param_idx["alpha0"]] = a0 # [mas] + test_samples[param_idx["delta0"]] = d0 # [mas] + + myGaia.ra = myHip.alpha0 + myGaia.mas2deg * a0 / np.cos(np.radians(myHip.delta0)) myGaia.dec = myHip.delta0 + myGaia.mas2deg * d0 - lnlike = myGaia.compute_lnlike( - raoff, deoff, test_samples, param_idx - ) + lnlike = myGaia.compute_lnlike(raoff, deoff, test_samples, param_idx) assert np.isclose(np.exp(lnlike), 1) # Case 3: only orbital motion explains Gaia-Hip offset - test_samples[param_idx['alpha0']] = 0 - test_samples[param_idx['delta0']] = 0 + test_samples[param_idx["alpha0"]] = 0 + test_samples[param_idx["delta0"]] = 0 mas2arcsec = 1e-3 deg2arcsec = 3600 @@ -175,49 +193,123 @@ def test_orbit_calculation(): myGaia.ra = myHip.alpha0 myGaia.dec = myHip.delta0 + 1 - sma = 2 * (myGaia.dec - myHip.delta0) * deg2arcsec * (plx * mas2arcsec) # [au] - per = 2 * (myGaia.gaia_epoch - myGaia.hipparcos_epoch) # [yr] + sma = 2 * (myGaia.dec - myHip.delta0) * deg2arcsec * (plx * mas2arcsec) # [au] + per = 2 * (myGaia.gaia_epoch - myGaia.hipparcos_epoch) # [yr] mtot = sma**3 / per**2 - test_samples[param_idx['sma1']] = sma - test_samples[param_idx['m0']] = mtot/2 - test_samples[param_idx['m1']] = mtot/2 + test_samples[param_idx["sma1"]] = sma + test_samples[param_idx["m0"]] = mtot / 2 + test_samples[param_idx["m1"]] = mtot / 2 # passes through peri (+sma decl for e=0 orbits) at Hipparcos epoch # -> @ Gaia epoch, primary should be at +sma decl tau = basis.tp_to_tau(myGaia.hipparcos_epoch, 58849, per) - test_samples[param_idx['tau1']] = tau + test_samples[param_idx["tau1"]] = tau # choose sma and mass so that Hipparcos/Gaia difference is only due to orbit. - deoff[1,:] = (myGaia.dec - myHip.delta0) / myGaia.mas2deg - deoff[0,:] = 0 - lnlike = myGaia.compute_lnlike( - raoff, deoff, test_samples, param_idx - ) + deoff[1, :] = (myGaia.dec - myHip.delta0) / myGaia.mas2deg + deoff[0, :] = 0 + lnlike = myGaia.compute_lnlike(raoff, deoff, test_samples, param_idx) assert np.isclose(np.exp(lnlike), 1) + +def test_hgca(): + """ + Tests that the HGCA module works for beta Pic b. + + Checks can download HGCA catalog if not available, and checks GRAVITY Collaboration et al. 2020 + orbital parameters against the data + """ + iad_filepath = os.path.join(DATADIR, "HIP027321.d") + gost_filepath = os.path.join(DATADIR, "gaia_edr3_betpic_epochs.csv") + astrometry_filepath = os.path.join(DATADIR, "betaPic.csv") + + hipparcos_lnprob = hipparcos.HipparcosLogProb(iad_filepath, 27321, 1) + hgca_lnprob = gaia.HGCALogProb(27321, hipparcos_lnprob, gost_filepath) + + # test a few things were read in correctly + assert np.all(np.isfinite(hgca_lnprob.hip_pm)) + assert len(hgca_lnprob.hipparcos_epoch) > 0 + assert len(hgca_lnprob.gaia_epoch) > 0 + + # Initialize System object which stores data & sets priors + data_table = read_input.read_file(astrometry_filepath) + this_system = system.System( + 1, + data_table, + 1.75, + 51.44, + mass_err=0.05, + plx_err=0.12, + tau_ref_epoch=55000, + fit_secondary_mass=True, + gaia=hgca_lnprob, + ) + + # create sampler just for lnlike function + this_sampler = sampler.MCMC(this_system, 1, 100, 1) + + # params from GRAVITY Collaboration et al. 2020 fit to HGCA DR2 + params = np.array( + [ + [ + 1.04e01, + 1.44e-01, + 1.5534e00, + 3.5325e00, + 5.5670e-01, + 1.80968e-01, + 5.1456e01, + 1.367e-02, + 1.783, + ] + ] + ) + + chi2 = this_sampler._logl(params) + + assert np.isfinite(chi2) + + # briefly run sampler and save result + this_sampler.run_sampler(100, 0) + this_sampler.results.save_results("hgca_test.hdf5") + + # check that the data stays the same + res = results.Results() + res.load_results("hgca_test.hdf5") + new_hgca_lnprob = res.system.gaia + assert isinstance(new_hgca_lnprob, gaia.HGCALogProb) + assert new_hgca_lnprob.hip_hg_dpm_radec_corr == hgca_lnprob.hip_hg_dpm_radec_corr + + os.remove("hgca_test.hdf5") + + def test_nointernet(): """ Test that the internet-less object setup works """ - hip_num = '027321' # beta Pic + hip_num = "027321" # beta Pic dr2_number = 4792774797545105664 num_secondary_bodies = 1 - path_to_iad_file = '{}HIP{}.d'.format(DATADIR, hip_num) + path_to_iad_file = "{}HIP{}.d".format(DATADIR, hip_num) myHip = hipparcos.HipparcosLogProb(path_to_iad_file, hip_num, num_secondary_bodies) - dr3Gaia = gaia.GaiaLogProb( - dr2_number, myHip, dr='dr2', query=False, gaia_data = {'ra':0, 'dec':0, 'ra_error':0, 'dec_error':0} + _ = gaia.GaiaLogProb( + dr2_number, + myHip, + dr="dr2", + query=False, + gaia_data={"ra": 0, "dec": 0, "ra_error": 0, "dec_error": 0}, ) - -if __name__ == '__main__': +if __name__ == "__main__": test_nointernet() # test_dr2_edr3() # test_system_setup() # test_valueerror() - # test_orbit_calculation() \ No newline at end of file + # test_orbit_calculation() + test_hgca()