{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Asteroid Detection Module\n", "\n", "**Lecturer:** Ashish Mahabal
\n", "**Jupyter Notebook Authors:** Quanzhi Ye, Ashish Mahabel, Dmitry Duev, & Cameron Hummels\n", "\n", "This is a Jupyter notebook lesson taken from the GROWTH Winter School 2018. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-astro-school-2018-resources.html\n", "\n", "## Objective\n", "Compute the detection rate of fast-moving asteroids, asteroids that pass within ~20 lunar distances from Earth, and leave trails/streaks on astronomical images. \n", "\n", "## Key steps\n", "- Figure out which known asteroids will show up as streaks in ZTF images, and\n", "- Examine the streak catalog to see how many of them are detected.\n", "\n", "## Required dependencies\n", "\n", "See GROWTH school webpage for detailed instructions on how to install these modules and packages. Nominally, you should be able to install the python modules with `pip install `. The external astromatic packages are easiest installed using package managers (e.g., `rpm`, `apt-get`).\n", "\n", "### Python modules\n", "* python 3\n", "* astropy\n", "* numpy\n", "* astroquery\n", "* matplotlib\n", "\n", "### External packages\n", "None" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import requests, numpy as np, glob, matplotlib.pyplot as plt, os, tarfile\n", "from astroquery.jplhorizons import Horizons\n", "from astropy.time import Time\n", "from IPython.display import Image, display\n", "import pdb\n", "import warnings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 1: fast-moving asteroids that fall on ZTF images\n", "\n", "We will loop over the catalog of known NEAs (Near-Earth Asteroids) and make a series of queries to JPL's Horizon service and IPAC's IRSA service, to determine if and when they will be visible at Palomar (where ZTF runs from)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's download the catalog of known NEAs from the Minor Planet Center (MPC)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "retrieving MPCORB...\n", "MPCORB file retrieved.\n" ] } ], "source": [ "os.chdir('data') # All relevant data are in data directory\n", "print('retrieving MPCORB...')\n", "# mpcorb = 'https://www.minorplanetcenter.net/iau/MPCORB/NEA.txt'\n", "# r = requests.get(mpcorb)\n", "# mpcorb = r.text.split(\"\\n\")\n", "mpcorb = open('NEA.txt','r').read().split(\"\\n\")\n", "print('MPCORB file retrieved.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let us loop over the MPC catalog to check the visibility for every NEA. We will do it for NEA 2018 CL (the first NEA found by ZTF), and it's up to you to write a loop for it.\n", "\n", "The format of the MPC catalog is described [here](https://www.minorplanetcenter.net/iau/info/MPOrbitFormat.html), it's best for machines but is somewhat human readable. First, we create an instance for the Horizons query. We only need to tell Horizons the name of the object we want to query, and the time window we want to query, and let Horizons do the heavy-lifting. Oh, and the three-letter-code for ZTF is I41, as described [here](https://www.minorplanetcenter.net/iau/lists/ObsCodesF.html). If the query returns zero entries, it means that the object is not observable at Palomar at the suggested time (interval)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, let us define a time window we would like to query for. Our default is 20180205, the day ZTF science observations started:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "start_date = '20180205'\n", "end_date = '20180206'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then let's do all the tasks described above. For now, we do our test with 2018 CL (the first NEA found by ZTF).\n", "(in your local copy of astropy - astroquery/jplhorizons/__init__.py - you may have to make a change.\n", "Replace http below with https:\n", "'http://ssd.jpl.nasa.gov/horizons_batch.cgi',)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def get_mpc_primary(inp):\n", " \"\"\" extract primary ID from an MPC-1 string parameter\n", " ---------\n", " input: MPC-1 data, strip of column 166-190\n", " \"\"\"\n", "\n", " if inp[0:1] == '(':\n", " # (433) Eros\n", " return inp[inp.find(\"(\")+1:inp.find(\")\")]\n", " else:\n", " # 2018 CL\n", " return inp\n", "\n", "for i, obj in enumerate(mpcorb):\n", " obj_name = get_mpc_primary(obj[166:190].strip())\n", " if not obj_name == '2018 CL':\n", " continue\n", " else:\n", " break\n", "\n", "obj = Horizons(id=obj_name, location='I41', \\\n", " epochs={'start': '{}-{}-{}'.format(start_date[0:4], start_date[4:6], start_date[6:8]), \\\n", " 'stop': '{}-{}-{}'.format(end_date[0:4], end_date[4:6], end_date[6:8]), \\\n", " 'step': '1h'})\n", "\n", "eph = obj.ephemerides(skip_daylight=True, airmass_lessthan=5.0)\n", "if len(eph) <= 0:\n", " print('object {} was not observable at Palomar'.format(obj_name))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We also want to calculate the motion of the asteroid to see if it will show up as a streak in our images." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "rate = np.sqrt(eph['RA_rate']**2+eph['DEC_rate']**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, since we just discovered the power of JPL Horizons... why don't we take a short detour and do something fun? Let's say, we want to see the predicted trajectory of 2018 CL..." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc4AAAEWCAYAAADvi3fyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXlYlcfZuO9hU4wIakDhIJEtBM4BXFBIjYoQNIkJJtTU2lptTWI+l5qaaG2/tvmZtmloG2PU2i81GsXaxtYlYKImrlnUIho1LsQEERSRmoBgEgEN+vz+eA8nRzgHjgsiZu7rei8488zyzLzL884z884oEUGj0Wg0Go1ruLW2AhqNRqPRtCW04dRoNBqN5grQhlOj0Wg0mitAG06NRqPRaK4AbTg1Go1Go7kCtOHUaDQajeYK0IYTUEq9opT6TWvrcT1RSj2ilCpRSn2llOrd2vq0NZRSA5RSBdb2e7i19dFoNDcPbd5wKqWKlVL3XkseIvI/IvK7a9Tjx0qp7deSx3XmRWCKiHQUkX32AqVUgFLqdaXUKaXUWaXUDqVUYoM4P1BKHVdKnVNKZSulutjJpiil9iilziulljYsWCn1PaXUx0qpL5VS+c0ZHqVUf6XUeqVUlVLqjFIqTyn1E6ssWSl18loa4ir5LfAXa/tl38iClVLtlFKLre3/pVJqv1Lq/gZxUpVSR5RS1UqpbUqpO+xk31NK7bTK3nWQf4pSaq9S6gul1DGl1IQGcofn/jrotVQpdcH6MlJ/uDtpg2u6n5RS45RSH1rreFIp9SellIedvItS6g1rHY8rpX5gJwtUSq213h+ilOrZIO8uSql/KaUqlFLlSql/KKU62cl7WutebW2Le+1kV62XVe6vlPqn9b6tVEr942rbSHP1tHnD2Rz2F+XNTAvoeQdw2ImsI7Ab6At0AbKAdUqpjlZdzMDfgB8B3YBq4K926U8Bvwdea5ixUsoELAeeBjoBM4B/KqUCHCmilLob2Aq8B0QAXYGJwP2O4t9AnLafMmjJe8cDKAEGA77Ar4F/1z/AlVK3A2uA32Ccvz3Av+zSnwFeBjId6O4JvIFxfn2BUcBLSql4q7ypc3+tegH8yfoyUn9cvMK2cZUOwM+A24FEIBWYbidfAFzAqOMPgf+z1h3gEvA28F0nef8e6AyEAuHWPGbZyV8H9mFcy78CViml/K+DXmC073+BECAA4wVZc6MRkTZ7AH/HuMhrgK+AnwM9AQEeA04A71vjrsS44M4C7wNmu3yWAr+3+/0gsB+oAnYCcXayHhgX7+dABfAXIBqoBS5a9aiyxvUFllnjHsd40LhZZT8GdgBzrPn8AeOBF2tXVgDGg8vfQd3drPkdBz6zluMLtLPqIMA5oNDFtvwC6Gv9/w/AP+1k4Rg3s0+DNL8HljYISwQ+axD2OXC3k3K3Awua0CsZOOliHYKAtdZ2PAo8YSebBfzb2k5fYhjFBCf5FDa4rtoB7wLPW89ZDYaRb668lRgvEV8CB4E7gV9az1cJMPQKrvUDwHet/08AdtrJbrPqdFeDNI8D7zYI62a9NjrYhe0GRl/Jub8avWhwnzVR1yu+n1zI82ngTTu9LgB3NniWZDZI42Ftq54NwjcAk+x+Twbesf5/J3Devr2AD4D/uVa9gKFAMeDu6nWjj5Y52nSPU0R+hGEcHxLj7fVPduLBGDfgMOvvDUAkhjHaCzh0cShjPPA14EmMN8a/AWutbip34C2Mm7YnYAJWiMjHwP8A/7Hq4WfNbj7GzR5m1Wcs8BO74hKBYxgPs98BK4AxdvLRwBYR+dyBqj+2HkOs+XfEcC2eF5GO1jjxIhLuqJ4N6twL8MJ4+AOYgY/q5SJSiPWGbi4vjF7Gx0qpdKWUuzLctOcxHrANy+0A3A2sciFfV1gBnMQwaCOBPyilUuzk6dY4fhgG7y+OMrG2mf11dd4q+hGGcfDBuAaaK+8hjAdfZ4weyDsYLzwmDFfw31yplFKqG0bb1/eAG56fcxjG3tw4daO6ncboEf3Een7uxuhd17tFXT73V6nXJGW44z9USjns0V3D/dQUg+z0vBOoE5FP7eQf4UL7WVkAPKiU6qyU6ozRM91glZmBYyLypYt5X4leScAnQJbVTbxbKTXYRZ0115E2bTibYZaInBORGgAReU1EvrQ+BGcB8UopXwfpJgB/E5FdInJRRLIwHvxJQH+Mh+QMa961IuJwHMZqZL8P/NJabjEwG+PhW88pEZkvInVWPbOA0UopZZX/COPB64gfAi+JyDER+QqjJ/P9K3X5Wsdm/g48JyJnrcEdMXrm9pzFMBhNIobrbRnwT4x2+yfwpPUh2pDOGNdg2ZXo7AilVA9gADDTel72A4swHq71bBeR9VYd/w7EX2ExS0XksIjUAd1dKO8DEXnHGn8l4I/Re/gaw+j2VEr50QRW1+o/gCwROWINvurzY+V14FmM8/MB8CsRKbmSvK9Sr3l88/L6G2CpUmqAKwq7eD85SzseSOAbt2ZHDA+LMz2bYy/Gi2aF9bjIN+5sl8/NVegVjNHr3IZx/c0Gcqwucs0N5FY2nPUPAqxv1plKqUKl1BcY7g4wxhkacgfwjDImqlQppaow3LNB1r/HrQ/C5rgd8MTomdRzHKO30UhHABHZheGaTVZK3YXhDlzrJP8gB3l7YPReXUIp5Q28CeSKyAt2oq8wxift6YThcmwuz3uBP2G4WL0wegaLrL3ahlRiuEQDXdW5CYKAMw3e9Bu293/t/q8G2l/hi4b9+XKlvNN2/9cA5fLNmF6N9W9HnGAdR/07Ro9vip3oWs7PXRhGeyzG+TEDP1dKDXc176vVS0T2ikiF9UVxPYbhzWhOZyuu3E+NsHo8XgDuF5FyV/R0gX8Dn2IYtE4YverlV5L3VepVAxSLyGIR+VpEVmBcky69fGiuH7eC4XS2vYt9+A+AEcC9GK6entZwRWNKgOdFxM/u6CAir1tlIU4etg31KAe+xjDE9YQApc3onoXhrv0RsEpEah1VDmOCTsO867j8Ye0UpVQ7IBvD1fhkA/Fh7HpjSqkwjHG+T2meXhjjyntE5JKI7AZ2YbT9ZYhINfAfnE/CuBJOAV2UUvZv9g3b+1qxP18tWp7V67AY40Xou9Zeaj0Nz89tGGORziaD2WMBPrX2hC+JyCfAOr6ZjNXkub/OegmO78F6mT2u3E+XoZS6D3gVw+V+0E70KeChlIq0C4tvQs+G9MLwSp2zenteAR6wyg4DYQ2ui8vyvga9DtC4XfT2Vq1Baw+yXusB5AIT7H73xLiYPOzCJmFM9umEMQD/V2ucCKt8KdZJCxiukxKM8UdljT8c4+3SHWPM4UVreHtggDXdfRg9WS+7cpdjzGD0wbjhjwCPW2U/xnAdNqxPD4zJJseBQU3U+3GgAGNmX0eMccLldnJb/Ryk9cToaWbbt5Od3IzhMhporedyjLHcermHte4vYPQ82tfng9HDLAd6WX/3xnBnOZwIA3wH4017BtDVGhZfXx7WyUHWMuwP5SCvDzDGLdsDcRgvEfdaZbMatE+j66RBXsX1aa2/360/d1dZ3r0YvQX7NhQg2En5r2Bc2x0dyPwxXHjftZb9RwyvQb3c3Rr+PxgT4doDnlZZuLW9UzCu73CMse0JLp77a9FrJMa16obhcvwSSHZS/yu6nxykT7Fedw7vIYxe9+vWOg6w6m0/YbC9VSZAFNDeTrYNY7zV23r8lcsnReViPCPaA49gTDL0v1a9MGYqVwLjrOd4JMaz4vareXbq4+qPVlfgmitg9CRPWC/O6Tg2nB2BHOuNehzDTeXQcFp/34cx07AKY/xtJdZZchhvudnWi78cmGcN98J4cz+D4ZIDYwxvOcYswBKMcSX7WbWNDKdVttn60GhkHOziuFnzK7HmvxzobCdvynAOtsqrMR6i9cdAuzg/sLbrOWvbdbGTzbKmtz9m2cmnYDyMv8SY/PRMM+ewP8bkirPW9tsFjLXKkh2U5bBuGGNAb1nzKMRuJiMtYzivpDyXDSeGURCMmaX25+eHDfI7guG+exe7mZ/Wa6they21k38POGQ9PycxDJybndzhub8Oen1gPcdfYLyAfr+Ja+KK7icH6bdheGDs9dxgJ++CcR+fs9b1Bw3SN7rm7GShGC+eFVb93gYiG1xb71rb4JMG19G16jUQY4b2VxgT8QY6a0N9tNyhrCfjW41SahlwVER+29q6ACilXsOYOPTr1tZFo9FoNJfTJhYHaEms45VRwKbW1gWMVUcwJkzoZfI0Go3mJuRWmBx0rfwXwyW7urUVUUr9DsOF9mcRKWptfTQajUbTGO2q1Wg0Go3mCtA9To1Go9ForoA2McZ5++23S8+ePVtbDY1Go2lTfPjhh+Ui4t98TM2V0CYMZ8+ePdmzZ09rq6HRaDRtCqXU8eZjaa4U7arVaDQajeYKuOUNZ0lJCUOGDCEmJgaz2czcuXNtsjNnzpCWlkZkZCRpaWlUVlYCcOTIEe6++27atWvHiy9evt3dnDlzMJvNWCwWRo8eTW1t4xXxiouL8fb2plevXrbjwoULTnVcunQpU6ZMcSqvR0SYOnUqERERxMXFsXfvXpvM3d3dVlZ6errD9CtXrsRsNuPm5nZZDz4vL8+WNj4+njfeeMNh+qKiIhITE4mIiGDUqFG2Op0/f55Ro0YRERFBYmIixcXFzdZFo9Fo2iq3vOH08PBg9uzZ5Ofnk5uby4IFC8jPzwcgMzOT1NRUCgoKSE1NJTPT2Pe3S5cuzJs3j+nTp1+WV2lpKfPmzWPPnj0cOnSIixcvsmLFCoflhoeHs3//ftvh5eV1zXXZsGEDBQUFFBQUsHDhQiZOnGiTeXt728pau9bxuvAWi4U1a9YwaNCgRuF79uxh//79vP322zz55JPU1TVex37mzJlMmzaNo0eP0rlzZxYvXgzA4sWL6dy5M0ePHmXatGnMnDnzmuuq0Wg0Nyu3vOEMDAykT58+APj4+BAdHU1pqbEudE5ODuPGjQNg3LhxZGdnAxAQEEC/fv3w9PRslF9dXR01NTXU1dVRXV1NUFCQy7qcO3eO8ePH079/f3r37k1OTo5NVlJSQnJyMpGRkTz33HMO0+fk5DB27FiUUiQlJVFVVUVZmes7ckVHRxMVFdUovEOHDnh4GMPdtbW1fLOr2TeICFu3bmXkyJHA5e1l344jR45ky5Yt6M+cNBrNrcotbzjtKS4uZt++fSQmJgJw+vRpAgONHa26d+/O6dNNbyxiMpmYPn06ISEhBAYG4uvry9ChQx3GLSwstLk/J0+eDMDzzz9PSkoKeXl5bNu2jRkzZnDunLFNZV5eHqtXr+bAgQOsXLnS4WSo0tJSevToYfsdHBxsewmora0lISGBpKQkm0G7Enbt2oXZbCY2NpZXXnnFZkgfeOABTp06RUVFBX5+frZw+7Lt9fLw8MDX15eKioor1kGj0WjaAm1iVu3VkL2vlD+/8wmnqmoI8vNmysBgXpw6mpdffplOnRpueQdKKYc9LXsqKyvJycmhqKgIPz8/Hn30UZYvX86YMWMaxa131dqzceNG1q5daxs3ra2t5cSJEwCkpaXRtWtXADIyMti+fTsJCQku1/f48eOYTCaOHTtGSkoKsbGxhIeHu5w+MTGRw4cP8/HHHzNu3Djuv/9+2rdvz/r16wEoLy9vJgeNRqP5dnBL9jiz95XyyzUHKa2qQYCTFV/yxNjRxCU/SEbGN/vmduvWzebqLCsrIyAgoMl8N2/eTGhoKP7+/nh6epKRkcHOnTvZtWuXrXfpbHwRDHfn6tWrbWORJ06cIDo6GqCR0VZKsWDBAlu+p06dwmQyUVLyzV7KJ0+exGQy9vGt/xsWFkZycjL79u1zvcHsiI6OpmPHjhw6dOiy8K5du1JVVWUb+2xYdr1edXV1nD171vYSoNFoNLcat6Th/PM7n1Dz9UXAMFYVG+bi1jmYT/0HXxYvPT2drKwsALKyshgxYkST+YaEhJCbm0t1dTUiwpYtW4iOjiYxMdFmDJ3NaAUYNmwY8+fPt43/2Ru3TZs2cebMGWpqasjOzmbAgAFMnjzZlm9QUBDp6eksW7YMESE3NxdfX18CAwOprKzk/PnzgNEz3LFjBzExMS63V1FRkc0gHj9+nCNHjtBwwQmlFEOGDGHVqlWN2su+HVetWkVKSkqzvXeNRqNps7T2vmauHH379pUroefMt+QO69Hth38UQDz9e4pnQKjEx8fLunXrRESkvLxcUlJSJCIiQlJTU6WiokJERMrKysRkMomPj4/4+vqKyWSSs2fPiojIs88+K1FRUWI2m2XMmDFSW1vbqPyioiIxm82Nwqurq2XChAlisVgkJiZGhg8fLiIiS5YskREjRkhycrJERETIrFmzHNbr0qVLMmnSJAkLCxOLxSK7d+8WEZEdO3aIxWKRuLg4sVgssmjRIofp16xZIyaTSby8vCQgIECGDh0qIiLLli2TmJgYiY+Pl969e8sbb7xhS3P//fdLaWmpiIgUFhZKv379JDw8XEaOHGmre01NjYwcOVLCw8OlX79+UlhY2MwZ0mg0NwJgj9wEz/Bb7WgTi7wnJCTIlawcNCBzK6VVNY3Cfdp5sPvX99Le0/16qqfRaDQ3JUqpD0XE9ckSGpe4JV21M4ZF4d3AOLorxZfn63hg7gfsOqZnfGo0Go3m6rglDefDvU28kBGLyc8bBZj8vJn9vXiWje/PhYuXGLUwl/994yBf1H7d2qpqNBqNpo1xS7pqm6L6Qh1zNn3K4u1F+Pu047cjLAwzd78ueWs0Gs3NhHbVtgy3ZI+zKTp4efCr4TFkTx5A5w5ePPn3D5m4/EM++6LxmrMajUaj0TTkW2c464kL9uPNn97DjGFRbDnyGakvvceKvBO0hR64RqPRaFqPb63hBPB0d2PykAjefmogMYGd+MWag4x+NZei8nOtrZpGo9FoblK+1YaznjD/jrz+RBKZGbEcPvUFw15+n7++e5SvL15qbdU0Go1Gc5OhDacVNzfF9/uHsOXpwaTeFcCf3v6E9L/s4MDJqtZWTaPRaDQ3EdpwNiCgU3v+b0xfXhnTl4qvzvPwgh38/q18qi803p9So9FoNN8+tOF0wn2W7mx6ejDf7x/Cou1FDJ3zPu9/+nlrq6XRaDSaVkYbzibw9fbkD4/E8q8JSXi5uzH2tTye/vd+Ks9daG3VNBqNRtNKaMPpAolhXVn/1EB+mhLB2v2nuPel98jZX6o/XdFoNJpvIdpwukh7T3eeGRrFW1PvIbhLB55asZ/xS3c7XExeo9FoNLcu2nBeIXd178Said/hNw/GkHvsDGkvvcfSHUVcvNR877OkpIQhQ4YQExOD2Wxm7ty5NtmZM2dIS0sjMjKStLQ0KisrAThy5Ah333037dq148UXX7wsvzlz5mA2m7FYLIwePZra2sarHxUXF+Pt7W3bELtXr15cuODc1bx06VKmTJnSbF1EhKlTpxIREUFcXBx79+61ydzd3W1lOdufdOXKlZjNZtzc3LBfTnHTpk307duX2NhY+vbty9atWx2md9ZeTeml0Wg01wNtOK8CdzfFY/eEsnHaIBJ6dmHWm/mMfGUnn57+ssl0Hh4ezJ49m/z8fHJzc1mwYAH5+fkAZGZmkpqaSkFBAampqWRmZgLQpUsX5s2bx/Tp0y/Lq7S0lHnz5rFnzx4OHTrExYsXWbFihcNyw8PDbRti79+/Hy8vr2tugw0bNlBQUEBBQQELFy5k4sSJNpm3t7etrLVr1zpMb7FYWLNmDYMGDbos/Pbbb+fNN9/k4MGDZGVl8aMf/chhemft1ZReGo1Gcz3QhvMa6NGlA1k/6cfLo3pRXH6O4fM+4KVNn3K+7qLD+IGBgfTp0wcAHx8foqOjKS0tBSAnJ4dx48YBMG7cOLKzswEICAigX79+eHp6Nsqvrq6Ompoa6urqqK6uJigoyGXdz507x/jx4+nfvz+9e/cmJyfHJispKSE5OZnIyEiee+45h+lzcnIYO3YsSimSkpKoqqqirKzM5fKjo6OJiopqFN67d29bPcxmMzU1NZw/f95h+Y7a61r10mg0mubQhvMaUUrxcG8Tm58ezINxQczbUsADcz9gT/GZJtMVFxezb98+EhMTATh9+jSBgYEAdO/endOnTzeZ3mQyMX36dEJCQggMDMTX15ehQ4c6jFtYWGhznU6ePBmA559/npSUFPLy8ti2bRszZszg3DljqcG8vDxWr17NgQMHWLlyJY52piktLaVHjx6238HBwbaXgNraWhISEkhKSrIZtKth9erV9OnTh3bt2gHw+OOP23Rx1l5N6aXRaDTXA4/WVuBWoWvHdswZ1YsRvYL41RuHGPnKf/hR0h38/L4otnz8GX9+5xNOVdUQ5OfNlIHBvDh1NC+//DKdOnVqlJdSCqVUk+VVVlaSk5NDUVERfn5+PProoyxfvpwxY8Y0ilvvqrVn48aNrF271jZuWltby4kTJwBIS0uja9euAGRkZLB9+3YSElzfmej48eOYTCaOHTtGSkoKsbGxhIeHu5we4PDhw8ycOZONGzfawhYtWuQwrivtpdFoNNcL3eO8ziRHBbBx2iDGDwhl+a7jDMjcys9XHaC0qgYBTlZ8yRNjRxOX/CAZGRm2dN26dbO5FMvKyggICGiynM2bNxMaGoq/vz+enp5kZGSwc+dOdu3aZetdOhtfBGMSzerVq21jkSdOnCA6OhqgkRFSSrFgwQJbvqdOncJkMlFSUmKLc/LkSUwmE4Dtb1hYGMnJyezbt8/1BrTm9cgjj7Bs2TKnBtdZezWll0aj0VwPWsxwKqVeU0p9ppQ6ZBfWSymVq5Tar5Tao5Tq31Lltya3tfPg2YdiWDPxO1RfuMgF62LxIkLFhrm4dQ7mU//Bl6VJT08nKysLgKysLEaMGNFkGSEhIeTm5lJdXY2IsGXLFqKjo0lMTLQZQ2czWgGGDRvG/Pnzbd+i2hu3TZs2cebMGWpqasjOzmbAgAFMnjzZlm9QUBDp6eksW7YMESE3NxdfX18CAwOprKy0jUmWl5ezY8cOYmJiXG67qqoqhg8fTmZmJgMGDHAaz1l7OdNLo9Forhsi0iIHMAjoAxyyC9sI3G/9/wHgXVfy6tu3r7RVes58S+6wHt1++EcBxNO/p3gGhEp8fLysW7dORETKy8slJSVFIiIiJDU1VSoqKkREpKysTEwmk/j4+Iivr6+YTCY5e/asiIg8++yzEhUVJWazWcaMGSO1tbWNyi8qKhKz2dwovLq6WiZMmCAWi0ViYmJk+PDhIiKyZMkSGTFihCQnJ0tERITMmjXLYb0uXbokkyZNkrCwMLFYLLJ7924REdmxY4dYLBaJi4sTi8UiixYtcph+zZo1YjKZxMvLSwICAmTo0KEiIvK73/1OOnToIPHx8bbj9OnTIiLy2GOP2cpx1l7O9NJovo0Ae6SFnvHf5kOJtNzqN0qpnsBbImKx/n4HeE1E/qWUGg08JCI/aC6fhIQEcTRBpS0wIHOrw0USAnzakfere1tBI41G821BKfWhiLg+QUHjEjd6jPNnwJ+VUiXAi8AvnUVUSk2wunP3fP55211cfcawKLw93RuFV3x1nvlbCrhQp/f81Gg0mrbEjTacE4FpItIDmAYsdhZRRBaKSIKIJPj7+98wBa83D/c28UJGLCY/bxRg8vPmtyNiuC82kNmbPuXB+R/w4fHK1lZTo9FoNC5yo121ZwE/ERFlTN08KyKNv8doQFt21TbFlo9P85vsQ5R9UcuPku5gxrAofNo3XuhAo9Forgbtqm0ZbnSP8xRQP500BSi4weXfVKRGd2Pj04P58Xd68vfc46S99D7vHP5va6ul0Wg0miZoyc9RXgf+A0QppU4qpR4DngBmK6U+Av4ATGip8tsKHdt58P8eMvPGpAH4dfDkyb9/yP/8/UNOf9F4wXaNRqPRtD4t6qq9XtyqrtqGfH3xEq9+cIy5mwvwcndj5v138YP+Ibi56VVxNBrNlaNdtS2DXjnoJsLT3Y1JyRG887NBxAb78uvsQ4xa+B+Oftb0risajUajuXFow3kT0vP22/jH44n8eWQcBZ99xf1zP2BOE7uuaDQajebGoQ3nTYpSikcTerD56cE8EBvIXOuuK3lFTe+6otFoNJqWRRvOm5zbO7Zj7vd7s+Qn/aj9+hLf+9t/+OWag5yt+bq1VdNoNJpvJdpwthGGRAWw6elBPDEwlH/tPkHaS++x4WAZbWFyl0aj0dxKaMPZhujg5cGvhseQM/ke/H3aMfEfe3li2YeUnW28Fq5Go9FoWgZtONsgscG+5EwewP8+cBfbj35O2kvvk7WzmIuXdO9To9FoWhptONsoHu5uTBgUzsafDaZ3iB//b+1hRr6yk0/+2/ynKyUlJQwZMoSYmBjMZjNz5861yc6cOUNaWhqRkZGkpaVRWWmso3vkyBHuvvtu2rVrx4svvnhZfnPmzMFsNmOxWBg9ejS1tY0XbyguLsbb29u2GXavXr24cOGCUx2XLl3KlClTmq2LiDB16lQiIiKIi4tj7969Npm7u7utLGd7k65cuRKz2Yybmxv23wpXVFQwZMgQOnbs2KQeztqrKb00Gk3bRhvONk5I1w4sG9+fOaPiOV5RzfB5H/DiO59Q+7XzT1c8PDyYPXs2+fn55ObmsmDBAvLz8wHIzMwkNTWVgoICUlNTyczMBKBLly7MmzeP6dOnX5ZXaWkp8+bNY8+ePRw6dIiLFy+yYsUKh+WGh4fbNsPev38/Xl5e11z/DRs2UFBQQEFBAQsXLmTixIk2mbe3t62stWvXOkxvsVhYs2YNgwYNuiy8ffv2/O53v2v0ktAQZ+3VlF4ajaZtow3nLYBSikd6B7P56cGk9wriL9uOcv/cD/hPYYXD+IGBgfTp0wcAHx8foqOjKS0tBSAnJ4dx48YBMG7cOLKzswEICAigX79+eHo2XoS+rq6Ompoa6urqqK6uJigoyGXdz507x/jx4+nfvz+9e/cmJyfHJispKSE5OZnIyEiee+45h+lzcnIYO3YsSimSkpKoqqqirKzM5fKjo6OJiopqFH7bbbdxzz330L59+ybTO2uva9VLo9HcvGjDeQvR5TYvXvpl/0m+AAAgAElEQVReL5Y/lsjFS8LoV3OZueoAVdXOXaLFxcXs27ePxMREAE6fPk1gYCAA3bt35/Tp002WaTKZmD59OiEhIQQGBuLr68vQoUMdxi0sLLS5TidPngzA888/T0pKCnl5eWzbto0ZM2Zw7tw5APLy8li9ejUHDhxg5cqVOFp2sbS0lB49eth+BwcH214CamtrSUhIICkpyWbQrgePP/64TRdn7dWUXhqNpm3j0doKaK4/90Tezjs/G8TcLQW8+sExthw5zf97yEzdxUu8uPFTTlXVEOTnzZSBwbw4dTQvv/wynTo13t1NKYWx+5tzKisrycnJoaioCD8/Px599FGWL1/OmDFjGsWtd9Xas3HjRtauXWtzidbW1nLixAkA0tLS6Nq1KwAZGRls376dhATXl908fvw4JpOJY8eOkZKSQmxsLOHh4S6nd8aiRYschrvSXhqNpu2je5y3KN5e7vzi/rtYO2UAQX7e/PT1fTyz8iNKq2oQ4GTFlzwxdjRxyQ+SkZFhS9etWzebS7GsrIyAgIAmy9m8eTOhoaH4+/vj6elJRkYGO3fuZNeuXbbepbPxRTAm0axevdo2FnnixAmio6MBGhkhpRQLFiyw5Xvq1ClMJhMlJSW2OCdPnsRkMgHY/oaFhZGcnMy+fftcb0AXcdZeTeml0WjaNtpw3uKYg3x5Y9IAfL09qP9aRUSo2DAXt87BfOo/+LL46enpZGVlAZCVlcWIESOazD8kJITc3Fyqq6sREbZs2UJ0dDSJiYk2Y+hsRivAsGHDmD9/vm0hB3vjtmnTJs6cOUNNTQ3Z2dkMGDCAyZMn2/INCgoiPT2dZcuWISLk5ubi6+tLYGAglZWVnD9/HoDy8nJ27NhBTEzMFbdfczhrL2d6aTSaWwARuemPvn37iuba6DnzLbnDenT74R8FEE//nuIZECrx8fGybt06EREpLy+XlJQUiYiIkNTUVKmoqBARkbKyMjGZTOLj4yO+vr5iMpnk7NmzIiLy7LPPSlRUlJjNZhkzZozU1tY2Kr+oqEjMZnOj8OrqapkwYYJYLBaJiYmR4cOHi4jIkiVLZMSIEZKcnCwREREya9Ysh/W6dOmSTJo0ScLCwsRiscju3btFRGTHjh1isVgkLi5OLBaLLFq0yGH6NWvWiMlkEi8vLwkICJChQ4faZHfccYd07txZbrvtNjGZTHL48GEREXnsscds5ThrL2d6aTQ3EmCP3ATP8Fvt0PtxfksYkLmV0qrGKwx1bOdB3q9S6eClh7s1mlsNvR9ny6Bdtd8SZgyLwtvT/bIwdzfFV+frGDrnfd779PNW0kyj0WjaFtpwfkt4uLeJFzJiMfl5owCTnzezH43n30/eTTsPN8a9lsdTK/ZR/tX51lZVo9Fobmq0q1bD+bqL/N+7hfx1WyHeXu786oFoHk0I1p9WaDRtHO2qbRl0j1NDOw93fnbvnax/6h6iuvnw89UHGP1qLsc+/6q1VdNoNJqbDm04NTYiAnxYMSGJFzJiyT/1BffN/YD5Wwq4UHeptVXTaDSamwZtODWX4eamGN0/hM3PDGZoTDdmb/qU4fM+YE/xmdZWTaPRaG4KtOHUOCTApz1/+UEflvy4H9UXLjLylf/wv28c5GzN162tmkaj0bQq2nBqmmTIXQFsnDaIx+8JZUXeCe596T3WHyyjLUwq02g0mpZAG05Ns9zWzoNfPxhDzuR76NapHZP+sZfHs/Y4XFBBo9FobnW04dS4TGywL9mTBvDr4dHsLKwg7aX3eG17ERcv6d6nRqP59tBihlMp9ZpS6jOl1KEG4T9VSh1RSh1WSv2ppcrXtAwe7m48PjCMjdMG0T+0C799K59H/rqDw6fOtrZqGo1Gc0NoyR7nUuA++wCl1BBgBBAvImbgxRYsX9OC9OjSgSU/7sf80b05VVVD+l928ML6j6m5cLG1VdNoNJoWpcUMp4i8DzT8hmEikCki561xPmup8jUtj1KKh+KD2PJ0Mt9LCOZv7x9j6MvvubTubUlJCUOGDCEmJgaz2czcuXNtsjNnzpCWlkZkZCRpaWlUVlYCcOTIEe6++27atWtn2/i6njlz5mA2m7FYLIwePZra2tpGZRYXF+Pt7W3bz7NXr15cuHDBqY5Lly5lypQpzdZFRJg6dSoRERHExcWxd+9em8zd3d1WlrPt1VauXInZbMbNzY2GK2S98MILREREEBUVxTvvvOMwfVFREYmJiURERDBq1Chbnc6fP8+oUaOIiIggMTGR4uLiZuui0Wia50aPcd4JDFRK7VJKvaeU6ucsolJqglJqj1Jqz+ef6wXIb2Z8O3jyQkYc/5qQhKe7a+veenh4MHv2bPLz88nNzWXBggXk5+cDkJmZSWpqKgUFBaSmppKZmQlAly5dmDdvHtOnT78sr9LSUubNm8eePXs4dOgQFy9eZMWKFQ7LDQ8Pt+3nuX//fry8vK65/hs2bKCgoICCggIWLlzIxIkTbTJvb29bWc429LZYLKxZs4ZBgwZdFp6fn8+KFSs4fPgwb7/9NpMmTeLixcY9+pkzZzJt2jSOHj1K586dWbx4MQCLFy+mc+fOHD16lGnTpjFz5sxrrqtGo7nxhtMD6AIkATOAfysnC6KKyEIRSRCRBH9//xupo+YqSQzryoanBvKzeyPZcPC/pM5+j3/vKXH46UpgYCB9+vQBwMfHh+joaEpLSwHIyclh3LhxAIwbN47s7GwAAgIC6NevH56eno3yq6uro6amhrq6OqqrqwkKCnJZ73PnzjF+/Hj69+9P7969ycnJsclKSkpITk4mMjKS5557zmH6nJwcxo4di1KKpKQkqqqqKCsrc7n86OhooqKiHOb7/e9/n3bt2hEaGkpERAR5eXmXxRERtm7dysiRI4HL28u+HUeOHMmWLVv0Z0QazXXgRhvOk8Aa6x6recAl4PYbrIOmBWm07u2q5te9LS4uZt++fSQmJgJw+vRpAgMDAejevTunT59uskyTycT06dMJCQkhMDAQX19fhg4d6jBuYWGhzXU6efJkAJ5//nlSUlLIy8tj27ZtzJgxg3PnzgGQl5fH6tWrOXDgACtXrmzkSgWjx9ujRw/b7+DgYNtLQG1tLQkJCSQlJdkMmqs0le8DDzzAqVOnqKiowM/PDw8Pj0Zx7NN7eHjg6+tLRUXFFemg0Wgac6N3L84GhgDblFJ3Al5A+Q3WQXMDqF/39l97Snhh/cfcN/cDfjokgkDf9szZXMCpqhqC/LyZMjCYF6eO5uWXX6ZTp06N8lFKNbtLS2VlJTk5ORQVFeHn58ejjz7K8uXLGTNmTKO49a5aezZu3MjatWtt46a1tbWcOHECgLS0NLp27QpARkYG27dvJyHB9c0mjh8/jslk4tixY6SkpBAbG0t4eLjL6Z2xfv16AMrL9e2j0dxoWvJzlNeB/wBRSqmTSqnHgNeAMOsnKiuAcaJ9R7csjta9nbHqAKVVNQhwsuJLnhg7mrjkB8nIyLCl69atm83VWVZWRkBAQJPlbN68mdDQUPz9/fH09CQjI4OdO3eya9cuW+/S2fgiGO7O1atX28YiT5w4QXR0NEAjo62UYsGCBbZ8T506hclkoqSkxBbn5MmTmEwmANvfsLAwkpOT2bdvn8vt11S+9XTt2pWqqirq6uocll2fvq6ujrNnz9peAjQazdXTkrNqR4tIoIh4ikiwiCwWkQsiMkZELCLSR0S2tlT5mpuH+nVvu97mRf1bkohQsWEubp2D+dR/8GXx09PTycrKAiArK4sRI0Y0mX9ISAi5ublUV1cjImzZsoXo6GgSExNtxtDZjFaAYcOGMX/+fNv4n71x27RpE2fOnKGmpobs7GwGDBjA5MmTbfkGBQWRnp7OsmXLEBFyc3Px9fUlMDCQyspKzp83JkiVl5ezY8cOYmJiXG639PR0VqxYwfnz5ykqKqKgoID+/ftfFkcpxZAhQ1i1alWj9rJvx1WrVpGSkqL3WNVorgcictMfffv2FU3bp+fMt+QO69Hth38UQDz9e4pnQKjEx8fLunXrRESkvLxcUlJSJCIiQlJTU6WiokJERMrKysRkMomPj4/4+vqKyWSSs2fPiojIs88+K1FRUWI2m2XMmDFSW1vbqPyioiIxm82Nwqurq2XChAlisVgkJiZGhg8fLiIiS5YskREjRkhycrJERETIrFmzHNbr0qVLMmnSJAkLCxOLxSK7d+8WEZEdO3aIxWKRuLg4sVgssmjRIofp16xZIyaTSby8vCQgIECGDh1qk/3+97+XsLAwufPOO2X9+vW28Pvvv19KS0tFRKSwsFD69esn4eHhMnLkSFvda2pqZOTIkRIeHi79+vWTwsLCJs6O5lYE2CM3wTP8VjuU0bY3NwkJCeJoUoambTEgc6vD9W3be7ixZXoyJj/vVtBKo7l1UUp9KCKuD8prXEKvVau5YcwYFoW3p/tlYZ7uiouXhLSX3mOxXvdWo9G0AbTh1NwwHu5t4oWMWEx+3ijA5OfNn0fGs3V6MomhXfiddd3bQ6V63VuNRnPzol21mpsCEWHdwTJmrc2nsvoC4wf0ZFranXTwutFfTGk0tw7aVdsy6B6n5qZAKcWDcUFseXow30vowasfFJH20vts+0QvZ6zRaG4utOHU3FQY697G8u8n78bby52fLNnNT1/fx+dfOl/3VqPRaG4k2nBqbkr6h3Zh3dR7mHbvnbxz6L+kzn6XFXknuKQnD2k0mlZGG07NTUs7D3eeujeSDT8bSHRgJ36x5iDfX5jL0c+cr3ur0Wg0LY02nJqbnnD/jqyYkMSfvhvHJ6e/5IG5HzBn06ecr9ObZms0mhuPNpyaNoFSiu/168Hmpwdzn6U7c7cU8MDcD9h1TO/2odFobixNGk6lVLBSarpSKkcptVsp9b5S6q9KqeFKKW10NTccf592zBvdm6U/6cf5ukuMWpjLL1Yf4Gz1162tmkaj+Zbg9DtOpdQSwAS8BewBPgPaA3dibA3WF/iFiLzf0krq7zg1jqi+UMfczQUs2l5E5w6ePPuQmYfiAvVC5hqNFf0dZ8vQlOG0iMghpwmV8gJCRORoSylXjzacmqY4fOosv1xzkAMnzzL4Tn9+/7CFHl06tLZaGk2row1ny+DU3dqU0bTKL9wIo6nRNIc5yJc3Jg3g2Qdj2F18hrQ577Hw/ULqLl5qMl1JSQlDhgwhJiYGs9nM3LlzbbIzZ86QlpZGZGQkaWlpVFZWAnDkyBHuvvtu2rVrZ9v4up45c+ZgNpuxWCyMHj2a2traRmUWFxfj7e1t28+zV69eXLhwwamOS5cuZcqUKc22gYgwdepUIiIiiIuLY+/evTaZu7u7rSxn26s5q29lZSWPPPIIcXFx9O/fn0OHHD8WioqKSExMJCIiglGjRtnqdP78eUaNGkVERASJiYkUFxc3WxeN5man2XFKpdRBpdSBBscHSqk5Sim9K67mpsDdTTH+nlA2PT2YeyJu5w/rj5D+lx0cOFnlNI2HhwezZ88mPz+f3NxcFixYQH5+PgCZmZmkpqZSUFBAamoqmZmZAHTp0oV58+Yxffr0y/IqLS1l3rx57Nmzh0OHDnHx4kVWrFjhsNzw8HDbfp779+/Hy8vrmuu/YcMGCgoKKCgoYOHChUycONEm8/b2tpXlbENvZ/X9wx/+QK9evThw4ADLli3jqaeecph+5syZTJs2jaNHj9K5c2cWL14MwOLFi+ncuTNHjx5l2rRpzJw585rrqtG0Nq5M8NkArAN+aD3exBjz/C+wtMU002iuApOfN6+OTeD/ftiH8q/O8/CCHfz2zXzOna9rFDcwMJA+ffoA4OPjQ3R0NKWlpQDk5OQwbtw4AMaNG0d2djYAAQEB9OvXD09Pz0b51dXVUVNTQ11dHdXV1QQFBbms97lz5xg/fjz9+/end+/e5OTk2GQlJSUkJycTGRnJc8895zB9Tk4OY8eORSlFUlISVVVVlJWVuVy+s/rm5+eTkpICwF133UVxcTGnT5++LK2IsHXrVkaOHNkovX2+I0eOZMuWLbSF9bE1mqZwxXDeKyK/FJGD1uNXwGAR+SPQs2XV02iuHKUU98cGsvmZwfwgMYTXdhSR9tJ7bM4/7TRNcXEx+/btIzExEYDTp08TGBgIQPfu3RsZi4aYTCamT59OSEgIgYGB+Pr6MnToUIdxCwsLba7TyZMnA/D888+TkpJCXl4e27ZtY8aMGZw7dw6AvLw8Vq9ezYEDB1i5ciWOxvtLS0vp0aOH7XdwcLDtJaC2tpaEhASSkpJsBq0hzuobHx/PmjVrbHocP36ckydPAvDAAw9w6tQpKioq8PPzw8PDo1HZ9np5eHjg6+tLRYX+hEjTtnFl6wl3pVR/EckDUEr1A+o3VWz8Gq/R3CR0au/J7x+O5ZHeJn655iCPL9vDA7HdSQrtwt/eL+JUVQ1Bft5MGRjMi1NH8/LLL9OpU6dG+Silmp2pW1lZSU5ODkVFRfj5+fHoo4+yfPlyxowZ0yhuvavWno0bN7J27VrbuGltbS0nTpwAIC0tja5djVGRjIwMtm/fTkKC6/M9jh8/jslk4tixY6SkpBAbG0t4eLjT+Pb1/cUvfsFTTz1Fr169iI2NpXfv3ri7G7f/+vXrASgvL3dZF43mVsAVw/k48JpSqqP195fA40qp24AXWkwzjeY60feOLrz104EsfL+QOZs+Zf3B/9pkJyu+5Imxo3n04QfJyMiwhXfr1o2ysjICAwMpKysjICCgyTI2b95MaGgo/v7+gGHgdu7cSWRkJE8++SQAv/3tb4mLi3OYXkRYvXo1UVFRl4Xv2rWrkdFWSrFgwQJeffVVwDBgJpOJkpKSb+p18iQmkwnA9jcsLIzk5GT27dvXyHA6q2+nTp1YsmSJTcfQ0FDCwsIuS9u1a1eqqqqoq6vDw8OjUdklJSUEBwdTV1fH2bNnbS8BGk1bpVlXrYjsFpFYoBfQS0TiRCRPRM6JyL9bXkWN5trx8nBjSkokt/u0s4WJCBUb5uLWOZhP/QdfFj89PZ2srCwAsrKyGDFiRJP5h4SEkJubS3V1NSLCli1biI6OJjEx0TYxx9mMVoBhw4Yxf/582/jfvn37bLJNmzZx5swZampqyM7OZsCAAUyePNmWb1BQEOnp6SxbtgwRITc3F19fXwIDA6msrOT8eWNnmfLycnbs2EFMTEyj8p3Vt6qqyjZDdtGiRQwaNKhRr1wpxZAhQ1i1alWj9Pb5rlq1ipSUFP2drabtIyJNHkA3YDGwwfo7BnisuXTX8+jbt69oNNeDnjPfkjusR7cf/lEA8fTvKZ4BoRIfHy/r1q0TEZHy8nJJSUmRiIgISU1NlYqKChERKSsrE5PJJD4+PuLr6ysmk0nOnj0rIiLPPvusREVFidlsljFjxkhtbW2j8ouKisRsNjcKr66ulgkTJojFYpGYmBgZPny4iIgsWbJERowYIcnJyRIRESGzZs1yWK9Lly7JpEmTJCwsTCwWi+zevVtERHbs2CEWi0Xi4uLEYrHIokWLHKZ3Vt+dO3dKZGSk3HnnnfLII4/ImTNnbGnuv/9+KS0tFRGRwsJC6devn4SHh8vIkSNtda+pqZGRI0dKeHi49OvXTwoLC5s5Q5rrCbBHbuCz+ttyOF0AoR6l1AZgCfArEYlXSnkA+8Tohd4Q9AIImuvFgMytlFbVNAr3cFMsfzyRpDDtRtTcOugFEFoGV2bV3i6GS/YSgIjUAXpbCk2bZMawKLw93S8L83J3o5O3B9/X695qNBoXcMVwnrMudCAASqkk4GyLaqXRtBAP9zbxQkYsJj9vFMZ3n38aGcf2mSk8OSiMlR+eJPWl93jzo1P6e0ONRuMQV1y1fYD5gAU4BPgDI0XkQDPpXgMeBD4TEUsD2TPAi4C/iDQ7l127ajU3ikOlxrq3B0vPknJXAL8dYSa4s173VtM20a7alsGVWbV7gcHAd4AnAXNzRtPKUuC+hoFKqR7AUODEFWmq0dwALCZf3pj0HX7zYAy5xyoYOud9Fm8v4uIl3fvUaDQGTr/jVEplOBHdqZRCRNY0lbGIvK+U6ulANAf4OZDjQKbRtDoe7m48dk8ow8zd+HX2IX73Vj45+0t5ISMWc5Bva6un0WhamaYWQHjI+jcAo7e51fp7CLATaNJwOkIpNQIoFZGP9Ldcmpud4M4dWPLjfrx1oIzn3jxM+l928PjAUH6WeifeXu7NZ6DRaG5JnBpOEfkJgFJqIxAjImXW34FcxeLuSqkOwP9iuGldiT8BmADGx+UaTWuglOKh+CAGRfrzwoaP+dt7x1h/sIznH45l0J3+ra2eRqNpBVyZVduj3mhaOQ1cjSULB0KBj5RSxUAwsFcp1d1RZBFZKCIJIpJQv4yZRtNa+HbwJPO7cayYkISnmxtjX8tj2r/2U/HV+dZWTaPR3GBcMZxblFLvKKV+rJT6McYWY5uvtCAxdlYJEJGeItITOAn0EZH/NpNUo7lpSArryvqnBjI1JYK3Dpzi3pfeY/WHJ/WnKxrNtwhXZtVOAV4B4q3HQhH5aXPplFKvA/8BopRSJ5VSj12rshrNzUB7T3eeHhrFuqkDCfPvyDMrP2LM4l0Ul59rbdU0Gs0NwOl3nEopJc28RrsS53qgv+PU3KxcuiT8M+8Ef9xwhAsXL/HUvZE8MTAMT3dXnDkaTcuiv+NsGZq6u7cppX6qlLpsPFMp5aWUSlFKZQHjWlY9jebmxs1NMSbpDjY/M5iUuwL409uf8ND87ewvqWo2bUlJCUOGDCEmJgaz2czcuXNtsjNnzpCWlkZkZCRpaWlUVlYCcOTIEe6++27atWtn27uznjlz5mA2m7FYLIwePZra2tpGZRYXF+Pt7W3bSLtXr1623U8csXTpUqZMmdJsXUSEqVOnEhERQVxcHHv37rXJ3N3dbWU52yHGWX3Pnj3LQw89RHx8PGaz2bbFWUM+/PBDYmNjiYiIYOrUqTbXubN8NZproSnDeR/GmrSvK6VOKaXylVJFQAEwGnhZRJbeAB01mpuebp3a839j+rLwR32pqv6aR/66g1lrD/PVeed7vXt4eDB79mzy8/PJzc1lwYIF5OfnA5CZmUlqaioFBQWkpqaSmZkJQJcuXZg3bx7Tp0+/LK/S0lLmzZvHnj17OHToEBcvXmTFihUOy63fSLv+8PLyuub6b9iwgYKCAgoKCli4cCETJ060yby9vW1lrV271mF6Z/VdsGABMTExfPTRR7z77rs888wzDg39xIkTefXVV206vP32203mq9FcC04Np4jUishfRWQAcAeQCvQWkTtE5AkR2ecsrUbzbWWouTubnh7E2KQ7yPpPMWkvvcfm/NMO4wYGBtKnTx8AfHx8iI6OprS0FICcnBzGjTMcOuPGjSM7OxuAgIAA+vXrh6enZ6P86urqqKmpoa6ujurqaoKCglzW+9y5c4wfP57+/fvTu3dvcnK+WZ+kpKSE5ORkIiMjee655xymz8nJYezYsSilSEpKoqqqirKyModxnaV3VF+lFF9++SUiwldffUWXLl3w8Lj8K7qysjK++OILkpKSUEoxduxYW3pn+Wo014JLAzEi8rWIlIlI8/4njeZbjk97T54bYWH1xO/Qqb0njy/bw+R/7OWzLxq7TuspLi5m3759JCYmAnD69GkCAwMB6N69O6dPOza+9ZhMJqZPn05ISAiBgYH4+voydKjjT6YLCwttrtPJkycD8Pzzz5OSkkJeXh7btm1jxowZnDtnTHbKy8tj9erVHDhwgJUrV+JovkFpaSk9evSw/Q4ODra9BNTW1pKQkEBSUpJTw+WsvlOmTOHjjz8mKCiI2NhY5s6di5ub8djq1auXrezg4GCHZV9pO2o0rtDUykEajeYa6BPSmbem3sPC948xd0sB7xd8zi/vj8bbw40XN33Kqaoagvy8mTIwmBenjubll1+mU6dOjfJRStHcSluVlZXk5ORQVFSEn58fjz76KMuXL2fMmDGN4ta7au3ZuHEja9eutY2b1tbWcuKEsZx0WloaXbsa+5RmZGSwfft2EhJcn29y/PhxTCYTx44dIyUlhdjYWMLDw53Gt6/vO++8Q69evdi6dSuFhYWkpaUxcOBAOnXq1KgOzeFKO2o0rqCn/mk0LYinuxuTh0Twzs8GEWvy5X/fOMgzqz6itKoGAU5WfMkTY0cTl/wgGRnfLA/drVs3m6uzrKyMgICAJsvZvHkzoaGh+Pv74+npSUZGBjt37mTXrl223qWz8UUwJvesXr3aNhZ54sQJoqOjARoZG6UUCxYssOV76tQpTCYTJSUltjgnT57EZDIB2P6GhYWRnJzMvn2NR3mc1XfJkiVkZGSglCIiIoLQ0FCOHDlyWVqTycTJkycdln2l7ajRuIJTw6mUilBKDXAQPkAp5fx1UaPRNCL09tv4x+OJ+HXwpH6jFRGhYsNc3DoH86n/4Mvip6enk5WVBUBWVhYjRoxoMv+QkBByc3Oprq5GRNiyZQvR0dEkJibajKGzGa0Aw4YNY/78+bbZqPbGbdOmTZw5c4aamhqys7MZMGAAkydPtuUbFBREeno6y5YtQ0TIzc3F19eXwMBAKisrOX/eWF2pvLycHTt2EBMT06h8Z/UNCQlhy5YtgOF2/eSTTwgLC7ssbWBgIJ06dSI3NxcRYdmyZbb0V9qOGo1LiIjDA3gLiHUQHgu86SxdSxx9+/YVjeZWoOfMt+QO69Hth38UQDz9e4pnQKjEx8fLunXrRESkvLxcUlJSJCIiQlJTU6WiokJERMrKysRkMomPj4/4+vqKyWSSs2fPiojIs88+K1FRUWI2m2XMmDFSW1vbqPyioiIxm82Nwqurq2XChAlisVgkJiZGhg8fLiIiS5YskREjRkhycrJERETIrFmzHNbr0qVLMmnSJAkLCxOLxSK7d+8WEZEdO3aIxWKRuLg4sVgssmjRIofpndW3tLRU0tLSxGKxiNlsltCUKYwAACAASURBVL///e+2NPHx8bb/d+/eLWazWcLCwmTy5Mly6dKlJvP9tgDskRv4rP62HE0tgLBbRPo5kR0UkdgWseQO0AsgaG4VBmRupbSqplH4bV7u7PrVvXRsp6cdaK4fegGElqGpMU6/JmTe11sRjebbwIxhUXh7Xr4lmbub4v+3d+/hUVVX48e/iyRAQEgACYThTiLNhXAnUBRD0oCCJRpAyiuF1lb8FSjVFqpvL1bb14pVy8VSq9UqvL5CCyhBAeWqCBgQBbmJQgAJIXIJCQpJwMD6/TEn45DMhATIDdbneeZxZp+z91lnHmRx9jmz15lz50n+63us/tSe+jSmpisrcW4RkftKNorIT4GPKi8kY65dd3Z38URqF1yhwQjgCg3mmZFdeWOC+6crP5mzhUmvfczxr63qijE1VVlTtS2AN4BzfJsoewF1gbu0Cqua2FStuR6cK7rA8+9l8OyafQTXDeC3Q6MY2bO1/YTCXDabqq0cfhOnZweRgUCs83GXqq6p9KhKsMRprif7jp3mN6/vYPPBk/SPaMaf7+pCu2YNqzssUwtZ4qwcZf0cJRFAVdcCb6nqs8VJU0RS/fUzxlyZiLAbmD++L4/fFcv2zFMMnrGO59/LoOj8heoOzRhD2fc4vUsvLCqx7XeVEIsxxlGnjnBPfDtW/vJWbolszhPL95AyewM7s05Vd2jGXPfKSpzi572vz8aYStAypD4v/LAnz93Tg2NfnyVl9gaeWPYpBefOV3doxly3ykqc6ue9r8/GmEoiItzeJZxVD97K3b1a8/y6/QyesY4N+05Ud2jGXJfKSpwdRWSJiLzp9b74c4cqis8Y4whpEMQTqXHMu68vAXWEe17cxNQFn5CX778QtTHm6ivr5yi3+tzgUNX3KiUiH+ypWmMuVvjNeWat3svz6/bTpEEQjw6LYWiXcPvpirmIPVVbOcoqZP1e8QvYDewu0WaMqSb1gwL49W3f4c1JN7tLk722lfvmbuGIj+X8SsrMzGTgwIFER0cTExPDzJkzPdtOnjxJcnIykZGRJCcnk5ubC8CePXvo168f9erV85QeKzZ9+nRiYmKIjY1l9OjRFBaWrjt68OBBgoODPRVVunXrxrlz/q+UX3nlFSZNmnTJc1FVJk+eTEREBHFxcXz88ceebQEBAZ5j+Vvg3t/5PvXUU56+sbGxBAQEcPLkyVL9Dxw4QHx8PBEREYwaNcpzTmfPnmXUqFFEREQQHx/PwYMHL3kuphbxt4gt7geA/gCcAE4CucBx4JGqXlDXFnk3xr9vis7rP9dl6Hd+t1yjf79c52w8oOfPX/C7/5EjR/Sjjz5SVdWvvvpKIyMjddeuXaqqOnXqVH3iiSdUVfWJJ57QX//616qqevToUd28ebP+5je/0aeeesoz1uHDh7V9+/aan5+vqqojR47Ul19+udQx/S0u78/LL7+sEydOvOR+S5cu1dtuu00vXLigH3zwgfbp08ezrWHDhpfs7+98vS1ZskQHDhzos//IkSN13rx5qqp6//3369///ndVVZ09e7bef//9qqo6b948vfvuuy8ZS2XAFnmvlFdZ9zgfBG4GeqtqU1VtAsQD/UXkwcpM5saY8gsMqMNPb+nIigcH0KNdEx5J28XI5z9g79Gvfe4fHh5Ojx49AGjUqBFRUVFkZWUBkJaWxrhx4wAYN24cixcvBiAsLIzevXsTFBRUaryioiIKCgooKioiPz+fVq1alTv2M2fOcO+999KnTx+6d+9OWlqaZ1tmZiYJCQlERkby2GOP+eyflpbG2LFjERH69u1LXl6ep/5mefg7X2/z5s1j9OjRpdpVlTVr1jBixIhS/b3HHTFiBKtXry6+IDHXgLIS5w+B0ap6oLhBVfcDY4CxlR2YMaZi2jRtwNx7+/DXu7uScfw0Q2a9z4xVn3O2yP9PVw4ePMjWrVuJj48H3DUvw8PDAWjZsiVHj5a96LzL5WLKlCm0bduW8PBwQkJCGDRokM99MzIyPNOfEydOBODxxx8nMTGRzZs3s3btWqZOncqZM2cA2Lx5M4sWLWL79u0sWLAAX885ZGVl0aZNG8/n1q1be/4RUFhYSK9evejbt6/PhFie883Pz+ftt99m+PDhnrYhQ4Zw5MgRcnJyCA0NJTAwsNSxveMKDAwkJCSEnJycMr9LU3uUVcMoSFVLPe+uqsdFpPQ/O40x1U5ESO3RmgE3NedPb+1mxqq9vLU9myeHdyHzZAFPvfMZR/IK3PdFb2nN05NHM2PGDBo3buxzrEs9bJSbm0taWhoHDhwgNDSUkSNH8uqrrzJmzJhS+3bq1Ilt27Zd1LZixQqWLFniuW9aWFjIoUOHAEhOTqZZs2YApKamsn79enr1Kv9zLl988QUul4v9+/eTmJhIly5d6NSpk9/9fZ3vm2++Sf/+/WnatKmnbdmyZYC7MLe5PpV1xVnWM+72/LsxNdiNN9Rj5g+68/KPe1Nw7jzDn/uAXy34hKy8AhQ4nPM1940dTVzCHaSmfruCZosWLTxTndnZ2YSFhZV5nFWrVtGhQweaN29OUFAQqampbNy4kU2bNnmuLpcsWeK3v6qyaNEitm3bxrZt2zh06BBRUVEApZKYiDB79mzPuEeOHMHlcpGZmenZ5/Dhw7hcLgDPfzt27EhCQgJbt24tdfxLne/8+fN9TtMCNGvWjLy8PIqKinweuziuoqIiTp065flHgKn9ykqcXUXkKx+vr4EqK2JtjLl8AzuHseLBATSsF8D5C+57bKpKzvKZ1GnSms+bX/yrs2HDhjFnzhwA5syZQ0pKSpnjt23blvT0dPLz81FVVq9eTVRUFPHx8Z5k6O+JVoDBgwfz7LPPeu7/eSe3lStXcvLkSQoKCli8eDH9+/dn4sSJnnFbtWrFsGHDmDt3LqpKeno6ISEhhIeHk5uby9mz7tJsJ06cYMOGDURHR5c6flnne+rUKd577z2/34GIMHDgQBYuXFiqv/e4CxcuJDEx0X4qdC2prKeOgH8Bx4CdXm1PAXuA7bhLloWWZyx7qtaYK9P+obe0nfNqcc+TCmhQ8/YaFNZBu3btqkuXLlVV1RMnTmhiYqJGRERoUlKS5uTkqKpqdna2ulwubdSokYaEhKjL5dJTp06pquojjzyinTt31piYGB0zZowWFhaWOr6/p2rz8/N1/PjxGhsbq9HR0Tp06FBVdT9Vm5KSogkJCRoREaGPPvqoz/O6cOGCTpgwQTt27KixsbH64Ycfqqrqhg0bNDY2VuPi4jQ2NlZffPFFn/39nW9xDKNGjSrV5/bbb9esrCxVVc3IyNDevXtrp06ddMSIEZ5zLygo0BEjRminTp20d+/empGR4fP4lQ17qrZSXpcsK3a5RGQAcBqYq6qxTtsgYI2qFonIk07ifuhSY9kCCMZcmf7T1pDl4zeeTRoE8fHvk+1q6BplCyBUjrKmaq+Iqq7D/ftP77YVqlrkfEwHWlfW8Y0x35o6uDPBQQEXtYlAbv43/HTOFrJPXXrhBGOMW6UlznK4F1jub6OIjBeRLSKy5fjx41UYljHXnju7u3gitQuu0GAEcIUG88yIrvxuaBQbMk4w6K/reG3TIfutoTHlUGlTtQAi0h53EezYEu2/BXoBqVqOAGyq1pjK80XOGR5etIMP9ufQr2Mzpg3vQrtmDas7LHMV2FRt5ajyK04R+RFwB3BPeZKmMaZytWvWkNfui+eJ1C7szDrF4BnrePH9/Z6ncI0xF6vSxCkitwG/Boapan5VHtsY45+IMLpPW1b+8lZujriR/1n6KanPbeSzL30v22fM9azSEqeIzAM+ADqLyGER+QnwN6ARsFJEtonIPyrr+MaYimsZUp9/ju3FrNHdyTyZzx3PupftO1d0obpDM6bGqNR7nFeL3eM0purlnD7LH9/aTdq2I3ynZSOeHB5H1zah1R2WqQC7x1k5qvOpWmNMDdbMWbbvxbG9yMv/hrv+voE/L/uUgnP+F4035npgidMYU6bvRbdgxS8HMKp3W15Yt5/bZ64jfb9V+jDXL0ucxphLalw/iCdSu/DaffEo8IMX0vntGzv4uvCb6g7NmCpnidMYU27f7XQjb/9iAPfd0oF5mw8xaPo61u45dsl+mZmZDBw4kOjoaGJiYpg5c6Zn28mTJ0lOTiYyMpLk5GRyc3MB2LNnD/369aNevXqesmPFpk+fTkxMDLGxsYwePZrCwsJSxzx48CDBwcGeairdunXj3Dn/hZ1eeeUVJk2adMlzUVUmT55MREQEcXFxfPzxx55tAQEBnmP5W9ze3/kCvPvuu3Tr1o2YmBhuvfVWn/0PHDhAfHw8ERERjBo1ynNOZ8+eZdSoUURERBAfH8/BgwcveS7m8ljiNMZUSHDdAH47NJrXJ/SnUf1AfvzKhzwwfysnz/hPSoGBgTzzzDPs3r2b9PR0Zs+eze7duwGYNm0aSUlJ7N27l6SkJKZNmwZA06ZNmTVrFlOmTLlorKysLGbNmsWWLVvYuXMn58+fZ/78+T6PW1wDtPhVt27dKz7/5cuXs3fvXvbu3csLL7zAz372M8+24OBgz7H8lVPzd755eXlMmDCBJUuWsGvXLhYsWOCz/0MPPcSDDz7Ivn37aNKkCS+99BIAL730Ek2aNGHfvn08+OCDPPTQJZcBN5fJEqcx5rJ0axPKWz+/hV8kRbJ0RzbJf32PNz854nPZvvDwcHr06AFAo0aNiIqKIisrC4C0tDTGjRsHwLhx41i8eDEAYWFh9O7dm6CgoFLjFRUVUVBQQFFREfn5+bRq1arccZ85c4Z7772XPn360L17d9LS0jzbMjMzSUhIIDIykscee8xn/7S0NMaOHYuI0LdvX/Ly8jw1PcvD3/m+9tprpKam0rZtW8/5l6SqrFmzhhEjRpTq7z3uiBEjWL16dbljMhVjidMYc9nqBtbhweSbePPnN9O6STA/n7eV8f/7EUe/Kj11WuzgwYNs3bqV+Ph4AI4ePUp4eDgALVu25OjRo2Ue0+VyMWXKFNq2bUt4eDghISEMGjTI574ZGRmeqdOJEycC8Pjjj5OYmMjmzZtZu3YtU6dO5cyZMwBs3ryZRYsWsX37dhYsWICvn8FlZWXRpk0bz+fWrVt7/hFQWFhIr1696Nu3ryehleTvfD///HNyc3NJSEigZ8+ezJ0719NnyJAhHDlyhJycHEJDQwkMDCx1bO+4AgMDCQkJAQgs88s0l8W+VGPMFftOy8a8PqE//1p/gKdXfMb3/voevxsaRd2AOjy94nOO5BXQKjSYSbe05unJo5kxYwaNGzcuNY6IXLLEWW5uLmlpaRw4cIDQ0FBGjhzJq6++ypgxY0rtWzxV623FihUsWbLEc9+0sLCQQ4cOAZCcnEyzZs0ASE1NZf369fTqVf6fQX7xxRe4XC72799PYmIiXbp0oVOnTn739z7foqIiPvroI1avXk1BQQH9+vWjb9++3HTTTSxbtgxwF+U21c+uOI0xV0VAHeG+AR1554EBRIc35qFFO/jVgk/IyitAgcM5X3Pf2NHEJdxBamqqp1+LFi08U53Z2dk+pyi9rVq1ig4dOtC8eXOCgoJITU1l48aNbNq0yXN16e/+IrinOxctWuS5F3no0CGioqIASiVtEWH27NmecY8cOYLL5SIzM9Ozz+HDh3G5XACe/3bs2JGEhAS2bt1a6vj+zrd169YMHjyYhg0bcuONNzJgwAA++eSTi/o2a9aMvLw8ioqKfB67OK6ioiJOnToFUIS56ixxGmOuqvY3NmTefX0JDQ6ieJ14VSVn+UzqNGnN580vflp02LBhzJkzB4A5c+aQkpJS5vht27YlPT2d/Px8VJXVq1cTFRVFfHy8Jxn6e6IVYPDgwTz77LOee7HeyW3lypWcPHmSgoICFi9eTP/+/Zk4caJn3FatWjFs2DDmzp2LqpKenk5ISAjh4eHk5uZy9uxZwH1luGHDBqKjo0sd39/5pqSksH79es99202bNnkSejERYeDAgSxcuLBUf+9xFy5cSGJiYpnfo7kCqlrjXz179lRjTO3S/qG3tJ3zanHPkwpoUPP2GhTWQbt27apLly5VVdUTJ05oYmKiRkREaFJSkubk5KiqanZ2trpcLm3UqJGGhISoy+XSU6dOqarqI488op07d9aYmBgdM2aMFhYWljr+gQMHNCYmplR7fn6+jh8/XmNjYzU6OlqHDh2qqqovv/yypqSkaEJCgkZEROijjz7q87wuXLigEyZM0I4dO2psbKx++OGHqqq6YcMGjY2N1bi4OI2NjdUXX3zRZ39/56uq+pe//EWjoqI0JiZGp0+f7mm//fbbNSsrS1VVMzIytHfv3tqpUycdMWKE59wLCgp0xIgR2qlTJ+3du7dmZGQosEVrwN/h19rL1qo1xlSK/tPWkJVXUKo9JDiQj38/iIA6Zd/LNFfO1qqtHDZVa4ypFFMHdyY4KOCitjoCpwqKuPv5D8g4frqaIjPmyljiNMZUiju7u3gitQuu0GAEcIUG88zIrswY1Y19x04zZOb7/HOdFcw2tY9N1Rpjqtyxrwr5zRs7WfXpUXq0DeWpkV3p1PyG6g7rmmNTtZXDrjiNMVUurHF9/jm2JzN/0I39J84wZOb7vLAuw64+Ta1gidMYUy1EhJRuLlY8OIABNzXnz8v2MOIfG9l3zO59mprNEqcxplqFNarPCz90X30eOHGGIbPe5/n37OrT1FyWOI0x1c776jPhpuY8sXwPw5/byL5jX1d3aMaUYonTGFNjhDWqz/PO1efBnDMMmbWef7yXQdH5C9UdmjEeljiNMTVK8dXnygdvZWDn5kxbvofh//jArj5NjWGJ0xhTIzVvVI9/jOnJs6O7c8i5+nzuXbv6NNXPEqcxpsYSEb7ftRUrHryVxM5hPPm2+97n3qNlX31mZmYycOBAoqOjiYmJYebMmZ5tJ0+eJDk5mcjISJKTk8nNzQVgz5499OvXj3r16nlKjhWbPn06MTExxMbGMnr0aAoLS9cbPXjwIMHBwZ5KKt26dePcuXN+Y3zllVeYNGnSJb8DVWXy5MlEREQQFxfHxx9/7NkWEBDgOZa/he1FpKmIrBSRvc5/mzjtCSJySkS2Oa9H/PTvICKbRGSfiPxbROo67fWcz/uc7e0veTLXCEucxpgar3mjejw3pof76vNkPkNnrefv7+7ze/UZGBjIM888w+7du0lPT2f27Nns3r0bgGnTppGUlMTevXtJSkpi2rRpADRt2pRZs2YxZcqUi8bKyspi1qxZbNmyhZ07d3L+/Hnmz5/v87jF9T+LX3Xr1r3ic1++fDl79+5l7969vPDCC/zsZz/zbAsODvYcq4xSag8Dq1U1EljtfC72vqp2c15/9NP/SWC6qkYAucBPnPafALlO+3Rnv+uCJU5jTK1QfPW58pe3khQVxl/e/ozU5zbyuY+rz/DwcHr06AFAo0aNiIqKIisrC4C0tDTGjRsHwLhx41i8eDEAYWFh9O7dm6CgoFLjFRUVUVBQ4Cn51apVq3LHfebMGe6991769OlD9+7dSUtL82zLzMwkISGByMhIHnvsMZ/909LSGDt2LCJC3759ycvL89TzLKcUYI7zfg5wZ3k7irtAaSKw0Ed/73EXAklyqSrk1whLnMaYWuXGG+rx3JiezP6vHhzOLeCOWeuZvdb/1efBgwfZunUr8fHxABw9epTw8HAAWrZsydGjR8s8nsvlYsqUKbRt25bw8HBCQkIYNGiQz30zMjI8U6cTJ04E4PHHHycxMZHNmzezdu1apk6dypkzZwDYvHkzixYtYvv27SxYsABfS4tmZWXRpk0bz+fWrVt7/hFQWFhIr1696Nu3r+cfAD60UNXiTPsl0MJrWz8R+URElotITHGjiCwTkVZAMyBPVYsLYh8GXMVfDZAJ4Gw/5ex/zQusrIFF5F/AHcAxVY112poC/wbaAweBu1U1t7JiMMZcu4bGhRPfsSmPpO3kqXc+4+2dXzKkS0teTT/EkbwCWoUGM+mW1jw9eTQzZsygcePGpcYQES51kZSbm0taWhoHDhwgNDSUkSNH8uqrrzJmzJhS+xZP1XpbsWIFS5Ys8dw3LSws5NChQwAkJyfTrJk716SmprJ+/Xp69Sr/0rJffPEFLpeL/fv3k5iYSJcuXejUqZPf/VVVRaR4ZYmPgXaqelpEhgCLgUhnvyEAInJjuYO5jlTmFecrwG0l2sqaazfGmAq58YZ6/P0e99Xn/uOnefLtz8jKK0CBwzlfc9/Y0cQl3EFqaqqnT4sWLTxTndnZ2YSFhZV5jFWrVtGhQweaN29OUFAQqampbNy4kU2bNnmuLsu4v4iqsmjRIs+9yEOHDhEVFQVQKmmLCLNnz/aMe+TIEVwuF5mZmZ59Dh8+jMvlvugr/m/Hjh1JSEhg69atvkI4KiLhzvjhwDEnrq9U9bTzfhkQ5CNR5gChIlJ8kdUayHLeZwFtnHEDgRBn/2tepSVOVV0HnCzRfNlz7cYY48/QuHAaBX97b1JVyVk+kzpNWvN581sv2nfYsGHMmeP+a2jOnDmkpKSUOXbbtm1JT08nPz8fVWX16tVERUURHx/vSYb+nmgFGDx4MM8++yzFlai8k9vKlSs5efIkBQUFLF68mP79+zNx4kTPuK1atWLYsGHMnTsXVSU9PZ2QkBDCw8PJzc3l7NmzAJw4cYINGzYQHR3tK4QlwDjn/TggDUBEWhbfkxSRPrjzwUWJT91BrwVGlOxfYtwRwBqtDeW2rgZVrbQX7inZnV6f87zei/dnH33HA1uALW3btlVjjClL+4fe0nbOq8U9TyqgQc3ba1BYB+3atasuXbpUVVVPnDihiYmJGhERoUlJSZqTk6OqqtnZ2epyubRRo0YaEhKiLpdLT506paqqjzzyiHbu3FljYmJ0zJgxWlhYWOr4Bw4c0JiYmFLt+fn5On78eI2NjdXo6GgdOnSoqqq+/PLLmpKSogkJCRoREaGPPvqoz/O6cOGCTpgwQTt27KixsbH64Ycfqqrqhg0bNDY2VuPi4jQ2NlZffPHFUn2dv0Ob4Z7h2wusApq6NzEJ2AV8AqQD39Vv//5dBrRy3ncENgP7gAVAPae9vvN5n7O9o15hzqgtr0qtx+n8ruct/fYeZ56qhnptz1XVJpcax+pxGmMupf+0NWTlFZRqrxtQh9W/upU2TRtUQ1TVy+pxVo6qfqrW51y7McZcqamDOxMcFHBRW1CAICi3zVjH/M2HqMwLBXP9qOrE6XOu3RhjrtSd3V08kdoFV2gwArhCg3lqRFdWT0mgS+sQHn59Bz+Zs4VjX5Ve9ceYiqi0qVoRmQckADcCR4E/4H7c+T9AW+AL3D9HKfkAUSk2VWuMuRIXLiivbDzIk2/vIbhuAP9zZyx3xJV/EYPayqZqK0el3uO8WixxGmOuhn3HTvOr/2zjk8On+H7XVvwpJYbQBle+LF5NZYmzctjKQcaY60ZE2A0s+tl3+VXyTSzfkc2g6etY+5k9amEqxhKnMea6EhhQh58nRbJ4Yn9CGwTx45c/5L9f387ps0WX7mwMljiNMdepWFcISybdzP0DOjL/w0xun7mOzQcu+ciFMZY4jTHXr/pBAfz3kCj+Pb4fgjDqhQ94fOluCr85X92hmRrMEqcx5rrXp0NTlv/iFkb3acs/3z/A959dz47Dpy7ZrzYUzBaRH4nI3y51LuI2yylMvV1EenhtO+9V8Nrnwrz+CmZ7be8tIkUiMsJP/54issM5/iyv5QDLHLc6WOI0xhigYb1A/nxXF175cW++KvyGu/6+gZmr9vKNn3JlcG0VzAZux10dJRL3kqfPeW0r0G8LXvtbmNdvEQ8RCcBd6HpFGcd/DrjPK4biIiE1rjiIJU5jjPGS0DmMdx4YwJAu4Uxf9TnDn9vIvmOli2VD7SmYDbQRkXedq7Y/+BkiBZjrLHObjrsqSni5Ayi7iMfPgUX4WS3OOU5jVU131tidi++C2TWiOIglTmOMKSG0QV1mje7O7P/qQebJfIbOWs9L6w9w4YL/373XxILZfPt3fB9gOBAHjBQRX7/t9BSmdngXra4vIltEJF1E/CUunwWzRcQF3MXFV7A424qLl7qc4/k6dlmFuKtFpRWyNsaY2m5oXDi9OzThvxft4E9v7Wbl7i95akTXUgvGnz59muHDh9e4gtlA8RzuSlXNcWJ5HbgZd+WU8mqnqlki0hFYIyI7VDXD386qFxXMngE8pKoXSn4HqtqtAjGUHLfaWOI0xpgyhDWqz4vjerFgy2Eee3MXt898n6FxLVm/9wRH8gpp2SiIwrf+zA/vucdnwezw8PAKF8wGPAWzIyMjuf/++wH44x//SFxcnM/+6hTM7ty5s6dNRIqfLiqZbFREJuK+pwgwBK/C1A5P0WpVLf7vfhF5F+gOlEycR0UkXFWzSxTx6AXMd5LmjcAQESlS1cVefbOc45U6dhnjVhubqjXGmEsQEe7u3Ya3HxhAi0b1+PeHh8nKK+SCKjvmTeOLC6F0HDjqoj41qWA2kOw8nRqM+x7hBlWd7fXAzxHcRTjGOk/X9gVOOcmqiYjUc76HG4H+wG4fIfgs4qGqHVS1vaq2BxYCE0okTZyp2K9EpK/zNO1YfBfMrhHFQeyK0xhjyqlN0wYUFH37G8+zWbs5s2st55q3556htxIZdgN//vOfGTJkCA8//DB33303L730Eu3ateM///kPAF9++SW9evXiq6++ok6dOsyYMYPdu3cTHx/PiBEj6NGjB4GBgXTv3p3x48eXO7bf//73PPDAA8TFxXHhwgU6dOjgvXkz7odzWgOvqqqvadpluK889wH5wI+d9ijgeRG5gPtia5qq+kqc04D/iMhPcIp4XCpmEdnmNV07AXgFCAaWO6/LGrey2SLvxhhTAR0eXlpq3hNAgAPThlZ1OGWyRd4rh03VGmNMBbQKDfbZXqeO8ElmXhVHY6qDJU5jnwqiewAADoxJREFUjKmAqYM7ExwUcFFb3cA63FAvgOHPbWT22n2cL+NnK6b2s8RpjDEVcGd3F0+kdsEVGowArtBg/jI8jnVTE7kttiVPvfMZo19I53BufnWHaiqJ3eM0xpirRFV5Y2sWj6TtQoA/3RnLnd1dl+xXWeweZ+WwK05jjLlKRITUHq1Z/otbuKllIx749zZ+MX8rpwq+qe7QzFVkidMYY66yNk0b8O/xffll8k28tT2bITPfZ9P+nOoOy1wlljiNMaYSBAbUYXJSJAv/Xz8CA4Qf/DOdv7y9h3NF/qutmNrBEqcxxlSi7m2bsGzyLdzdsw1/fzeD4c9tJOP46eoOy1wBS5zGGFPJGtYL5MkRcfxjTA8yc/O5Y9Z6Xtt0iNrwcKYpzRKnMcZUkdtiw3n7FwPo2a4Jv3ljB+P/9yNyTp8ts09mZiYDBw4kOjqamJgYZs6c6dl28uRJkpOTiYyMJDk5mdzcXAD27NlDv379AHqIyEUVs0XkQRHZJSI7RWSeiNQveUwRaS8iBSKyzevlt1q2iPxIRP52qfN31sGdJSL7RGS7iPTw2nbe61hL/PRvKiIrnbqiK0WkidOe4oy3zSl/drOf/j1FZIdz/FnOurh+x/XHEqcxxlShliH1mXtvH343NIr3PjvObTPf573Pj/vdPzAwkGeeeYbdu3eTnp7O7Nmz2b3bvVTstGnTSEpKYu/evSQlJTFt2jQAmjZtyqxZswAuKgLq1MacDPRS1VggAPiBn0NneC0C301Vz13puQO3A5HOazwX1+gs8DqWv9XsHwZWq2oksNr5jPO+q7Pu7b3Ai376P4e7IkxxDLddYlyfLHEaY0wVq1NH+OktHUmb1J8mDYIY96/NPLpkF4XfnC+1b3h4OD16uC/MGjVqRFRUFFlZ7opbaWlpjBvnLhwybtw4Fi92Fx0JCwujd+/eULqcGLiLewSLSCDQADhS3rhFpKGI/EtENovIVhHxLvnSRkTeda7a/uBniBRgrrqlA6FOqbDySgHmOO/n4K70gqqe1m/nvRvi47yd4zRW1XRn37nF/f2N60+1JM7yTBUYY8y1Liq8MUsm3cyPvtueVzYeZNjf1vNp9ld+9z948CBbt24lPj4egKNHjxIe7s47LVu25OjRo377gqeu5tPAISAbd+mwFX527+Q1dTrbafstsEZV+wADgadEpKGzrQ8wHIgDRoqIr4UXXECm1+fDThtAfWeaNV1E/CWuFk4JMoAvgRbFG0TkLhHZAyzFfdVZ3F5c8dvlHM/Xsf2O60uVlxXzmiqIVtUCEfkP7qmCV6o6FmOMqW71gwJ4dFgMA78TxpQFn5Dytw38+rbONGtQl6dXfs6RvAJahQYz6ZbWPD15NDNmzKBx48alxhERnFt2fjn37lKADkAesEBExqjqqz52z/Aq+VVsEDDM675pfaCt836lquY4x3kduBmoyJJv7VQ1S0Q6AmtEZIeqliyW7aGqKiLq9fkN4A0RGQD8Cfie017yHMpUclxfqmuq9rKnCowx5lp0603NefsXt3Br5+b8z9JP+dXCT8jKK0CBwzlfc9/Y0cQl3EFqaqqnT4sWLcjOdl8oZWdnExYWdqnDfA84oKrHVfUb4HXguyIS73V16b9atrt62nCve5FtVfVTZ1vJZKMiMtFr3FZAFtDGa5/WTlvx1TCquh94F+ju4/hHi6d2nf8eK7mDqq4DOjpFt71lOccrdezyjOutyhNneacKRGS8c9m+5fhx/zfOjTHmWtHshnq88MOehAYHUVxgRVXJWT6TOk1a83nzWy/af9iwYcyZ4741N2fOHFJSUkoOWdIhoK+INHCeKE0CPlXVTV7J0OcTrY53gJ97PY3qndySnadTg3HfI9ygqrO9xj0CLAHGOk/X9sX993+2iDQRkXrOmDcC/QFfxbKXAOOc9+OANKdPhFdMPYB6wEVLNTlTsV+JSF9n37HF/f2N60+VJ84SUwWtgIYiMqbkfqr6gqr2UtVezZs3r+owjTGmWojIRWvbns3azZldayk8tJ0Pp/+Ubt26sWzZMgAefvhhVq5cSWRkJKtWreLhh90Pg3755Ze0bt0a3Pfqficih0WksapuAhYCHwM7cOeAFyoQ3p+AIGC7iOxyPhfbDCwCtgOLVNXXNO0yYD+wD/gnMMFpjwK2iMgnwFpgmqr6SpzTcCfovbivnqc57cOBnc79zNnAqOKHhbzuceIc70Xn+BnA8kuM61OVV0cRkZHAbar6E+fzWKCvqk7w18eqoxhjrif9p60hK6+gVLsrNJgNDyeWexyx6iiVojrucfqcKqiGOIwxpkbyVSw7OCiAqYM7V1NExluVP1WrqptEpHiqoAjYSsWmCowx5ppWXMPzqXc+8zxVO3Vw52qt7Wm+ZYWsjTHmGmVTtZXDVg4yxhhjKsASpzHGGFMBljiNMcaYCrDEaYwxxlSAJU5jjDGmAmrFU7Uichz44ioMdSNw4iqMU5Us5qphMVcNi7lqFMfcTlVt6bWrrFYkzqtFRLbUtkezLeaqYTFXDYu5atTGmGsTm6o1xhhjKsASpzHGGFMB11virI1L+1nMVcNirhoWc9WojTHXGtfVPU5jjDHmSl1vV5zGGGPMFbHEaYwxxlTANZM4ReRfInJMRHZ6tT0lIntEZLuIvCEioU57exEpEJFtzusfNT1mZ1uciHwgIrtEZIeI1K/JMYvIPV7f8TYRuSAi3Wp4zEEiMsf5fj8Vkf+u6ngvI+a6IvKyE/MnIpJQHTGXEfefnJi3icgKEWnltIuIzBKRfc72HrUg5u84/w+eFZEp1RHvZcR8j9O+Q0Q2ikjX6or7mqGq18QLGAD0AHZ6tQ0CAp33TwJPOu/be+9XS2IOBLYDXZ3PzYCAmhxziX5dgIxa8D3/FzDfed8AOAi0r+ExTwRedt6HAR8BdWrQd93Y6/1k4B/O+yHAckCAvsCmWhBzGNAbeByYUh3xXkbM3wWaOO9vr67v+Vp6XTNXnKq6DjhZom2FqhY5H9OB1lUeWBkqGPMgYLuqfuLsl6Oq56ss2G/ju9zveTQwv5LD86mCMSvQUEQCgWDgHPBVVcXqFV9FYo4G1jj7HAPygGr58bufuL2/v4a4v2OAFGCuuqUDoSISXjWRXhRfuWNW1WOq+iHwTdVFWFoFY96oqrlOe437e7A2umYSZznci/tft8U6iMhWEXlPRG6prqAuwTvmmwAVkXdE5GMR+XU1xlWWkt9zsVHAvCqOpby8Y14InAGygUPA06p60l/HauQd8yfAMBEJFJEOQE+gTbVF5oOIPC4imcA9wCNOswvI9NrtsNNWI/iJuUYrR8w/wff/n6YCrovEKSK/BYqA/3OasoG2qtod+CXwmog0rq74fPERcyBwM+7/IW4G7hKRpGoKzycfMRe3xwP5qrrTZ8dq5CPmPsB5oBXQAfiViHSspvB88hHzv3AnnS3ADGAj7nOoMVT1t6raBnfMk6o7nvK41mIWkYG4E+dD1RHbteSaT5wi8iPgDuAedSb5VfWsquY47z8CMnBf0dUIvmLG/RfjOlU9oar5wDLc9zhqBD8xF/sBNfBq00/M/wW8rarfONOeG6imaU9f/Px5LlLVB1W1m6qmAKHA59UYZln+DxjuvM/i4ivj1k5bTeMdc21xUcwiEge8CKQU/91nLt81nThF5Dbg18AwJ9kUtzcXkQDnfUcgEthfPVFezF/MwDtAFxFp4Nx/uxXYXR0xllRGzIhIHeBuqun+pj9lxHwISHT2aYj7oZU9VR9haWX8eW7gxIqIJANFqloj/mwAiEik18cUvv0+lwBjnadr+wKnVDW7ygP0oYyYayx/MYtIW+B14IeqWlP/QVW7VPfTSVfrhfuKJhv3TfvDuKck9uG+h7LNeRU/ZTYc2OW0fQx8v6bH7Ow/xol7J/CXWhJzApBei/5s3AAscL7n3cDUWhBze+Az4FNgFe5SUjXpu17k/JndDrwJuJx9BZiNe8ZnB9CrFsTc0tnnK9wPYR3G62nWGhrzi0Cu15+bLdX15+NaedmSe8YYY0wFXNNTtcYYY8zVZonTGGOMqQBLnMYYY0wFWOI0xhhjKsASpzHGGFMBljjNdUVEzjvVI3aKyJviVX3G2f6AiBSKSEgZY4SLyFt+tr0rIpe1YIKI3CEif7ycvsaYqmOJ01xvCtS9wk4s7kWyJ5bYPhr4EEgtY4xfAv+shNiWAt8XkQaVMLYx5iqxxGmuZx/gtai4iHTCvQDC73AnUH+GA287fYJFZL64a3e+gbuiSvF4g5zajR+LyAIRucFpHyLuupofOfUo3wJQ94+q38W9pJ4xpoayxGmuS86Si0m4l30r9gPcSwO+D3QWkRY++nUAclX1rNP0M9wL2EcBf8BdmQQRuRF3Av6eqvbAvQD7L8VdfPx54HZV7Qk0L3GILUBNrdZjjMESp7n+BIvINuBLoAWw0mvbaNxFrC/gXr5spI/+4cBxr88DgFcBVHU77uXOwL3GbTSwwTneOKAd8B1gv6oecPYrufj9MdyVWYwxNVRgdQdgTBUrUNVuzn3Ed3Df45wlIl1wL/a/UkQA6gIHgL+V7A/UL8dxBFipqhdN+YpIt0v0q+8cwxhTQ9kVp7kuqbu6yGTc9TYDcV9tPqqq7Z1XK6CViLQr0fVz3IuqF1uHuxQZIhILxDnt6UB/EYlwtjUUkZtwL8beUUSKxxhVYvybcC/UbYypoSxxmuuWqm7FPbU6Gvf9zTdK7PKG0+7d5wyQUZwQgeeAG0TkU+CPwEfOfseBHwHzRGQ77geRvqOqBcAE4G0R+Qj4GjjldYiBuJ+uNcbUUFYdxZgKEpG7gJ6q+rvL7H+Dqp4W95zwbGCvqk53HkZ6TVWTrma8xpiry644jakgVX0DOHgFQ9znPDC0CwjB/ZQtQFvgV1cWnTGmstkVpzHGGFMBdsVpjDHGVIAlTmOMMaYCLHEaY4wxFWCJ0xhjjKkAS5zGGGNMBfx/+VFuzqiL4WQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "\n", "ax.plot(eph['RA'], eph['DEC'], 'o-')\n", "ax.set_title('trajectory of 2018 CL on from %s to %s' % (start_date, end_date))\n", "ax.set_xlabel('RA (deg)')\n", "ax.set_ylabel('DEC (deg)')\n", "\n", "for i, xy in enumerate(zip(eph['RA'], eph['DEC'])):\n", " ax.annotate('%s' % eph['datetime_str'][i], xy=xy, textcoords='data')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What about the change of azimuth and elevation angle?" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdIAAAEWCAYAAADSGRaUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4VuXZwH93BiGMBAgrJETIAEJCCAgERDSADFFBI4goFSe2gta2oNS2Vq0WLCoiYr9WLeACFRQQUFFwgoAsUZasBJJAgCwgA0hyf3+ck5c3iwzGm/H8rutc73vOs+5n3s86zxFVxWAwGAwGQ9Vwc7UABoPBYDDUZIwiNRgMBoPhAjCK1GAwGAyGC8AoUoPBYDAYLgCjSA0Gg8FguACMIjUYDAaD4QK4JIpURIJE5JSIuF8K/0sJb66IPFvd/KpkuE+IyBuXO9wL4VKllYi0EpFvReSkiLx4sf2v7YiIt4h8IiKZIvKhq+UxGGo7l0SRqupBVW2kqvkX228RuVtEvq+Cuz4isvZiy1MVRCRWRBKdn6nqP1X1flfJVM0YDxwHfFT1T8UNRWSyiPxiK9oDIjK5mHk7EflKRLJFZJeIXOdkFikin4vIcREp8RK17XaFiKSLyBEReVVEPMoSVER8RORlETlodx732ffNbfN45/AvEyOBVoCfqo66zGEjIuNEZJOInBCRRBH5l3MaikgzEflYRLJEJEFE7nAy8xeRpSKSLCIqIu2K+d1MRN4XkVQ7D98VER8n8/Pl/YXIFSsiBXYeF17jzpMGKiKhVUy/DiKyRESOiUiaXV47FrPzB7t8nhCR/4mIl5PZP0TkZxHJE5GnSvH/YbvenBCRjSJytZOZiMjzdvqm2v/lYshlm//eDjtLRHaKSIeqpFF1oy5N7d4ArHC1EIYKcQWwQ8s+LUSAu4CmwFBgoojc7mQ+H9gC+AF/ARaKSAvb7CzwAXBfGX6/BhwF/IFo4FrgoVKFEKkHrAIibDl8gD5AKtCr3FheOq4AflXVvNIMz9cxuEg0AB4FmgMxwEBgkpP5bOAMlrK/E/i3iETYZgXAZ8CtZfj9LFa+twdCbD+ecjI/X95fiFwAyfYAofCaV25KVI0mwFKgoy3LBmBJoaGIDAGm2PJfAQQDTzu53ws8Biwv7rGIxADTsDpbvsCbwMdybvZwPHAz0BWIAm4CHrwYconI/Vj17gagEXAjVoe55qOq572wEmYfcBLYAdziZPYTcMrpUiAWaGf/97DtfY1VAdba9j7BKujvAieAH4F2tt0ibp3c3w+EA7lAvu1Phm0+F6sSLLflXA+EFIvHZqC7/b8T8AWQBuwGbnOyNxd41un+RmArkGHLH2U/fxxYWCyMmcAr9v97gJ22PPuBB+3nDYEcrAajMN3aYDUG7zj5NRzYbof7NRDuZBaP1QBsAzKB94H6ZeRfCLAaq3E/bqd5k4r6hVUhDwPJdh4oEFqZtCpDrqvsfM+0f69y8vMsVoN2CriuAmX0FWCW/b8DcBpo7GT+HfDbYm5CAS3Fr53AMKf76cB/ygj3fiAFaHQe2eIrGAc34K9AApYifwvwLVYnxgEH7Xz8Sxn+PG2n3Vk7/e4D7gbWADPscvBsBcO7BzgEpAO/BXra5SQDeLW8ODnJ9EfgE6fyfwbo4GT+NjCtmBsPW4Z2xZ5/CjzkdD8B+LwyeV8VubDatcQKxvdbW/YsOw9G288fwFJyaVgKqU0F/Wtm++dn378H/NPJfCBwpBR37wBPFXs2GtjgdN/Q9tvfvl8LjHcyvw9Yd6Fy2eXtEDCwouWmJl0VycRRWA29m50JWYWJXszeeGAXVq+8sCI6K9K9WI26L5ZC/hW4zq4wbwFzbLtF3Dq5v9/+fzfwfbGw53JuFOCBpSwWOJn7A0lYI5mGdobeY9vthtUwdXby61n7fzesRiYGcMdqyOIBL6weVzZ2pbXNDwO97fsb7PgK1qgmm3OKPJZilRInRYrVIGQBgwBPLGW2F6hnm8dj9Qbb2IV5J2U3FqG2P15AC6xK/rKTeZl+YY2yjmCNuBpgVcxSFen50qqMCpgO/MbOgzH2vV9xfytQPgVrBFIo8y3AzmJ2XsVWtMXSRUvx70Gs8tgACAB+wanzWMzuAmBeOfLFUzFFeq+dx8FYvfWPgLeL1YnXAW+s0cJpnDpXZZUlpzqTBzxsp7d3BcP7P6A+MBirA7sYaGmny1Hg2grm0WLOKaRuQHYx80nYCs3pWVmK9EasmaWm9rUaeLQyeV8VubDq7BmsjtMBrE5Jw/PE2VFP7PsBWO1Md6y6OAv4toLpdzNw2On+J2zlbN83x0mhOT0vTZH6AJs4V08fxqo/YptnAjFO9nsAJy9ULiDI/v97rPb3AFanz60iaVDdr3KndlX1Q1VNVtUCVX0f2EOxaSt7jv1ZYLiqnijDqzmquk9VM7F6lftU9Uu1pp8+xCrIF8LHqrrB9u9drGm5QoYBn6mVuzcC8ao6R1XzVHULsAirw1Cc8VijkfWqmq/WVM5pLGWZgDXKvcW2OwCrIq4DUNXldnxVVb8BVgL9KhiX0cByVf1CVc8CL2A1flc52XnFzpc0rBF+dCn+oKp7bX9Oq+ox4CUsxe5MWX7dhpVv21U1m6JTaMUpM61KsXsDsEdV37bzYD5WJ+ym8/hfFk9hdfLm2PeNsBoDZzKBxhX071usjsMJIBHYiNXgloYfVufpYnAn8JKq7lfVU8CfgduLTcM+rao5qvoTVqPVtRL+J6vqLDu9cyoY3j9UNVdVV2J17Oar6lFVTcIa6ZVbZ0XkXqzG+AX7USOstHWmMvmzGaiH1XFOxZqdes3J7wrlfRXk2oVVL/yx6vqVWHWpotwJ/E9VN6vqaaz07lN8DbgUOQOxZtv+6PS4eDwL/1ckDU9itXffY9XPv2ONQPU8fjcqXCe9ALkC7f+DgS5Af6wOdFlLLDWKchWpiNwlIltFJENEMoBIrJ5GoXlbrDWncar663m8SnH6n1PKfaNKSV6SI07/s4v5N4xz66NXADGF8bHjdCfQuhQ/rwD+VMxuW6zRG1hTGWPs/3fY9wCIyPUiss5elM+wZWhOxWiDNeUGgKoWYPXiAioYXwdi7YBdICJJInICq5daXI6y/Gpjh1uI8//ilJdWzhSJn00CReNXLiIyEWut9Aa7cQJrKs2nmFUfrAakPP/csNbnPsKauWiONep5vgwnqVgN68WgeJokYI3KWjk9q1Cel0HxvKtIeBdUZ0XkZmAqcL2qFq6FVTl/bD7Ams1qbLvbh1WmK+x3VeRS1SOqusMeUBzAmiUqax23NIrX6VNY5afMMm+v7a4EXrM7m4UUl7Xwf0XS8D6s2bgIrA7JWGCZiBTW09L8PuWkaKsqV479/1+qmqGq8cB/sNrFGs95FamIXIE1nTQRa9qgCdZUV+EuLm+s3vrLqvrpRZIpy/5t4PTMWckplUBEPLFGYF/Yjw4B36hqE6erkar+rhTnh4Dnitlt4FR4PgRi7d7ZLdiK1N6ptgirt9vKTrcV2OlWgTgkYymmwjgIllJKqkzcbf5ph9dFVX2wKo6c34mDw5zrSWLLUBblpZUzReJnE0Ql4mePKKZgrbk474DeDgSLiHPvvKv9vDya2XK8ao/gU7FGumVV9i+BISLSsKJyn4fiaRKENR2bUrr1SlO8zF3S8ERkKFbbcZOq/uxk9CvgISJhTs8qmj9gjQr/o6pZtjL6P87lT7l5fxHlUiq3WbN4nW6INaNRapkXkaZYymqpqj5XzHg7RWcjugIpdnktj2hgmar+ancKPsOq54WzXaX57Zx+VZVrN9bUuHM5rFRbXp0pryAULkQfAxCRe7BGpIX8D9ilqv+6WALZ049JwFgRcbcbzBAnKylAoL1jsiJcDWxzmnJeBnQQkd+IiKd99RSR8FLcvg78VkRi7G3hDUXkhsKKasv6NVZje0BVd9ru6mGtgxwD8kTkeqwpDec4+ImIbxkyfwDcICID7Y7An7CmYary+k5jrJ5ipogEAJPLsV9cjntEJFxEGgB/O4/d86ZVMVZg5cEdIuIhIqOBzlh5Uy4icidWB2GQqu53NrNnRbYCfxeR+iJyC9buw0W2WxGR+lh5hG3Hy3Z7HGvt5ne2XE2w1nq3lSHK21gdiEUi0klE3ETET6x3gp2Vr6cdTuFV2q7Z+cAfRKS9iDSy4/e+lrHz9iJwycITkQFYyyu3quoGZzNVzcIa8T9jl5G+wAistCx0Xx+r/gB42feF/AjcL9a7st5YSwrbbL/Ly/sqyyUi/UXkCrv8tMXa+bqEsknBWn8uZD5WXYq2y9s/gfX2yKx4+vkAnwNrVHVKKX6/BdwnIp3tMvpXrH0Fhe497TRzw+oc1Jdzu3J/xGpbgu24DMLak/GLk99/FJEAe5T6p0K/L0Que2nofeAxEWlsDz7GU8E6X+3R8he6n8PaZXYca03gG85t/FGsKSbnnbv9KH2z0f1Ofj4LzHW6vw7Y63R/PVaDlgG8WCzMeli7c9OA4/azuRTdPRqLvZkHa1Q4qVicOtp+HMOaXlkNRJfh11CswpeB1XP7kKK7An9jx3VysTAmYFWmDKzKuKCYv/+zw86g9F27t2Btysq04x/hZBaP0waW4m6LyRGBtbngFFYj8yecNjqV5xfWWs4RrB717+y4tq1KWhWT62pbrkz792onsyL+luL2AOd2pRZe/+dk3g6rzOVg9YSvK2amxa54J/No2206Vpn/AGtWoSxZfIGXsRTqKaypxpc4t3EqvpTwSsQNq9F70vbnGNZ0ZdNiMpe6Aa8Uv4rn4d2U3KBX2fASgVin+3eAv5YR/ldYo1vn/PnUybwZ1kxWFtYu5DuKuS+eXupk1h5rHT8Vqw34DAirYN5XWS6stcAkrPbuENZO8VLLtm3/t1h1IAP7rQD72T5b7mVAYBlux1F012/hFVRMnhSsdd05OG3qw6o/xdPwbttMgGfs+J3E2lz4Gye3AvzLljHN/i8XSS4frHbwpJ2GTxb6XdOvwgSqtYjIDmCkqu5wtSw1HXvU/gtW5bhUIyWDwWCoUdTqAxns6d+3jBKtOiJyi4h42Wsjz2O9DmCUqMFgMNjU+hGp4cIQkc+wTuvJx5pifkhVL9YrHwaDwVDjMYrUYDAYDIYLoFZP7RoMBoPBcKm51IdXOxDrKwHvOz0Kxtq19Zb9vB3WDsfbVDX9fH41b95c27Vrd0nkNBgMhtrKpk2bjqtqi/JtGiqDS6Z27XeakrDOe5wApKnqNBGZgrUF//Hzue/Ro4du3LjxMkhqMBgMtQcR2aSqPVwtR23DVVO7A7HO2k3Aeum58HNE87AOQjYYDAaDoUbgKkV6O9ZJH2C97F64C/QIRc/6dCAi48X6CO3GY8eOXQ4ZDdWIQ4cO0b9/fzp37kxERAQzZ850mKWlpTFo0CDCwsIYNGgQ6enWysCuXbvo06cPXl5evPDCC0X8mzFjBhEREURGRjJmzBhyc3NLhBkfH4+3tzfR0dGO68yZM2XKOHfuXCZOnFhuXFSVRx55hNDQUKKioti8ebPDzN3d3RHW8OHDS3X/888/06xZM+rVq0ejRo2YOnUqANOnTycyMpLGjRvj5eWFiLB///4SafHEE08QExNDaGgoo0ePZvr06URERBAREUFQUBDBwcHExMQQHx9f7dOirLzPzMzkpptuomvXrkRERDBnzpxS3W/atIkuXboQGhrKI488UnhwQJn+GgylcrlPgMA6meg49mkx2N8UdTJPL8+PK6+8Ug11i+TkZN20aZOqqp44cULDwsJ0+/btqqo6efJknTp1qqqqTp06VR977DFVVU1JSdENGzboE088odOnT3f4lZiYqO3atdPs7GxVVR01apTOmTOnRJgHDhzQiIiICss4Z84cnTBhQrn2li9frkOHDtWCggL94YcftFevXg6zhg0bluv+d7/7nU6cOFFVVZ966ilt2rRpibRYunSpBgcHl5oWUVFROn/+fFVVHTt2rPr5+Wl2drbOnj1bg4ODdc6cOTp//ny97bbbqn1alJX3zz33nOP/0aNHtWnTpnr69OkS7nv27Kk//PCDFhQU6NChQ3XFihXn9bemA2zUanASUG27XDEivR7YrKqFh2OniIg/gP171AUyGao5/v7+dO/eHYDGjRsTHh5OUpJ13veSJUsYN24cAOPGjWPxYuurZy1btqRnz554enqW8C8vL4+cnBzy8vLIzs6mTZvSPlJTOllZWdx777306tWLbt26sWTJuSNXDx06RGxsLGFhYTz99NOlul+yZAl33XUXIkLv3r3JyMjg8OGKv5q7atUqnnjiCQDGjx/PmTNnSqTF/Pnz+e1vf1siLTw8PNi7dy8jR44EYNSoUZw6dYqcnBw+/vhjWrduTZs2bRg5ciSrVq1yjNCqa1qUlfciwsmTJ1FVTp06RbNmzfDwKLq38vDhw5w4cYLevXsjItx1110O92X5azCUhisU6RjOTeuC9aX4cfb/cZz/IGiDgfj4eLZs2UJMTAwAKSkp+PtbXzNr3bo1KSnn/4BJQEAAkyZNIigoCH9/f3x9fRk8eHCpdvft2+eYXpwwYQIAzz33HAMGDGDDhg189dVXTJ48maws66NFGzZsYNGiRWzbto0PP/yQ0jbFJSUl0bbtuQ/pBAYGOhRhbm4uPXr0oHfv3mU23s7xzc3NJTs7u0ha+Pr68tlnn3Hvvfc60mLYsGEkJyeTk5ODt7e3Q6l069aNpk2bEhQUxFdffUXLli0ZPHgwHh4e+Pr6kpp67oMi1T0tnPN+4sSJ7Ny5kzZt2tClSxdmzpyJm5vV3EVHRzvCDgw893Ej57ArW6YMdZvL9voLOD4dNAh40OnxNOADEbkP63t9t11OmQzVm8Vbkpj++W6SM3Jo08Sbif0CeeGRMbz88sv4+BT/fKQ1EhE5/1fi0tPTWbJkCQcOHKBJkyaMGjWKd955h7Fjx5awGxISwtatW4s8W7lyJUuXLnWsu+bm5nLw4EEABg0ahJ+fHwBxcXF8//339OhR8U2SCQkJBAQEsH//fgYMGECXLl0ICQkpkhYncvPoO221Iy28vb2LpMUnn3xC37598fPzc6TFihUrSg0vIyODkydPcuDAAWJjY8nOzq4RaVEaznn/+eefEx0dzerVq9m3bx+DBg2iX79++Pj4lIhDeVSkTBnqNpd1RKrWNwT9VDXT6Vmqqg5U1TBVvU5V0y6nTIbqy+ItSfz5o59JyshBgcTUkzxw1xiiYm8kLi7OYa9Vq1aO6cDDhw/TsmXL8/r75Zdf0r59e1q0aIGnpydxcXGsXbuW9evXO0ZcS5cuLdO9qrJo0SK2bt3K1q1bOXjwIOHh1lf4ije4IsLs2bMd/iYnJxMQEMChQ+e+s52YmEhAgPV958Lf4OBgYmNj2bJlS4m0cG/YhPiEBB64awyhPa4tMqpq1aoVc+fOZcyYMaWmhbe3t2NKG2Dp0qU0adKEFi1aEBgYSExMDGvXrmXNmjUkJCRw3XXXVbu0cKasvJ8zZw5xcXGICKGhobRv355du3YVcRsQEEBi4rlP2TqHXdkyZajbXNYRqcFQGaZ/vpucs/mA1WCnfjoTt6aBbGx0FS+t3O2w1ybqau7/6wsMHD2eVe//lzZd+/HSF786zH/Yl0o97xzc7GcJh2H5qm+Z9slPeHrVZ8HbH9O2QyRrTzRl3PQPAOtbVz9+t5/UU2eY4eSXCDTr0JP7H/sHcRP/hoiQuHcHgaGd2bD9CMuXfcpzH23A06s+b7zzAWP+9E+COnbhnhesqeMPd2ShQVfyzIx/k9L8ShJ2/kSuePHB9lPkrN9IPS9vPOrV41RmGsu//Jo214wmedUeXv9uvyMtvEN6cezjZ6kfGMm6hJOMHjHCId+QIUN44403WLhwIbNmzWKEk5klvxASEsLChQu5/fbb2bhxI/n5+WRnZ3PTTTcxa9YsJk6cyKFDh4iLi+ODD6z0cN7B68yQIUOYNWsWs2bNQkTYsmUL3bp1A+CLL74gLS0Nb29vFi9ezP/+9z969OjhmBYGGD58OK+++iq3334769evx9fXF39/f9LT02nQoAFeXl4cP36cNWvW8Nhjj5UIf/jw4cybN48pU6Ywb948R3yDgoJYtWoV/fr1IyUlhd27dxMcHFzErb+/Pz4+Pqxbt46YmBjeeustHn744fP6azCURo08a9ccyFA3aD9lOYWlMzdxOynvPo5ni3aWNgOaXnsXDUJ6kp9zgmOLp5F34hgePi1pPmIK7t6NyT+VzuF5j1JwJhvEDTfP+rS5/9+4eTUg47t3ydr1HeLmRr1WIfgNfQTxKLopKS8zhaMLn6bNfa8VeV5w9jTpq17ndNJOUMWjSStajvw7p37+kuw96yg4nUX+yeM07NyfJlffUSJeqkraF/9H7oFNiIcXfsMexcs/jNzEnaR9/qoVP1Ua9xhB464l126z923g2MJnwN0DcfMgomMozz//PMOGDePVV1/lH//4Bz4+PlxxxRV88MEHNGvWjIEDB7Jjxw7H+mVubi6BgYFceeWVdOjQgUWLFuHu7s6pU6fw8PDAz8+PBQsWOJRPfHw8N954I7/88ksRWXJycnj00UdZu3YtBQUFtG/fnmXLljF37lwWL15MZmYmiYmJjB07lr///e+lpsXEiRP57LPPaNCgAXPmzKFHjx6sXbuWBx98EDc3NwoKCnj00Ue57777SrhPTU3ltttu4+DBg0Xim5yczN13383hw4dRVaZMmeKYro6OjnZM727cuJG7776bnJwcrr/+ekeHoCx/azrmQIZLg1GkhmpH7tl83vhuPy+u/JXSSmdAE2/WTBlw2eWqTF2pTLWqiNV+/1pNckbJd13rubvx/oO96RbUtOIBGuosRpFeGsyh9YZqg6qy9KdkBr74DS+s/JXIAB+8PIoWUW9PdyYP6egS+Qo3nVTkcnOr+OVegeuxIZ3w9nQvIo+nu+DlIdzy2loemb+FxPRsl6SLwVDXMWukhmrB1kMZ/GPZDjYlpNPZ34cXRnWlT4hfiV27k4d05OZuAa4W97JTGOfiaTGocyv+880+/vvdfj7bfoT7rm7PQ7EhNK5f8t1Zg8FwaTBTuwaXcjgzh399tpuPtyTRvJEXjw3pyK1XBuLuZl43qAzJGTm88PluPtqShF/DevxxcAdG92iLh7uZdDKcw0ztXhqMIjW4hOwzefznm/3859t9FCg80K89v4sNpZGXmSS5ELYlZvDs8p1sOJBGWMtG/OWGcGI7mlc3DBZGkV4ajCI1XFYKCpTFW5P412e7OXIilxuj/Hl8aCfaNmvgatFqDarK59tTmPrpThJSs7mmQwv+Miycjq0bu1o0g4sxivTSYBSp4bKxMT6NfyzbwU+JmXQN9OVvN3amR7ua/0pBdeVMXgFvr0tg5pe/cup0HqN7BvHHQR1o0djL1aIZXIRRpJcGo0gNl5xDadlM+2wXy7cdppWPF48P7cTN0QG4mXXQy0J61hleWb2Ht39IwMvDjYf6h3Lf1e2pX2wXsKH2YxTppcEoUsMl49TpPP799V5e/+4AbgIPXhPCg9cG06CeWQd1BfuPnWLap7tYuSOFgCbePDa0IzdFtTEdmjqEUaSXBqNIDRed/AJl0aZEpq/czbGTp7mlWwCPDe2Iv6+3q0UzYB2Z+OzyHWxPPkHXtk342w3hZoq9jmAU6aXBKFLDReWHfan8Y9kOdhw+QfegJjx5UwTRbZu4WixDMQoKlI+2JDH9812knDjNsC6tmTI0nCA/s+mrNmMU6aXBKFLDRSEhNYt/rtjJ59utacPHr+/ETVH+5vNT1ZzsM3m8/u0B/u+bfeQXKHf3bceE/qH4epsDHWojRpFeGowiNVwQJ3LPMnv1XuasicfDXXgoNoT7+wWbjSw1jJQTuby4cjcfbkqkibcnj17XgTtigvA0BzrUKowivTQYRWqoEnn5BSz48RAzvviVtOwzjOweyOQhHWnpU9/VohkugO3JmTy3fCdr96US3KIhfxkWzoBOLc3MQi3BKNJLg1Gkhkrz3Z5jPLtsJ7tTTtKrfTOevLEzkQG+rhbLcJFQVVbtPMo/P93J/mNZXBXix19uCCeijcnjmo5RpJcGo0gNFWbfsVP8c/lOVu06Sttm3jxxfThDI1ub0Uot5Wx+Ae+tP8jLX/5KRs5ZRl0ZyJ8Gd6SVmXWosRhFemkwitRQLhnZZ5i5ynqhv76nOw8PCOXuvu3w8jDroHWBzJyzzP5qL3PWHMDDzY3fXhvCA9e0N+8D10CMIr00GEVqKJOz+QW8uy6Bl1ft4UTOWXPEXB0nITWL5z/bxYqfj9Dapz6Th3Tklm7mhKqahFGklwajSA2l8tXuozy7bAf7jmXRN9SPv97QmXB/H1eLZagG/BifxrP2mcmRAT78ZVhn+oT4uVosQwUwivTSYBSpoQi/ppzk2eU7+fbXY7Rv3pAnhoVzXbjZtWkoSkGB8sm2ZJ7/dBfJmbkM7tyKPw8Lp33zhq4WzXAejCK9NJiXxOoAhw4don///nTu3JmIiAhmzpzpMEtLS2PQoEGEhIYSHN2HwdNWsPVgOvd38SJn4Z+5sVsQL774YhH/ZsyYQUREBJGRkYwZM4bc3NwSYcbHx+Pt7U10dLTjOnPmTJkyzp07l4kTJ5YbF1XlkUceITQ0lKioKDZv3uwwc3d3d4Q1fPjwUt0XxjcsLIxBgwaRnp7uMPv666+Jjo4mIiKCa6+9tlT3Bw4cICYmhtDQUEaPHu2I0+nTpxk9ejShoaHExMQQHx9fblxqMm5uwojoAFZPimXykI6s2XucQS99w9OfbCcju+x8LqQiZbJ4Hu3atYs+ffrg5eXFCy+8UMS/2lgm09PTueWWW4iKiqJXr1788ssvpbo3ZbIaoKo17rryyivVUHGSk5N106ZNqqp64sQJDQsL0+3bt6uq6h//NEnjxk/WyL9/ps1ix2mfW+7R1FOnNSUlRTds2KBPPPGETp8+3eFXYmKitmvXTrOzs1VVddSoUTpnzpwSYR44cEAjIiIqLOOcOXN0woQJ5dpbvny5Dh06VAsKCvSHH37QXr16OcwaNmxYrvvJkyfr1KlTVVV16tSp+thjj6mqanp6uoaHh2tCQoIbfhgwAAAgAElEQVSqqqakpJTqftSoUTp//nxVVX3wwQf1tddeU1XV2bNn64MPPqiqqvPnz9fbbrutXFlqE0dP5OqURdu0/ZRlGvXU5/rGd/v19Nn8Mu2fr0yWlUd1rUxOmjRJn3rqKVVV3blzpw4YMKBU95Upk8BGrQZteG27zIi0DuDv70/37t0BaNy4MeHh4SQmJrJy+xH+PW8B6z270D2oKUtn/oXU7Wto1rAeLVu2pGfPnnh6ljwqLi8vj5ycHPLy8sjOzqZNmzYVliUrK4t7772XXr160a1bN5YsWeIwO3ToELGxsYSFhfH000+X6n7JkiXcddddiAi9e/cmIyODw4cPVzj8JUuWMG7cOADGjRvH4sWLAXjvvfeIi4sjKCgIgJYtW5Zwq6qsXr2akSNHlnDv7O/IkSNZtWoVqjVv2aSqtGjsxdS4Lnz6+2uICvTlH8t2MHjGN3z2y5FS06G0MpmUlASUnUd1rUzu2LGDAQMGANCpUyfi4+NJSUkp4rayZdJwaTCKtI4RHx/Pho2beX2XO+Pf3sTZU+m8/fAQ5t3bi6u6hJaoqMUJCAhg0qRJBAUF4e/vj6+vL4MHDy7V7r59+xzTWhMmTADgueeeY8CAAWzYsIGvvvqKyZMnk5WVBcCGDRtYtGgR27Zt48MPP6S0dfCkpCTatm3ruA8MDHQ0wLm5ufTo0YPevXs7GpPipKSk4O/vD0Dr1q0d8f31119JT08nNjaWK6+8krfeesvhZtiwYSQnJ5OamkqTJk3w8PAoEbazXB4eHvj6+pKamnretKyNdGzdmLfvi2HuPT3xdHfjt+9sYvR/1/FzYmaZbuLj49myZQsxMTFA2XlUFrW1THbt2pWPPvrIIUdCQgKJiYlA1cskYN5ZugRc1kQVkSbAG0AkoMC9wG7gfaAdEA/cpqrpZXhhqASLtyQx/fPdJGfk0KaJN+N6tuLJ8aMoiLmLPRn5PDMigkf+z4P+Ha3Rl4iUu6koPT2dJUuWcODAAZo0acKoUaN45513GDt2bAm7ISEhbN26tcizlStXsnTpUscaV25uLgcPHgRg0KBB+PlZuz/j4uL4/vvv6dGj4vsiEhISCAgIYP/+/QwYMIAuXboQEhJSpn3n+Obl5bFp0yZWrVpFTk4Offr0oXfv3nTo0IEVK1YAcPz48QrLUteJ7diSq0ObO46RvOnV74nrFkBUW19e//aAo0xO7BfIC4+M4eWXX8bHp+Su8LpcJqdMmcLvf/97oqOj6dKlC926dcPd3Xp325TJ6sXlHpHOBD5T1U5AV2AnMAVYpaphwCr73nCBLN6SxJ8/+pmkjBwUSEw9ye/v/w1nrriKh++9k28m9eeuPu1o1aqVYxrq8OHDpU5pOvPll1/Svn17WrRogaenJ3Fxcaxdu5b169c7evpLly4t072qsmjRIrZu3crWrVs5ePAg4eHhACUaTBFh9uzZDn+Tk5MJCAjg0KFDDjuJiYkEBAQAOH6Dg4OJjY1ly5YtJcIvK76BgYEMGTKEhg0b0rx5c6655hp++umnIm79/PzIyMggLy+v1LAL5crLyyMzM9PRANdVPNzdGNv7Cr6eHMvvYkNYsjWJp5buKFImH7hrDFGxNxIXF+dwZ8qkFV8fHx/mzJnD1q1beeuttzh27BjBwcFF3Fa2TAJ5501MQ5W4bIpURHyBa4A3AVT1jKpmACOAeba1ecDNl0um2sz0z3eTczYfsBqK1E9n4unXltCBt/O3Gzvj28BaZxo+fDjz5lnJP2/ePEaMGHFef4OCgli3bh3Z2dmoKqtWrSI8PJyYmBhHQ1TW7kSAIUOGMGvWLMe6mXPD8sUXX5CWlkZOTg6LFy+mb9++TJgwweFvmzZtGD58OG+99Raqyrp16/D19cXf35/09HROnz4NWL30NWvW0Llz5xLhlxXfESNG8P333zvW2NavX+9oTAsREfr378/ChQtLuHf2d+HChQwYMMC8MmTTuL4njw/tRHOngzwKy6Rb00B+bVF0h7Qpk1Z8MzIyHDtw33jjDa655poSo/bKlknDJeJy7WoCooENwFxgC9YUb0Mgw8mOON8Xcz8e2AhsDAoKUsP5aff4Mr3Cvlrd+bwC6tminXq2bK9du3bV5cuXq6rq8ePHdcCAARoaGqoDBw7U1NRUVVU9fPiwBgQEaOPGjdXX11cDAgI0MzNTVVWffPJJ7dixo0ZEROjYsWM1Nze3RPhl7ZDMzs7W8ePHa2RkpHbu3FlvuOEGVbV2SI4YMUJjY2M1NDTUsVuxOAUFBfrQQw9pcHCwRkZG6o8//qiqqmvWrNHIyEiNiorSyMhIfeONN0p1X1Z8VVX/9a9/aXh4uEZEROiMGTMcz6+//npNSkpSVdV9+/Zpz549NSQkREeOHOmIe05Ojo4cOVJDQkK0Z8+eum/fvvPkTt3ElMnKlcm1a9dqWFiYdujQQW+55RZNS0tzuKlqmcTs2r0k12U7kEFEegDrgL6qul5EZgIngIdVtYmTvXRVbXo+v8yBDOVz1bRVJGeUfJcuoIk3a6aYnqnh8tN32mqSMnJKPDdl8vJhDmS4NFzONdJEIFFV19v3C4HuQIqI+APYv0cvo0y1lt7BJdfnvD3dmTykowukMRhg8pCOeJfywfcR3Sr+qorBUB25bIpUVY8Ah0SksCUfCOwAlgLj7GfjgCWlODdUgr1HT7F822E6tW5MQJP6CFavf2pcF27uFuBq8Qx1lJu7BTA1rgsBTbwRwN+3Pi0a12PRpkTSsso/DclgqK5c1rN2RSQaa220HrAfuAdLmX8ABAEJWK+/pJ3PHzO1WzZn8wuIe20th9KzWfnoNbQ03440VGN+Scok7rW1XNOhOa/f1cNs0LrEmKndS8NlfY9UVbcCpWXiwMspR21m1qo9/JyUyb/v7G6UqKHaExngy5TrO/HMsh289UMC465q52qRDIZKY042qkVsPpjOq1/tJa57ANd38Xe1OAZDhbinbzsGdGrJcyt2siP5hKvFMRgqjVGktYSs03n88f2t+Pt689TwCFeLYzBUGBFh+sgofL09eXj+ZnLO5LtaJIOhUhhFWkt4bsVOEtKyefG2rvjUL3mot8FQnfFr5MXLo6PZfzyLZ5btcLU4BkOlMIq0FrB6VwrvrT/I+H7Bpb72YjDUBPqGNufBa0KYv+EgK36u+NdTDAZXYxRpDSf11GkeW/gznVo35o+DO7haHIPhgvjT4A50bduEKYu2kZie7WpxDIYKYRRpDUZV+fNHP3Mi5ywzRkfj5VHyZXeDoSbh6e7GK7dHU6Dw6IKt5OUXuFokg6FcjCKtwSzclMjKHSlMGtKBcP+Sn6AyGGoiV/g15LlbItmYkM4rq/e6WhyDoVyMIq2hHErL5ulPdhDTvhn3XR1cvgODoQYxIjqAW7sH8urqPazbX/c+kG6oWRhFWgPJL1D+9IH1rcwXb+uKu5s5DcZQ+3h6RARX+DXkD+9vJd0cIWioxhhFWgN5/bv9bIhP4+nhEQQ2beBqcQyGS0IjLw9eub0bx0+d5vFF27icx5kaDJXBKNIaxo7kE7y4cjfXR7Ymrrs5gN5Qu+kS6MvjQzuxckcK76w/6GpxDIZSMYq0BpF7Np8/vL+VJg3q8dwtXcwB34Y6wb1923Nthxb8Y9kOdh0xRwgaqh9GkdYgXly5m90pJ/nXyCiaNaznanEMhsuCm5vwwijrxK5H5m8h96w5QtBQvTCKtIawdt9x3vj+AGN7B9G/Y0tXi2MwXFZaNPbipdu68mvKKZ5dbo4QNFQvjCKtAZzIPcukD36inV9DnhgW7mpxDAaXcE2HFoy/Jph31h3ks1+OuFocg8GBUaQ1gKeWbCfl5Gleuq0rDepd1k/IGgzVikmDO9IlwJfHF20jOSPH1eIYDIBRpNWe5dsO89GWJCb2D6VbUFNXi2MwuJR6Hm68MqYbefkFPLpgK/kF5pUYg+sxirQak3Iil78s/pmugb5MHBDqanEMhmpB++YN+cfNkWyIT+NVc4SgoRpgFGk1RVWZvHAbuWfzeWl0NJ7uJqsMhkLiugdyS7cAZq76lR/j01wtjqGOY1rnaso76xL49tdj/GVYOCEtGrlaHIOh2vHMiAjaNmvA7+dvITP7rKvFMdRhjCKthuw7dornVuzkmg4tGNv7CleLYzBUSxrX9+SV27tx9ORppnxkjhA0uA6jSKsZZ/ML+OP7W6nv6c70kVHm9CKD4Tx0bduEyUM68ukvR5i/4ZCrxTHUUYwirWa8unovPyVm8s9butDKp76rxTEYqj0P9AumX1hznlm2nT0pJ10tjqEOYhRpNWLLwXRe/Wovcd0CGNbF39XiGAw1Ajc34cXbutKwngcPmyMEDS7AKNJqQvaZPP74wU+09qnPUyMiXC2OwVCjaNm4Pi/c1pVdR04ydcVOV4tjqGMYRVpNeG75TuJTs3jxNutwboPBUDn6d2zJfVe3Z94PCXyxI8XV4hjqEJdVkYpIvIj8LCJbRWSj/ayZiHwhInvs3zp3fM9Xu47y7vqDPNAvmN7Bfq4Wx2CosTw2tCMRbXyYvPAnjmTmFjE7dOgQ/fv3p3PnzkRERDBz5kyHWVpaGoMGDSIsLIxBgwaRnp4OwK5du+jTpw9eXl688MILRfybMWMGERERREZGMmbMGHJzi4YHEB8fj7e3N9HR0Y7rzJkzZco/d+5cJk6cWG48VZVHHnmE0NBQoqKi2Lx5s8PM3d3dEdbw4cNLdS8io0Rku4gUiEiPYmZ/FpG9IrJbRIaU4b69iKy37b0vIvXs5172/V7bvF25kakFuGJE2l9Vo1W1MPOmAKtUNQxYZd/XGdKyzjB54TY6tW7MnwZ3cLU4BkONxsvDnVljunEmr4BH399S5AhBDw8PXnzxRXbs2MG6deuYPXs2O3ZYX5KZNm0aAwcOZM+ePQwcOJBp06YB0KxZM1555RUmTZpUJJykpCReeeUVNm7cyC+//EJ+fj4LFiwoVaaQkBC2bt3quOrVu/BPIH766afs2bOHPXv28N///pff/e53DjNvb29HWEuXLi3Li1+AOOBb54ci0hm4HYgAhgKviYh7Ke6fB2aoaiiQDtxnP78PSLefz7Dt1Xqqw9TuCGCe/X8ecLMLZbksOPeMQzt2IuGbD5kxOhovD/c61zMuK75ff/01vr6+DvfPPPNMqe4PHDhATEwMoaGhjB492hGn06dPM3r0aEJDQ4mJiSE+Pr7cuBhqB8EtGvH08AjW7U/j31+fO0LQ39+f7t27A9C4cWPCw8NJSkoCYMmSJYwbNw6AcePGsXjxYgBatmxJz5498fQsudySl5dHTk4OeXl5ZGdn06ZNmwrLmJWVxb333kuvXr3o1q0bS5YscZgdOnSI2NhYwsLCePrpp0t1v2TJEu666y5EhN69e5ORkcHhw4crHL6q7lTV3aUYjQAWqOppVT0A7AV6OVsQ6528AcBC+5Fzu+3cni8EBkodeIev0opURNxEpJuI3CAiA0SkMh/HVGCliGwSkfH2s1aqWlgCjgCtygh3vIhsFJGNx44dq6zY1YrCnvEzb6+k0ajncdu1Ek1PBOpez7is+AL069fP4f7JJ58s1f3jjz/OH/7wB/bu3UvTpk158803AXjzzTdp2rQpe/fu5Q9/+AOPP/74BcfVUHMYeWUgw7u2YcaXe9iUUPIIwfj4eLZs2UJMTAwAKSkp+PtbO+Vbt25NSsr511gDAgKYNGkSQUFB+Pv74+vry+DBg0u1u2/fPkeHcMKECQA899xzDBgwgA0bNvDVV18xefJksrKyANiwYQOLFi1i27ZtfPjhh2zcuLGEn0lJSbRt29ZxHxgY6OgU5Obm0qNHD3r37u3oEFSCAMD5hdxE+xkiskJE2gB+QIaq5hW34+zeNs+07ddqKqxIRSRERP6L1UOZBowBHgK+FJF1InKPiJTn39Wq2h24HpggItc4G6p1NEmpx5Oo6n9VtYeq9mjRokVFxa6W+Pv706JdJ55aup3enQKJ6RZVZ3vGZcW3Iqgqq1evZuTIkSXcO/s7cuRIVq1aZU6+qUOICM/eEkmbJvW5b+5G+kxdRfspy+k7bTXz1/zKrbfeyssvv4yPj0+pbssbRKWnp7NkyRIOHDhAcnIyWVlZvPPOO6Xade7Azp49G4CVK1cybdo0oqOjiY2NJTc3l4MHDwIwaNAg/Pz88Pb2Ji4uju+//75ScU9ISGDjxo289957PProo+zbt69S7stCVYepavJF8ayWUZkR6bPAO0CIqg5R1bGqOlJVo4DhgC/wm/N5oKpJ9u9R4GOsKYMUEfEHsH+PVj4aNYv8AuVPH/wEwKMxTdi6te72jM8X3x9++IGuXbty/fXXs337dsfzYcOGkZycTGpqKk2aNMHDw6NE2M5yeXh44OvrS2pq6nnT0lC78Knvya3dA8nIOcvhzFwUSEw9yQN3jSEq9kbi4uIcdlu1auXoAB4+fJiWLc8/0fbll1/Svn17WrRogaenJ3Fxcaxdu5b169c76th51idRVRYtWuRQsAcPHiQ8PByghBIXEWbPnu3wNzk5mYCAAA4dOjdwTExMJCDAGhQW/gYHBxMbG8uWLVsqnmiQBLR1ug+0nzmTCjQREY9S7Djc2+a+tv1aTYUVqaqOUdVvtZRuvaoeVdWXVXVeaW4BRKShiDQu/A8MxlrwXgqMs62NA5aU7kPNZ/GWJPpOW03IEyvYEJ/GkA6+PHzfWNMzLiW+3bt3JyEhgZ9++omHH36Ym28+t3S+YsWKSo26DXWXDzcmOv6rKqmfzsStaSC/tri2iL3hw4czb57VfM2bN48RI0ac19+goCDWrVtHdnY2qsqqVasIDw8nJibGUcfK2hcAMGTIEGbNmuWYJXFWdl988QVpaWnk5OSwePFi+vbty4QJExz+tmnThuHDh/PWW2+hqqxbtw5fX1/8/f1JT0/n9OnTABw/fpw1a9bQuXPnyiTZUuB2e/dteyAM2OBswdYBXwEj7UfO7bZzez4SWF2azqhtVGWNNK6Ua2AF1kpbAd+LyE9YGbNcVT/DmiYeJCJ7gOvs+1rH4i1J/Pmjn0nKyAFA8/P4z5MT6nzPuKz4+vj40KiR9dWbYcOGcfbsWY4fP17ErZ+fHxkZGeTl5ZUadqFceXl5ZGZm4udX65dqDMVItusbwOmkHWRt/4rcg9v4ccb9REdHs2LFCgCmTJnCF198QVhYGF9++SVTplgvDxw5coTAwEBeeuklnn32WQIDAzlx4gQxMTGMHDmS7t2706VLFwoKChg/fnypMpTG3/72N86ePUtUVBQRERH87W9/c5j16tWLW2+9laioKG699VZ69OhRwv2wYcMIDg4mNDSUBx54gNdeew2AnTt30qNHD7p27Ur//v2ZMmVKqYpURG4RkUSgD7BcRD4HUNXtwAfADuAzYIKq5ttuCtdIAR4H/igie7HWQN+0n78J+NnP/0hdeQtDVSt1AcuBNGCRfaUCK4E9wG8q619VriuvvFJrGldNXaVXPL5Mr3h8mQY99ok2jOivja8crldNXVXE3qRJk3Tq1Kmqqjp16lSdPHlyEfO///3vOn36dMf9unXrtHPnzpqVlaUFBQV611136SuvvFIi/AMHDmhERESJ53/+8591woQJWlBQoKqqmzdvVlXVOXPmqL+/v6ampmp2drZ26dJFf/zxxxLuly1bpkOHDtWCggL94YcftGfPnqqqmpaWprm5uaqqeuzYMQ0NDdXt27eXcF9WfA8fPuyQaf369dq2bVvHvTMjR47U+fPnq6rqgw8+qLNnz1ZV1VdffVUffPBBVVWdP3++jho1qoRbQ+3Hud45X8XrXV0B2KiXoY2ua1dVFOnnWDttC+9b2c+aAb9cDqFroiJt51SJW935vALq2aKderZsr127dtXly5erqurx48d1wIABGhoaqgMHDtTU1FRVtRRLQECANm7cWH19fTUgIEAzMzNVVfXJJ5/Ujh07akREhI4dO9ahwJwpS5FmZ2fr+PHjNTIyUjt37qw33HCDqlqKdMSIERobG6uhoaH61FNPlRqvgoICfeihhzQ4OFgjIyMdynbNmjUaGRmpUVFRGhkZqW+88Uap7suK76xZs7Rz584aFRWlMTExumbNGoeb66+/XpOSklRVdd++fdqzZ08NCQnRkSNHOuKek5OjI0eO1JCQEO3Zs6fu27evnBwy1EY+3pyonf76aREl2umvK/TjzYmuFs0lGEV6aS6x0rbiiMgOVe3sdC/AdlXtLCJbVLXbhYyQK0KPHj20tI0v1Zm+01Y7pnWdCWjizZopA1wgkcFQN1i8JYnpn+921L97rrqCvw+PdLFUrkFENum5w3AMF4mqHMjwtYgsE5FxIjIOa3H5a3sDUcbFFa/28OC1wSWeeXu6M3lIRxdIYzDUHW7uFsCaKQPY989htG3mzbakE64WyVDLqIoinQDMAaLtax7WgnSWqva/mMLVJk6dtjbEtGzshWCNRKfGdeHmbgHnd2gwGC4K7m7CvX3bsykhnc0H010tjqEW4VG+laKoqtoHzmeq6pci0gBoBJgv6pZBQYHy3vqD9An2Y/743q4Wx2Cos9zWoy0vffErb353gO531rnvYxguEVV5/eUBrDMU/2M/CgAqfQ5VXeKbPcdITM/hzt5BrhbFYKjTNPTy4I6YID795TCH0rJdLY6hllDVqd2+wAkAVd0DVOa83TrHu+sO0rxRPQZ3bu1qUQyGOs/dV7XDTYQ5a+JdLYqhllAVRXpaVR2fDbGPgar1J1dUleSMHFbvSuG2Hm2p51EdPrZjMNRt/H29uTHKn/d/PMiJ3LOuFsdQC6hKy/6NiDwBeIvIIOBD4JOLK1btYcGPh1BgTC8zrWswVBfu7xdM1pl8Fmw46GpRDLWAqijSKcAx4GfgQWAF8NeLKVRt4Wx+AQs2HCS2QwvaNmvganEMBoNNZIAvvYObMXdNPGfzC1wtjqGGU2lFqqoFqvq6qo5S6+svr2tlT3WoI6zamcLRk6e5M+YKV4tiMBiKcf/VwSRn5vLpL0dcLYqhhlPh119E5GfOsxaq1ufUDE68u/4gbXzr07+T2YtlMFQ3BnRqSXDzhrzx3X5uivIv90tLBkNZVGZEeiNwE9YXAT4D7rSvT7Gmdw1OxB/P4rs9x7m9VxDubqaCGgzVDTc34d6r27MtMZMf480BDYaqU5nvkSaoagIwSFUfU9Wf7etxrG+LGpyYv+Eg7m7C6J5ty7dsMBhcwq3dA2nawJPXv9vvalEMNZiqbDYSEenrdHNVFf2ptZzOy+eDjYcY3LkVrXzqu1ocg8FQBt713Bnb+wq+3JnCgeNZrhbHUEOpigK8D3hNROJFJB54Dbj3okpVw/n05yOkZ581m4wMhhrAb/pcgaebG3PWHHC1KIYaSlV27W5S1a5AV6Crqkar6uaLL1rN5d31CbTza8BVIX6uFsVgMJRDy8b1GR7dhg83JpKRfaZ8BwZDMSqsSEVkrIg47KtqpqpmOpmHiMjVF1vAmsbuIyf5MT6dO2KCcDObjAyGGsH9/dqTczafd9ebAxoMlacyX3/xA7aIyCZgE9ahDPWBUOBa4DjWYQ11mvfWJ1DPw42RV5pNRgZDTaFTax/6hTVn3tp4HugXbI7zNFSKyuzanQl0B+YDLYCB9n0S8BtVvdU+wL7Okn0mj482J3FDF3+aNaznanEMBkMluL9fMEdPnuaTn5JdLYqhhlGp75Gqaj7whX0ZirF0azInT+dxZ4w5V9dgqGlcE9acsJaNeOP7A8R1DzAHNBgqjJm/uIi8u/4gHVs15sorzAeDDYaahohwf7/27Dx8grX7Ul0tjqEGYRTpRWJbYgY/J2VyZ+8g05M1GGooI6IDaN6oHm+YAxoMlcAo0ovEu+sO4u3pzs3dAlwtisFgqCL1Pd35Te92fLX7GHuPnnS1OIYaQqUVqYh4icgdIvKEiDxZeF0K4WoKmTlnWfpTMjd3a4NPfU9Xi2MwGC6Asb2D8PJw483vzQENhopRlRHpEmAEkAdkOV11lo83J5JzNp87epmTjAyGmo5fIy/iugeyaHMSqadOu1ocQw2gUrt2bQJVdWhVAxQRd2AjkKSqN4pIe2AB1nuqm7Bepakxx4uoKu+uP0jXQF+6BPq6WhyDwXARuO/q9szfcJC31yXw6HUdXC2OoZpTlRHpWhHpcgFh/h7Y6XT/PDBDVUOBdKyzfGsMP8ans+foKXOursFQiwht2YgBnVry9g8J5J7Nd7U4hmpOVRTp1cAmEdktIttE5GcR2VYRhyISCNwAvGHfCzAAWGhbmQfcXAWZXMa76xNoXN+DG7v6u1oUg8FwEbn/6vakZp1hydYkV4tiqOZUZWr3+gsI72XgMaCxfe8HZKhqnn2fCNSYba+pp07z6c9HuCMmiAb1qpKUBoOhutInxI9wfx/e+O4At/Voa15rM5RJVb7+kgA0AW6yryb2s/MiIjcCR1V1U6WltNyPF5GNIrLx2LFjVfHiovPhpkTO5BeYk4wMhlqIiPBAv/bsOXqKb36tHm2OoXpSlddffg+8C7S0r3dE5OEKOO0LDLe/YboAa0p3JtBERAqHc4FYZ/eWQFX/q6o9VLVHixYtKiv2RaegQHlv/UF6tW9GWKvG5TswGAw1jhuj2tDKx4s3vjOvwhjKpqof9o5R1SdV9UmgN/BAeY5U9c+qGqiq7YDbgdWqeifwFTDStjYO6/Waas/3e49zMC3bjEYNhlpMPQ83xl3Vju/3Hmfn4ROuFsdQTamKIhXAeRtbvv2sqjwO/FFE9mKtmb55AX5dNt5dn4Bfw3oMjWztalEMBsMl5I5eQXh7upsDGgxlUpUdMnOA9SLysX1/M5VUfqr6NfC1/X8/0KsKcriMI5m5fLnzKA/0C8bLw93V4hgMhktIkwb1GNUjkPkbDvLYkI609KnvapEM1YyqbDZ6CbgHSLOve1T15YstWHVmwY8HyaARVMcAACAASURBVC9Q7uhlpnUNhrrAvX3bk1egvPWDta/y0KFD9O/fn86dOxMREcHMmTMddtPS0hg0aBBhYWEMGjSI9PR0AHbt2kWfPn3w8vLihRdeKOL/jBkziIiIIDIykjFjxpCbm1tChvj4eLy9vYmOjnZcZ86UfXbN3LlzmThxYrlxE4tXRGSv/UpjdyezfBHZal9Ly3A/SkS2i0iBiPRweu4nIl+JyCkRefU84TcTkS9EZI/927Q8uaobFVakIuJj/zYD4oF37CvBflYnyMsvYMGGQ1zToQVBfg1cLY7BYLgMtGvekEHhrXhnfQLZZ/Lw8PDgxRdfZMeOHaxbt47Zs2ezY8cOAKZNm8bAgQPZs2cPAwcOZNq0aQA0a9aMV155hUmTJhXxOykpiVdeeYWNGzfyyy+/kJ+fz4IFC0qVIyQkhK1btzquevXqXYzoXQ+E2dd44N9OZjmqGm1fw8tw/wsQB3xb7Hku8DdgUgkXRZkCrFLVMGCVfV+eXNWKyoxI37N/N2Ed8Vd4Fd7XCVbvOsqRE7lmk5HBUMd44JpgMrLPsmhzEv7+/nTvbg2QGjduTHh4OElJ1gsHS5YsYdy4cQCMGzeOxYsXA9CyZUt69uyJp2fJD1vk5eWRk5NDXl4e2dnZtGnTpsJyZWVlce+999KrVy+6devGkiXn9mseOnSI2NhYwsLCePrpp8vyYgTwllqsw3qTosInzKjqTlXdXcrzLFX9Hkuhno8RWIfxQNFDeS5IrstJhRWpqt5o/7ZX1WCnq72qBl86EasHhVM5tw3qw9E5E9i+8lyPsaZN5agqjzzyCKGhoURFRbF582aHmbu7uyOs4cNL74CWFd9C/r+9O4+vqj4TP/55sgARIWwBwg0QNoEkhKBhUaxCIuKIDRjEqS0jdtrSVlrH6U8qbadqO9OCVcet1NalFoZaq6LAC1BZxFrQsEjYEVmFLAQIyJ4AyfP745xcL+Te7Mm9SZ7363VeOfd7tuccTvJwvud7v9/169cTERHBW2+95Xf7Tz/9lEGDBtG3b18eeOABVLVK+zUmmFJ7tmdwXDR/Xr2f0lL1lh84cIDs7GyGDx8OQEFBAbGxzt/7rl27UlBQUOF+PR4PDz30ED169CA2Npbo6GhuvfVWv+vu3bvX+/s5bdo0AH7zm9+QlpbGunXrWLVqFdOnT+fsWWcckXXr1jF//ny2bNnCm2++CeCvGs0DHPL57NsxTiv3+/tZIlJnvc6JyMs+1cBdVDXfnT8MdKlCXCGlJt8jXVmVsqYmIiKCh375P7Sf8nseeelt/vjCHxptVc67777L7t272b17Ny+++CI//OEPvcuioqK8x1q0yO8rkYDnC1BSUsLDDz8c8A8BwA9/+ENeeuklbwzvvfdepfs1JthEhOTu7dh/7Cy9f76UkbM+4G9rPmfixIk888wztG3b1u82lfWIdOLECRYuXMj+/fvJy8vj7NmzzJs3z++6vn8PZs+eDcCyZcuYNWsWKSkpjBo1iqKiIg4ePAjAmDFj6NixI1FRUWRmZgJcXc3T7qmqqcA3gWdEpE81t/dLVb+rquVqMtX5X7X62SSkVecdaSv3XWgnEWnvviDuICLxhOj/EupSbGwsW4s7EibCfTcnNOqqnIULF3LvvfciIowYMYIvv/yS/Px8v+sG2t7f+QI8//zzTJw4kc6dO/vdNj8/n1OnTjFixAhEhHvvvde7fUX7NSbYFmTn8uaGrx6QcgpP87177yF51B1lSQqALl26eH+f8vPzA/4ulFmxYgW9evUiJiaGyMhIMjMz+fjjj1m7dq336TPQf2rBqWGaP3++N8EePHiQgQMHAvhN4iIyzacBUTecTnC6+6zi7RhHVct+7sP5psWQCk+mZgrKqmzdn0fc8oBxhZrqPJF+H+d96AD3Z9m0EAjYIqupKL5UwpsbDpE+oDNFJw43mqqcDRvKv77Ozc2le/ev7s+4uDjvfwqKiopITU1lxIgRARNZoPPNzc3lnXfeuewJt0xKSop3nbi4OL/Hru51NKYhPfH+LooulgJO8ip891nC2sfxeczNl62XkZHBnDnOK785c+Ywfvz4Cvfbo0cPsrKyOHfuHKrKypUrGThwIMOHD/cmx0CvWQDGjh3L888/731Fkp2d7V22fPlyjh8/zvnz58t+n8+o6myfBkR5wCLgXreV7AjgpKrmuw9MLQFEpBNO73Q7qnPNqmgRTmc8cHmnPH7jqofj11qVv0eqqs8Cz4rIj1X1+XqMKaQsyM7lifd3kfvleQC6XU2dVuW0a9eOSZMmMW/ePCZPnlxu3bKqHF/Lli1j0aJF3veu/qpyADIzM1m9ejWpqalU1RdffIHH42Hfvn2kpaUxaNAg+vQJXJvje74PPvggjz/+OGFh5f9/duU5VKYq19GYhpTn/g0AKM7dwdntq4iMiWf9098l5fW2/Pa3v+X2229nxowZ3H333bzyyiv07NmTN954A4DDhw+TmprKqVOnCAsL45lnnmHHjh0MHz6cu+66i2uvvZaIiAiGDBnC1KlTqxzXL3/5Sx588EGSk5MpLS2lV69eLF68GIBhw4YxceJEcnJymDx5Mlu3bj3nZxdLgduBPcA5nK83AgwE/iQipTgPXbNUtVwiFZE7geeBGGCJiGxS1bHusgNAW6CF+471VlXdISIvA390q3dnAW+IyHeAL4C7K4kr5FS7QwZVfV5EkoAEoJVP+dy6DCwULMjO5Wdvb+W8Ox6hllziyenfZ9IE/1U5sbGx1a7KAbxVOf369eP73/8+AL/+9a9JTk72u31ZVU7//v0vK1+7dm255CMizJ49m5deegmApUuX4vF4OHTIp4oqJwePx6mdL/vZu3dvRo0aRXZ2drlEGuh8N2zYwDe+8Q0Ajh07xtKlS4mIiGDChK/aKHg8HnJycvweu7rX0ZiG1K1dlPc/1K3iEun5sJOsPO2iWDMjzbtex44dWbmyfLORrl27Xnbv+/rVr35VUataAOLj49m2bVu58qioKP70pz+VK7/vvvu47777Lit77LHHyq3nvpec5qf8Y6DSsadV9R3gnQDL4gOUf9dnvhBIr2pcoagmjY0exfnfx/PAaOB3QOB6h0bsifd3fZVEG2lVzsiRI5k2bZp3v926dSMjI4O5c+eiqmRlZREdHU1sbCwnTpyguLgYcBLhmjVrSEhIKHf8QOe7f/9+Dhw4wIEDB7jrrrv4wx/+cFkSBeddc9u2bcnKykJVmTt3rnf76l5HYxrS9LH9iYq8vCezFuFhTB/bP8AWptlQ1WpNwFacBLzZ/dwFWF7d/dRmuu6667QhxD+8WHu6U5dvPa6ARsbEa2TnXjp48GBdsmSJqqoeO3ZM09LStG/fvpqenq6FhYWqqpqfn68ej0fbtGmj0dHR6vF49OTJk6qq+sgjj2j//v01MTFRJ0+erEVFReWOv3//fk1MTCxXfu7cOZ06daomJSVpQkKCjhs3TlVVX331VR0/fryOGjVK+/btq4899pjf8yotLdX7779fe/furUlJSbp+/XpVVV2zZo0mJSVpcnKyJiUl6csvv+x3+0Dn62vKlCn65ptvej8PHjzYO79+/XpNTEzU3r1767Rp07S0tLTK+zUmmN7ZmKM3zFyp8Q8v1viHF+ukF9YEO6RqATZoA/6tbi6TONe26kRknaoOE5FPcZ5ITwM7VXVAHeb3CqWmpqq/RjR1beSsD7xVOb6urMoxxjQ//++NzSzbfpj1/3ULrSIbR5/bIvKpOl9nMXWoJqO/bBCRdsBLOK12NwKf1GlUIWL62P60irj8EkVFhltVjjGGCUO6cbr4Eqs+O1L5yqZJq0mn9fer6peq+kdgDDBFVUO2NVVtTBji4espzvc7BedJdGbmICYMafJfmzXGVOL63h3pdHVLFmwKya82mgZU7Va77ggArwMLVfVAnUcUYo6fuUD3DlF8NH20fR3DGOMVER7G1wfH8tesg5w8d5Hoq8p3vGKah5pU7T4F3AjsEJG3ROQuEWmSA/Sdv1DC6j3HSB/QxZKoMaacCSkeLpSU8u62kOwnwDSQmlTt/kNV7wd6A3/C+fJsk3xJsHrPMYovlXLLwC6Vr2yMaXaS46Lp1ak1CzflBTsUE0Q1eSJFRKKAicAPgKF8NQROk7JyZwFtWkYwrFezGW7VGFMNIsL4lG5k7S/k8MnKRgszTVVNOmR4A9gJpOH0sdtHVX9c14EFW2mpsmLnEW7qH0OLiBr9f8MY0wyMT/GgCos2W6Oj5qomGeIVnOT5A1VdpaqldR1UKNiSe5JjZ4q5ZaB1U2eMCaxXp9YMjotmQbZV7zZXNUmk/wR+JiIvAohIPxG5o27DCr6VOwsIExh1jSVSY0zFxqd42JF/it0Fp4MdigmCmiTSV4ELwA3u51zgf+osohCxfEcBqfEdaN+69gNpG2OatjsGxxImWKOjZqomibSPqv4OuAigqudw+itoMnJOnOOzw6etWtcYUyWd27RiZN9OLNycS3W7XTWNX00S6QW31a4CiEgfoLhOowqyD9wuv9Ltay/GmCqakOLh0PHzbDx4ItihmAZWk0T6KPAe0F1E/gqsBH5ap1EF2YqdR+jdqTV9Yq4OdijGmEbi1sQutIwIs0ZHzVBNOmRYDmQC9wF/A1JV9cO6DSt4zhRfImtvIelWrWuMqYY2rSK5JaELS7bmc7GkSX6ZwQRQ5UQqIteWTUBPIB/IA3q4ZZVt30pE1onIZhHZLiK/cst7ichaEdkjIn8XkaC27vnn50e5UFJq1brGmGqbkOLh+NkLrN59LNihmAZUnU7rn6pgmeJ00FCRYiBNVc+ISCSwWkTeBX4CPK2qr4vIH4HvAC9UI646tWLnEaKjIknt2T5YIRhjGqmbr4mh3VWRLNiUy+gBVqvVXFQ5karq6NocyB2d/Yz7MdKdyhLwN93yOcBjBCmRlpQqq3YdYXT/GCLCrTcjY0z1tIgI4/ZBsbyzMZezxZdo3bLaA2yZRqg6Vbs/9ZmfdMWy31ZxH+Eisgmnk/vlwF7gS1W95K6SAwRtsM/sgyc4fvaCVesaY2ps/OBunL9YwvIdBcEOxTSQ6jx2fcNn/mdXLLutKjtQ1RJVTQHigGHAgKoeXESmisgGEdlw9OjRqm5WLSt2HiEiTLi5f0y97N8Y0/QNje9At+hWNuB3M1KdRCoB5v19rpCqfgmsAq4H2olIWf1HHE5PSf62eVFVU1U1NSamfhLdyp0FDO/dgbatbIBeY0zNhIUJGSke/rn7GIVnmtRX7E0A1UmkGmDe3+dyRCRGRNq581HAGJxRZFYBd7mrTQEWViOmOvNF4Vl2HzlD+gCr1jXG1M6EId0oKVWWbLUBv5uD6iTSwSJySkROA8nufNnnQVXYPhZYJSJbgPXAclVdDDwM/ERE9gAdcUaXaXArdjq9Gdkg3saY2hrQtS0DurZhQbZV7zYH1Wm1G16bA6nqFmCIn/J9OO9Lg2rlzgL6db6aHh2vCnYoxpgmICOlG797bxcHC8/Z35Umzr7jAZw8f5F1+49zS4I9jRpj6kbG4G4ALLRGR02eJVLgH58f5VKp2mgvxpg6E9f+KobFd2DBJhsRpqmzRIpTrduhdQtSultvRsaYujN+SDf2Hj3L9rxTwQ7F1KNmn0gvlpSy6rMjpA3oTHhYkxpW1RgTZOMGxRIZLla928Q1+0S64cAJThVdsmpdY0yda3dVC26+pjOLNudRUmrVu01Vs0+kK3cW0CI8jK/1s96MjDF1b3xKNwpOFbN2X2GwQzH1pFknUlVlxc4CRvTpaJ1LG2PqxS0Du9C6RTgLN9mA301Vs06ke4+e5UDhOcZYta4xpp5EtQhnbFJXlm7Lp+hiSbDDMfWgWSfSlTud0RnSrDcjY0w9mpDi4XTRJT7cdSTYoZh60MwT6REGxrbF0y4q2KEYY5qwG/p0pNPVLVmQ7VTvHjp0iNGjR5OQkEBiYiLPPvusd93jx48zZswY+vXrx5gxYzhx4gQAn332Gddffz0tW7bkySefvGz/Tz/9NImJiSQlJXHPPfdQVFRULoYDBw4AXCsim3ymFoFiFpH7ROT3lZ2bOJ4TkT0iskVErvVZVuJzrEUBtp8kIttFpFREUn3Kx4jIpyKy1f2ZFmD7DiKyXER2uz/bVxZXXWu2ifTE2Qts+OK4VesaY+pdRHgYdyTH8sFnRzh5/iIRERE89dRT7Nixg6ysLGbPns2OHTsAmDVrFunp6ezevZv09HRmzZoFQIcOHXjuued46KGHLtt3bm4uzz33HBs2bGDbtm2UlJTw+uuvBwqlWFVTfKYLdXB6/wL0c6epwAs+y877HCsjwPbbgEzgoyvKjwFfV9VBOAOa/F+A7WcAK1W1H7DS/VxZXHWq2SbSVbuOUKrYIN7GmAYxYYiHCyWlvL/tMLGxsVx7rfOA1KZNGwYOHEhurvNd04ULFzJlyhQApkyZwoIFCwDo3LkzQ4cOJTKy/DCPly5d4vz581y6dIlz587RrVu3KsclIq1F5M8isk5EskVkvM/i7iLyofu092iAXYwH5qojC2dozNiqHl9Vd6rqLj/l2apa1kJrOxAlIi0DHH+OOz8HmFAXcVVHs02kK3ceIaZNSwZ5ooMdijGmGRgcF018x6vKDfh94MABsrOzGT58OAAFBQXExjp/77t27UpBQUGF+/V4PDz00EP06NGD2NhYoqOjufXWWwOt3tKnqnW2W/YL4ANVHQaMBp4QkdbusmHARCAZmORb9eobAnDI53OOWwbQSkQ2iEiWiEwov2mVTQQ2qmoxgIi87BNLF1UtG6/uMFD2dFRRXHWqWSbSC5dK+cfnR0kf0Jkw683IGNMARIT+Xdrw8d5Ces1YwshZH/C3NZ8zceJEnnnmGdq2bet3G5GK/0adOHGChQsXsn//fvLy8jh79izz5s0LtLpv1e40t+xWYIaIbAI+BFoBPdxly1W1UFXPA28DN1bztHuqairwTeAZEelTze0RkUTgceD7ZWWq+l1V3XDluup0atzgPV80y0S6dn8hZ4ov2dijxpgGsyA7lw8/Pwo4f+lzCk/zvXvvIXnUHWRmZnrX69KlC/n5zgNWfn4+nTtX3I5jxYoV9OrVi5iYGCIjI8nMzOTjjz9m7dq1pKSkkJKSwqJFftv5lBFgok+C7aGqO91lVyYlFZFpPk+13YBcoLvPOnFuGapa9nMfTpIuN5RmhYGJxAHvAPeq6t4AqxWUVdm6P8uaRgeMq641u0R66NAh7s0cR/4rP+SBu9IarLVcVFSU96ZOSUnhwoXA7/j/8pe/8KMf/ajSc1FVHnjgAfr27UtycjIbN270LgsPD/ceKyPD/zv+QOe7cOFCkpOTSUlJITU1ldWrV/vd/tNPP2XQoEH07duXBx54wDvCRaD9GtOcPfH+LoovlQLO727hu88S1j6Oz2Nuvmy9jIwM5sxxXvnNmTOH8ePHl9uXrx49epCVlcW5c+dQVVauXMnAgQMZPnw4mzZtYtOmTQH/BrjeB34s7qOviPgmuzFuq9gonHePa1R1tk/SzQMWAfe6rWRHACdVNV9E2pe90xSRTsBIYEcVLxci0g5YAsxQ1TUVrLoIpzES7s+FPuXl4qrq8atFVRvddN1112lN5ebm6uAfv6D//uo6PXXqlPbr10+3b9+uqqrTp0/XmTNnqqrqzJkz9ac//amqqhYUFOi6dev05z//uT7xxBPefeXk5Gh8fLyeO3dOVVUnTZqkr776arlj7t+/XxMTE6sc46uvvqrTpk2rdL0lS5bobbfdpqWlpfrJJ5/osGHDvMtat25d6faBzvf06dNaWlqqqqqbN2/W/v37+91+6NCh+sknn2hpaanedtttunTp0gr3a0xzFv/wYu3pTl2+9bgCGhkTr5Gde+ngwYN1yZIlqqp67NgxTUtL0759+2p6eroWFhaqqmp+fr56PB5t06aNRkdHq8fj0ZMnT6qq6iOPPKL9+/fXxMREnTx5shYVFZU7/v79+xWnFe1lf0+BKOBPwFacRj2L3fL7gAXAKmA38OiV27rrCTAb2OvuI9Utv8H9vNn9+Z0A29+J8/6yGCgA3nfL/ws4C2zymTq7y172OU5HnNa6u4EVQIeK4qqPKehJsSZTbRLpzvyT2vPhxfra2i9UVTUjI0OXLVumqqrXXHON5uXlqapqXl6eXnPNNZdt++ijj5ZLpHFxcVpYWKgXL17UcePG6fvvv1/umIES6ZkzZ/Tb3/62Dh06VFNSUnTBggWq6iTSjIwMvfnmm7Vv37762GOP+T2XqVOn6muvveb97Bt/VRJpZeerqvrxxx/rgAEDypXn5eVdlmBfe+01nTp1apX3a0xzc8PMld5E6jvdMHNlg8UAbNAQ+Bve1KZmV7W7cqdTfZ4+oHODtpbbu3evt6p12jTnHf9vfvMb0tLSWLduHatWrWL69OmcPXsWgHXr1jF//ny2bNnCm2++yYYN5d6rk5ubS/fuX70CiIuL8zahLyoqIjU1lREjRnibz1+povN95513GDBgAOPGjePPf/6ztzwlJcV77Li4OL/Hru51NKY5mD62P1GR4ZeVRUWGM31s/yBFZOpKs+mpfUF2Lk+8v4vcL88TGS6s3HqQJx+4p85ay7Vr145JkyYxb948Jk+eXG7dPn36sGnTpsvKli1bxqJFi7zvXYuKijh48CAAY8aMoWPHjgBkZmayevVqUlP9tTz374svvsDj8bBv3z7S0tIYNGgQffoEbjB35fneeeed3HnnnXz00Uf88pe/ZMWKFQDlzqEyVbmOxjQHE4Y437wo+zsE8OCYft5y03g1iyfSBdm5/Oztrd6b98KFiyHRWk5VmT9/vrdBwMGDBxk4cCBAueQjIsyePdu737y8PDweD4cOffU1qZycHDwe55ey7Gfv3r0ZNWoU2dnZ5Y5flfO96aab2LdvH8eOHbus3OPxkJOT4/fY1b2OxjQXE4Z4WDMjjayfpQNQfLE0yBGZutAsEukT7+/ivDvqgmrotJYbO3Yszz//PM6rCy5LdsuXL+f48eOcP3+eBQsWMHLkSKZNm+bdb7du3cjIyGDu3LmoKllZWURHRxMbG8uJEycoLi4G4NixY6xZs4aEhIRyxw90vnv27PHGtHHjRoqLi71Px2ViY2Np27YtWVlZqCpz5871bl/d62hMc9M1uhVD49uzZEv9NCI1DSzYL2lrMlW3sVEotJbz19jo3LlzOnXqVE1KStKEhAQdN26cqjqNjcaPH6+jRo2qsLFRaWmp3n///dq7d29NSkrS9evXq6rqmjVrNCkpSZOTkzUpKUlffvllv9sHOt9Zs2ZpQkKCDh48WEeMGKH//Oc/vdsMHjzYO79+/XpNTEzU3r1767Rp07wtfQPt1xjzlb+s2a89H16snx8+1WDHxBob1cskzrVtXFJTU9Vf45tARs76wFut68vTLoo1M/wOKGCMMfXqyKkihs9cyX+k9+PBW65pkGOKyKfq9DRk6lCzqNq11nLGmFDTuW0rhsV3sOrdJqBZJNIJQzzMzByEp10UgvMkOjNzkLWWM8YE1R3Jsew+coZdh08HOxRTCw329RcR6Q7MxemZX4EXVfVZEekA/B2IBw4Ad6tqnfcpN2GIxxKnMSak3JYUy6OLtrNkSx79u1oNWWPVkE+kl4D/p6oJwAhgmogkEHhQVmOMadJi2rRkRO+OLN6aT2Nsr2IcDZZIVTVfVTe686eBnThjwwUalNUYY5q8ccmx7Dt6ls+serfRCso7UhGJxxlOZy2BB2U1xpgm77bEroQJ1uioEWvwRCoiVwPzgQdV9ZTvMvd7Tn7rN0RkqjvS+oajR482QKTGGFP/Ol7dkhv6dGLxljyr3m2kGjSRikgkThL9q6q+7RYHGpT1Mqr6oqqmqmpqTExMwwRsjDENYFxyLAcKz7E971TlK5uQ02CJ1B009hVgp6r+r8+iQIOyGmNMs3BbYlfCw4QlW616tzFqyCfSkcC/AWkissmdbgdm4YzCvhu4xf1sjDHNRvvWLRjZtxNLtljr3caowb5HqqqrcUYs9ye9oeIwxphQdMegWH46fwvbck8xKC462OGYamgWPRsZY0youzWxCxFhwuItecEOxVSTJVJjjAkB7a5qwY39OrHYqncbHUukxhgTIu5I7kbul+fZnHMy2KGYarBEaowxIWJMQhciw4UlVr3bqFgiNcaYEBEdFclN/WKs9W4jY4nUGGNCyLjkWPJOFrHx4JfBDsVUkSVSY4wJIbckdKFFeJj1vduIWCI1xpgQ0rZVJDddE8PSrfmUllr1bmNgidQYY0LM1wfHcvhUERsPngh2KKYKLJEaY0yISR/YhRYRYSy26t1GwRKpMcaEmKtbRjC6v1O9W2LVuyHPEqkxxoSgccndOHK6mA0Hjgc7FFMJS6TGGBOC0gd0plVkmA2t1ghYIjXGmBDUumUEaQM6s3TrYaveDXGWSI0xJkSNG9SNY2eKWbffqndDmSVSY4wJUaMHxBAVGW5Dq4U4S6TGGBOirmoRQdrAzry37TCXSkqDHY4JwBKpMcaEsDsGxVJ49gJr3erdQ4cOMXr0aBISEkhMTOTZZ5/1rnv8+HHGjBlDv379GDNmDCdOOB06fPbZZ1x//fUA14rIQ777F5H/FJHtIrJNRP4mIq2ujEFE4kXkvIhs8plaBIpZRO4Tkd9Xdm7ieE5E9ojIFhG51mdZic+xFgXYfpIbe6mIpPqUD/PZdrOI3Blg+14istY9/t/LzklEWrqf97jL4ys6D0ukxhgTwkYP6MxVLcK9nTNERETw1FNPsWPHDrKyspg9ezY7duwAYNasWaSnp7N7927S09OZNWsWAB06dOC5554DKPDdt4h4gAeAVFVNAsKBbwQIZa+qpvhMF+rg9P4F6OdOU4EXfJad9zlWRoDttwGZwEd+ylNVNQW4DfiTiET42f5x4GlV7QucAL7jln8HOOGWP+2ujx9gAQAADDZJREFUF5AlUmOMCWGtIsO5ZWAX3tuWz6WSUmJjY7n2WufBrU2bNgwcOJDc3FwAFi5cyJQpUwCYMmUKCxYsAKBz584MHToUwF/z3wggyk00VwFVfiErIq1F5M8isk5EskVkvM/i7iLyoYjsFpFHA+xiPDBXHVlAOxGJrerxVXWnqu7yU35OVS+5H1vh57xFRIA04C23aA4wwSeuOe78W0C6u75flkiNMSbEjUuO5cS5i3y8t/Cy8gMHDpCdnc3w4cMBKCgoIDbWyUNdu3aloKCg3L58qWou8CRwEMgHTqrqsgCr9/GpLp3tlv0C+EBVhwGjgSdEpLW7bBgwEUgGJvlWvfrwAId8Pue4ZQCtRGSDiGSJyITym1ZMRIaLyHZgK/CDssQqIktFpBvQEfjSJ+H6Htsbl7v8pLu+X/4edY0xxoSQm6+JoUW48IN5n3L+Qgnd2kXxo6/F8eQD9/DMM8/Qtm3bctuICBU8RJWt0x7n6asX8CXwpohMVtV5flbf61aV+roVyPB579oK6OHOL1fVQvc4bwM3AhuqdsYA9FTVXBHpDXwgIltVdW9VN1bVtUCiiAwE5ojIu6papKq3uzF1qkYsFbInUmOMCXHvbTtMSSmcu1CCAjmFp/nevfeQPOoOMjMzvet16dKF/HznXWp+fj6dO3eubNe3APtV9aiqXgTeBm5wn+bKnj4DvZ8EEGCiz7vMHqq60112ZXWqisg0n/12A3KB7j7rxLllZU/LqOo+4ENgSGUn448bzxkg6YpFhThVyWUPlN5j+8blLo921/fLEqkxxoS4J97fRYk6eUlVKXz3WcLax/F5zM2XrZeRkcGcOc6rvTlz5jB+/Phy+7rCQWCEiFzlvgNMB3aq6lqf5Oi3xazrfeDHZe8PRcQ32Y0RkQ4iEoXz7nGNqs722W8esAi41229OwKnajlfRNqLSEt3n52AkcCOyq+Uw22NG+HO9wQGAAd811FVBVYBd7lFU4CF7vwi9zPu8g/c9f2yRGqMMSEu78vz3vni3B2c3b6KooNbWP/0d0lJSWHp0qUAzJgxg+XLl9OvXz9WrFjBjBkzADh8+DBxcXEAXYD/EpEcEWnrVn++BWzEeZcYBrxYjdD+G4gEtrjvI//bZ9k6YD6wBZivqv6qdZcC+4A9wEvA/W75QGCDiGzGSXazVLVcIhWRO0UkB7geWCIi77uLbgQ2i8gm4B3gflU95m5T9o4U4GHgJyKyB+cd6Ctu+StAR7f8J8CMii6CVJBkQ1Zqaqpu2FCdqnZjjGm8Rs76gFyfZFrG0y6KNTPSqrwfEflUVf01+jG1YE+kxhgT4qaP7U9UZPhlZVGR4Uwf2z9IERlfDZpI3e8bHRGRbT5lHURkuftdo+VuKzJjjDGuCUM8zMwchKddFILzJDozcxAThngq3dbUvwat2hWRm3BaT811e9FARH4HHFfVWSIyA2ivqg9XtB+r2jXGmOqzqt360aBPpKr6EXDleEC+PUj49ixhjDHGhLxQeEfaRVXLhoA/jNOqrBwRmer2crHh6NGjDRedMcYYU4FQSKRe7vd0/NY1q+qLqpqqqqkxMTENHJkxxhjjXygk0oKyTordn0eCHI8xxhhTZaGQSH17kPDtWcIYY4wJeQ3davdvwCigE864eI8CC4A3cDo6/gK4W1WvbJB05X6OuuvWh07AsXrad12w+GrH4qsdi692gh1fT1W1d2N1rFH2bFSfRGRDKDcPt/hqx+KrHYuvdkI9PlMzoVC1a4wxxjRalkiNMcaYWrBEWl51Rj4IBouvdiy+2rH4aifU4zM1YO9IjTHGmFqwJ1JjjDGmFiyRGmOMMbXQ7BJpdYZyE8dzIrJHRLaIyLVBiO0JEfnMPf47ItLOLY8XkfMissmd/lifsVUQ32MikusTx+0+y37mXrtdIjI2SPH93Se2AyKyyS0PxvXrLiKrRGSHiGwXkf9wy0Pl/gsUX0jcgxXEFxL3YAXxhcw9aOqJqjarCbgJuBbY5lP2O2CGOz8DeNydvx14FxBgBLA2CLHdCkS484/7xBbvu14Qr91jwEN+1k0ANgMtgV7AXiC8oeO7YvlTwCNBvH6xwLXufBvgc/c6hcr9Fyi+kLgHK4gvJO7BQPGF0j1oU/1Mze6JVKs3lNt4nLFTVVWzgHZl/QI3VGyqukxVL7kfs4C4+jp+ZQJcu0DGA6+rarGq7gf2AMPqLTgqjk9EBLgb+Ft9xlARVc1X1Y3u/GlgJ+AhdO4/v/GFyj1YwfULpEHvwcriC4V70NSPZpdIAwg0lJsHOOSzXg4V/+LWt3/HeUIp00tEskXkHyLytWAFBfzIrfb7c1m1JKF37b4GFKjqbp+yoF0/EYkHhgBrCcH774r4fIXEPegnvpC6BwNcv5C6B03dsUR6BVUNOJRbMInIL4BLwF/donygh6oOAX4CvCYibYMQ2gtAHyDFjempIMRQFfdw+ZNA0K6fiFwNzAceVNVTvstC4f4LFF+o3IN+4gupe7CCf9+QuQdN3bJE6gg0lFsu0N1nvTi3rEGJyH3AHcC33D+0uNVVhe78pzjvf65p6NhUtUBVS1S1FHiJr6rOQuLaAYhIBJAJ/L2sLFjXT0Qicf7I/lVV33aLQ+b+CxBfyNyD/uILpXuwgusXMvegqXuWSB2BhnJbBNzrtp4cAZz0qYJrECJyG/BTIENVz/mUx4hIuDvfG+gH7GvI2Nxj+76zuxMoazG7CPiGiLQUkV5ufOsaOj7XLcBnqppTVhCM6+e+I3sF2Kmq/+uzKCTuv0Dxhco9WEF8IXEPVvDvCyFyD5p6EuzWTg094VSt5AMXcd6ZfAfoCKwEdgMrgA7uugLMxvmf4lYgNQix7cF5z7PJnf7orjsR2O6WbQS+HqRr93/utdmC84cr1mf9X7jXbhfwL8GIzy3/C/CDK9YNxvW7EafadovPv+ftIXT/BYovJO7BCuILiXswUHyhdA/aVD+TdRFojDHG1IJV7RpjjDG1YInUGGOMqQVLpMYYY0wtWCI1xhhjasESqTHGGFMLlkhNkyEiE0RERWRAFdb9uI6OGS8i3/T5fJ+I/L6K277lfn/wyvIq78PPti1E5CO3AwBjTAOwRGqaknuA1e7PCqnqDXV0zHjgm5WtdCURScQZiaROv4CvqhdwvpP6r3W5X2NMYJZITZPg9m96I04nEd/wKf+1z3iPuSLyqlt+xv05yu0wfKGI7BORWSLyLRFZJyJbRaSPu95fROQun/2ecWdnAV9z9/+fblk3EXlPnPFFfxcg5G/xVQ9GiMi3ReRzEVkHjPQpjxGR+SKy3p1G+pQvF2fcy5dF5AsR6eRutsDdvzGmAVgiNU3FeOA9Vf0cKBSR6wBU9RFVTQFG4Qyx5q/KdDDwA2Ag8G/ANao6DHgZ+HElx50B/FNVU1T1abcsBeeJcBDwryLS3c92I4FPwdvF3a/cshtxxtEs8yzwtKoOxekJ52W3/FHgA1VNBN4Cevhssw0YWkncxpg6Yu9RTFNxD07SAXjd/VyWqASYB/yvOp2DX2m9un3YisheYJlbvhUYXYNYVqrqSXd/O4CeXD6cFziDQB9154cDH6rqUXebv/NV5+W3AAnOKQDQ1ufp+04AVX1PRE6UraCqJSJyQUTaqDMupjGmHlkiNY2eiHQA0oBBIqJAOKAiMl2dPjAfA3JU9dUAuyj2mS/1+VzKV78jl3BrcEQkDGhRQUi++yvB/+/ZeaBVBfsoEwaMUNUi30KfxBpIS6CospWMMbVnVbumKbgL+D9V7amq8araHdiP8+7y6zhPdQ/U8hgHgOvc+Qwg0p0/DbSpwf52An3d+bXAzSLS0R2Ga5LPesvwqV4WkRR3dg1wt1t2K9DeZ52OwDFVvViDuIwx1WSJ1DQF9wDvXFE23y3/CeAB1rkNgn5dw2O8hJPsNgPXA2fd8i1AiYhs9mlsVBVLcN7b4lYrPwZ8gpMgd/qs9wCQKiJb3GriH7jlvwJuFZFtOIn3ME5SB6c6ekm1zs4YU2M2+osxQSAiUcAqYKSqltRg+5ZAiapeEpHrgRfcRlWIyNvADLfhlTGmntk7UmOCQFXPi8ijOE/LB2uwix7AG+772gvA98DpkAFYYEnUmIZjT6TGGGNMLdg7UmOMMaYWLJEaY4wxtWCJ1BhjjKkFS6TGGGNMLVgiNcYYY2rh/wO5Zqt/+AIAigAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "\n", "ax.plot(eph['AZ'], eph['EL'], 'o-')\n", "ax.set_title('azimuth/elevation angle of 2018 CL on from %s to %s' % (start_date, end_date))\n", "ax.set_xlabel('Azimuth (deg)')\n", "ax.set_ylabel('Elevation (deg)')\n", "\n", "for i, xy in enumerate(zip(eph['AZ'], eph['EL'])):\n", " ax.annotate('%s' % eph['datetime_str'][i], xy=xy, textcoords='data')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The complete lists of the values supported by `astroquery.jplhorizons` can be found [here](https://astroquery.readthedocs.io/en/latest/api/astroquery.jplhorizons.HorizonsClass.html#astroquery.jplhorizons.HorizonsClass.ephemerides). Play with it and be creative! For example, can you plot the change of RA/DEC rate over time? And what about the change of predicted V magnitude?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's come back to our main journey. How fast does the asteroid need to be in order to show up as a streak? Let's take anything longer than 5 full-width-half-maximum (FWHM) as a streak. We know the exposure time of ZTF is 30 seconds, and the typical FWHM is 2''. Therefore, the minimum rate of motion is not too difficult to calculate:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "minimum rate of motion for a streak: 0.33 arcsec/sec.\n" ] } ], "source": [ "min_rate = 2*5/30\n", "\n", "print('minimum rate of motion for a streak: %.2f arcsec/sec.' % min_rate)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the object will be a fast-mover at some point within the time window being queried, we will use IPAC's Moving Object Search Tool (MOST) to get a list of the images that the object is a fast-mover on them." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "if max(rate) > min_rate:\n", " params = {'catalog': 'ztf', 'input_type': 'name_input', 'obj_name': '{}'.format(obj_name), \\\n", " 'obs_begin': '{}-{}-{}'.format(start_date[0:4], start_date[4:6], start_date[6:8]), \\\n", " 'obs_end': '{}-{}-{}'.format(end_date[0:4], end_date[4:6], end_date[6:8]), \\\n", " 'output_mode': 'Brief'}\n", " irsa_return = requests.get('https://irsa.ipac.caltech.edu/cgi-bin/MOST/nph-most', params=params)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "...and MOST did give us something, right?" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\\output_url = \"https://irsa.ipac.caltech.edu:443/workspace/TMP_17w7Dr_2536/MOST/pid2536\"\n", "\\catalog = \"ztf\"\n", "\\user_input_begin_time = \"2018-02-05\"\n", "\\user_input_end_time = \"2018-02-06\"\n", "\\object_name = \"(2018 CL)\"\n", "\\element_epoch = \"58154.0000 (2018-Feb-05 00:00:00.0000)\"\n", "\\eccentricity = \" 0.242496961324341\"\n", "\\inclination = \"11.987721695251381\"\n", "\\argument_perihelion = \"142.423969683856910\"\n", "\\ascending_node = \"136.338359646885010\"\n", "\\semimajor_axis = \" 0.851979508669249\"\n", "\\mean_anomaly = \"236.573041656847607\"\n", "\\magnitude_parameters = \"25.50 0.00\" \n", "\\job_time_stamp = \"Mon Apr 1 22:42:14 2019\"\n", "\\description = ZTF processed science image\n", "\\identifier = ivo://irsa.ipac/ztf.sci\n", "\\notes = ZTF science images\n", "\\datatype = catalog\n", "\\archive = IRSA\n", "\\set = ZTF\n", "\\ generated at Mon Apr 1 08:00:01 PDT 2019\n", "\\ \n", "\\fixlen = T\n", "\\ \n", "\\ ra (deg) \n", "\\ ___ Right ascension of image center\n", "\\ dec (deg) \n", "\\ ___ Declination of image center\n", "\\ filefracday \n", "\\ ___ Integer (YYYYMMDDdddddd) used in product file names\n", "\\ field \n", "\\ ___ ZTF field number\n", "\\ ccdid \n", "\\ ___ CCD number (1..16)\n", "\\ qid \n", "\\ ___ Quadrant ID (1..4)\n", "\\ rcid \n", "\\ ___ Readout-channel ID (0..63)\n", "\\ fid \n", "\\ ___ Filter ID\n", "\\ filtercode \n", "\\ ___ Filter code (abbreviated name)\n", "\\ pid \n", "\\ ___ Science product ID\n", "\\ nid \n", "\\ ___ Day/night ID\n", "\\ expid \n", "\\ ___ Exposure ID\n", "\\ itid \n", "\\ ___ Image type ID\n", "\\ imgtypecode \n", "\\ ___ Single letter image type code\n", "\\ obsdate \n", "\\ ___ Observation UT date/time\n", "\\ obsjd (d) \n", "\\ ___ Observation Julian day\n", "\\ ra1 (deg) \n", "\\ ___ Right ascension of first image corner\n", "\\ dec1 (deg) \n", "\\ ___ Declination of first image corner\n", "\\ ra2 (deg) \n", "\\ ___ Right ascension of second image corner\n", "\\ dec2 (deg) \n", "\\ ___ Declination of second image corner\n", "\\ ra3 (deg) \n", "\\ ___ Right ascension of third image corner\n", "\\ dec3 (deg) \n", "\\ ___ Declination of third image corner\n", "\\ ra4 (deg) \n", "\\ ___ Right ascension of fourth image corner\n", "\\ dec4 (deg) \n", "\\ ___ Declination of fourth image corner\n", "\\ \n", "| ra_obj| dec_obj|sun_dist|geo_dist| dist_ctr| phase| vmag|match| ra| dec| filefracday| field| ccdid| qid| rcid| fid| filtercode| pid| nid| expid| itid| imgtypecode| obsdate| obsjd| ra1| dec1| ra2| dec2| ra3| dec3| ra4| dec4|\n", "| double| double| double| double| double| double|double| int| double| double| long| int| int| int| int| int| char| long| int| int| int| char| char| double| double| double| double| double| double| double| double| double|\n", " 130.243140 10.591412 0.9936 0.0077 0.3695 9.6662 15.69 1 129.946762820000 10.364530940000 20180205231528 517 6 4 23 2 zr 400231532315 400 40023153 1 o 2018-02-05 05:33:25 2458154.731539400 130.387829180000 10.796160510000 129.508055050000 10.798784690000 129.506884090000 9.932328820000 130.383994010000 9.929949140000 \n", " 129.946506 11.101031 0.9935 0.0077 0.1305 9.6479 15.67 1 129.949891860000 11.231928510000 20180205251296 517 6 1 20 2 zr 400251292015 400 40025129 1 o 2018-02-05 06:01:52 2458154.751296300 130.392207100000 11.663277090000 129.509976810000 11.666108690000 129.508681130000 10.799672070000 130.388448090000 10.797038270000 \n", " 129.625980 11.645721 0.9934 0.0076 0.5217 9.6703 15.65 1 129.950695400000 11.232704280000 20180205272106 517 6 1 20 2 zr 400272102015 400 40027210 1 o 2018-02-05 06:31:50 2458154.772106500 130.392985830000 11.664048400000 129.510776420000 11.666846170000 129.509495230000 10.800441100000 130.389245740000 10.797862500000 \n", " 129.618211 11.658854 0.9934 0.0076 0.3954 9.6713 15.65 1 129.648394110000 12.053413850000 20180205272593 1563 2 3 6 2 zr 400272600615 400 40027260 1 o 2018-02-05 06:32:33 2458154.772604200 130.091191430000 12.485962890000 129.206230010000 12.486622890000 129.207327970000 11.620350100000 130.089107560000 11.619537940000 \n", " 129.304417 12.186967 0.9934 0.0075 0.2490 9.7350 15.64 1 129.070327660000 12.285688820000 20180205292488 517 10 3 38 2 zr 400292483815 400 40029248 1 o 2018-02-05 07:01:11 2458154.792488400 129.513283440000 12.718624570000 128.627359850000 12.718548800000 128.629102050000 11.852312030000 129.511895640000 11.852222940000 \n", " 129.296470 12.200286 0.9934 0.0075 0.3748 9.7371 15.64 1 129.649259050000 12.054213010000 20180205292986 1563 2 3 6 2 zr 400292980615 400 40029298 1 o 2018-02-05 07:01:54 2458154.792986100 130.092073220000 12.486767920000 129.207114190000 12.487391270000 129.208189100000 11.621136410000 130.089970620000 11.620309610000 \n", " 128.976880 12.733908 0.9933 0.0075 0.4280 9.8423 15.62 1 129.071104030000 13.152379000000 20180205312789 517 10 2 37 2 zr 400312803715 400 40031280 1 o 2018-02-05 07:30:26 2458154.812800900 129.515614710000 13.585218500000 128.626757470000 13.584987180000 128.628208840000 12.718695140000 129.514139550000 12.718795740000 \n", "\n" ] } ], "source": [ "print(irsa_return.text)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is a loop over the query result from MOST to find out which images will contain the target as a fast-mover. Note that we want the predicted magnitude of the object to be <20 since that's the typical depth of ZTF." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "limiting_magnitude = 20" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We also make another query to Horizons to ensure the object is really a fast-mover in these dates. We do this as a double-check since Horizons produces most precise ephemerides for NEAs. Horizons also tell us what the expected positional uncertainty will be. Obviously, we don't want the positional uncertainty to be too big -- say, over 20\"." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "max_unc = 20" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And then we print out all the metadata that are potentially useful. We will just return the first entry (or it will be a very long entry!)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "hit: object 2018 CL\n", "object name: 2018 CL\n", "image date (UT): 20180205\n", "fracday: 231528\n", "filter code: zr\n", "distance to the Sun: 0.9936 AU\n", "distance to the Earth: 0.0077 AU\n", "predicted V magnitude: 15.69\n", "ecliptic latitude: -7.4574757\n", "positional uncertainty: 2.31 arcsec\n" ] } ], "source": [ "for ztf_entry in irsa_return.text.split(\"\\n\"):\n", " try:\n", " ztf_entry_vmag = float(ztf_entry[62:67])\n", " except:\n", " continue\n", "\n", " if ztf_entry_vmag < limiting_magnitude:\n", " ztf_entry_jd = float(ztf_entry[270:287])\n", " ztf_entry_pid = int(ztf_entry[185:199].strip())\n", " ztf_entry_fracday = int(ztf_entry[120:126].strip())\n", " ztf_entry_filefracday = int(ztf_entry[112:126].strip())\n", " ztf_entry_field = int(ztf_entry[135:139].strip())\n", " ztf_entry_filtercode = str(ztf_entry[173:177]).strip()\n", " ztf_entry_ccdid = int(ztf_entry[144:146].strip())\n", " ztf_entry_imgtypecode = str(ztf_entry[235:238]).strip()\n", " ztf_entry_qid = int(ztf_entry[150:152].strip())\n", " ztf_entry_sundist = float(ztf_entry[24:31])\n", " ztf_entry_geodist = float(ztf_entry[34:40])\n", " ztf_entry_dateUT = str(ztf_entry[112:120]).strip()\n", " \n", " obj_this = Horizons(id=obj_name, location='I41', epochs=ztf_entry_jd)\n", " eph = obj_this.ephemerides()\n", " ztf_entry_rate = np.sqrt(eph['RA_rate'][0]**2+eph['DEC_rate'][0]**2)\n", "\n", " if ztf_entry_rate > min_rate:\n", " ztf_entry_unc = np.sqrt(eph['RA_3sigma'][0]**2+eph['DEC_3sigma'][0]**2)\n", "\n", " if ztf_entry_unc < max_unc:\n", " print('hit: object {}'.format(obj_name))\n", " \n", " print('object name:', obj_name)\n", " print('image date (UT):', ztf_entry_dateUT)\n", " print('fracday:', ztf_entry_fracday)\n", " print('filter code:', ztf_entry_filtercode)\n", " print('distance to the Sun:', ztf_entry_sundist, 'AU')\n", " print('distance to the Earth:', ztf_entry_geodist, 'AU')\n", " print('predicted V magnitude:', ztf_entry_vmag)\n", " print('ecliptic latitude:', eph['ObsEclLat'][0])\n", " print('positional uncertainty: %.2f arcsec' % ztf_entry_unc)\n", " \n", " break" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voila! You made it. Note that it takes quite a bit of time (a day or two) to loop over the entire NEA catalog -- this is because Horizons has to calculate the ephemerides for each of these ~15,000 NEAs and this is slow! Therefore, we provide a pre-generated result file for you to proceed to step 2. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But do take a few known asteroids and see what you get for them. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 2: compare the catalog derived from step 1 to the source catalog generated by ZTF streak detection pipeline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's read in the result of step 1... (OK, we cheated, it is pre-generated)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/chummels/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:4: VisibleDeprecationWarning: Reading unicode strings without specifying the encoding argument is deprecated. Set the encoding, use None for the system default.\n", " after removing the cwd from sys.path.\n" ] } ], "source": [ "check_fmo = np.genfromtxt('check_fmo.txt', dtype=None, \\\n", " delimiter=(12, 10, 6, 10, 10, 10, 5, 10, 10, 10), \\\n", " names=[\"object\", \"dateUT\", \"filter\", \"rate\", \"sundist\", \\\n", " \"geodist\", \"vmag\", \"ecllat\", \"unc\", \"fracday\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And also we read the source files from ZTF... essentially, these are a large number of \"streak_qa\" files, each file contains candidate streaks in one image. Apparently, we want to loop over them and see if there is any candidate that coincides the position of a known streaking NEA. We also record the rate and mag of detections and non-detections for diagnostic purpose.\n", "\n", "For demonstrative purpose I am picking one of them to show you how this whole thing works. Let us start with our old friend, 2018 CL." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "for cfi in check_fmo:\n", " if cfi[0].strip().decode('UTF8') == '2018 CL':\n", " break" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What is the uncertainty range for this entry?" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.3094 arcsec\n" ] } ], "source": [ "print(cfi[8], \"arcsec\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is very precise! Well, we know it should be small because the orbit from JPL has included ZTF observations.\n", "\n", "Now, let's obtain the predicted coordinate for 2018 CL at this time, using the knowledge we just learned in Step 1." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2458154.731528" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "expTimeJD = Time('{}-{}-{}'.format(str(cfi[1])[0:4], str(cfi[1])[4:6], str(cfi[1])[6:8]), format='iso').jd + float(cfi[9])/1000000\n", "\n", "expTimeJD" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table masked=True length=1\n", "\n", "\n", "\n", "\n", "\n", "
targetnamedatetime_strdatetime_jdHGsolar_presenceflagsRADECRA_appDEC_appRA_rateDEC_rateAZELAZ_rateEL_ratesat_Xsat_Ysat_PANGsiderealtimeairmassmagextinctVilluminationillum_defectsat_sepsat_visang_widthPDObsLonPDObsLatPDSunLonPDSunLatSubSol_angSubSol_distNPole_angNPole_distEclLonEclLatrr_ratedeltadelta_ratelighttimevel_sunvel_obselongelongFlagalphalunar_elonglunar_illumsat_alphasunTargetPAvelocityPAOrbPlaneAngconstellationTDB-UTObsEclLonObsEclLatNPole_RANPole_DECGlxLonGlxLatsolartimeearth_lighttimeRA_3sigmaDEC_3sigmaSMAA_3sigmaSMIA_3sigmaTheta_3sigmaArea_3sigmaRSS_3sigmar_3sigmar_rate_3sigmaSBand_3sigmaXBand_3sigmaDoppDelay_3sigmatrue_anomhour_anglealpha_truePABLonPABLat
------dmag---------degdegdegdegarcsec / harcsec / hdegdegarcsec / minarcsec / minarcsecarcsecdeg------magmag%arcsecarcsec---arcsecdegdegdegdegdegarcsecdegarcsecdegdegAUkm / sAUkm / sminkm / skm / sdeg---degdeg%degdegdegdeg---sdegdegdegdegdegdeg---minarcsecarcsecarcsecarcsecdegarcsec2arcseckmkm / sHzHzsdeg---degdegdeg
str9str24float64float64float64str1str1float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64int64float64str1int64int64int64int64int64float64float64int64int64float64float64float64float64float64float64float64float64float64float64str2float64float64float64float64float64float64float64str3float64float64float64int64int64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64
(2018 CL)2018-Feb-05 05:33:24.0192458154.73152799825.50.15130.2430410.59115130.4917110.52376-2179.933843.943124.290355.1746612.91693.91594314.595242.96122.7336.78755582461.2170.2115.5599.291--612957.6*----------235.370.0----136.0615-0.05880.993552869965-4.89240910.00768971056396-5.55238570.06395327.287358.81488170.266/T9.658769.673.20.071455.334274.9696.00716Cnc69.18491130.0804119-7.4577146----215.72888129.09429621.53361001240.0003541.321.8952.0980.965-61.0812.722.3091487.72610.0070343111.58405.530.009925217.3296-1.9118911799.6624132.954-3.7651
" ], "text/plain": [ "\n", "targetname datetime_str datetime_jd ... PABLon PABLat\n", " --- --- d ... deg deg \n", " str9 str24 float64 ... float64 float64\n", "---------- ------------------------ ----------------- ... ------- -------\n", " (2018 CL) 2018-Feb-05 05:33:24.019 2458154.731527998 ... 132.954 -3.7651" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obj_this = Horizons(id=obj_name, location='I41', epochs=expTimeJD)\n", "eph = obj_this.ephemerides()\n", "\n", "eph" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A simplistic approach is to loop through each small folder generated by findStreaks and see if there is any candidate matching the predicted RA and DEC:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "detectionRadius = 5/3600" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "YES! Streak found in run_merged/ztf_20180205231528_000517_zr_c06_o_q4/ztf_20180205231528_000517_zr_c06_o_q4_streaksqa.txt\n" ] } ], "source": [ "d = sorted(glob.glob('{}/ztf*'.format('run_merged')))\n", "targetRA = eph['RA'][0]\n", "targetDec = eph['DEC'][0]\n", "\n", "for di in d:\n", " fileName = os.path.basename(di)\n", " streaksQA_path = '{}/{}_streaksqa.txt'.format(di, fileName)\n", " \n", " if os.stat(streaksQA_path).st_size > 500:\n", " try:\n", " streaksQA = np.loadtxt(streaksQA_path, delimiter=',', comments='\\\\')\n", " for streaksQA_item in streaksQA:\n", " if (abs(streaksQA_item[1]-targetRA) < detectionRadius and \\\n", " abs(streaksQA_item[2]-targetDec) < detectionRadius) or \\\n", " (abs(streaksQA_item[3]-targetRA) < detectionRadius and \\\n", " abs(streaksQA_item[4]-targetDec) < detectionRadius):\n", " print('YES! Streak found in', streaksQA_path)\n", " except:\n", " continue" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wonderful! We have one hit! We can even take a look at the cutout to see what it looks like:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "image/jpeg": "/9j/4AAQSkZJRgABAQEASABIAAD/2wBDAAEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQEBAQH/wAALCACdAJ0BAREA/8QAHQAAAwEBAQEBAQEAAAAAAAAABAUGAwIHAQAICv/EACgQAAMBAQEAAwADAAMAAgMBAAMEBQYCAQcUFRITFgAIESQlFyImI//aAAgBAQAAPwD/AEz6H6zVGAgxEG4JlpWdRMpQpzV6sdk2L5LnDr0n1qUysqtPn3pzgXFM6stnrIMpTDavBVsAWJMsgwUa9L+6o3S/mnPPKpRfZl+RxQmKS1grTO5Sgcre0E+QN0PKoJ1zjT+UWX1Gj1qlkx4hT6bVktNUgt8sUSGmpAHYyVnpqZoM+Rpe4oHR1VpPHP8AnOTRtNT1iGN8tBnk0CNmbQtJbOmhD4iLhmLWmlFtj51YkmkprvLqL6B2UFSmeTxWzYH3oZJqwc3oSS6lGK4sN+jP08SqDkN/Itiflj/sozUqXmIjWSy1CaPuRNK/ODUzOIiZ33W6GIfr2jX7mcHosU5qU5vLRAZxToEaUHY6ko1Hk7F+FBYcBAv09jJ8tG8eojlrbbToQpFFb23ZrVzAQkcx0/VVMxpAWabnvOhBNMxm6OofO1bCkJudG9C3Io1M6uuRfFl4gg77XlrkIOLC0LBqjDx6M/zn9eh3RwKbFdCVr1IbEA7UuwcLFejXAG1G1E1yJwaul49abhFrsncU94La7kMTMpFn9SlEdda746JsS86CmFy1NOZQpVIqPGdRk05sWxeUle0EI078vcdTVkukeSPiJ1bKm6QdgMseud2lfKamfL4xC3dbOS3iVnKjWpg7kzI/AGrdIC/MSOZrnjNPeUp6FMfdV9ilNZccOmzbA+QcjQN6WNz7bk6FYzB0iu/y8PhmMvLPm0Wqceqi4SNcppfe5PZOOvee6qgPmz5+774vo2HGWMvKXnokBFopzk/P4155HihHfHz/AGChnhbNNJb0KMuBG49Z+Q+c3Fr/AFSvT0uoWkzGj5lWg9Gi59WDOQmr3EvKvrNVReLlp2m99ck+S8iZxzCWuo7lOx10qTkbIbk6gEdrxV5iJvEKnog1Fpns2Vpaf+Jidp9QXY61O305H4ytSM+87P1LRMyZTkEkS+Mz7DHKcwOfcjN0K1HkURtHzsTy32+zCll4jd9tEKLCUuLXuaLQ1A+86k79L8EGaWcpjlRPsa68z6CnJn3+B1/NQerUldEptTvaCxrbhaQld85nzvToxmTFpvpZ1zHJFVXs4Ts0aepJn52T6QD7tH6JtisB5mdHA/641NVXz1NRaO0lH/csz5YfAPU03GtTyO85JiTED1qtbQ1gQaHghT55u54pUL3gVivHD0WSzPlaOLPpzej/AKI+T74hi1j7EmDY6ETyrXGSz15kUR9Z7PchjdR4wXBPG4xR8/Fh8+u1Ibv8EbT0Kj/GVWSqFoi9V9Mpzf6lY2pPL0Q1+4lz6FIUS4/l807bWVBPGP2ryFJvTU478atBo6bT6qd3xznGp4kW23ndBopjA9pZ77ChU8cWs2YLYualTSHHMr+l8qdJ0ZAeRTush+0/5mFl1m/ATPNF4iOfEy7lS2mKkeF/CO5/LpRW6Sm3oWtFTbVL7HHFao6md1TlMz709BRUkrTqs/8Aiecg02nEM2e7WogJ7M8Z0TGa8RUi52ceFOu3don2XNdahFNYPfQVYewZ3dJS2abb5y4o7IOa8FjVzWXps7CM8CpYJo2fsDUkaLZ90s7z3CCFUmVTLXxXcTUcsUm5icuBT6qmy3pksSmG35Pk6DinokZLnPEkt1JO0vVyWq69Ttg/YOhBlJTGJ4QO5OwdhJmiRr49oeUcyNpHzhstOSO6xWqzAd5OivURoUH6wL57AgwjM283Bl5hBLNIdLTK02po71ZZ2DNgrsB0NtmYh4OPHnIqMMhS6QTAGF1RlCfnURenSzKzDEn5D0LXTXVSpEj02L3kmBPU78qoF5qc3bIJznNJ8HvaMqrVoxnaksvCkdasxH44K9EbqaR6VSdVotB9jxHdJPH5yE9Dn7kWxE55q0s6YFgn9dw+gi4orxvdGunUfL0pjxxTK3ZIk7G5APfJwyaneR5z+O0rhPtC785Vne6tDuY67ofDctk0cqZNc6aGpJoLlO/1eTVn2aN2wP44j92bznWSpbiPiefX4F7T6nmvNHrflPLrs18rnUtHNa5lZ3fKWs78ccH9mxs6z/PTGMPt8I8ZP1fiKScvn5jENzhDj0lxINDF9saFKc8zyhxArZrqI0xQpeLr6mnldyNloKCT/wDxFBlWaOiVfmS2ajDKYWBcJgeY10eh9Tyr55JVpqqyG/wrrud1GbjR+OXaIcIo15Sd0ctQXuvPvNzU3WkiXNK3ReqI2sTKSkjUljITSqqa/HUHO2gTOC5du/8AhqTGJrFLXZGX7FSVUYWSlYAFp/VfTTtB7Pi2QPDyQFRowRLbLlzm7Y98scZmlSVVzyPyT3YH9aZ31W6DzGgxx3vWMkzQ++pOaK9R0z3lUXZZuhc6fZ5cnGQlWQTOZ3yh7lO0ZTUScttqUNimWcz29Q4aCsyRGutnoHHnhGIsSuoXTcqeUKPa2qYp0Mw+CpEuk5z67nlnTlHyO9OpGVu2FiOqBc7lWiswZ7jPPSH3Qi3Uw7oqHvs26HPSKvfuYFwV1GM6dKtFTdCUBLWPr1kLlKThVaKlRu9j/F495SYoLVx+uZ3fkxOpVmeanUcRUfUAHP8Aptqkl5jHVXvy+TMMfTamIT4GkJG6AG7Qab87qLVXrDHqRpiTvs5/uRShe0Nbc7h9RZWQuIdaGJEvTBe53zJr0aFSiMNcuf2WVQZpO42S3sq5BWqc7V45lFKtGvoR8/1g338Xrm1UCehlNStGwF+snnkGw8NL9HK9ME1UYb8FL99geLKNwll17Sm6zxrLOIam8zTaX2NPyEThBpFXF/IjOW6nOn+2DpG8vLm2n36Y8ulji8s0Fpac5Yw81GCVaAbleCUa4mFnRKHPtI1/+o78vWxfAUpji+g4Iu94ci5hhpUDIY95ZIxwuUXrK0nQ6C/bbNwjTVeISkuWtLtTfa39HWgr5rkx5NNHkeOZ0DsWSsg0gr3OAEeoEr5aQf4abxKAXpDHjkuKfOmtU7zmjQdbveWZkh7rw1n+4DnnavVPTbaQSN487WcHnUr4nwXEHVXm493vmVB6SioVHXEHqk7x35QmxldNmJGRR4j1m01MxWvq8p5hbLMg4d7pT/76We0AqU2SvRRgzE+POGwyREeVysqnIpMg/j5XR5+T1oM1ofFZKvtBGdqdRdnT8ik/q63TTSsdPNTVbQKFN6x90d7la52+0VakA5CLSNZgkvRXtE7ak6Ja9Azb1D2AQCv6LUbu2OmHS8PyUqMx2vQXZ0XUySnU0B4vtiWOFEgejAZk8j7Ssd9SaIa1O2LlCRRhwhHDOjgu8oMrROIufmURgz19b32PGZkJ1fXhf1N6Rz1LYaMmcmT6kb/eS/CLsciTGsuJwTrjL9Iv3+RUxte+gblROC5q7HS+oHQXqcgOcemRbraxqD4U08nYtISaolB53DC3kYUeK/nh6/pMPWJTwnyFLhM0stovkT+2nlXsI9SN/IZPkqXBOl3EbPqArRneOI97aKedWx9oTOdDRuLGZz+3PNt9ohnyPGGfKlSgrouDPAd7UpUv2SUx6ht48/ZnwqYIIFtcil6xw6CzYXy1XOS+k6LF+hKDSFC87APOjFYmaTqNKeeXhI1EPq1XXGA/8be8JVIgtCxVoRLsKu/nWqM+MexFrqJ0gOc5ykeRU8laDPSEiztp/eWkXPa+LfjLHeVz5V3pB3b+zKeRGYmS9UrEVormNGi66lmGSoHg03fdRD0DfbzTugvRUOJq3ppK4aU2t5LRrSlKLUge3nyr1dD2pD5byjGa6RJCzrUvQwY8T9Ta8UgrXKhp3GU6nfXNwX1xd0863naWaAuxPYlUDTcodOurInntew+c9NYcLb0dGZr3G6NVANRhaBPLn9FRheuRU9WZVbQzasTioW11sEXzvNrmbG44Fxo7zPeux7f5J2U/Y8+K9oilnIHoaEynjl387QtP89NOBsAkkY1Ins/z4/T9s1u6p/BwjpQ+AyvzKRK1lf8A1U+rytawyjGf95qty3iT/CFHoWmpwq/Nqeh5Cq/mXkuIMHmYksXTY4M8RHltyETOd8iR0A04HvK8bXuuExLFW62lMYYi/q2eu/OvHURrWp5Xcg4DN04Ew0wQdVnJDVl9xZaOdt6XO0FGa4rdnXypMfTGU72GNuSG06PBe/YThssYVP2rmCaeT11kl1KSieUnZS5R2Mj78i/G1HK48HsjJTNRt7sG1n0v3tGU0x4Gb+R0aoyLA9npLJafIZhwUtGdUS4cWtdAmMktC9JZm0ZhpnNkvlFm2KJ4OM6eGBPT1MzL8dsQDIzqSemLSbtY5zga7suavQWo5um7bVqsvrxPe81mc9bCh3na20rruRIH2VC6QyDFXqvMv+CR6x+mTm1cZr+Z3Lb8zRuaDX2Ks0NRX601FQW1E7uXkNS607N0dlpvyvXSjjhUUbQo0GLxxOiOf3hlKI5vhKZ/GjRZVpod2OXkSdXOibx8vVtQ5JECtAWPS7NV+60tdHbkUy+8KLrOzqvJ3l55pTkimFxad1C07XmjGnMoMjZ4mVlaEzQXuF7LH5ujqSsuZrjTcq9R7JVrn1nBu1u+4i6RhNZZjP8AUKI5Rf8ASlinLQfPXF4fS0DvM5NcZms/XeDmarrjdqb7l6JnPNAtnb3lvkcwQJ1rnSx6VXtNP++lbSVXmoU2WP5ttUvQva1ye95ys2tV+S1VBsenDGhlbzzmqzMy/U4H36nJ5q/rqFh0ZqN+KdSjHb/rSVowXWJLih3V5VHlzm5MlDjCpvaNjgrz3r4Mt2rSqUU2YmgHQfsp62pWe+1OKErGrC34Xz+aEymTPmPUn5sMkEx+lyPSSkSf/wAHFeXh57RHiMmQ74YjS9EpWmWgotHLUoTLO0y94youso46LVF+fomgYFIA5kLuY+KDXt1scki1sqtWPfjD1NOfN71R6cwzSOjf0dQubqKzHW3ekPVtAgB2nsq/h+aqdyc5DpQIGd8y99BfvMvahHP0UbCeU7XbDdXJn46zzAQmizzbNGvNlpt45piLY/4YwgyvNStuUK+bBiewdHBnVBMRXtHPoo+UxqG5HF05Trq9qkKtQzk8DDGg/wBNsl+Vmpjj7B6+I4aKK+dqg0nsd6VaSvccbJnnHXaHY6IPqHRUCCVfEtF5sVHzLx4sp0pc/ZWnvaPr/k2jmIw4afTRN+1Bc/DttlmHEBrQqLCBn5MpjTdZD8icqv4cwZOFZUFNqTpVX3QvVLNKHIhfZV0OZXc9Yyyd0o6knpHFcjqROOJFRbO8EjrSK82ZyZJLIMy/kCy0bTBu9cUeGa5+/ZM0vRU59NLIYqvLR1GfA+0k94lNic1dB8lVX5EbRxUby9Of080vGEGGKfoW+BQV5ItUnIn3mZUyjAo17fj7Nf2MVzQnC340Q2tuERQbhQlO1XbW4lTQcs+KWkNtdRFNnM1VilJXoCn5PVzOrE2Sgjn/ALB/aIs5VJBeYYqPXmf1+k4D+lmOdVVpC7UwsJ4z6+2vtMkjNxKNKl2jfmQgL9J1aZBd7X9UH7FMPPx5p9dRNn18G3W0Kun5rw3SvU/j9NYWW9iEs3T0Yi4Ki5IwGtpBBVk3mfV/qJpT975T+fvnjMsR/mTKfHWFT+NGfij5Ch/E3zh8lfF/yPDura6Fq1MV8m5j5/0WDpv5H5CwmqI9h8f8X/FC+Dw+o+MSVB6KjpiI5ENGPpPsBGqcjvdLUL5R3Egu1DGXuLD4pq/WVPX9FFoUgaRot/2HySQoRxlWeHzRa2LCz5EeXQim5fePSzwq+QFmC1M1TlVlrsIE5qs3Wo58NhGahV7r/IWe46dKge1orOG/qxfIa+3X0TZswhWUp9C2O0ECDXWgDKLSz/jHkyXtcwzppRXRkcRzkel14sGpGd/tzc9utT9Mqa7d0ZZzHpBVHQwEzFnZ9r2nMlIpJaVuEu1Y7al6KZI86pPvUn9HQYUIWqahWflpSXeolhT0UpPo/ToM69V/fMpoSe63xuBNpLyajFSFoi7GBqrnHhpibU4kTpOn1nBKCiizR7ExRoXGa/bz/XMaG+ZI6TdEE7gOrN4s932iIcpwvbzTa8P/ANpceD5gSvSto6iwBGOmvpF1IUNx5cFlTmtI5ukZ0zn3PQGXacsQe41GjQAgrlToXfRUmHnjLuOZphwEbMaYXD3UOXcXS8T08JMo6lXDOzLguAGX6Mw7Zzsczl8Iq2dkLMox8zEmMwEeU05My0xCS6C3ZUWg+x2k2UWF1CI96U026t34h31P7C0gkfriQUcwKaKE9eFKQkQDu+T0V1NDJz6Ohn5yfSz153tCLOZqyM2Gzn3mjYslIA7cWUsdudoSLS3Lcamg43IU9zvHydC0A49MwI+em6D86XWGsJmxYZWkTM93qu59+DTWInhXNK1Jfof19Zew9w7SJlXgTe7c8cjlinLMKhzQa/UmKs5nzzQuLKRVF9OZ2pRSPauSV2iXpKs+ZCszbZI/1o87JjpUfHHc22kvyqnzTYVUai1azL9bK0fTyzMd1v2IJ8/JrqWNYFBrQ85iof8AkCB35QkyFpEJ7xJzigz7kHhTnY1RjwVRJZmbMl0wzOOdG3d6UzN3RmiPgdfnV1cnQP8AmcWcQHT1lmBcy87RrOtLztBKhVdAlHqN+0m+04rbyFVHRCSQhvNTOOJK3CdUb3yBCeIwyHI13lufKGiDcqqFzrTEelebr+Kv+af1hSCMr4hdvQZUNufTTesemopIwFaB+DNIPeLc3Rvjc0dUyq/FZuO6CW3LjC4aS+QHxlsBQHT66GelazHtDQNy9Vo4/MXP6R0uVBBq6B3RcLlZ+bB8DOKgYDsWhl5EwXDmbnOT1QqZvj44VYApwBu/raaujdgu8oKdZ3qUnLlT9Hx3kXM1pZmj+WpFJzfznRbLd7GXEFKDjcoqVmzLCKhKmEbJQFJYGDY5+Dai1T//AGKytS7+YysbLBROdo53l/mn1oVab0EhqcvggaL7bCBdPOSWzqNPtGXAFK9a+ujaj0jHrzXXbtGBRg/YY1wuPPref+N8dttHrsr8Eza17ZqzZ2z0HPxblEr9B/LUdC1MSqaA/wD0u+Utjo5EpDVe8Zgel+VtWhARaYlZidn0eGva1jBnTBWM85lSlUUOYNRtWbpD9mNH5ZzVLrRjkEeJ6raPTi5PUoSVNGw/k785Jjjn/wDo0CMIyaBmhTzlWs3Pa50IN2STUrt0ylpDShPRAdzCgXmVJkuRpI+fmR1pw+VSrG9kd8pQWRVs6RBodpT5k1uiZvj0dAFNpjMuEa0Y8pnf0A53aNvxokNiZ7ruNGqp+EavQm8nltemjcry/dBuFVhxdWoSkboui0wv7Jy7rtelbRTtUbs2Poejd0Z0Cygh3LiNNtR+Xwzv7uMz7/ESWhWpDKIs5CTyGTyXmu1TRkO2QpaX/wBWvxhZWVM98YqyuzSbTOSckK2Owpvyb+QBnIxCpVMxcQDVZoVpnuksz5ia/TyBKMu5mdDeofbW6M3KJEkqKyej+aQz8+ijUkssc6etrv8AOWXBVu7h1HmWP17YdbbxVHvVqE5k9uaACHqehzquuYAwLwK789M6NXY2OVdF4pRKWbZ1MNqz1Mz+U0FtJT4w8tqKLS7xS+6nnik2zZbaZJapKZr3UMyRy5XI0ltAFhX5FcYQ/o0/DrTOjm/YmrnEw7NUlb9DCLQ/y5ojMSeuy1rzLfNI5s57oFU+M500z7B4mK81ik0fegWDzodtaZ75a1ryRYtVosVGnQA+kFG0wpra7bKmjxhST7hrGma00VMUthXzqCbEcNZjN6Wz/fU09uBFRW1irIKlzjhXiW2VhSiZaXsJ71ZfUR88Gt9M1oJT6dCoYzueMUhNOy0o8zSUp6OlIpnq/Huk80FieStOjFaoxfYiHaC1nzVdXEjr9uy5yUhatHgVmM2OFCAsTqrn86GpqZ1fj/4WcyyeaI6jxoJucBI50tGcbR50ui5NBilaAG4GmT6mn88u/wBr33G281GVd5tntItxSy6K10XrCtl+L5zWayj2mkypYN9fRTx1Zi5s/wAcIP39BmimQVWOO7CEF8Boll9yd1046KTar/1tJW1bhcl5ovNE3TZiyNNExwU53rdkqItRqU5H2Z7FyYx+09czH/8Arr1FIEoEul+bVv5nTGuK/wBHn7vcfYT5XZqbtUMHjVyeL8Px8MwpHve+GLOnnbLSUJG51sW6zrFKKnvEollLOTthIbTlCJ0F2xm6uoKvHENhz47zUpW3QGZbq55jYV61XgDK5pimonyev1QsZiY9lg7Geq2uj4VEoyZ9jLljz/e/NQ9o4wnLzFs5OheKcq0vPGXbMHWx0XE87StFaXU8ieh4mlfROxotHQmUUZNBp2ga19mepnRZ+quR73HCiDnsnqmpxqvadT2kt3dp+6qmFROzdTvl7SO2ZZ8B6zrpm5rMi1LmZWa1eFJIBbWtNesnXEBP0fsNPiGHqfEmvNmniabpAL3VzsByO9H9o1K/kkanqHk5zy3PJN6IDu5dJeeQnsah2p75PDpbmnd74CH1nUr/ANz0KWirLcaWVJKD0JsPTVACyy8bpMmngR88tWs0E/vaL/TJ/J0iO2ax5VGEbmg5m9STqBjWlPyzgH1kdE5VQ1LGgZRWrU5aSVD1FKTI4pijGGTgs8ziCoVrnBa2ch0gUASD7BqcK+ciFF++hWKjGV0GmozApVhHWYIyj6tDqMryKKdGs04F9d5ntxWSX/Hp3oMINVUmv4kWcSECVtaxKx3Ird9aV4QMi7nGGWbiPl2DzM6rTqAPZyqt5ptN3SJWGjpeuQ+SWKs3l/BTIFjSesh42GxqR3GbUkZhxM+T3ntXleHqs80ePP8AF9PnEVsaKGj5nGUXVnx1FbGg0JuetDpY2jVgy/8AP8L06Wihc0OkGN0lJNrq38qJcrLGMSzDlpWc3xoCJ1SlGpp6nNu0mnbgcwsBTpZlqEb60ytxoFYVfn2seWGkinoXmbtZPTjViM3wnc0s+H1x2V/fSWQnRhcDACp5QTNzsWIN2jtQXepk/Ru3Cz0k/LAYXGwJJWL9/Pt2fI3dpYK16dBTTYiWTzLn8bU801ih/b1cjBz0WAxNpZj8TovozcDBPee18IEMflKcjUyRoGS2nH2Vs/SmdR1qEdm1X/noJFZNkvRmvxe+QoecMOL7gVWQxYa7FpuDv5j73tVFcQojGlTG8hO1ly71e54xmlY4RDVpRSaBXplk23qyFZSdN982dt0FZ1a48fqYpqtDPxv9dyX8dodEkc6F2ECSxYUoT/RRoiiZ9vqiwwFXIp7B79LnL6dSr7T/AErNyAcbFDyrxBXZznnoRMBPre1LK9Xuq4I8NgrE6PugA4YQV8vAZzcjRR0+aYkpg3niLo8WG+m+M+83xHdFCC/7Ph2Z/TZailZVP1qXoMxovSfZWY1NFDJNWpcUtObDp3OUl++22tSi6iPXab9FjmapRafLdRm8t5oK9tqTnUSvtjSaSpeLl4NgcigK5oWI/c2LIVUsU7vZ71aOWGtfnsty/eemkhwqvs3xvMIRDKeZRufsdKo1DpPA9cZ+TUakMeuf4bpPubcuW1k1GWqPkI+lmoxgcG5QLlFQq0Ycq44ebSi/Iqc6zHM0upFX5p/8V8eByzycSXXDx9Kx1Ike7/8A1EOd6yky7cVk6JWSLwWku1QTdhQqDEz5++CLbmLfzu6k6Vhl3zJmd/YkR/tBdpdWGdCWmhBIkGvnGCPRXXu1bNie0ZkdszeyWleUe2oi9G7GTuRjush3EZ5H19Arxg8zPUmTVLwbb0Yap42sqiXQjqawK3jLTa8DBmD3YRS9QY7oPJUvF8ohKcsvMbMvVO3W7lD/ABbT8+HaFyD5SPmDehtBcdPDd/qSLnCtyKvtOsxIhlH3lmNQaEdXuUtAsWZnh+PYkn+RNF67IdA9Q9qx5059zkTa5auy7UKtH/pSLP0amHznNRr1ssV/3PppwKVCNo0mOxl8/bl95rHu95Psvn69CM+HNcDjpWK+GIqx3fFWdnG+rrRWVR1GQymZkP2tbvB9sWXcbUTQmZYmw/a8r9VBJWbnISPkXhX1tNdYcZo8JOW676WqbO3YPpKqZw35lxBOw17IZdqDkXa63N0KKi0z69DP+XealhtMn4s45Y9tkaoy42gtTGrSyU+LzSzys2mF9rR+Xxy87ZGReTfhy6UIIuQJw1EJTiwrUPlHImV0dbdVOs/e0LObrup/yZb+vINKyUzQVXJ484HPOqJ8TT/w0Cptyaetj8xe5+wQwGmhoSvKE44ejX2LERSCXbTFY05V6jAnwjSli8mzw1LyQqNN5WUKlOWPLadPz95b9+lfrVJoZC9lqsO1WuzMzrPQHYHN8q+1xY6h+taamscpoNWVtJ1i5FRyjo9O/XsxFapzf3UXaSXfiQMVRMkiTIaBufydmWrUiddItSc9tECSEAUnO2JUVid0vS9b7aq0lxNNHsULnXPsafSVzzgNSIVZhGZfb0o/jifoVWfdalXYM1cHNWNRTsBSrWZnGZgCX/ppU5k5A6cLQ+6gFf0sGw915A4Qh0UCeUnrQNTJ9Qt+wE/VYtWCr95CCPP5pxValRo5ApyvT0cl5TliEQixe3ANsrwBQ25ik3ngyaMhi1h5y+eLJJYfo0p6zgW+6wKCCdCdKuTGgYlzl/F6qV534w937RJjlVangVXoqMHl9NTsHjYiJjej/wBwEpegX0GhxtG6k1ARrswV9SSSOXm4X9tCcw0LcKduq7NXuM/rc1WWViTO6Gi2tauGh7IcI/MpqrBvNFBPtUUO2gSHGO7/AE2KvkRgoh/Eft/nQ6L5VVG2Gil5E7Ic5K89p5kKOusqjHlEmV+YWw9mfwzaetC/79wIo9Yom5ndKWmZzqAJUlULHGZ5dB2xpdxbvmUocIVotvET2HWv2azSkq4Hr49MtxND+7LQFXdJQT4JWryAec17FGn6u2Vkb1wNFHew8i3mmPkTqStV80ZB6quoq6/JoGFoc/HzQ2VF+437SWj03qO0r+dyeb7BYTRczxO0UiXon+s0RzXF6ufvVn4nHXebuKfG3AqJUs23IowkKT7bNFDt8shW0014FP69lP3Gq+VM24mMkxZxPWUKsB1/hYNPT5zSagHUAOfB4u3mtIhpGKkj19s3kMTeiXnsYrnRrweTqrNyaiVrPR7uovZeKun7yyhOJLGJuTrjTOhZfxJde0or2R2vM0z8Asy+GC1hu1G5vf2ZWQLih5BnVnrlz47urFuX4316RGEoKGmz/wCatqXsqigxi0l5mr9Yka0fKHtOcFGzMnJjieTyR15124jV6mTePZmrqEt+Oeo4tXubJlZxVNsTEwqZvXyTLFe+OixloU0IfUizHhqzOJeckIwHBPAuOzajNbX2Cn0xeaS8jxfyY7MucQIekOmv5KowG1Mr2NrS6xeMkagHZ93s1Omq6U/xw8watV9mCmw/Ur1NCcsw9pq4Gf3xUe7c2ZrxqFvdM+PvKY6Je3d1DiAPjQZdYctaGfvOp5uv1ln4bFOwpR0q01Q0ZlfRaoWgDexsULaudZ2uhlAJPn3BLx+F1Zer0M9G5UlL0ES/1Rmc7zYlykq6NilrMlTlSuJKUliKD0DOiBxZB03aHFz4h8T2GlGE56slxDMaGm7Jm4yfLn3GqISOSu5T2qP1ihBl51ovugpTVMoPPZDTg4z86VYnqxku5pjh6kKXOalR67aydzLET7Gz30VWoyg7JnC0+NI9Zei86rSVWGoCwdnqa0qrOtJR6mieh3NJlXWFTMzbQQCHDfbMtzYkCnT7dtqXY8rStOxuKnsG8NGFQdrgeLPUnIKIQvKxTNz0fW4uW2dyk4cPL0j47o28Sd+rFyQfFJ2im8w3YDNswLM/rmMarRQj9RuvY/0aVX9RcwwX+/T9WFJ6lihIrSwNJiJ0sXSyqtdRl/HCUpP5utIFini64E3U+PjhRv1tmo7EaHO/h2eaXRD4BTvvU9NH4F3LkUOMo55ptJzmArxbb87KVYYjSphUKacdYrga4O63DUVlwC7l2elYyxv9zIXhVkRj4R7CxCQkqXBP8zUZblZrXPKjqus9v1G7EnRVhlNBoTIBNEPcM3ftf3TGV5fXcUJB2lSpAuB7c1WnkTUPILNfqsWKn2nKAMcj0dgnAdFmA2GpEdh2Qf1EI6eanLJ24EzUcSbZlJylM9GebpAY0aSzXB9CrJBn59HTvqAYejH9FTTFMWoLoB0ILPIHeZs3W0kHHFYXby55UfLod8t6VAlc+mkxtO2SQB2twgKxZIp41M5zoqik2wfome4swhORnnBVXmOuZ9cXksb3VTz0QD2dW7xPiXl+mVk9CQhLBJRmXM/JlzgX2/Kpo3QKN2g3eTs1aPVpSvqG1O9oDST5cjTUvo1E3Zz/AJy4HG28+CxpVluXOF64k9NDLv40r720ErQ88yBpQMWPVm0dTGemgif5Gr3onuTc5fKB8SPUt9+PSrL3fBhCq2ERshYF1VQ06ebz39Kein3a8nQz9V+hlQM3e77ZA/x0Glo26tCn5mXGonG2BpU3laVROQ3Z+lO/F0f1GNC3MOy1Ky4ZvSXupv8ACGkQNKEN3RF00uG7c0KVF0ADWMamKecb3vjVKK9krHEu57DY0lPSw7EOdd5F1Vg7RYy/PP8AGNQzfyLRJB57AuIUkVH420KPbQ0U1U0y+0iEnhdWbDMXXs8aw7zYq4LL9mc/i/M7SOcZc/PT6oScakLVV8gJLqfnYdp7cRz6bT0KKfx9qx+inZFD5BW8eUlw3vRKtRWGU6flwHbunaE5pNfzF06xiKd2lg95nPl8rHEn6jXm5MtW79BmiYceY5PmcUEcv4VuenptzPa8lvsEhxGZo8/6O5Nnu9MOc8+/atAlJjTpSzZ42d6pIjkgAGovGChSryFsxCQPqPz3060gLfmnPQn6cmohIaFjRo0dxfKsnvvdHmenzaJtcflSfQzFOBHXUiM8Hx5cVOuHMPv9exzBpp2aT1873mLhVaH6ViroVKQmV4Nujnq9rPS8yUTbXsi0bipPIwCV7fDW807OPe1MdizUkR0/tVVONMq7mkczJiWewQ66zpGZE0cTwepnqLMcI0loUXxHPWM6Kd2uzrLtQ71mn5LjeSX8Ca8Rm95qem6xmoHM6ZA0YB8cz9AGBwLPk1oU3uuRV9VMa1JiaRBXupCmNU8u0my7RNTK0eiJHWgKds9qoWessm8IWclwUFE/cnctwp/thUaJ9UXQ8U8W/tcYSmlSUyIhfgGK+9n/AM1R8fqJ6becpUdy2qvWoarIJe6nhKdSPH2fSaSaTMZtgM9m1/pvlfCoPN/xKm1rc1o406HTT1pJpv5yyn/cnI2Pmf5m/wCvmNmbmB8kfGHyVAz+mz93AtbfN5iv8+fI3/XOD/1wtaj5HyHzj/2VndRNql/3E/68lUqZz5NzbfeX/wCsH/Z2vG+N/iuG0CTgfSs6tV+RdDjX8JJxXsLP/HElea42OSt8kaD/AGivw9oM5Q2dLU087QQqcYLcxaNGanTUpydIfVOZxojrHlR31KQNkK4HIFGifScWMseW2tHB3P8AR3l1X+atNmY5DX51nHSc9WLfTqNoZbPVNFllNRFyrqpp6vimy5lc91DGkcc6NnpdKF+qEtYeY0MpOsyCBmWSmqzRTNALrLS2YlHm90rZ8o8TYMZdoXfOnleDmO6R160nLyEytRc5h8Zy1Sco6WPM7HIYzDDfEsqIgVL3Y7D4Sms6YrOhK5Im5x6yy57QyVpxCGZRW3CSgRYrzwHmvPeLMtR1+zLOh5/YXRgDeSU6zwknZehCRdZTQNV6lfxrtRKyzzYqU4VvWp7VAyA6sv1ZqmeMlK60aEixM01NoX6/teZWnUBVUVjQBq9H2juepSqaVjM/s8tcqUksy/KLl8iAi7K4grtoaFWl9cRKV7y2Ff08uQCgc3vVJF2rJ1Tmv2MoPePlcd6Wi44Poh/6Ly+OZI9nmVkZ6Uv7dP77pXvUZ7AUzdKkcqmV1jVMQpQKpl+QoUspN4Onq4v2E+Zo8Gpv0qnsLRttq6JuhOeGq13NzXRgE6sT+ubbM+klHPs8dmVfqtt+m5VUerZ4y7LPedlF/aEM1iXBw5K8pO2QSa2MLeTQjcdRYLLE14GV/fRBSXoMOr9A89Cv4kukykXsDf58mhRpJVGDpzOc7qHzq3EZ5qzY/T9v0uNLRU0TlRFZ3+9/US9jSoeo0kuE+aUQ7PIlRufZCfO0w3LUZPte8F6m1Lj+0o79lfo02avEer0g6Vrm1rCoJR8vWsT/AEPnvoMpOmH70dR1MvFPRh8VQCGVvIvKvKZ3IrCWpKrUDB76TZRtLQvj+FNpvec1M+b5AybkBdFuZ77do/H+laApd4RzGPoQUmZ/tar9YFZeASLPcN7oRpurCZ19B2zLf4lTaFeLDCieSCaiEL8J4acrTFy5EfCSzlg6vTdBMgJ30OXMjl/red+d6BiV6xVmqBfC7PTY4piUjTzarhvW0pkNViqPgU2ZqUFNNRuQr0fiT43wXxR8rOipZLRWN5v8dB9s/HOSsE3cCf7R+Wfkbva4daBsfkWRp7+0lq0PhXcydDUzpPMm7j6UOrGafhff85Xq6mWi0W74gyozHdBShA5o1h+Tp9kokmwas16exfLoH2Wlg6MU27JsV5U1XP8ADpfEe3zR5aF1ns067jdhYW19Il+G/Qca9t8aKrRSTFWWNUNgZczZo9qzh3PbN/QSLMbrjRjqVgxz01fqGbvdnyvHNEE2OUZAGfF7eKRxCJ5eFdZ1fdn3QpdTG5tZFGFFscKUe3l/JUsiTrSEGevAi0gxxqUn6Qpamk1L6KaNILWV2WgQUY2l6F7X09ilozN7iZ8S/GOvqv3tF1qE/lLcaPa1Lr3QLMtTxPpzYERnN0X1EyTpvtgOFzc46eqDx3LZF+lr4/qFF+fxlpafpzsIq8F/UTEBea0l8mBQcUjka0P9XA9it1oM+/N7ieOvkz9vKkZI8uzEVp1O6Ghf8yrxZEl/L3LMxiUdOBFeUz/KiLHiMlvJ6ACk1eg3udt5Tkd6J+v8hctF/tcRJ/oiKX1YHKLP96sZhWSpFjuoSpnXtGQ4tL9TW2KgX1p8rzVP8OUFlh0IdktgBuZdFieFhigb+V6xyWku0rMBJ0WjJQSDzAyOo60gHo7ebqI+R1mnEMDPjogLyhRzrMvhGd0+4WMrOje84EjjCnsIOg0jVaNrWab+c5PU84b64CxnQ6Fv1tlxeFmwjsUwcpK0Jlyc5NRSSLOpbX9K0l10Oko1ljWc720v3n7ea0xqvDM2UGU2TUfNpqGiamUXtSzedPKaddWViZKpZq3QfmKOqC/ISlTjmp6SXHYhti1ybNdoCSviOb07rAVrIIqeJHhzQt/ev08/12FyHEhbMWa4oWjS3p3NJPDJ/l+AdO12WjTnJZBjpqIV0OUF0NUKWLKZirxUM4uD72tnrG5jaDXJL2LCTyEaaVyzJ0CZeOeofTnbSBM/7lqajLWije1U6Uhqho9aj0cka5Xp8t++er52fF32uPBgJFPNhUaHvxn+eZQ+h+nQ9C3ZV5QtLSRPwSHL5cOSOF3N+7O6aLj+O88S/EdcNwjKor6GXXpA8ib6KJRDP6dQoGdObR1mS8mYX8z9uWipRu6WKF3I5+odtP2tOxHvEutTC3MJNUSWLc42NOxySihMImidXYUfbiee5dollsF5cpzPIdAxWf2j8p5+Cx51FJYD6OT5r3JbKx63calJwATpZfgfRrAmYOkBezqll7xSCBRFdl1tTMurjCHnT+tSm7J0GuWE/eZqpFXV3a/SSKamZezfLOhZg1A2fZLiNu70i9PoTNkvwhLbWzdbqBRnQcZlap36/M/0QJ1zXz6sQkJg3KwdDrai2aqWJjUjQocTCWE6k/x+X+QHRTKbbNGUp7rtokpQaeTREv8AZZFkNG/PVanOWf6pGh8jRFzqvZ0/CaGbLogxwZu6Qifie/wsOj3qogLYqRqyEdd8qYTtJlUs0oPC6VLqmld4f8i5budxT19wJXFL8vMaHsGfz/qqf9PmKE3qBdSKiqj6kN+TGeioxJBZHfV9wKE+l4oFWD3TknmcL4EZqKtmg2GutMtXlYKPPhP7fHlKSlnSt0o6LSDXzyqtxRz65ZnemZZcII7L3Mx4hX8vo6LKtbEtPdTfZ5mlWNg0CbZZKy2j5l4kxKlmDLKd1YdMzM6UdvFn1dPQyMTQV25nPNuQoWVR0c+DxzeuWnpXoye+/wANhWmzJh71RwLtSwzyzyyzlxJgyzTKlbmYXqI1+SeukWaikEEgjJq9FfOX2G2UreT/AGk+pl87JjwG9PYbvRFIUPP6f1T9PhibLem5iWOnxMjsS6LfraR1gTQ300SaO8m0n9JOhscwOnclAUJ2e6qATEwB4l7lUZojXCZZvgc76qnCmztDVY8HoJcaVT2kmgxIpdTXuXYfNEcBjzlk/Og8nOfqvz2avbMPpg6EnQm1KR3WaA5y6eLlS2zRpXiyCubFV7oxiBVn0ImnsUlpsxSZYrtQeDv6NyVHz1C+x+lz013QVMySzKSi0qKEQv2JAxmp8p50QedWxEBnXUJzbGr+OtB+o/P6hBi0En2vWpjfu7p1pfci3/ZQ1fnOXXLT0AYTCZjSGGqOqNKKCZKNV8Kdi8VXOPiYLl/UOUitF0ubSzw2k7fJRFh1XjIka8pt+cq5Rm9aD2r1na1wXtWSGzT74nItRrJ6KWESHPx7PZQ1p2Wkd5/9OEb1dY2KbExcMTyE4ddyjq88hcYXz2hkoqNCgyXeDX1s5JnDLkrcO+tRd7Uiu1F7YaCTekS53fgzfk7aUcVaez4/EXYmVAY/OFl5ukQSYMzY4iTPf6eVxdaqk/sNXjosqYyBMPMrOD4WpTjT0FhuNV3inPfVc98nvUakXPd0c8hu9zHJjR8I9vqwXH4ecv8AdYyZ69BHgdW0epIEhyjTUVsfqHe1rF5lyG+LSTX86TUr2TZ2VT1nQXYKKpH3vJnZjrdxF6W5tvufYH9YtKMuZ6Zpa+soigrb8WmRNwVrufzbSiMX0v78zSPO1J2g5llmP0ZNUY1JvmQOBCsS90/P5N+A5fUXpIJgL1SWZuba5o3bVwDkbiLxDl99eEQvXeTDuBp7xzPNuK1DxrSs0SdOU96B2IdGuJfxT3wxpC+czuWdx9qznVxTWqUpbrQ5JluPSgFRn0caRdXQ1pCkVQNy7/W1Krf3V6PXI6ootUONlCBIF/8AZUn10GYg1Iux/wD27hVVNNQ7e7kUK2d0i9Kezc6nPaSg5nYr4Hql5xb3QgyWezMVdtjSVKrAfkRXQMpI96Wg3WURi6NHQcTHLXlPqx9OHVN53Jc94PzSeRlUb3ed9eDIZz3bGhosw7PgENnYbfEiWDQTvuVJkr/FCcQ46eLKNOaj1wtUFTto5+nNXo9T+mPa/vH6EGrKgwNBO0mg03hqeE1pqbTJQ1evFKPN6S01YAJM3spP1X30CcbyQsELkKKlRipq0l3YrN/Fp1EyD65PP/5vyEHLnBlaD5slCW1mOyfMU4LEM5IrCmWzekzxCpkYPT+v5P5LErVV85Rvw86tBAwKdwH1R5JrZvP/ANkwOMjIKZXQ+eD9nSyREdTG0EpNftALo5YBVH3tfpPkp3KocxbtDpWsK17gbR6yWmHxWYzwmZMubHlkoUwxJEW31KmsTM2Hxp7JgLR5Zodm7hEGkOK1UMOrS1yq2scYm3oWdZQ9OUHmGSEQRrPJ59QNPnSvIgS8eg2odJ33Sv1fNdSG9nJafWh0Ps4wngZ8tWFdQg0JymXh/wCtYSQ58nPTDicMvRdNFv1s9xM8W/pqtTHKv+d0Ff8APMlLyMhPMvRlFrg+qZKmZNmQpRUlqHPUJFnxUfUisamieChdtztCx17Uy02jsZ+ftaqWLRODviTDclOm8ygpLCFJd4rA+FdLC7nJuGdtUjn4JwoB+fJY0sLKZV4HM/pVerH5nBve6CmWqyrPfPO5R60M2gZ2i9BvUGeIg/s1kl7upUfsdrzYU7W42zDYHnZpxNaKWzFoeyQJSpOq0PiyBX+n1FcsZCqScnUVn1+rPIaP+jCt4iRzkAKquckJuyFkbmZWK5YtjRSfB4PKErhzCa3Q6VBdXGwdBK+1OqeT0KTuhzrcoDH9URlb9KjSn3A9WnrgNbfogvUzTKdpA2H0yYi4lu3/AOfhrNSqPdjmCVv60CFanWFxt6AXcRCM7jpebE3pnYZV1FldAL1vlozrXAa3cn4/gY492l6ylxy396i2nToyc1zM+QXqmLu+JTKfVdZFTSozyJFfg13T4vg9NuktJzsaMZcqT8HKdTqDOhlVo06o05LCh10itQq/ncc8UCglEvejafkjZF2PwDyWu9Ypu6CImsP1dT+QD6HJS2uELTXtHy8ke6utsZXqiioZEk4qWKXP3Q+nYRYpSmNhVv8AHaZCSaKpFpKrLkbUpx824syjm7Mnlf3g2sCt3E9B/LjL59iLoHh8KvS89xpKhwYbtJt97xz1zMUJ6Gs2Kdfmlzpn1KKa3k9zxOCn1EncyKHIZWjlAjeHqWRz5YG0LNHOZtBTssVk9I/nKdjPOOOaFCIzdE2ZVKzLDHiNcKDFoVqRPDZDvUvUPL1NWgEPsTYVEQ/wj3ekc/rp9wM5IpAKaTNc/ltyMzYpJTtLp52erPZie97Sb8d80cnP6AzAp3vXOUl613wPEzml/muKVxmScrvF12UYPtKNVyn9Ob/JBFYYgXQ0lmH5ufn0wRZMOpOKM0miDiY/YJSf0acZdzqJOsWJlmZx7JbzOjzXadNHxsEjvPrbF3L5ugFHjqpzXm0akrqzXXY9+vwbqitpZLjoPU44eOxbVYjRH+jrySzbcVVL85+l6kGn/pWBS39KPS9fozUG8g6s4mD6fq86HqMdMUanF5lLTu2GfZpOUc1Bo3Y9u78b55NETNDSalTc1kZApcia3YhZeepNlXrAMvSo09UZdBDmxCrLQPxGs7oFw4My/aK/q8Zm5f0SUbLTx1uNBmfM9TS/QoAaqgcVBK4Wudv0ZjmmiarO0UOm12nwqv02+2HVB/LM0kkHGShrKvBnmz7Lf6CCtDM197SBnqa7/rICorjXoVLoK9XyrY7XToUNnX8sYzH6Duzq5mTyBqFvGD+03VVMW1EyUJKhq57Ixy8R3jZVRRAgtKxpMd/fds11pplbr96M05LnqALLrBZ5P/SWR2X0l+R8Ue7k2A7zdSu/03bDJlCdqBpz8+4wgHBEYNfS0LGeI5pSq0ZyL9bnnolOylaOwgzUYYJOeayHScnv7uosDzrl3r2ihk6nAYwH/Sq5TO2s7ZI17QMWdRl03dR05r6bHOSoHH5ZXQaTvmAtWrXOegv6Z0cCeroZHYDpMz0yaKdGVoL0O9b7Grpc/l6AsfNIA0ZoiULIfM1E/AggzOY71sU6lWUxr+0WyKwJXePhQi2Fx1Ui0iJxumaXoTn/ADJA1aZOJZJ0AS3o2XzPcaP1Vlsd9+PvafequDlP9SWF82HCu+/oLKm0YW4PudF+Lf0CKLH6ljkougBzvTfYFatGpoLihqnLuty+qB1Hccz77aEVeC3C71MwzlZ5ufb4gaFyStCqrJrSIA6aHfVKGHVQuuWVkN5BJ+LbNOnsu3/7YaFBV95IHI4pp87fueH9NOREXR5KfJA6h5Ieol0hNOKJHBLf4XYDuc6ZrlDPz8oXdLxwmq9lXGmsikzpD9m6L4UOjz6JbNFWcmXT8Ajg57vAdqSW3sXWx5eWLhhsOZ197vpcKHUDRFz7vDlyt7w2krOcsM6jQI41LSK24vLtBih1dqTpyOh+P6VaQh0lIzr2ejqa12v63RfExZrq5qYKPUTq8tDuya6uovgUn9OczKS9r/N6kTEn1KUZpDsKSsNHP101BL5A037jBJ6OfSgvaZbaRzvap+opd+8Scrpayt1ofjD3hGlCMzAAkckOKml1q/KJILDCZonUkiRrZU7VRvxR0NUJc28S+RA+RhFLNz/UDO6X/QMGyHCHM93NUulWDraIhyJdsomplS0Ggp/oZZZziS0qAcdcnT7vHcvxwFp4EWl6h/ms+9Re4b/SHJQp9g8BV7i9czOpvjPLtlclb7N1rRFmrbe9zcdZ12dZsrNajm/oLv5MiWiv1PM1UU7q9dLvljSgweSEnmRd2FNnSLpdSgKYOiTFlE6OBji2PZtjL0160gSSM6r4SHlqLkzpZhBKDHvrnFKcgrTRUeqOSHAmqJV4EVaQ1oS5tUiC4psvOKejaqMUUOPP26ziClVRan7Zls6dLU0QSMxlp7MxzyXYUmU7lu1YKQm1H1EB4Px/XzaKdbPdNgjOyrvJs1tH8zfrnq5tVCpQBRx+b/BVhvO/1pOhEWuxOsUG+P8Aj6OXw93P0Kiw16B70ViR4t/mV7LJtJU7zmWP1kaFFgmNh5XnMphpsvUZzJNqSVy1q6D9/TeT0eems3gQfyVV8/8AHZipg6Qyeim0yI8TddY46QWWGQ34bqXT2did6dw7K7GLjT8Ve6mn84izUWHTjR6Qdb07RTseaPZ6a1wSmrA7ou9TNQ5oOYqCVK/xCkSl2poYy0tR6Kejl8xM4Wy+HfTlTnJNDdmo2K4sPaoL0clWoLudMvNS25WNs5uSxRZlEpaLvOWU7JNLDzCXyE4Plk7mmWYkd483UUHLdqfFDpCa6pQfJMyygNM3LVpsrWJ9VZzWqV9Dp5jlGtrHJ8bOfIXESj5TbXTXGFO+nR4AkxxUaajGDQn5/RcXZ08cs5x5/Q8ktWlncuvO/T1Dz0XauvvIfn22c7GPN9Lqg2kM26mdx+oGq2l86dCfGNZj/Me2gM9Qp+p4KRBqNN5ztbQWvO3AId9ZsHOgzGdnhjGGFbPBDOo/jg7FrJqfC6mK9PkNOk0Oep2LlUIkAlOXQu1aEftKlH20CcjnlkEdHzl4kkdacNSXliR0cYynzn9ZLpq7NE8wlRahQChMqtTl3JvnSylG2qttsl/sX73ZLibM/iH5HnL+xQGrhmZKirMHc529t5l8rltim3M4Ri9Q85Sofi6pecTld2Pw1mYdDiFPY9EMtyH7IuKwumJKy0nqZBXIk3+9ZLQLSUgfIFhtO3l2V5L1KKH7YePeEJH+cVv6VfIfVZcV5bu1jL81KOg0HtusnTaor9IryxD64bkGmHW8/pYdid/q5X0s+HmTPaaSCQd1jtvRB7niKnOzj9RfLreHzpxDBsvuxPRP1KUtzm0XzCKwGtRJHIh6fwNeZqEHlXXIk25rsurclcZ00ylMHX6mc0JVJztsFOdppzX9fpmuT1K+o/QsaHx5Dys1Bd1DOf8AdNy1/D66+p7FeY5KV1iLDQKRpPOaj/QeeloyfA0xVl/wO/Vkde56Ojai00HKURNXqf1zomVkLbApXNOV6dCplFGO8x7yNuJ2FvL2GYXasn+plxgx4E0vfCk1bUMKL56c6mdue1Cl2cY+D2vFyr+tbWWC9dbCJbqc6FzSTe6/ViJIjzZSzHoL8Gi4gmAKp+ZVVzUX1jSw6TtLQd8xotldu3K1UW0ONReoFJojE+wOlp3VUjqUfqgLWBx5XaaMSv2l+rcqnI6OTTdmylG5nIKSfGQb0bGS0jSd6blRoSF9FU6cfow7hJqzOcfgeDny/YkyHmpMAkKcGfUvPewZ7BNMZSW/SfE6zTePV9rq8O3Z9rzye/MtIEJxdRWB9TZBMGppG4RDvQfPqzEkTvmny/f1ef2YMtdtWX6khkQ8jXfCwUNeV7RUZZziWbtL1F8439sSWtVcvIa7HzmSk5/27x1JVNDyyY7a2dzeTpZtF1qLLa6lFDmj19JYmumdh6lL5IhSjkvi+w0aalz3nvwpm2JQb3He4kcaMNkR58/RfxLCzw6i4QM6CdLcCk1a5IQ1M1tBidNb0hZQoXS+Zq0G4sjlBfevRaZzY/v3wH5aC/8AmYOPs59SYuM1/wDvyn/47XdydfliVOm0NNEh45mgy9N65nc211bTK3C3gH7HHrS+mw2o4QoqdnuqoaAi2lpuvtdRb9Loc72uedLXxkODpYyDLBxdeucO6fIu+L09JyhzikvWKlTU+7UY2NTjUEbsZ+Wo5Z+2zGuIfHb4hptGLFjKQ4qKHPHTM/vYZWC1JqNeNomd0Ysd12XXpLEmoaYNV3PAMavQm2bussTqCT0+TUA+Er2dzKPo1k7ZRtn8890vhpdNb/S6ZlAdxHhLzizaIdMIbFPM19PXgeZI1CMmwWU5Q7oY9pcZjQeTuOvJidD7Y9YR0YvNOpSfBVKQoZdu2nW0S6XMzMz2Kpl7F8P6cLLqis7HXHfWVGOQk3CJXnvsYSo/32cVnHqrWlk3J6s3RsmHemQpgTTOZxq4el57OZY5bVp1YurK1kjLUlEyOzdDoHKcniklKUOdFLYvzjryn1b7PV9hvRPzZtr/ACu17UV0fsGLJ7mv6RnuRrJeoi5KnLt//A1KfqXlhK0Qwav+agLTGG3VtUB64VXqbGm9Ny4ks3ZFRaIyNN46yjkVPiTCVlJDsWENAO9K7krpdzTgOBtbzxnyv4y8f1s1CzoKagFsgyLSqnRp9Ts5+RmocB2pBMzSXnOr0E/X8/MbvSYmmNk46g+dNkk6j8uq/wAUSKuPhakqectQLkumsvPziBYWq/MpU5aiUmJ7JpSmo9FwLxEJoydqe59CmFdy5C6EAKyi9esI7Abu2Vey393cNfIf/W8nVyJ2E9Onc6QkSCjLQlqrIWt/6h6VH4+yVWRKcnf6X6+Vt03EmLZY55P60h6ZouxFQkWpUztQAlLimmaRuKduFXoTEAVmfGGpf9zfjkBgnbcXiLEs5NrqmpcD0h+Tx5OVRt5Nas+PMMyxdA8m9jTHNgVtD+RCh0fpc06LHHalel45Uuj79FxsDVHiZb3p93ETAoodtdxbLtaKbaDc6pxUKGdvzgUbwynO+jmD2uLZHO4EViov+Zo0mOpuIsY9xH2oQQc5mOKezW5ZUn8bQvgEmjKdAZdL2nwGEpPTm9dHRYxjv0NMkx5w9PaelWFmEq+itOscWA26dWgg/wDV65PkV/6XGk/WJZua/wAksduHj7+TpmkaKwoXUyv2pHz6N9bR6dsFT+YH8+LPNY4TqyoDuhCxQi3enl807z6svQUpo3ic2c3PJpkerviPJ0AxdhUF92pjnoEahXRHYkwe6QYs/Rnx6Io1oZnec27xockaaWvjZ3R4os8JaVm+7WdEkj1d1DXWo7bO3JbD7JrUOkCp9x4Hwk5y2J3qH7+W81lJLP5vjv3aHExn31GMRe4isyhZtuW3s96pokf6agcQUUpWmaiZWLJMkNWrfZqzTSErQrx4iBaKPjNYj+QscrmuPCTdOVreTVK8dTU0XQtANwxn5obq/UTJPGHvm/EIrsOKdzhfL2rTM+OsjbMdO9h2ZfHkuPk/weQdvc2qDmHGKdapaEh8vQDlXavnivNXprETbQhJyIFjgCwwNWco9NUJ6sUXM3Ha+govRm0wOLrCn50E1o1iuwC0zIaJOYT0LklzRqmzu4ZSn0OgiwKT5pc98gLGatPr/HRznBRrPU1HeZfFTHZzQcMW4fPmgnJhqkccPdUu56pOX+6O5N583FAhbCCy6kK/qBWJnHHZpiNNCPbz4JNButHNoKJZVDq9nl1Vplu1DUp9yJmWXn29Aq0H19HOmzGUEScJKVBL+dMJ0haxeBH6vofhBCsz59Iv6NiC9DSnpd04woWjTQnnH7QqEs38fVNozMT4wa+mrwfcXYjSaeWozaa6LP710us7m0lM/S5rRlHJub0jCFmgdfYk8i5W9aXSymaXUyMRLpBjTqZ4uYlx9E7b5njZY0bV3M4tZqXLOEMqFao5fxpADsGXUWVUAmwxGdyajVvKCrrhkTZ1FjesZvRaARXrrCuWjM+3RrG6Rh1G2e6sv/xw0kADGR+gl9Dw4xDAZXtapSu287cZBlWMfA70bDN6HO10GtJPGngJDoD1OfcemUIabOgsZ+6rzUrM9+uD4ckJeouZ3sGi4r+dhtUYjPk5mpY9lX7uEUSVLc9++yv56xWcUbtoZ5SnPPK5EvtvLbNIyjHSaidkGm8hBCb8f3+wgHnE4MAJpbmZpZusvRTAZsbMdfpBTP5Dgkv48C1os3ZRc8dVnz+8bS1Baj3TiFTt5wZHoUzl5/PU2SZXY6StDRbn2YSWnKmJWeKe/wCynPoNatn66Xhs7Ft2M+ac8Edf79+TnWx9JIig9JKP85mto5kPQ8LoVP8A7j1thcJ3FKFe41WtNua5snqhAzDK+uUqOqGbxlxhdK9KkQsJsZ1LvSyODJ0aTnKehle++OdaUefYpAzwY1+vYj0M0PDzPLaubgP7IcXVw9B7WOxDYIt7bl4gajB1e0aKAu/WnZNONfocgWTetYtdq33RCyg4DiyswrnZHXvJfsVNFEV0f6AvDvKOY8g/XjK+cLzG50+01VZz2bghB2gsM9mIwMq+e1b07vZLPT3NPkXEHJFJnxHPsnkKUSCoeL14XmrmUrzECkRY7xlXZUasrxbfPSvWNVE7jp2m4UnMGNWdOtrkDC5LQJArdK3E0c48sgBpoLflMOjqgdqPFkJaTC3OP7QMczVRLVZ/F7kalgcKSxTyMUn95xKBrsu2s2OeezQnvbGr1MmfmNNSbwFugvrozZvzme+zaUPiVZhCjHCSnpYJImb48QmQ9NkIsEacV2k3q4XRnhc+PKg0D5fON6cqRG1e5OSiULULiDUe4Un5rpytVNC16yLCc9PgGOWzixXE+T3K2fzbrrfKs6R4qS6qk34xnb/5XFCaykKPWC8tnH1C/tyZtmlPcBUVh1qpIkhTvh1EE5lWYyDPlXi02i6Ios2CaHG0zFLQK577w5V8NHWVsg7Oso2eOzngZwTLU0H6fXme+9zLW4VkdwX7XMO57DHmnOHL40Kc73Nz8wrbDTTWv3LtJOJTpFiyefkbyMCQLQmiPLPTGOLl5O+1Qg1FIcitKZT7uvNkozlIviXnpb3vucpza1jZtiZXrWCcs31hG0ak/MuA6btc/wCRXM2BDjlPRpux0LcMwqIsC9/xfNUh0lZnf7ObLLsdR/TPlUnJW87U/wANuFzLo2l2Yr2L5zdOhp0bOUajns+86bvKssXPKiSk5pWkcZ/mqd2Pn/cqJXF0AXPSI1k6/Fi1eBFsUV3OS2kc6xKN8gaAz/h3e5cLMOzF+CywHzVfHir44LL0u8wCT8bXnFt6zBrlzWR9jWH8+VPm7pZ5F4nRpNPT4byRNZci93Z0Bmk4WW9nZpSsenfFT/xvLyox6DP7Svb6Z5+4/GyHdRoQu505xaK/HGXniMvn/HSyprbB2690YGG7ZuaS5wi8hkh60lfN0uxLl7o4G+2s0rorGn90ebkRP6JyQ+qlH+xaceo50QPPFQaqdYfvegAypdhr/wDD9F25O2cX1QbWgVH/AO8TNAO80s35Vbje0tHQx7uaAbMSmQUFEWGoFnv/ANfFNlKyyS4ZCyaQtogZ0wz4L8+ZpU2YynFZEp+c6rx1TOLxN2WXhp/2oIt6rZzVJjr7OPYorROTmnRkZtErXdMoj1GdvqZbw4blgx6UAk7zxS1ktAH17xhWvSVEv/5Pr5oegUpKdbPHYKcpMzbBctJUDNTh8qwI/icXKNtMZVvD2Q+2/JYZ3yI93Llwc04qJJHm5Q6z/hubifviCkRqNX0OX75AkSD0zS79dXiootE93SQdW1UgzLum5olYV4i2dP4HM5fqsnThwNd6CscM3QVaC+bq2uPF32NU8GchAxmxhudy2FFvAV4PWc9T0uc9te5F5zLwss1W58mAJ5McY5BTmU1RsiF3keeVeYI18rn68IJvMvg6ehR8n2Or0uMwEAxQJHETRxZdFPSZI72iSbeViE0L915SdB5sVtDnFl0jbtW5oYq42/n2ogKv3XDGn2aGoqBOOPqZ2noaBdrhkVdRWHjERXq8uXQXm1UpNL7wrzEJeOyvPKwHtddkwW++Yi/2+FskxPTz0l9omd75zSL2Ts8RGJddU8i3TFfk5gH8EXZg1IeuUcm6NbJ1be0XVNRDlul6LkrQ5ZatVi5UC/JCt0U3m8rJHOyNDOu9fffD/oXrjNlTxetao9tQHurVxcgakXPIWJ6ChbGebW0N6W/C8MrnrzzdZB3O9Pl+QS5bmtG59jWU+dcbs8K/dfXmNV+v1B6BQlCnMda9ujR/iF5haHAn/wBDBC/kPcWa7idu2i1JsshmTTdZ3jAMXs/OZzYBz53CIEbLwg595hvnMf6ZlFHYzvXlFyEZI8uOID26lfs6Kcmiyr4+daNWmMz06GNlJT6HGhQq/wA1+l0dFls8OrjNKhR9N73Pr5xM16xLpZdrNH75i+tvfqRlNlHVYZ6GWVehqMozfLBUyJzJftNjw6FtMLuPL1MGennP3kLoqGioaiTA54nWWSRHF7UMQGAXnWk9WRHjRPxdRSYz8Zadyb1T9ah6utbnfYaPorsh/jUPoUq3tmQ4N0CDk+mglBfbEjpWZQVLGdeuB/zwgvIfWkyhzyZ2hd5FkJYbBTUE2p1PjekztJQM2LoNStWybfslK5WO/Vvp19G+w5Hpz9j32hbyDeZZzFCPplP9yf2dnAzmh0pyKqC+du5XgJJTghGhoK7SDd33ROA0DMxtE8WrLd9Xclsy6FReBWnWIwTGlzosC45bUdRCY1YOEo2g4s/Gh070WKP9mzasEZcqRYkufTQn35U+Cxos+jz9KhLoO5nO8nNboPjlqu0aay7CfV492ULE9JBVF/WtIcZi++Rk6FZXSS8lGvMS1xk8H4hZq+6NySs4BSsn9vLMdiksrkQzjakYbXx5XHLF92yUL7bJDDI14zXRVrMl8J0p/Sfl5zCpq2p76zkOjHcIktIluKJ0xLurXI3Y3LiAGBd3pcC7OWMyhFtG64HP5e6m8FN6MPq1pkncp1qqkF4z7cryVy/6ADFKbwnWUzDPi9gLZ/kfL9EtiPWTW+rrbGLRZVj0XGp5Bq0DFv8AXVMdO82PrmK5ozhCu6I/LqG+RrSjOiZ/9paYb1cLc4XkkU1qXuvj2tFJ4GZ2mchJHE06CvSLskJBPdsNLnc8YO4pva7mdA0l9mKi+YGYXbMn4ZpJBh7MZpemx2VYBPRuI26qrjlGfb5s8NcVXUqhaqfoxc0kzxlC7u0M6+7ELjUMlWk8cH6PJDVQqKUbNRZEXqdOfS0IGJ8phxG8Aisuf0iPwqbAVkltFcyO56yXLXXrjl7FlNpuALiqc0tXPz1BugEH9faPJuq76dh9ooDWatFH+1izxxzHHDz0oUPjy9guDLsXeBY35J0A+xN+Z467M2+Kjb8U9kL8KjJoqazDvR2E2/Y3pw8Sxi6AwV4yYU7IEbXJAKWpk6heyjq0yR7znzfF01ekun303Pap009QTTnFquy1VqNngtgr77TVYB5UfS0SaQNrSoSe6bGISy1Rnn9LtYNzp5+uvfBxx2q33C5udNvNudyyc9maYBy/+ijMlJJUvegA7evuT5/U3QJTT0/apmh1k+DWp00PPkuM8l6vF/M8YfYQ+kfwfTZlOnhOJJeoMs/b5AMuKUOGLIcDJ0GxKKkyBoFWdax84k/LqCaRZ/GUkF07386PvdN18XRuOfpGbYY65+Tjp4V2xIaHWv8Aa+iXLRoe3WZbFd298c7PcfeMIIGVlSryuSZ7tdQQ0PaBzamcpJc76R9p86EWtnin/wAmILQNobHr0YnSoffUsJKZjdHenMpszaBq/Pi4v/DLc+JSU1ESdUKCq1sMZrJ50sGle/jFH7nvf8715Ll0I9Ok6LKVtKhYZrJ3v5e+rCPr4rKRVGP1VLObfpNsWsRPpPuM7n5r+wXT0nHdKnqdszIcrThpI+DzZH/k5KXCEjaT04kgJLq+moGimje2aDB3CiV55WWW35m/f0FLO+/n+I56jZOuYshL1swluMrIkLdcKcoywfhF1YWIzqssNScGMHlN4FWjStnnvuv/AJ+fo0OZDBLzK6PKqM1iZPjq5qjJRhKylPKTXnCs6RFyagPuEceG9CNaA+F+l/JL074r+HhfNcMKaV7jGgyKCpA+/wCel6Mjq92lcWXU6G17OnzxzOc56QXk1BYZPaREwhTjzIkud//Z\n", "text/plain": [ "" ] }, "metadata": { "image/jpeg": { "height": 157, "width": 157 } }, "output_type": "display_data" } ], "source": [ "for di in d:\n", " fileName = os.path.basename(di)\n", " \n", " streaksQA_path = '{}/{}_streaksqa.txt'.format(di, fileName)\n", " \n", " if os.stat(streaksQA_path).st_size > 500:\n", " try:\n", " streaksQA = np.loadtxt(streaksQA_path, delimiter=',', comments='\\\\')\n", " for streaksQA_item in streaksQA:\n", " if (abs(streaksQA_item[1]-targetRA) < detectionRadius and \\\n", " abs(streaksQA_item[2]-targetDec) < detectionRadius) or \\\n", " (abs(streaksQA_item[3]-targetRA) < detectionRadius and \\\n", " abs(streaksQA_item[4]-targetDec) < detectionRadius):\n", " folder_name = di.split('/')[-1]\n", " jpgs = glob.glob('run_merged/' + folder_name + '/' + folder_name + '_cutouts/*.jpg')\n", " display(Image(filename = jpgs[0], width=157, height=157))\n", " except:\n", " continue" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voila! Now, homework for you: can you use the example shown above to work out the detection efficiency of ZTF?" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.1" } }, "nbformat": 4, "nbformat_minor": 2 }