{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# ***WARNING***: This document is obsolete.\n", "\n", "## Please refer to the [LIGO-Virgo Public Alerts User Guide, https://emfollow.docs.ligo.org/userguide/](https://emfollow.docs.ligo.org/userguide/).\n", "\n", "---\n", "\n", "![LSC/Virgo Logo](http://www.ligo.org/images/logo02.gif)\n", "\n", "# LIGO-Virgo EM Follow-Up Tutorial\n", "\n", "by Leo P. Singer (NASA/GSFC) \n", "\n", "This document is LIGO-G1500442-v10.\n", "\n", "## Abstract\n", "\n", "This document explains how to receive, filter, and process gravitational-wave (GW) detection candidate alerts from Advanced LIGO and Virgo. We provide sample code in Python and document alternatives for users of other programming environments. You can download this document and run the code samples in [IPython Notebook](http://ipython.org/notebook.html).\n", "\n", "\n", "## Table of Contents\n", "\n", "[Introduction](#Introduction)\n", "\n", "1. [Sign up for GCN/TAN network](#1.-Sign-up-for-GCN/TAN-network)\n", "1. [Sign up for a GraceDb robot password](#2.-Sign-up-for-a-GraceDb-robot-password)\n", "1. [Install some dependencies](#3.-Install-some-dependencies)\n", "1. [Write GCN handler script](#4.-Write-GCN-handler-script)\n", "1. [Working with probability sky maps](#5.-Working-with-probability-sky-maps)\n", "1. [Basic observability calculations with Astropy](#6.-Basic-observability-calculations-with-Astropy)\n", "1. [Submitting observation coordinates to GraceDB](#7.-Submitting-observation-coordinates-to-GraceDB)\n", "\n", "[Appendix: Full example code](#Appendix:-Full-example-code)\n", "\n", "## Introduction\n", "\n", "The LIGO-Virgo data are analyzed in near real-time to search for GW transients due to compact binary coalescence (CBC) events or unmodeled \"burst\" sources. For each detection candidate, a series of alerts are produced and distributed via the [Gamma-ray Coordinates Network/Transient Astronomy Network (GCN/TAN)](http://gcn.gsfc.nasa.gov).\n", "\n", "\"GCN/TAN\n", "\n", "The GCN/TAN system may already be familiar to some, as it is has been in use since the early 1990s to transmit times and coordinates of gamma-ray bursts (GRBs) detected by gamma-ray space missions to ground-based observers. However, there are three peculiariaties of GW observations that veteran GCN users should be aware of:\n", "\n", "1. Unlike GRBs, GW sources cannot (always) be localized to a unique sky location; GW position reconstructions are (sometimes) multimodal and non-Gaussian. As a consequence, LIGO-Virgo GCN notices do *not* contain a RA, Dec, or error radius; instead they contain a pointer URL to a FITS file containing a probability sky map in the [HEALPix](http://healpix.sourceforge.net) all-sky projection.\n", "1. GCN notices will produce LIGO-Virgo candidates that span a range of significances. The significance of a candidate is reported as the false alarm rate (FAR). The FAR of a detection candidate measures approximately how frequently an event of similar strenght occurs due to chance noise fluctuations or instrumental glitches. A FAR of 1/century ($\\sim 3 \\times 10^{-10}$ Hz) is generally regarded as sufficient for a \"confident\" detection. However, alerts will produced for all events with FAR≥1/month ($\\sim 4\\times10^{-7}$ Hz), such that there should be at least one alert per month on average.\n", "1. The first few GW alerts will be distributed privately to LIGO MOU partners, over a purpose-built GCN connection.\n", "\n", "## 1. Sign up for GCN/TAN network\n", "\n", "The first step is to sign up for the GCN network. Signing up for LIGO/Virgo GCN notices is slightly different from the [standard signup process](http://gcn.gsfc.nasa.gov/invitation.html). However, if you are already receiving GCN notices (for example, for GRB follow-up), then you can reuse your existing GCN configuration and add the LIGO/Virgo notices.\n", "\n", "There are [several distribution methods for GCN notices](http://gcn.gsfc.nasa.gov/gcn_describe.html#tc7). For the purpose of this tutorial, we will focus on VOEvent over VOEvent Transport Protocol, which is among the more convenient methods for autonomous operations. However, you can use any other distribution method of your choice.\n", "\n", "To receive VOEvents, you will need a computer with a static IP address and you will need to register the IP address from which you will connect to the GCN network. Do the following two steps to submit a new GCN site configuration:\n", "\n", "1. Go to http://gcn.gsfc.nasa.gov/lvc_config_builder.html, fill out the form, and then click the \"Send LVC Notice request to GCN Admin\" button.\n", "2. If you are registering as a new GCN recipient (as opposed to adding LIGO/Virgo notices to an existing site configuration), then click the [\"New Site or New User\"](http://gcn.gsfc.nasa.gov/cfg.html?) button and fill out the resulting form. Select the `VOEVENT2.0` option for the `Distribution Method` field.\n", "\n", "## 2. Sign up for a GraceDb robot password\n", "\n", "You will need a *robot password* to download files that are linked from the GCN notices. This is an automatically generated password that you can use to download files from GraceDb in a script. It is different from your log in for the GraceDb web site, but you use the GraceDb web site to manage your robot password. Here is how to obtain a robot password.\n", "\n", "1. Got to https://gracedb.ligo.org/ and log in.\n", "1. Click the **Options** tab (see screenshot below).\n", "![GraceDb screenshot 1](gracedb-screenshot-1.png)\n", "1. Click the **Manage Password** link (see screenshot below).\n", "![GraceDb screenshot 1](gracedb-screenshot-2.png)\n", "1. Click **Get me a password!** (See screenshot below.) ***Warning: if you had an old robot password, this will reset it and generate a new one.***\n", "![GraceDb screenshot 1](gracedb-screenshot-3.png)\n", "1. Record the robot user name (which is the same as your user name for the GraceDb web site) and the generated robot password (see screenshot below).\n", "![GraceDb screenshot 1](gracedb-screenshot-4.png)\n", "1. **Optional, but recommended:** For convenience, we suggset that you add the following line to your [~/.netrc](https://www.gnu.org/software/inetutils/manual/html_node/The-_002enetrc-file.html) file:\n", "\n", " `machine gracedb.ligo.org login` *albert.einstein@LIGO.ORG* `password` *ABCDEabcde0123456789*\n", "\n", "replacing *albert.einstein@LIGO.ORG* with the robot user name and *ABCDEabcde0123456789* with the robot password. This will make it so that you can call [`curl`](http://linux.die.net/man/1/curl) or other HTTP downloader programs without explicitly providing the user name and password. To test your user name and password, you can run the following command:\n", "\n", " curl --netrc https://gracedb.ligo.org/apibasic/\n", "\n", "## 3. Install some dependencies\n", "\n", "You will need to install a few third-party Python packages to run the example code in this tutorial. These include:\n", "\n", "* [pygcn](https://pypi.python.org/pypi/pygcn) for connecting to GCN (alternatives: [comet](https://pypi.python.org/pypi/Comet))\n", "* [requests](https://pypi.python.org/pypi/requests) for easy HTTP downloads in Python (many [alternatives](https://docs.python.org/2/library/urllib.html) in Python standard library)\n", "* [healpy](https://pypi.python.org/pypi/healpy), for decoding HEALPix images (alternatives: [DS9](http://ds9.si.edu), [Aladin](http://aladin.u-strasbg.fr), [HEALPix bindings for C/C++/Fortran/Java/IDL](http://healpix.sourceforge.net), [reprojection with the Python reproject package](http://reproject.readthedocs.org/en/latest/healpix.html))\n", "* [astropy](https://pypi.python.org/pypi/astropy) version 1.0.1 or newer (optional, for computing observability windows, etc.)\n", "* [numpy](http://www.numpy.org) and [matplotlib](http://matplotlib.org), popular math and plotting packages for Python\n", "\n", "If you are on a Mac and use the [MacPorts](http://macports.org) package manager, you can install all of the above with the following command:\n", "\n", " $ sudo port install py27-gcn py27-lxml py27-requests py27-healpy py27-astropy\n", "\n", "Otherwise, the fastest way to install the dependencies is with [`pip`](https://pip.pypa.io/en/latest/quickstart.html), a package manager that comes with most Python distributions. To install these packages with `pip`, run the following command:\n", "\n", " $ pip install pygcn lxml requests healpy astropy\n", "\n", "## 4. Write GCN handler script" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we'll write a GCN handler script. First, some imports..." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Python standard library imports\n", "import tempfile\n", "import shutil\n", "import sys\n", "import glob\n", "\n", "# Third-party imports\n", "import gcn\n", "import gcn.handlers\n", "import gcn.notice_types\n", "import requests\n", "import healpy as hp\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*IPython Notebook **only***: If you are following along and running commands in the IPython Notebook, you should also run the following command so that plots are [rendered inline inside the notebook](http://ipython.org/ipython-doc/stable/interactive/magics.html?highlight=matplotlib%20inline#magic-matplotlib):" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we'll write a function that we want to get called every time a GCN notice is received. We will use the `@gcn.handlers.include_notice_types` [function decorator](https://docs.python.org/3/glossary.html#term-decorator) to specify that we only want to process certain notice types. There are three notice types:\n", "\n", "1. `LVC_PRELIMINARY`: Provides the time, significance, and basic parameters about a GW detection candidate. No localization information. Sent with a latency of a minute or so.\n", "1. `LVC_INITIAL`: A rapid sky localization is available. Sent with a latency of a few minutes.\n", "1. `LVC_UDPATE`: A refined sky localization is availaable. Sent with a latency of hours or more.\n", "\n", "In the following example, we will process only the last two types (`LVC_INITIAL` and `LVC_UPDATE`), which contain links to sky map FITS files. The following handler function will parse out the URL of the FITS file, download it, and extract the probability sky map.\n", "\n", "Events come in two very general flavors: 'CBC' or compact binary coalescence candidates detected by matched filtering, and generic 'Burst' candidates detected by model-independent methods. Most users will want to receive only 'CBC' or only 'Burst' events. In this example code, we are going to keep only 'CBC' events.\n", "\n", "***Very important:*** Note that mock or 'test' observations are denoted by the `role=\"test\"` VOEvent attribute. Alerts resulting from real LIGO/Virgo science data will always have `role=\"observation\"`. The sample code below will respond **only** to 'test' events. When preparing for actual observations, you **must remember to switch to 'observation' events**." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def get_skymap(root):\n", " \"\"\"\n", " Look up URL of sky map in VOEvent XML document,\n", " download sky map, and parse FITS file.\n", " \"\"\"\n", " # Read out URL of sky map.\n", " # This will be something like\n", " # https://gracedb.ligo.org/apibasic/events/M131141/files/bayestar.fits.gz\n", " skymap_url = root.find(\n", " \"./What/Param[@name='SKYMAP_URL_FITS_BASIC']\").attrib['value']\n", "\n", " # Send HTTP request for sky map\n", " response = requests.get(skymap_url, stream=True)\n", "\n", " # Uncomment to save VOEvent payload to file\n", " # open('example.xml', 'w').write(payload)\n", "\n", " # Raise an exception unless the download succeeded (HTTP 200 OK)\n", " response.raise_for_status()\n", "\n", " # Create a temporary file to store the downloaded FITS file\n", " with tempfile.NamedTemporaryFile() as tmpfile:\n", " # Save the FITS file to the temporary file\n", " shutil.copyfileobj(response.raw, tmpfile)\n", " tmpfile.flush()\n", "\n", " # Uncomment to save FITS payload to file\n", " # shutil.copyfileobj(reponse.raw, open('example.fits.gz', 'wb'))\n", "\n", " # Read HEALPix data from the temporary file\n", " skymap, header = hp.read_map(tmpfile.name, h=True, verbose=False)\n", " header = dict(header)\n", "\n", " # Done!\n", " return skymap, header\n", "\n", "\n", "# Function to call every time a GCN is received.\n", "# Run only for notices of type LVC_INITIAL or LVC_UPDATE.\n", "@gcn.handlers.include_notice_types(\n", " gcn.notice_types.LVC_INITIAL,\n", " gcn.notice_types.LVC_UPDATE)\n", "def process_gcn(payload, root):\n", " # Print the alert\n", " print('Got VOEvent:')\n", " print(payload)\n", "\n", " # Respond only to 'test' events.\n", " # VERY IMPORTANT! Replce with the following line of code\n", " # to respond to only real 'observation' events.\n", " # if root.attrib['role'] != 'observation': return\n", " if root.attrib['role'] != 'test': return\n", "\n", " # Respond only to 'CBC' events. Change 'CBC' to \"Burst' to respond to only\n", " # unmodeled burst events.\n", " if root.find(\"./What/Param[@name='Group']\").attrib['value'] != 'CBC': return\n", "\n", " # Read out integer notice type (note: not doing anythin with this right now)\n", " notice_type = int(root.find(\"./What/Param[@name='Packet_Type']\").attrib['value'])\n", "\n", " # Read sky map\n", " skymap, header = get_skymap(root)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we will listen for GCNs. You need to tell the `gcn.listen` function on what port you want to connect to the GCN network; this will be port 8096. You also need to tell `gcn.listen` which function to call whenever it receives an GCN; this will be the `process_gcn` function that we just defined.\n", "\n", "When you run the following code snippet, the `gcn.listen` function will continue until you interrupt the program (by pressing the stop button in IPython Notebook, typing Control-C in the terminal, or sending a kill signal to your Python script).\n", "\n", "Note: `gcn.listen` will automatically reconnect to the GCN network if the network connection is broken." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Listen for GCNs until the program is interrupted\n", "# (killed or interrupted with control-C).\n", "gcn.listen(port=8096, handler=process_gcn)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That was fun! Now kill the listener: if you are running the example code in an IPython Notebook, press the stop button. If you copied the example code into a standalone Python script, then kill the script by typing control-C.\n", "\n", "## 5. Working with probability sky maps\n", "\n", "Let's take a look at what is inside one of the LIGO/Virgo probability sky maps. They are FITS image files and can be manipulated and viewed with many commonplace FITS tools. However, they are a little unusual in two regards. First, since they are all-sky images, they are stored in the [HEALPix](http://healpix.sourceforge.net) projection, a format that is used for [Planck](http://www.cosmos.esa.int/web/planck) all-sky CMB maps and by Aladin for [archival all-sky survey images](http://aladin.u-strasbg.fr/java/nph-aladin-new.pl?frame=aladinHpxList). Second, the value stored at each pixel is the probability that the gravitational-wave source is within that pixel.\n", "\n", "![HEALPix projection](http://healpix.jpl.nasa.gov/images/healpixGridRefinement.jpg)\n", "\n", "Let's download one of them with `curl`:\n", "\n", " $ curl --netrc -O https://gracedb.ligo.org/apibasic/events/T125706/files/bayestar.fits.gz\n", "\n", "or with `wget`:\n", "\n", " $ wget --auth-no-challenge https://gracedb.ligo.org/apibasic/events/T125706/files/bayestar.fits.gz\n", "\n", "We can look at the metadata inside the FITS file by printing its header with tools like\n", "[funhead](http://linux.die.net/man/1/funhead) from [Funtools](https://github.com/ericmandel/funtools), [imhead](http://linux.die.net/man/1/imhead) from [WCSTools](http://tdc-www.harvard.edu/wcstools/), or [fitsheader](http://astropy.readthedocs.org/en/latest/io/fits/usage/scripts.html?highlight=fitsheader#module-astropy.io.fits.scripts.fitsheader) from [Astropy](http://www.astropy.org):\n", "\n", " $ funhead -a bayestar.fits.gz\n", " SIMPLE = T / conforms to FITS standard \n", " BITPIX = 8 / array data type \n", " NAXIS = 0 / number of array dimensions \n", " EXTEND = T \n", " END \n", " \t\tExtension: XTENSION\n", "\n", " XTENSION= 'BINTABLE' / binary table extension \n", " BITPIX = 8 / array data type \n", " NAXIS = 2 / number of array dimensions \n", " NAXIS1 = 4096 / length of dimension 1 \n", " NAXIS2 = 3072 / length of dimension 2 \n", " PCOUNT = 0 / number of group parameters \n", " GCOUNT = 1 / number of groups \n", " TFIELDS = 1 / number of table fields \n", " TTYPE1 = 'PROB ' \n", " TFORM1 = '1024E ' \n", " TUNIT1 = 'pix-1 ' \n", " PIXTYPE = 'HEALPIX ' / HEALPIX pixelisation \n", " ORDERING= 'NESTED ' / Pixel ordering scheme, either RING or NESTED \n", " COORDSYS= 'C ' / Ecliptic, Galactic or Celestial (equatorial) \n", " EXTNAME = 'XTENSION' / name of this binary table extension \n", " NSIDE = 512 / Resolution parameter of HEALPIX \n", " FIRSTPIX= 0 / First pixel # (0 based) \n", " LASTPIX = 3145727 / Last pixel # (0 based) \n", " INDXSCHM= 'IMPLICIT' / Indexing: IMPLICIT or EXPLICIT \n", " OBJECT = 'T125706 ' / Unique identifier for this event \n", " REFERENC= 'https://gracedb.ligo.org/events/T125706' / URL of this event \n", " INSTRUME= 'H1,L1 ' / Instruments that triggered this event \n", " DATE-OBS= '2010-08-30T10:11:34.322149' / UTC date of the observation \n", " MJD-OBS = 55438.4247028027 / modified Julian date of the observation \n", " DATE = '2014-08-13T19:21:12' / UTC date of file creation \n", " CREATOR = 'bayestar_localize_lvalert' / Program that created this file \n", " ORIGIN = 'LIGO/Virgo' / Organization responsible for this FITS file \n", " COMMENT \n", " COMMENT This simulated detection candidate is a copy of event 11178 \n", " COMMENT from the 2015 scenario of the \"First Two Years\" paper \n", " COMMENT (http://www.ligo.org/scientists/first2years/). \n", " END \n", "\n", "There are several useful pieces of information here:\n", "\n", "* `COORDSYS=C`, telling you that the HEALPix projection is in the Celestial (equatorial) frame (as all LIGO/Virgo probability sky maps will be).\n", "* `OBJECT`, the unique LIGO/Virgo identifier for the event.\n", "* `REFERENC`, a link to the candidate page in the [GraceDb](https://gracedb.ligo.org) gravitational-wave candidate event database.\n", "* `INSTRUME`, a list of gravitational-wave sites that triggered on the event (`H1` represents LIGO Hanford, `L1` is LIGO Livingston, and `V1` is Virgo).\n", "* `DATE-OBS`, the UTC time of the event. In the case of a compact binary coalescence candidate, this is the time that the signal from the merger passed through the geocenter.\n", "* `MJD-OBS`, same as `DATE-OBS`, but given as a modified Julian day.\n", "\n", "You can view the sky map in many common FITS image viewers such as [Aladin](http://aladin.u-strasbg.fr):\n", "\n", "![Aladin screenshot](aladin-screenshot.png)\n", "\n", "or [DS9](http://ds9.si.edu/doc/ref/file.html#FITSHEALPIXTable) (although DS9 shows HEALPix sky maps in an unusual orientation; see [Figure 4](http://mnras.oxfordjournals.org/content/381/2/865/F4.expansion.html) of [Calabretta & Roukema 2007](http://adsabs.harvard.edu/abs/2007MNRAS.381..865C) for information):\n", "\n", "![DS9 screenshot](ds9-screenshot.png)\n", "\n", "Now, let's go through some examples of manipulating HEALPix sky maps programmatically. The [HEALPix project](http://healpix.sourceforge.net) provides official libraries for many languages, including [C, C++, Fortran, IDL, Java, and Python](http://healpix.sourceforge.net/documentation.php). However, since this is a Python tutorial, we are going to demonstrate how to manipulate HEALPix maps with the official Python library, [Healpy](http://healpy.readthedocs.org/).\n", "\n", "First, if you have not already downloaded an example sky map, you can do so now by having Python call `curl` on the command line:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Download sky map\n", "import subprocess\n", "subprocess.check_call([\n", " 'curl', '-O', '--netrc',\n", " 'https://gracedb.ligo.org/apibasic/events/T125706/files/bayestar.fits.gz'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we need to read in the file with Healpy:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "NSIDE = 512\n", "ORDERING = NESTED in fits file\n", "Ordering converted to RING\n" ] } ], "source": [ "hpx = hp.read_map('bayestar.fits.gz')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can suppress printing informational messages while loading the file by passing the keyword argument `verbose=False`. You can read both the HEALPix image data and the FITS header by passing the `h=True` keyword argument:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "hpx, header = hp.read_map('bayestar.fits.gz', h=True, verbose=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The image data is a 1D array of values:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 3.08905283e-22, 3.04303258e-22, 1.21389563e-22, ...,\n", " 9.39631031e-22, 3.19591130e-22, 2.66302276e-22])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hpx" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Healpy has [several useful plotting routines](http://healpy.readthedocs.org/en/latest/healpy_visu.html) including `mollview` for plotting a Mollweide-projection all-sky map:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAFvCAYAAAAsUj00AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xt8XGWB//HPQzO0CUmISWmqTSWFVlqkUEqXokVbFeQi\nCgIKCmpdcNWfN1RWXeS1i6vsiuIdV3fVtWrxCoiIFKViq6AtC6VQbIsNkNoEm9LGkIQmkJTn98fz\nPD0n00mayyRz5sz3/XpNz8y5ZdIkZ77nuRprLSIiIiJS3A4p9BsQERERkbFTqBMRERFJAYU6ERER\nkRRQqBMRERFJAYU6ERERkRRQqBMRERFJAYU6ESkaxpjlxpg/xF4/b4w5Ko/nv8QY8+shtq8xxlyW\nr68XO2+XMaYx3+cVkdKiUCciE8IY02yMedYYU5e1/kEfzl5cqPcWWGtvtNaeMdQu/pHvr1tlrW3O\n93lFpLQo1InIRLHA48BbwgpjzHygnHEISiIipUahTkQm0krg7bHX7wC+D5iwwhhzuDHm+8aYXb50\n75PGGJN9ojhjzCxjzN9jr79ljGmLvf6BMeZDsfN/xxjzpDGmxRjzaWPMIX5bdvXu6caYrcaYDmPM\n1/z7jL/XfzTGbDbGtBtj7hystNEYs8oY876sdQ8ZY87zz/dXIxtjJhtjrjfGbDfG7DTGfMMYM8Vv\nW2uMOd8/X+KPO9u/fo0x5sGh/p9EJN0U6kRkIq0Dqo0xc40xk4CLcEEv7mtAFTALWIoLge8c6qTW\n2ieATmPMiX7VK4EuY8zc2Os1/vkK4DngaOBE4LXA5dnnNMZMBW4GrgLqgMeAJfhSRWPMucC/AG8E\npgJ/AH40yFv8IQNLKI8FXgz8Kse+nwVmAyf45QzgX/22NcAy/3wpruTzlbHX4XsUkRKkUCciE+0H\nuKB2OrAZaA0bYkHvX6y1z1hrtwNfAN42jPOuBZYZY6bjgtdNwFJjzCyg2lr7kDGmHjgL+LC1tsda\n+xTwZeDiHOc7G3jEWnuLtXaftfbLwM7Y9vcA/2mtfdRa+zzwn8ACY8zMHOe6NWvbJcDN1tq++E6+\nRPJdwEestR3W2m5/3vD+fo8LbwCv8NvC66X+/0BESlRZod+AiJQUiwt1f8CVxA2oesWVeGWA7bF1\nf8WVVh3MWuANQAsu/KzFhcFe/xrgSH/+v8VqdA/xXyPbi/y54nbEnh8JfMUY84WsfWZk7Ye1tssY\n8ytcad3ncCHtgNJB4AigAngg9v4M0Q34n4CXGGOmAQv89/sp3/nkH2Lfp4iUIJXUiciEstb+FVdt\neBZwS9bm3UAf0Bhb92IODFe5rMWVXi3DVUPeg6sujZdg7QCeBeqstS/wj8OttfNznO9JYH+pmy9F\ni5fC/RX4p9h5XmCtPcxau26Q9/cj4C3GmJcBU6y1v8uxz26gBzg2ds4aa201gLV2L/AAcAWwyZf0\n/RH4KNBkrW0f8n9IRFJNoU5ECuEy4NXW2p74SmvtPuCnwLXGmEpjzJHAhzmw3d0BrLVNuFK5S4G1\n1touYBdwAT7UWWv/BvwG+KIxpsoYc4gx5mhjzCtznPIO4KXGmDcaY8qADwLTY9u/CVzl28eFDhhv\nGuIt3oEr3fsU8ONBvofngW8BXzbGHOHPO8MY89rYbmuB9xEF1TXA+1HVq0jJU6gTkQlnrX3cWrsh\nvir2/APAM7jSvD8ANwLfje1nBzkOXMDZba1tjb0GiH+ttwOH4trztQM/Iwpr+89vrd0NvAnXcWE3\nrtPCPbHv4VbgOuDHxpingU3AoGPcWWufw5VMvgbXcWLA5tjzjwNNwDp/3ruAl8S2rwUqiapafw8c\nhqpeRUqesVbDQ4mIiIgUO5XUiYiIiKSAQp2IiIhICijUiYiIiKRA0sapO5todHQRERERcX6P60U/\nqKSFumXAPxf6TYhI8hhzTaHfwoSw9ppCvwURSaZDKLJQJyIlolRC2kiN9P9FIVBEAoU6ERkRhbFk\nydfPQ+FQpPgp1IlITgpvpWWwn7fCnkjxUKgTKXEKbzIUhT2R4qFQJ5JCCmoy3kbyO6YAKDIxFOpE\nipwCnCRdrt9RBT2R/FOoEykyCnGSBtm/xwp5ImOnUCeSMAptUooO9nuv0CdycAp1IgWkACcyPKrC\nFTk4hTqRcabgJjI+hvrbUuCTUqRQJ5JnCnEihac2e1KKFOpE8kBBTiTZ4n+jCniSVgp1IiOg8CZS\n/DSgsqSVQp3IQSjIiZQGleZJsVOoE8lBQU6ktCngSTFSqJOSpvAmIgej6lopFgp1UnIU5EQkH1Sa\nJ0mjUCclQUFORMaTAp4kgUKdpI4CXC4Zv+yLPSe2TkTyRbNfSKEo1ElqlG6YyxAFs3KgJ/a8P2td\ndey4Hr+tHOiM7ZcrACr4iYxFuD4p3Ml4UqiTopf+MBdCW3Z46489LyMKZeW48BYCWpXfr9Yv2/1y\nRux5I9Dm1+2K7d8V+xrhfAp4IqOlcCfjSaFOikZphDeI/ixDKVtYV561X1wIbuWx7SHsQRTouoDZ\nuDBX6x/tQD0urNX7R3PW81pc2IuX7CngiYyWetTKeFCok8RLX5jLVb1ZlmN7fFtZ1jpi27LXx89V\nHnueIQp39bH9ZvhlI1FJ3zxc2CuLHTMX2I4LkPEq3vBcRMZKJXkyFgp1kjjpCnG5StXKspbZ++UK\naiGc9RG1i+shCljlQAWwF6jz58iufs1kHdcV20ZsW8afr9F/vS6i8FYH7CEKeq1+WztqfyeSP9nX\nQYU8GQ6FOkmM0glzufbJLl3rJwpc8WrYeJVqvV9W+HUhwMWrXzOxbUHYJ95pIsheF9rn9RGVGO7J\ncUw1LuCVoVAnkn8qwZPhUKiTgktPmMsV5GDwMFfOwI4PZD2PV8PmOkfo2ABRCVszLuxV40rYav22\nEAzjbe6ySwSrYvuG1+BK4cJxDf78rbh2daEjRfx8oZ2eiOSbwp0MRaFOJlT6A1yQHZgGaxMXqjo7\niUrJ2v3zXQzssRrC0jRc+Grz5w29VWf4822JHRdK1br89k5ctSr+a8RDHkDFwCZ8ZQB1LueFB9W4\ncNfpzx/vLducdU6V2omMB42FJ7kcUug3IKUjHYEuw9AlcuER368Ml4ayj+vEhZ7O2OtOv28o6WrP\net6JKyVr8ufbFTu+FRf0wHVo6AK2EpWmhf3u8+vC+bbgQmKsw0Ovf9uV/jHdLxv8cyCqdj0ZV6oX\nQiK4ThUwsARSRMZTOq6xMhYqqZNxVRoXmcFK4YJQnRk6NHTG1ndl7RvfHlQTDQ4czkPWfvHn2e8j\n3pnhSL9uC1H1ap3fttDt2x+qb2PKcIGugyjo1QAdddAd3mMXsMSfW0QKQdOVlTZjrS30e4j7HPDP\nhX4TMjbpCXLD6ewALpxlrw9BLonVj9mdJiCqrt0FLMWVBC4D6l2A6wYW+V1q/LIXF/LKgJ3+sb8a\n+Ta/U7NfhtJGjW0nUggKeKnweeBjQ+2gkjrJm9IMc/Hn8WFEkh5cwvvrii1DSWCTX+5x67r3AI3Q\nZNy4xaGGOVTL7sYFvF6gptrnuCXAhtg54zNViMhEUweL0qA2dTJmxlyTkkA3WHu50CYOBvYYiAe3\nPqJq0WIJdLnWxQcSXg3cy/7g1w1s9MuQYUO17AJcCV54UE005Ep4LSKFlp7rteSikjoZlXRcFIbq\nwZo912ofQ4ehYhV//2Hqr9B7NvSwbQM2QH8jsBDur3b9ILqB04Ap/lGDC3tzgdszUHM8tMwDbol9\nnZAENQuFSCGp7V06KdTJsKUjyMGBYS4+XVd26Vv2urQL32sYm66ZaAgUgDXQ+wbYuBemVrh2dFP9\npgaiAs3duG0tGVydbRicGDS1mEiyKOClh0KdDCldQS679C2ulILbcIXg1UwUyqYB3wIWu+C28Xi3\nSyNwJi7ghZK73biOFJwELXs5sJ2diCSNAl5xU5s6ySl97S5KsdQtn7JL1lYDndC8zQW3XqKOrmW4\nkDeVKOhRwcBx7GaP79sVkTFL3+dA+inUyQH0RywHCm0KWxnYY7bThbp11rWxa8aV0FXi2tZNBY4D\nTgUX6N7AwJK6g83MISKFps+E4qHqVwH0RysjEUrtNuAGM94CLIaVc1wv2EZc+7pQUrcbaAE65sAj\nAPP8MbNxM18UwxAwIqVN1bLFQaGuxCnMyejFO1W0AxY2tkN/nRv65DSiKcX6/etmoOYkaGmKnacM\nhTqR4qEx75JL1a8lTIFORic+vEuYKWIN+8NdB+6xEZjSB7N7YZF1Ae9CfBu71+OqYXVfKVKs9BmS\nPLqilhD9AUr+haFP+nHj0WWg5V64fbnrAXtcBhoyVDY8RffFNW6Ik5XA1grorcfV1a7BjYvXmvMr\niEhyZX+uqPSusFRSVyIU6GT8hMHpYp0emp9wbenWAN3wbO9kKqd2uFK6c3DLhsW4uWc1xIlIWuiz\nprAU6lJOXdJl4nUC7dEAxKuhr/dQ6g7bw5Tp7a4adqp/UOePCdOIqTesSLHT507hqPo1hfTHJBMr\ntK9rJppebBfwE1hzEdwETJnCk5e+iGPqHoWTHuORz/4D3A88Ug391f6YEOjUG1YkDdRjduKppE5E\n8ih0nNgA9ESzSvwY+jqqaKOeCvbScMY2N+zJlQDzcV1jA91rioiMhkJdiqjIW5IhtLErh46fAHvd\nwMS7DV1PV9JFFTV0cPjFO90gxVTjSvhAgU4knfT5NDF0BS1y+iOR5AlVp1uAeqDJDTp8/fH0nlJL\n2xU9LJp0PxWT9/LEJ5/jqUdeDLcD3WcC9/pjs6clE5E0UJXs+FJJXRFToJPk62H/PLE3dcJOaF85\nA4AX8iRT2Q2X4ztNhM4S6g0rUgr0GZZ/CnVFSMXYUjx6/KMV2AO3Au+H1W2uDd1JPMDhp+6Eb4Mb\ns+6DROFORNJOn2f5pVBXRPTLL8VrC9DuhjjphudbDmMjJwJw2uTV0GiBjC+xA5hRmLcpIgWhz7f8\nUKgrEvpll+IV2thtgu4VQCfcAG1PT2MT8zmU56htfBI+i58r9mRcFWwGjVsnUlr0WTc2CnUJFu5c\n9EsuxSs+3lynf2yDW6H3ylp27JvJoTzLhZNuggXAxeAGJK7D9eNSXy6RUqPPvtFTqEso/TJLOvmO\nEx1AC7RfP4NuqjiUZzn5jLWuWd3yatz0YSJS6vRZODIKdQmkX2JJp/gwJZ1w523QDH+nhueYzDR2\nMe+SDW5QYqqBs1BPWBHRZ+LwqW4jIfRLK+kXnwpsjVvcBHefeg5/PKedKw+/3q37NGxZtxBW7wXW\n42ac2DTRb1ZEEkTj2w2PSupEZIKFdnZt+6cR622uZR9l1NPGC3nSl9ZVFO4tiogUIYW6AlNjUCkt\nsUAHwFfhy0AT9FBOI81U0AM1+OFNamPHqiesiOhzcygKdQWiX0opbT3Roxu4EG7kErqo4uX8kaO+\n9Gd4D3DKBbi5ZBcX8s2KSALpc/RACnUFoF9CkXb2d5zYeTOwl6cem8kalgEwkx1kruiMdZooR02A\nRSQXfaZGdJWcQPrFE4mL94ZdBcsu4L51i3l2xmTeyo3sq5vEppXzefr+s6C5BXgYF+7aC/N2RSSx\nwudrqXeiUEndBFGgEzmIFuATU+ighjL2cQa/pmZyh2tfR73fqXbw40Wk5JX6Z61C3QQo9V8ykaG1\n+8deaIYdbTPpZxIV7OVF+3vCZlCgE5HhKOXPXIW6caRGnCLD0Qp0AVvgEXh+0WHczIUAXMhNTFnZ\n7oPdElynCRGRoZXq569C3TgpxV8mkTHruAt2whPM4o+8HICZh++Am8DNBxumD9NMEyJycKX2WayO\nEnlWar9AIvnRB6wGZkN/Hx17athRNxOAF9BBZnYnfVTjesIuATYU7q2KSFEppU4UKqnLIwU6kdHq\nIapavY2+c6q5789L+QkXcQE3sbhuPZwDcHpsP5XWicjwlcJntEJdnpTCL4vI+NsK9MO6TrgV7nt0\nKQALeBBOA84k/ONplgkRGb60f1Yr1OVB2n9JRArizujpi/iby3K94ErqZqNOEyIyGmn+zFaoG4NS\n7V0jMj7CvLBNwF1wz154D9zB69hDHUcd82f4DHBchijYiYiMXFo/vxXqRimNvwwiyZKBNX38oe0V\nAJzIgzQs2Qanguss0RftJyIyCmn7LFeoG4W0/RKIJEe8tO42oIzn7zyMzRzL0TzGsWzGTQ/bAywE\nZhTmbYpIaqTpM12hboTS9MMXSaZ4sLsT1sCqX53PJPbRSDMnXLQOGucAW3Dj1vUNdiIRkWFJy2e7\nQt0IpOWHLlJUVgKrYQczqWMPx/AonAdMvcDvoLZ1IjJ2afiMV6gbhrQ2qBRJrh6/fBj618NKWPnY\n5XRQw9E8xiGfeMa3rbvA76tqWBEZu2L/vFeoO4hi/uGKpEYlZGq6eJSXsI9JHFP/KJyCD3ZBbYHe\nnIikTbF+9ivUDaFYf6gi6dCD69l6LzRvpm91NRv3ncheytlLhQt1xwGcRVSyJyKSH8WYARTqBlGM\nP0yR9OnETVG9Gj4L7Z+dwRpexTJ+xwlL17lQ19iAG+KkHzc3rIhIfhRbFlCoy6HYfogi6VcOG/vg\nHtjcdiwV9DCTHbAI6ADowg1xop6wIpJfxZQJFOqyFNMPT6Q0lPtlBu6E59ccxnMcSh27mTK33W9b\nRlQFW37AGURExqJYsoFCXUyx/NBESkurf6yEqcD7YT2L6aGCVx2+Bi4HFgCcDMxDc8KKyHgohoyg\nUOcVww9LpLT1wO4HYKqrgn2SF1FFlwt0jeDa34mIjJ+kZwWFOpL/QxKRoBy2tvH8xsNo4mi3agGu\nw8TsWX6f+QV6byJSCpKcGUo+1CX5hyMiQQbYBWwC1sMNsPPGo2jiaOa89CG4FJgLbniT9iHOIyIy\ndknNDiUd6pL6QxGRbH24oU2agXJYA9wPjz07m6nsoXZ2qw91IdBpIGIRGV9JzBAlG+ry98PI+IeI\njK8wGPEm6N4MK+HpH0/nUJ6jcdITcCEwvZ6BY9Zp3DoRGT9JC3YlGeqG/iGMJKApzIkURjXsboGN\n0EUlk3mOOYsfgtOAsjmx/TTThIiMryQFu5ILdfkJdCqdEymMTtw4dL8EamEj7GEq+5jEZJ6D2UC/\nxT1RoBORiZGUYFdyoS63eEgbakT67DDXF3uIyMRoA6qAW2AjbP/tXJ7A93w9FVhgcOPVzQCORDdg\nIlIqSirUHZikh1PilhlkPwU5kcLqh44noBmeemwmHdRAo/UDEc/BdawA/a2KyERIQmldyYS63IEu\nLleJW/Y++pAQSYbt7J8ObCOw1dDxTA2V03e7mlcRkQIodLAriVA38D85lLqVxdb15dge9imL7duP\nAp1IUjQDm2AFcCd033kElYd1wfJeH+xOBxYC9agXrIhMlEIGu5IIdQPFQ1q2zCD7aC5JkWQJf59b\noPsBuBPYCvsoo7Kmy7Wtq6kHuog6TCjYiUi6pT7UHTwxh5K3eKCLv+6P7adSOpFkiAe1WvdnuhWe\ne/ZQt7oBPxjx4hzHiIiMr0KV1g1WZFX0Bv6HlmdtjZe8ZVfFxnvBqrpVJLl6cKHuXmieBR3w9Nbp\nZBo6XaDr9o9HTgPWo1AnIhMp5BBrr5mwr5n6kroopGWXvOUSwl9f1lJEkqnNLx92ixbo66iC6bjS\nut3xfWdM6DsTEZloKQ518WFIBgt0oe1cOQcGOt3VixSVjjZoAnYaqAGmAOeAG7NOY9WJSPqlLtQZ\ncy3GXEsU1Mo5sD1cdpjrjz16UKATKSbbcTNMPAzX44Y4mYLrAXscMD0DnIybjaK2UG9SREqUMddM\nWBu71IW6watMQyldaD+X3RNOYU6keJWz/2+6GegFKnHhbje40rrFuL/x7Da2IiLpkKpQ55JwKH2L\nDxTcT1Q6V00U6DJEQx6od6tIcQp/t1ugpRPuB7bi/twb8KV1Bld0F9rVqTpWRCbWRJTWpSrURXfr\nZQzsuZpdOhded/qlwpxIcQt/yw/DI8BOoANXUrfMbyqr80/UYUJE0ik1oc6Y64bYGi+dU6ATSadm\noB52b3OBLvR8rQEagX4b21d/9yIy8YbOKmOXilBnzDeIglt8OJJQcldFVC3b5R+hU4SIFL/Q232L\ne6zDdZjYiQt1i/BVsIuBdtRhQkQmXgbo9505x0cqQp0T2syFDg+5GkP3ZC1FJB3CzVwTUOtCXTeu\ntG4KLtjNxv+zBHWYEJGJN/5t94s+1BnzVYaeGKPKL7MDnapfRNIl1mGiAWghCnWVxC4T7bFj1GFC\nRCZan88u+VfUoc6YFbjq1X6iatVQzVrrHxm/vRN3MVcvV5H08m1lt94Da3C9YMEFu0UABjgd11lC\nHSZEpBDc+Lguw+RXUYc6J7sqNR7mwl14mEpIYU4k/XYBtdC9x7Wr240rpZsCLICB14GqAw8XERl3\n45NHijbURQm3nNxt5EJdSwh0akcnUjo2RE+7/TIMRNygKlcRKZRQe1iLq4ZdkdezF22oczqJpv6p\nBabhqmPLcVWtbX6pQCdSWnYBW6Bps6uCDaV1i3DhjvOBer+vOkyIyEQJTcDC1KT5LbEr4lAXglq4\nINfiAl1nbFvXRL8pESm4WIcJ6l2g6/CrpoR92nFTh2mGCREphJBTsqcsHZuhuo0mkjEricJaCHTh\nwhz+k1oZjwQsIsWiE9dergnW1UVhLgxt0tTgVzTjrhvbJ/wdikipc5nFmJ8AfVh76ZjPWMQldRBV\nu8LAQAcKdCKlzlfBdne6Kthuog4T08EFv3BDqA4TIjJRMriMEvJL/srXijDUhVK6UN1aTRTomv1S\nbehESluoTvWl+ZUMrIKdClG1Rxi3TrNMiMhEy25KNjZFGOqm+ccM3H9CBtchQp0iRCQI0wS2Antg\n4xNueJN+XKgrwwU9qoE5RDPRiIiMtz6iAqkeBg7BNjZFFeqM+VbsVfYYdLtQY2cRifiBiLkXqHcz\niPXiAl0/rn0ddbhGdo1+X11DRGQihJvIaYRQNzDjjE5RhTqXbOuJhiII7ed2FebtiEjChWDX5Ba9\n/mVNfNnJwOuKgp2ITIR47UCoXRibogl1rtdrPdE0YC24UrpmonFf1DlCRLK14q4TLa4KthsX5ir9\ng3rc3TK4a4uuIyIynkJHiXDdKSPcTLqsM3pFEeqMuQt34a3CfeNdRKV0GjhURA5mC9DsCuzivWC7\nwV1bunDj1kG+x40SEYmEmoD4UGwDawfGEuyKItRFU2qEKcGa/fp2op5rIiKDqQZqXa/XJqKOElPB\nXVBnE1XVHniRFRHJj1Cr2MXACRLKiAqp6kZ99iIJdbOJ7p63+aVmixCR4QjtVjbBTtw9YegFW4kf\nIiq01z0SmIuqYEVkfMWHMomqX93rvaM+a5GEujKi6b92odI5ERmZcM3wN4WhCnb/mJ8ZonHq1PFK\nRMZTCHDTYutCsMswliYgRRLqenDpNX6x1Z20iAxHuIA2AQ9Db6crsYNoiBMsro1LuMbUIyIyfuI9\nXzNZr0c/ZmbiQ50xLbiq1rW4epNOorYvIiIHE9qwlONK7JpdSV03rk3dFKDMAAY4Ddfcowd1mBCR\n8dGHu4kM7erKcNebcANa77PPyCU+1LlvvhlXNdI29K4iIoMKN4Pr3SJMG1YZtlu/DEOcqMOEiIyX\ncD1qzFpfzVjGrEt0qDPmYeA+XLVJE7rAisjYtAIZaN7shrrsxoW6fnAldWW4oZPAdZpQMw8Rybcw\nkkcV0VBtfbjrT5jXHowZea1kokOdqyoJgw2DLrAiMjYZ3EW03pXU7cZVv04P2+twF9uFuOtPbY5z\niIiMVgZ3Fznbv+4nGs4khLv4EG4jk9hQZ8wD/pmf3keBTkTyYiv7x7oMVa/7Z5cowwU7DZkkIhOh\nhajXa5iYGqIZJh7Ifdggyg6+S6Fs8MseFOhEJD/CXXCTWz5SDwsqXGldP0RNPBqJ2vB2xY4VERmL\nPlz1ajtwPC6GharX0DmrmiiejazDVmJL6tyFV+NFiUi+hTHr7gMq3LAm3fgq2DB8Uq4Apza9IjJW\nIaTV4zpLxKtaQ6/YMONEOfDwiM6e4FAX5vIREcmneDjb7NrW7Z8Htgp3QY3vEwYIVUmdiIxVD1Fb\nuXr/vJ/omtOPawJShQt987JPMKQEhzpQ1auI5F8YJqkZ6HEDEfcDNeB6wIY76Xm4O+VwxywiMlZ9\nuJ71nbjgFjqDhurXqti+bYy0xjLBoW70IyqLiBxcLfvb7oZgVwZRFWwPA2eWULATkXwI7XT7ical\nA3e9qYjtV8v+Tl3DlMhQZ8y1hX4LIpJq8ZvGTtfztRffWaKaaAypHqI2LiIiYxVuFOcR1QpUcOBU\nYe2EUrqRZKJEhjoRkfHXirtoNsPOtmjaMCBqzxs6Tcxg4DQ+IiIjVc7AzhGh/W7o/TrNT1sI7g6z\njZH2LUhoqNNdsYiMtxDQ2oBa6O5zpXWVEFW1luGGN9GwJiIyVmEMulAiFx+6pM7PQY0Ldvu3tzKS\nYKfupSJSosKYdVuI5l+c46+K8XYt5URVsOFiPPLpe0Sk1M2IPQ+dIvpx15hq11mrH98Tfy/u2lTN\nSOa9T1xJnTHXFPotiEjJaMe1cdnA/inBesO2TGwZny6sB1XDisjIhTEy5+GuKSb3bv3gxtGMrjvD\nzUaJC3UiIoXR6u6Q99dfxMNbqI7VXLAiMhrhWjKDA0v6y12V6/7xMuPas1cMSdWvIlLiWv2yFtgL\n/X7asN4Q5MIFOMwJG8aRGtnFVkRK3QzcdWMeA5t4VESdtPqtf9LKaKikTkQEcONBtcaqX+N31GEQ\nO5XUichoxDtZxdvWZbLa0hkGjk03svExFepEpMRlV4W0+TYtYXaJ/gOOiAYoFhEZjnpcKd3C3Jtr\nwpO9fhlqBEbWKUvVryIidPqHn5exPwwQGqbvifd6DaV124nGmBIRGUwYn24hA0vpvBpck49+cKV0\nm1D1q4jImJTjhhAIE2mDa0cXL5Fr8EuNWyciwxHvKX880AVlFdFYdFP9ox/o3kM048ToagIU6kRE\n9msDmmKv46EtVMPW+aXa14nIcPTjSuh8aKv0q0OV6/460zAf9fDHpcumUCciAkRVrO3AKv88QzTp\ndlAGzPZ9TWlgAAAgAElEQVT7a+owERlKmGawH6YcG00DNgWYjiul68DfM95IFOh6DjzVMCjUiYjs\nVx97HqpgZzCwY0TGbzuS0V54RaRUVLO/g0QvrnSuzC97caGuDNgaP2b0tQAKdSIi+8WrXvf4ZTUw\nDRfg+vyjGteuLjR6VmmdiGTL4K4bs4H2aCw6cCVzU4nCHW1E15Hto/6KCnUiIvuV43qdtQK7cMML\nhAYvR+KuxKFtXS2uxK4KdZgQkQP1AXP9cplrSxdmjjgON+V0GdACrtdrM+7aM/qBSRTqRET2i48J\nVYW7yPbgqmW7Ytv6idrKgMasE5ED1eNuDs8HrMtqZUThrgx3Kel9AlcbEIy+WYfGqRMRGSBcUO8D\n5hOVzFUxcMy6ClwVbS0as05EBgo3evMB40rmQgeJGv9oBh4BNy4d5KONrkrqRERyKidqY1eGC3Hx\nErl4gNPwJiIShLZxjW4xG9duLlS9VuLuFXfn/ysr1ImIDNCHC3HNuA4RoUq2CteLLS4+1ImCnYiE\nQBfGpZsT9XoNgW4Kbl0ZsNv6/TeRjw5XCnUiIgcIQ5hswwW50EFinl/vpxOjL+u1iJS2PqJ5W88H\n6l1JXQ2ut+ts3Ph09wMbO3G1AU2468jYm28o1ImI5NSOa+i8hYEdKOqJqmFDezt1mBARiErsl7hF\nA1GniPC8H1cRAMC9fhm/xoyeQp2ISE4Z3B10I26YgVB6l33xjQ9xEo4TkdIT/9tvh8oKF+TCYMMh\n3LWEfbb5ZX4CHSjUiYgcxL24UrtO3EW7GhfgwkDEcSqpEylt1UA9NJzk7gen+sd0XKhr9g/WEI1L\nlz8KdSIiOYX2cuBurTuz1gWhLV151lJESksYa26hW4S2dJW4Ert+4B78MCaQNTdYXiQu1Fl7TaHf\ngoiIF2aMiCvH3Y1XE1W9BuF1NSJSivqAOW5cutDTNVTBbvXL3j3Awxx4bRnccLORBh8WERlUBjci\n/FxcVUloNxd6qlXhZprowV1Ow/Y2ol6xIpJu4SZuCZQd66papxPN7TodNybdaqBjG7Aed93I/2Dl\niSupExFJjjBmXRiEOHucujB1WPz+uIeoFE+dJkTSLUN08zbPLRpwQa6SaKDhbqAju0NE/q8PKqkT\nERlSD64ErhU3nElYNvrn7UQ9Y8uIBinONVixiKTPXLeYbuBUXKBrJOrxuhq4FdzwSGv8MeMzpWAi\nS+rUrk5EkiWEs2YGDjZcS9QuJtyt98deq9OESLrFZpI5DhfkGnAldFNxM0fsBFqeADYwmva2I8lE\niQx1IiLJE6pat/hlNa7tXFyujhKqghVJp3DTthB4gwt1jUS9XcFVu95ObAWMZ1tbVb+KiAzbLuBI\nXK9YPx4VPUSdJeKqUfWrSFqFQHcyMAtOIyqpm4q7v7sfWAf0tuHGMmlnvKpdA5XUiYgMS7gYbyfq\nuXaweV81y4RI+mRwqa0RmAcLiNrQhVDXjat2vRVcoMv/mHS5KNSJiIxIObGJG3HBLQw6mh3uQru6\n8b07F5GJNtstFlTAKbhSugaiMenWAStwIY8mv3L8rwOqfhURGbZOXEirwrWtm4O7a+8iCm/xGSaq\nYts0Zp1I8Qul768HKlxv19OAqZCZ3Unf7mqX4T4LsAf4qd9/Yv7+E1tSpx6wIpJMYayp+ti6GX6Z\nq5q1fIhtIlI8snqzLwAWATVwyOxn2Nc/yVW7rg473OuXow90I81CiQ11IiLJEy7q8fHpQulcuHCH\n8NZD1Bs2ezoxESlOPcBCmF7hAt0CYHYvFZV7eb75MFft2gJwm99/YjtLJbr61dprMOaaQr8NEREv\nBLdqoh6woa3MNFwniiAEwDCdWH/stYgUj3CjNgNYDNTDJ4C5MP2Ex5nMc+ylgu4VuI4RO+/BNc8Y\n283caGosVVInIjJioQq2xS/r/LI2x76hpxwo0IkUoz7cUEbA9DmwqBqmwyELnqGKLp7lUJ564MU+\n0G3OcezESXRJHai0TkSSKIMbs242rnqlnmgOyFBCF0r1yvwjbFOHCZHiEW7U2oHL4WLgFDjqoj9T\nw9/Z/PSx9N5eC3cCO78RO25sf+ej7VegkjoRkRHrwwW1Zv86BLZcQ5uEUrpy/1wdJkSKR/hbPgsu\nrIBGOGTZM9TwdybzHL1ba11v1zUQdZjKnmlm4ijUiYiMSbiAx8MbDF3tomAnkmwZBkwDVjYHzgPO\nhLr6PdTQwT4mwY+Bm4CW0CSjsLPIJL76VUQkmXqIqmZCgAthLYxZlx3swiVXvWFFkq8ROBkaZ8E5\n0HDJNhpp5kU8yf0s4vFrX+oD3RpgPe7vurDtZlVSJyIyau24ABdK6+KDDw9F99MiyeZnjJg9C5YD\nF8ICNtLom1w8fvNLXbVry2ZcoINCBzooklCngYhFJJlCh4lyoh6x02LbRKT4hBL4N8MyYC4ctfTP\n+6tcm2l0U4CtAFiLG9oof8aSeYrmdlG9YEUkefpwga4VV1XTj7usVjGwxC67erbPPy/8nb2IBLEZ\nI467CBbAnG89xPFs4oU8STOz+PWeM+i7qRpuv9nvmN9OEWMtxCqaUCcikkyh52srUe83iNrVDUXB\nTiR55rnFcpjNY8zDjT13+6Nvcm3obi3YGzuooqh+DVQNKyLJFErgwn1yuOPPrpbJ1XlC1bQihZXx\nj9n+cTpcAczt5STup4IedjATHgGuB+4H16BuK/n8+81HxlFJnYjImHUysH1deHT6ZZffL37JDc81\nGLFIYfUB893T8y6AM+G9l32R2TTxKMfwB17BlrUL4cIncH/ja2LHJUtRldSJiCRTfCiTbD0MPZRJ\nrmNEZGLM8I92aLzAzef67seZyQ6msYsdzGTLjQvhaoBf4nq6JndmGIU6EZG8Cxf84VTNaMw6kcJb\nDlfCIVc8w4ls5BgepYtKVj10PqwG1sX3LewAw0MpuurXUOesnrAikhyhGiZc7Gtj27Lv6kOIi19+\n1WFCZOLES9YXwoKT4HI4930/4hJ+SB17uJ4rWbX+fDhlD66H6wbGK8zls79A0YU6EZHk68R1kigj\nGr5ksBK5siG2icj4ORKmngSnwQnvW8eF3MxMdvAEjay6+XxYCa79XBeud3vyKdSJiORNKJEL49SV\nEQW67PCWq8RORMbfkcBimDuHhVvu4Qx+wyv5PZPYx428lRt+/TG4HOhYj+vhCsVSkl60beo0vImI\niIiMzDwgA8fNgWvglfyBi/gx89nEf/Nubmj9INwOdGzDldKF0vbxke8sU7ShTkQkucIQJsMpjQv7\naLw6kfFV7xZTL4KL4dSL7uKt/JCZ7KCJ2dy89lJ4zxS4oRPXOwKS2st1MEVd7q+pw0QkWeLTgYV2\ndcTWHaztnDpMiORf+Du8BE6p5s1/+h5ncwev4Pd0UcV/8x6u+sWX4D3Azp/4ffM7/Vcu41HjWNSh\nDhTsRCSJ+nCX1zCFWDCcThEKdiL5EwLdElhUDdfDh/kSJz39EJlOuGjmZ/npY2+Ha4gFui3j/q7G\nqwmZql9FRMZFqE4N1TdlsWWu6cJyHSsio5fB/f3Nh+mL4Qq4bMkNLHjmITK7YMPMefz0V++Ayw1s\n3OuPGf9AN56KvqQOVFonIkkTPkyyZ4uIt5/LDnYhyKmUTmRswt/SkcBFcGGG6T97nDs5kxPu3sa2\nVzdw6ZyV3PedpXD5elxTifsm7N2NZ0fP1JTUqTesiCRHCGahlC6EuaGmBIu3xxORsZkLLIOrM1Su\neIpv8l5OuHsbzIIP8RXu+/BS+Ca4Hq4b/DHjf0M13lklFSV1IiLJFUrs4j1hB2tXFwYqFpHRm+sX\nDbzs03fzNr7Pq/b9DibBf866glX/fT58eQ8DZ4hIRwl5qkKdqmFFJDniJW8jCWvp+HARmVi17C8J\nr7kA7oFzX/ojbv35W+FJ4DgoX7CH3qpa6F6Bq3KFiRyyZCJqFFMV6kREkif0hA00JZjI+OgBLoeV\n8LWXXs7b930fHgfOh2/Mege9c0Kgg/Gax7XQUhfqQhJWiZ2IJEsogcu+7KpkTmT0fOncgg/CaXDO\n53/GLze+Ga4EFsDij67hvq8vhRVA07eAXUz039xEtvlPXagTEUmeXD1hQYFOZLRCc4aFQDVcAwvP\nvYdbnn4zfB1YCv936XHc966l8G2An+ACXbqlpveriEiyqdpVJD9i7VOnnA4Ni7n63Ku4jXPJfBp4\nLXzm0o9y5r47faC7DWjyB6T7Riq1JXWqhhWR5Bisk4RmjxAZmQxu/LnzAXhZz91cxbWcc/7dUAfb\nv3UEjQ/sgiOA3Q8AK3E3VBP/d1aIodZSG+oC9YgVkcILvV/7OfCyq2AnMnwzgNPh1Aq4FO5+5jVM\n+XdgCfBmWMrvYdE2v+9qvyyNQAclEOpERJJhsGCnQCcyPPPcYm49U25v532H/xdTXgM8A3euW8pZ\na9fAPwJ8NXbMxA1ZkgQlEepUWiciyaL2dSLD1wiUw3EXwWfAHm3gncAjwG3QOHcL281c4Ft+/8IO\nV1LIGa5KpqOEphETkcLLVSqnGSREBtfoFg0XcdSmP/Plc98NlwP3Ax+CZXNXsX1OPNAVVqGzRkmU\n1AUqsRORwtNUYCLDU++Xy6nc+hS/55XMuLLdjUzyFfjAG69j7QvOhI6VuCGDmgY/1QQodKCDEiqp\nExFJjr5BnouIm/KrHrics2w1v7Uvp/PF03hBZTvMgaMffwTzH5YbzJug41u4MFfYQJcUJVVSByqt\nE5EkUaATiVT75XK3uCnDHesvgI8CZVDxMzj71TfzuJkFtAFrgdZCvNEDJKGUDkow1IGCnYgkgQKd\nSKQ6enpqNUf+YSuf5H0wAzY/CbX2cF64/Sl4WwY3O0Qf0FyYt5olKYEOSjTUgYKdiIhI4YUw93qY\nO4fMPZ0893PDnjLYtA/4Pzhv0Ua2mRNwM0P04apak3FTlKRAByXepi5pPwwREZHSEeZDPg2umMM/\nbfkKW+uO4YF3QfM+WPavYB6zbJsZAt0mYCsKdIMr6VAHyfyhiIiIpFe5eyz6OCz/CBfYh7GbDa8z\nV9BtdjLDHs6irRZzh4WL10DLdbhA14cC3dBKtvo1TlWxIiIi4y2UzC0DyuDLsGbJYpZ+7j5W/AbO\nnwzVj4P5cAfcA9z/ALDeH5OMMAfJDXSgkrr9kvxDEhERKW7l0dPjFsPqk7CHG5aefx83fxyWvw++\n3fte6upb4MsP+0C3GjfNV3Km+kp6VlBJXYxK7ERERPIpVjp33GK4HuxRBl4D63fA4hvguVvOxZx4\nK5htuJ6tAJ2FebtDSHqgA4W6AyjYiYiI5EM1rtp0GRy3mKWb7uRKPs9PjItsL7SvZgHX8tB3ToGN\nTwC/9Mcp0I2WQl0OCnYiIiKjFYYpOQ3ogwUnserBZZz59bU8/H6YMQnq/hXMT34LNwD3tAA/98co\n0I2FQt0gFOxERERGohboB14PDXM4d8ePuJLrOfV1G7jGwInA8WvA7LZw4R7gq/649kK94YMqpkAH\n6igxpGL7YYqIiBRGrV++Ea6YA6vh1n97KwsrN3DtHXDNVbDVnozpCIHuRr+/Al0+KdQdRDH+UEVE\nRCZOvV++Ho6bxT996StsOaaRr/47fPsZ+OQSWHbtKpadvh7OU6AbT6p+HYbww1V1rIiISFAPXA70\n0GDbeDdf5+qvf4FrDOwGGuxZXPC9O/jQCsB8FfiLP05hbrwo1I2A2tmJiIg0+uVyuBC4OMOOHYfT\nNx+ufRquWQ7/+923cNnLfgjrWoBb/P7JDXNQ/IEOVP06Ymn4oYuIiAxfxj/mAY0wfTl8Yjlr7GKs\nMdx0oaFvPuxrhUPsFZj7LZeZT8G6L+ICXTsKdBNDJXWjoBI7EREpDRm/nO0WK5Zz+MU72Tp5LtvN\n01wLvBH4UMcX+Paey+k7rhqav+iPSd7wJLmkJdCBSupGLU2/BCIiIgcKge5IaLwI5l7Ede/4AB3f\nfiH95ml6gE9+DA6xR/KND3yEvqnN0LzZH6NAVwgqqRsDdaAQEZF0CQMHzwdOh/Ng+s8f529rDbwT\nbjbAdXCrvYwPzPw2/BD43D1uJZCkeVqHkrYwFyjU5YGqY0VEpPjV40LZEqhZDLfDfy1ZzuVPf48V\ny2AhcIqtpeqZrXRfewS03AVs8McWR5iD9AY6UPVr3qT5l0RERNKsnv1jzc39CMxdzJTmdiyG99Z/\nj9tqYPm7ocmeRcMH9tBdeQRcvQoFuuRRSV0eqcRORESSrxw3nZfv/MBFcA5kVnTyXLuB82FNDXAp\n3N32Mv7E2Vz4sqvBdALXxs7TN8Hve2zSHuhAoS7v1M5ORESSKwS6+UAPNFwEN8EXFv8/3soPWT8V\npgHL3g0zv/kXWn49B84B+lfEzqEwl1QKdeNEpXYiIpIs5X45H1gCy+s44rt/ZdfGI+Ft8PBKWPwx\n+Nl153AtZ9ByxBzY/Q3cvK7NBXvXY1FKgQ4U6saVSu1ERKTw6oFLgDa4cg7zPr+Bu3gJM/6hnTYD\n3AA3/OAyVv/gNfxi/lvgc224QYO/4Y6hrYDvfXRKLcwFCnUTQKV2IiIycUKJ3DygFbgErqyGBdX8\n9JLX86YVt8M3gX6ofxDMET3wvSnwWWDrF2PnKL4wB6Ub6EChbsIo2ImIyPiLBzqg7L1UdjzFgsPu\n4he8gdplvbANmAufWvcxfseroHwK9N5GNEZdJ8UyeHC2Ug50oFA3oVQdKyIi+ZPBdVqo9a+XAfNg\nWQauhptfczZnP3MuU84AJgNvgP+35gt84zsfgWsA8wRwL1GP1uLqABFX6mEuUKgrAJXaiYjI2IRA\nNxtoB5bBccfD1XDWRbdwBr/m/H9bBeuA98Djl07nTfyMDeZUaABaVgJNhXv7eaRAF1GoKxCV2omI\nyMhU4wb6nQ90wvRLYRE0/HIby1jDpziDoz6+E1YDc+HDn/oPfs4b2T5rLrxtG25DM7QozKWVQl2B\nKdyJiMjgMsAMXGkcwFVAM1w5C06F1577C379xHmuFnUL2I/D5tqjeB2/YnvVXGgEmq8jai9X/IFO\nYW5wCnUiIiKJk/HLGX75eqicA8cBV8/i6tddxbFs5i3rfwG/BebAhmvn8R9cxc3XXgobge5V8MjD\n/vji7MkqI6NQlxDxOw+V2omIlKpq4EhgHpQdC4uAc+DIT25lESv5Bu/liLu74TZofUMtaxefzLKO\n9bAc+CCwcwXwE1yJXPF2fMim0rnhUahLIHWkEBEpFeVEpXG7gIXusbzalcotgo8u/Qzv5Lu85OnH\nyWwEZsLH51zD5/78b3A78AmAFbgQ10rUiSIdFOiGT6EuoRTsRETSLIwnFwLdQrdYcCy8B1797tuZ\nyQ5OZCMf2vU/8AxQBp9Z+lF+xzLufts5sHIbbraIW/w5Wv1Sga5UKdQlmKpkRaQwQnuu9ISDwgr/\nn2E8uYX++Rw4FVfStqCX6TOe5BO8mzP4NTOf3cGjk4+hgxqOnvYIj3/4pW54knWrgCrgq/5c7aSN\ngtzoKdQVCZXciUhhZLJeK+gdXKj+DD1Oj/TLhVA2C/r3wvIKN8TchXDZMTdwNI/xIp7kFfyeXdSz\nevJp/Dv/ylNrXwwXArv3Ar8Emv25eibyG5owCnRjo1BXRDT8iYhMrOxAF1+vcBcJ/x/lQL9fFwLd\nbL/NB7qLgdkVHP6JnRw9uYlXsYbX80sq2AvAd3knP+DtbP/KXLgCN1Dw7hX+3M0ozMlQFOqKkMKd\niOTHYOGsj6EDXfzY7GXaZQe4abgq0Nm4wLUUaIW5p0ILrlfqIphyXjtvPPx/qaeNOvZwNE38hWPY\nwUyW3bve5bUrgN1hkOD7gSZoSff/q8JcfinUFTFVyYrI6A0W2oKhgl32OQ62XxqE77HMP0L7uAZc\nqdx8mH6sWzV1lgtzNcCZvZw14w6m0caxbKaMfTzJi/gaH+RP61/tstv7w9dYD2zwz8MgwQp0MnwK\ndUVOnSlEZHRCaBtJIBvqI6M/61zFHEbC91GOK4kLVZ71uN6qJ7mXc3EdHSqB6cCZcMj0Zzi2fjPz\n2UQVXUxiH49xNA+wiO99571wA7CxD9iGG08O0t5OLk5Bbnwp1KWIqmVFZHz1k/tjI+O3hWE6+vx+\n/SS3ajZe8gYuUNUCXUTTcs3D1aEeDyx2uy0AdgLnAVNxoa4GDl+wk8bJzdTTxqE8xw5msoc6tjy6\nEO4BrgQ6+oiC3PbY11eYk/xQqBMRkTGIl2qB+1gJbc66iCahrwY6hzhHvtrnDedc2YGu3D+qiTo4\nvBmYAZXG1bDuBM7EVanWAKe45eGn7GTm5B37Ozo8yYvooIaWn8yBrbjBgZuAjidwE7TiV4T3lbSw\nK8VMoS6FVCUrIqMTSteyl0OJt72Lh6T4shYX7Opxwa4aVxJWhQt+xM6RHQKzl9lC6VpchmhQ3yq/\nbI19/Ua/PN4fX+d2afS7TsHN5lDml5W4YNcANPQypXIvk6c8x7O9h/L0PdN5unc63ISbb3UjwM24\ntnHN/v8g11hy6Q9zKp2beAp1KacqWREZXK5ANthyOOeJV89mt9cLAS+UhNX646pi+/QPsm/2Mvt9\nhXNklxqG9dW4RFbuv24tUOE2NfhdanDt4ipx1apT/euG+PpeAHo31tK7E7gTVxr3CNDdhqtnDZr9\nMn2DAx+MwlzhKNSVCPWUFZHRGUuVaK6SvlA6V4ULPLVZ23fhStlacyzDIL6tRAEvBLh2orZwjcAe\nYI7/OmFfXyI33b8MValTcaVzDX5ZiRuhJHSAKMNVvzZNgTX+W1oNNFlgE1G1aghwpdFOLhcFusJS\nqCsh2X9sCnkiEgW2Hg4stYsPWVLGwI+M+P5BfHvoLJHrXLV+/QyisBiqXOtxqSkEuSW4Uq/TcD1G\n5+GqTUNwqo8dW8b+Erigxq8OYW0KLqhNISqdq4ltB+jFTcnVgQtxzRZqDHSETg5luHZx1UAbB0p/\n1WqgEJcsCnUlTKV3IuL0xZbxatShxqoLbd/KY8vhhpkQ+OJfs8KvyxCVvjX65Ty/PN4vy/1+1bH9\nywHjNocQB64UDqJq1VD6VhnbB1xu7MW1iWvxD57A9VJ92AW8/W33QqAszdK4QIEueRTqSpza3InI\nyGW3XSvPWh/kGgKlhygUlXOgeOAj9rycwUsQY4EuXuJWQxTmpvhlOEWHf3vNQDeuerUFopK30D6u\n1S9DoFOPVYW55FKoE0A9ZkUkyFVqF6bEio/pVsXAsenigStMZp9dkpX9kRNK2kLQy8TOW0vU67Wa\nqDQufG3jglp4lBGVvvX75yG4bcQV+t3eB5UZ6H4A14DuNr982H+91tj3U9qlcNkU5IqDQp0cQNWy\nIuLEA16oXh1qBor4sCRhGe9kUY3rwFCPKxELw5FMY2CQa4wts6pWyUQhLt5kr9I/b8GVyK0DsLjO\nC+2wuxFYBd0zcMONbIi20UNUIqcwl02Brngo1ElOqpYVkQOFqcD6cGEojDcXwll2T1YYWEVahgt0\nELWXm0bULi60k4OoE0UZA9rK9ePavsVNxw0tUgl074HdnbjQVgdsARbieqi2+UcokVOAG4rCXPFR\nqJMhqVpWRJzQ27WTKBCFqtMQklpxIS20S8sefLgO2OuXYTDg7UTzq5b59Y3+6zT742cATdDf6Pdr\njp23GXZWu+O6j8Slu1oG9kpdFXuuErmhKMgVN4U6GTYFPJFSF6pjQyDq5MDq2O1EY9GFUNeEC2Zb\ncGPNbcANBLceN5bcFqIq2XpcaNtCVBLY41/Hv25YdvlHJy7Q9ZF7mBEZjIJceijUyago4ImIk92x\nAqKeovFqzl244LUJF9Y24ErU1vt9w8C9Yfy3zti6DX65haEHQS7tXqkjoSCXTgp1MmZqfyciTq5Q\nlV26Fn8en0KrPcf24X4NGS6FuXRTqJO8UbgTEUkmhbnSoFAneaeqWRGRwlOQKz0KdTKuFPBERCaO\nglxpO6TQb0BKhy42IiLjR9dYUUmdTKhcFx2V4ImIjIwCnOSiUCcFpw4WIiLDozAnQ1Gok8RQuBMR\nyU1hToZDoU4SJ/vipZAnIqVGIU5GQ6FOEk8leCJSKhTmZCwU6qRoDHaxU9gTkWKj8CbjQaFOip5K\n8kSkWCjMyXhSqJPUULgTkaRSmJOJoFAnqaOx8ESkkBTgpFAU6qQkaLoyERlPCnKSBAp1UnIU8EQk\nHxTkJGkU6qSkqUetiByMwpsUC4U6kRxUmidS2hTkpBgp1IkchAKeSGlQkJNip1AnMgKqrhUpfgpv\nklYKdSJ5oNI8kWRTkJNSoFAnkmfZHx4KeSITTyFOSpFCncg4G+rDRYFPZPQU3EQGUqgTKSDNfiEy\nPApwIgenUCeSMAf78FLokzRSaBMZO4U6kSKjNnuSBgpxIvmnUCdS5FSFK0mnACcyMRTqRFJoJB+i\nCoAyGgpqIsmjUCdS4jSgsgxF4U2keCjUiUhOCnulReFNpPgp1InIiOTrw1/hMD8UxkQkUKgTkYIY\naRgplRCokCYio5W0ULcGeL7Qb0JEkkdhR0RK3O8PtoOx1k7EGxERERGRcXRIod+AiIiIiIydQp2I\niIhICijUiYiIiKSAQp2IiIhICijUiYiIiKSAQp2IiIhICijUiYiIiKSAQp2IiIhICijUiYiIiKSA\nQp2IiIhICijUiYiIiKSAQp2IiIhICijUiYiIiKSAQp2IiIhICijUiYiIiKSAQp2IiIhICijUiYiI\niKRAWaHfgEi+GWNsod+DiEictdYU+j1I+inUSUpd55dlQGaQJf75wbYPdo7BtmcxwKTYLuHU2evi\nr8mxLv6aQY4byblzHZu9fbDzlw2xz4DtFsr2Qdk+DpnU7zZl9jGpLDz6KSvbx6RD9jGJff6t72MS\n/ZSxzz+PXkfb9+Xcz23P57GDHX/gfkN97aG/t5Ecm3X8vn1M6vf/r/ueZ1I/TOoHsw/ojz0A4uv2\nZW3LXpe9Pft8Q20fydceznvLw9fu64f+fujb55f90eY+/+jPsWSQbX2x/7bhHHsNIhND1a8iIiIi\nKZA75zQAAAcFSURBVKBQJyIiIpICCnUiIiIiKaBQJyIiIpICCnUiIiIiKaBQJyIiIpICCnUiIiIi\nKaBQJyIiIpICCnUiIiIiKaBQJyIiIpICxlpNkynporlfRSRpNPerTASFOhEREZEUUPWriIiISAoo\n1ImIiIikgEKdiIiISAoo1Mm4MMb8izHmz8aYTcaYHxpjJmdtv8QY85Ax5mFjzL3GmOOztk8yxjxo\njPllbN3njTFb/HG3GGMO9+sPNcZ8159rozFmaeyYNcaYrf5cDxpjpo7w+7jRH7/JGPMdY0xZ1vZ/\nMMb0G2POH8l5RWTiGWOa/XXiQWPMfTm2Xxm7Vmzyf9s1xphjYusfNMY8bYz5oD+m1hhzlzHmL8aY\n3xhjamLrf2eM6TLGfC3r67zTn/8hY8wqY0zdCL+PWcaY9caYbcaYHxtjMn79Mv/ewvu8evT/W1KM\nFOok74wxjcC7gIXW2vnAJODirN0eB15prT0e+DTwP1nbPwRsBuI9eX4DvNRaewLwF+Bf/Pp3Ac/7\nc50OfCF2jAXeaq090T92j/DbWWmtneu/j3Lg8tj3OQm4DrgTUM82keSzwDJ/LTj5gI3WXh+uFbjr\nyxprbYe19tHY+pOAvcDP/WGfAO6y1r4E+K1/DdALXA1cGf8axphDgeuBpf5a9jDw/hF+H9cBX7DW\nzgH+DlwW27Y2dr37zAjPK0VOoU7GQyfQB1T4kq0KoDW+g7X2T9bap/3L9UBD2GaMaQDOBr5NLCxZ\na++y1j6f45h5wO/8Pk8BHcaYRbEvd0DgMsYcYYy5yRhzn3+8PNc3Yq1dFXv5f/H3CXwAuAl4Ktex\nIpJIw70BeyvwoxzrTwMes9bu8K/fAHzPP/8ecB6AtXavtfZe4Nms4/txQazSGGOAw/HXx+Fcl/wx\nr8JdewZ8zRF+f5JCCnWSd9badlxp2V+BJ4EOa+3qIQ65DLgj9vpLwD8Dz+feHYB/jB3zEPAGX2U7\nC3cnPTO27/dyVEV8BfiSv1u/EBcgB+WrNy4FVvnXM4BzgW/4XTQ2kEjyWWC1MeZ+Y8y7BtvJGFMB\nnAHcnGPzxcAPY6/rrbVt/nkbUJ/ja0Yv3I3ph4BHcGFuHvAdv3k416U63DU1XB9bgRmx7S/31bp3\nGGOOHex7lHQqO/guIiNjjDkauAJoBJ4GfmaMucRae2OOfV+FC2hL/OtzgF3W2geNMcsGOf8ngees\nteHC+r+4C+P9wHbgj8A+v+0Sa+2TxphK4GZjzNustT/A3W3Pcze9AFQZYyqstXsH+bb+C1etca9/\n/WXgE9Za6++cdXcsknxLrLV/M8YcAdxljNlqrf1Djv1eD9xjre2Ir/RVp68HPp7r5P56MOQNnjGm\nGvgqcIK19gnf3u4q4FpGfl3K9gAw01q71xhzFnAr8JJhHispoFAn42ER8Edr7R4AY8wtwMuBAaHO\nd474FnCmtfbvfvXLcaVuZwNTgGpjzPettW/3xyzHVc2+JpzHWrsP+EjsvPfi2txhrX3SL7uNMT8E\nTgZ+gAthi621z2W9pztxd9r/Z639J7/u34A6a238zv4k4Mf+4jsVOMsY02etvW3k/10iMhGstX/z\ny6eMMT/HXQ9yhbqLyV31ehbwgG/mEbQZY6Zba3caY14I7DrI25gHPGGtfcK//hlRSBzWdcl33jjE\nl9Y14KtvrbVdse91lTHmv4wxtb72REqAql9lPGwFTjHGlPtSrNNwnR72M8a8GLgFuNRa2xTWW2uv\nstbOtNbOwl1Y744FujNx1bLnWmt7Y+cqN8Yc5p+fDvRZa7f66tipfn0Gd4e9yR/2G+CDsXMs8F//\nTN/AOAS6y4HX4trX7GetPcpaO8u/z5uA9yrQiSSXMabCGFPlnx+G+7velGO/w4FXAr/IcZq3cGDY\nuw14h3/+Dlzp2IBTZr1+HJhrop74pxNdH4d1XcK1IX5T9tc0xtT7ay7GmJNxs0Yp0JUQTRMm48IY\n8zHcxeZ5YAOuh+o7Aay1/22M+TbwRly7O3BB7OSscyyF/9++HapEEEVxGP8OrBaxrF3YhxGb0SCi\ngkFcmz6ET2A2CYugxSKYTPoOlgHtGsSy4RjuHZTFsgbBy/drM8wwzA1/zsy5h5PM3KjHT8Ai0IfU\nQ2aO67TtbX3WC7Cfmc81uO+BBcoE7h1wXFskK8AZ5at5QGmtjn94jynQAe/11NXsRFlEnAM3mXk9\n7zpJ+ht1v20/sToALjLzNCIOoORSvW4XWM/MrZn7lyjbO0bf/4hFxBC4BFYpWbHZt20jogOWKbn1\nBqzVD84dvvYNd8BeZr7OkUsjYAIMKfm6nZnTiDgCDinDGB+UvHv89aLp37GokyRJaoDtV0mSpAZY\n1EmSJDXAok6SJKkBFnWSJEkNsKiTJElqgEWdJElSAyzqJEmSGvAJYnXzEX9OKFEAAAAASUVORK5C\nYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "hp.mollview(hpx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each entry in the array represents the probability contained within a quadrilateral pixel whose position on the sky is uniquely specified by the index in the array and the array's length. Because HEALPix pixels are equal area, we can find the number of pixels per square degree just from the length of the HEALPix array:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.013113963206424481" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "npix = len(hpx)\n", "sky_area = 4 * 180**2 / np.pi\n", "sky_area / npix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The [`pix2ang`](http://healpy.readthedocs.org/en/latest/generated/healpy.pixelfunc.pix2ang.html) function converts from pixel index to spherical polar coordinates; the function [`ang2pix`](http://healpy.readthedocs.org/en/latest/generated/healpy.pixelfunc.ang2pix.html) does the reverse.\n", "\n", "Both `pix2ang` and `ang2pix` take, as their first argument, `nside`, the lateral resolution fo the HEALPix map. You can find `nside` from the length of the image array by calling [`npix2nside`](http://healpy.readthedocs.org/en/latest/generated/healpy.pixelfunc.npix2nside.html#healpy.pixelfunc.npix2nside):" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "512" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nside = hp.npix2nside(npix)\n", "nside" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's look up the right ascension and declination of pixel number 123. We'll call `pix2ang` to get the spherical polar coordinates $(\\theta, \\phi)$ in radians, and then use Numpy's [`rad2deg`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.rad2deg.html) function to convert these to right ascension and declination in degrees." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(129.375, 89.269029291573901)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ipix = 123\n", "theta, phi = hp.pix2ang(nside, ipix)\n", "ra = np.rad2deg(phi)\n", "dec = np.rad2deg(0.5 * np.pi - theta)\n", "ra, dec" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's find which pixel contains the point RA=194.95, Dec=27.98." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "833621" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ra = 194.95\n", "dec = 27.98\n", "theta = 0.5 * np.pi - np.deg2rad(dec)\n", "phi = np.deg2rad(ra)\n", "ipix = hp.ang2pix(nside, theta, phi)\n", "ipix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's find the highest probability pixel. What is the probability inside it?" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "5.7702549383975565e-05" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ipix_max = np.argmax(hpx)\n", "hpx[ipix_max]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Where is the highest probability pixel on the sky? Use `pix2ang`." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(213.22265625, -37.450292350169015)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "theta, phi = hp.pix2ang(nside, ipix_max)\n", "ra = np.rad2deg(phi)\n", "dec = np.rad2deg(0.5 * np.pi - theta)\n", "ra, dec" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How do we find the probability that the source is contained within a circle on the sky? First we find the pixels that are contained within the circle using `query_disc`. Note that `query_disc` takes as its arguments the Cartesian coordinates of the center of the circle, and its radius in radians.\n", "\n", "Then, we sum the values of the HEALPix image array contained at those pixels." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.056844148994230181" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# RA, Dec, and radius of circle in degrees\n", "ra = 213.22\n", "dec = -37.45\n", "radius = 3.1\n", "\n", "# Spherical polar coordinates and radius of circle in radians\n", "theta = 0.5 * np.pi - np.deg2rad(dec)\n", "phi = np.deg2rad(ra)\n", "radius = np.deg2rad(radius)\n", "\n", "# Cartesian coordinates of center of circle\n", "xyz = hp.ang2vec(theta, phi)\n", "\n", "# Array of indices of pixels inside circle\n", "ipix_disc = hp.query_disc(nside, xyz, radius)\n", "\n", "# Probability that source is within circle\n", "hpx[ipix_disc].sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, we can use the `query_polygon` function to look up the indices of the pixels within a polygon (defined by the Cartesian coordinates of its vertices), and then compute the probability that the source is inside that polygon by summing the values of the pixels." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.0011913816551896161" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Vertices of polygon\n", "xyz = [[-0.69601758, -0.41315628, -0.58724902],\n", " [-0.68590811, -0.40679797, -0.60336181],\n", " [-0.69106913, -0.39820114, -0.60320752],\n", " [-0.7011786 , -0.40455945, -0.58709473]]\n", "\n", "# Array of indices of pixels inside polygon\n", "ipix_poly = hp.query_polygon(nside, xyz)\n", "\n", "# Probability that source is within polygon\n", "hpx[ipix_poly].sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These are all of the HEALPix functions from Healpy that we will need for the remainder of the this tutorial.\n", "\n", "Other useful Healpy functions include [`ud_grade`](http://healpy.readthedocs.org/en/latest/generated/healpy.pixelfunc.ud_grade.html) for upsampling or downsampling a sky map, and [`get_interp_val`](http://healpy.readthedocs.org/en/latest/generated/healpy.pixelfunc.get_interp_val.html#healpy.pixelfunc.get_interp_val) for performing bilinear interpolation between pixels. See the [Healpy tutorial](http://healpy.readthedocs.org/en/latest/tutorial.html) for other useful operations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. Basic observability calculations with Astropy\n", "\n", "Now we are going to teach our GCN handler how to determine whether a gravitational-wave event is observable. We are goint to use the [astropy.coordinates](http://astropy.readthedocs.org/en/stable/coordinates/index.html) module of the [Astropy package](http://astropy.org) to do [observation planning in Python](http://astropy.readthedocs.org/en/stable/coordinates/observing-example.html). First, we will need to import a few extra Python modules:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import astropy.coordinates\n", "import astropy.time\n", "import astropy.units as u" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The LIGO/Virgo probability sky maps are always in equatorial coordinates. Once we have looked up the coordinates of the HEALPix pixels, we will use [the positional astronomy features of Astropy](http://astropy.readthedocs.org/en/stable/coordinates/observing-example.html) to transform those coordinates to an alt/az frame for a particular site on the Earth at a particular time. Then we can quickly determine which pixels are visible from that site at that time, and integrate (sum) the probability contained in those pixels.\n", "\n", "Note: users may want to do something more sophisticated like determine how much of the probability is visible for at least a certain length of time. This example will illustrate one key function of HEALPix (looking up coordinates of the grid with `hp.pix2ang`) and some of the key positional astronomy functions with Astropy." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def prob_observable(m, header):\n", " \"\"\"\n", " Determine the integrated probability contained in a gravitational-wave\n", " sky map that is observable from a particular ground-based site at a\n", " particular time.\n", "\n", " Bonus: make a plot of probability versus UTC time!\n", " \"\"\"\n", "\n", " # Determine resolution of sky map\n", " npix = len(m)\n", " nside = hp.npix2nside(npix)\n", "\n", " # Get time now\n", " time = astropy.time.Time.now()\n", " # Or at the time of the gravitational-wave event...\n", " # time = astropy.time.Time(header['MJD-OBS'], format='mjd')\n", " # Or at a particular time...\n", " # time = astropy.time.Time('2015-03-01 13:55:27')\n", "\n", " # Geodetic coordinates of observatory (example here: Mount Wilson)\n", " observatory = astropy.coordinates.EarthLocation(\n", " lat=34.2247*u.deg, lon=-118.0572*u.deg, height=1742*u.m)\n", "\n", " # Alt/az reference frame at observatory, now\n", " frame = astropy.coordinates.AltAz(obstime=time, location=observatory)\n", "\n", " # Look up (celestial) spherical polar coordinates of HEALPix grid.\n", " theta, phi = hp.pix2ang(nside, np.arange(npix))\n", " # Convert to RA, Dec.\n", " radecs = astropy.coordinates.SkyCoord(\n", " ra=phi*u.rad, dec=(0.5*np.pi - theta)*u.rad)\n", "\n", " # Transform grid to alt/az coordinates at observatory, now\n", " altaz = radecs.transform_to(frame)\n", "\n", " # Where is the sun, now?\n", " sun_altaz = astropy.coordinates.get_sun(time).transform_to(altaz)\n", "\n", " # How likely is it that the (true, unknown) location of the source\n", " # is within the area that is visible, now? Demand that sun is at\n", " # least 18 degrees below the horizon and that the airmass\n", " # (secant of zenith angle approximation) is at most 2.5.\n", " prob = m[(sun_altaz.alt <= -18*u.deg) & (altaz.secz <= 2.5)].sum()\n", "\n", " # Done!\n", " return prob" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we need to update our GCN handler to call this function." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Function to call every time a GCN is received.\n", "# Run only for notices of type LVC_INITIAL or LVC_UPDATE.\n", "@gcn.handlers.include_notice_types(\n", " gcn.notice_types.LVC_INITIAL,\n", " gcn.notice_types.LVC_UPDATE)\n", "def process_gcn(payload, root):\n", " # Print the alert\n", " print('Got VOEvent:')\n", " print(payload)\n", "\n", " # Respond only to 'test' events.\n", " # VERY IMPORTANT! Replce with the following line of code\n", " # to respond to only real 'observation' events.\n", " # if root.attrib['role'] != 'observation': return\n", " if root.attrib['role'] != 'test': return\n", "\n", " # Respond only to 'CBC' events. Change 'CBC' to \"Burst' to respond to only\n", " # unmodeled burst events.\n", " if root.find(\"./What/Param[@name='Group']\").attrib['value'] != 'CBC': return\n", "\n", " skymap, header = get_skymap(root)\n", " prob = prob_observable(skymap, header)\n", " print('Source has a %d%% chance of being observable now' % round(100 * prob))\n", " if prob > 0.5:\n", " pass # FIXME: perform some action" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's run the new GCN handler now..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Listen for GCNs until the program is interrupted\n", "# (killed or interrupted with control-C).\n", "gcn.listen(port=8096, handler=process_gcn)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 7. Submitting observation coordinates to GraceDB\n", "\n", "Suppose you have performed some EM observations to follow up on a candidate GW event, and you now want to supply the coordinates of those observations to the LV-EM group. The first step is to obtain a robotic access password from GraceDB and add it to a netrc file (see Section 2). One can then use the Python GraceDB client in order to submit a list of coordinates to GraceDB. To install the GraceDB client package:\n", "\n", " $ pip install ligo-gracedb \n", "\n", "Note: It is highly recommended to install the GraceDB client (as well as other packages described used in this tutorial) inside a [virtual environment](https://virtualenv.pypa.io/en/latest/). This isolates the packages required for interacting with GraceDB from other packages on the system.\n", "Now that the GraceDB client is installed, one can use a script such as this to submit observation records to GraceDB:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from ligo.gracedb.rest import GraceDbBasic, HTTPError\n", "\n", "service = 'https://gracedb.ligo.org/apibasic/'\n", "g = GraceDbBasic(service)\n", "\n", "graceid = 'T125706'\n", "raList = [45.0, 47.0, 49.0]\n", "raWidthList = 2.0\n", "decList = [45.0, 47.0, 49.0]\n", "decWidthList = 2.0\n", "startTimeList = ['2015-05-01T12:30:10.95', '2015-05-01T12:31:10.95', '2015-05-01T12:32:10.95']\n", "durationList = 100.0\n", "comment = 'Some text comment goes here. This is optional.'\n", "\n", "g.writeEMObservation(graceid, 'Test', raList, raWidthList,\n", " decList, decWidthList, startTimeList, durationList, comment)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the above example, the lists of ra and dec values were chosen arbitrarily, as well as the start times and durations (or exposure times). Since the comma-separated lists have three elements, there will be three footprints associated with this observation. And since only one value was provided for the raWidth and decWidth, these widths will be assumed to apply to all three footprints.\n", "\n", "It is possible to create the same record using curl instead:\n", "\n", " $ curl -X POST --netrc --data \"group=Test&raList=45,47,49&raWidthList=2&decList=45,47,49&decWidthList=2&startTimeList=2015-05-01T12:30:10.95,2015-05-01T12:31:10.95,2015-05-01T12:32:10.95&durationList=100,100,100\" https://gracedb.ligo.org/apibasic/events/T125706/emobservation/ \n", "\n", "However, the data argument in the above curl command is rather unwieldy--and this is only with three simple footprints. So using Python (as in the above code example) or your own favorite language and HTTP library is much preferred. We are also planning an email submission mechanism which may be eaiser for some users. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Appendix: Full example code" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is a complete, working GCN processing script. Copy it into a `.py` file and customize as needed." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Python standard library imports\n", "import tempfile\n", "import shutil\n", "import sys\n", "import glob\n", "\n", "\n", "# Third-party imports\n", "import gcn\n", "import gcn.handlers\n", "import gcn.notice_types\n", "import requests\n", "import healpy as hp\n", "import numpy as np\n", "import astropy.coordinates\n", "import astropy.time\n", "import astropy.units as u\n", "\n", "\n", "def get_skymap(root):\n", " \"\"\"\n", " Look up URL of sky map in VOEvent XML document,\n", " download sky map, and parse FITS file.\n", " \"\"\"\n", " # Read out URL of sky map.\n", " # This will be something like\n", " # https://gracedb.ligo.org/apibasic/events/M131141/files/bayestar.fits.gz\n", " skymap_url = root.find(\n", " \"./What/Param[@name='SKYMAP_URL_FITS_BASIC']\").attrib['value']\n", "\n", " # Send HTTP request for sky map\n", " response = requests.get(skymap_url, stream=True)\n", "\n", " # Uncomment to save VOEvent payload to file\n", " # open('example.xml', 'w').write(payload)\n", "\n", " # Raise an exception unless the download succeeded (HTTP 200 OK)\n", " response.raise_for_status()\n", "\n", " # Create a temporary file to store the downloaded FITS file\n", " with tempfile.NamedTemporaryFile() as tmpfile:\n", " # Save the FITS file to the temporary file\n", " shutil.copyfileobj(response.raw, tmpfile)\n", " tmpfile.flush()\n", "\n", " # Uncomment to save FITS payload to file\n", " # shutil.copyfileobj(reponse.raw, open('example.fits.gz', 'wb'))\n", "\n", " # Read HEALPix data from the temporary file\n", " skymap, header = hp.read_map(tmpfile.name, h=True, verbose=False)\n", " header = dict(header)\n", "\n", " # Done!\n", " return skymap, header\n", "\n", "\n", "def prob_observable(m, header):\n", " \"\"\"\n", " Determine the integrated probability contained in a gravitational-wave\n", " sky map that is observable from a particular ground-based site at a\n", " particular time.\n", "\n", " Bonus: make a plot of probability versus UTC time!\n", " \"\"\"\n", "\n", " # Determine resolution of sky map\n", " npix = len(m)\n", " nside = hp.npix2nside(npix)\n", "\n", " # Get time now\n", " time = astropy.time.Time.now()\n", " # Or at the time of the gravitational-wave event...\n", " # time = astropy.time.Time(header['MJD-OBS'], format='mjd')\n", " # Or at a particular time...\n", " # time = astropy.time.Time('2015-03-01 13:55:27')\n", "\n", " # Geodetic coordinates of observatory (example here: Mount Wilson)\n", " observatory = astropy.coordinates.EarthLocation(\n", " lat=34.2247*u.deg, lon=-118.0572*u.deg, height=1742*u.m)\n", "\n", " # Alt/az reference frame at observatory, now\n", " frame = astropy.coordinates.AltAz(obstime=time, location=observatory)\n", "\n", " # Look up (celestial) spherical polar coordinates of HEALPix grid.\n", " theta, phi = hp.pix2ang(nside, np.arange(npix))\n", " # Convert to RA, Dec.\n", " radecs = astropy.coordinates.SkyCoord(\n", " ra=phi*u.rad, dec=(0.5*np.pi - theta)*u.rad)\n", "\n", " # Transform grid to alt/az coordinates at observatory, now\n", " altaz = radecs.transform_to(frame)\n", "\n", " # Where is the sun, now?\n", " sun_altaz = astropy.coordinates.get_sun(time).transform_to(altaz)\n", "\n", " # How likely is it that the (true, unknown) location of the source\n", " # is within the area that is visible, now? Demand that sun is at\n", " # least 18 degrees below the horizon and that the airmass\n", " # (secant of zenith angle approximation) is at most 2.5.\n", " prob = m[(sun_altaz.alt <= -18*u.deg) & (altaz.secz <= 2.5)].sum()\n", "\n", " # Done!\n", " return prob\n", "\n", "\n", "# Function to call every time a GCN is received.\n", "# Run only for notices of type LVC_INITIAL or LVC_UPDATE.\n", "@gcn.handlers.include_notice_types(\n", " gcn.notice_types.LVC_INITIAL,\n", " gcn.notice_types.LVC_UPDATE)\n", "def process_gcn(payload, root):\n", " # Print the alert\n", " print('Got VOEvent:')\n", " print(payload)\n", "\n", " # Respond only to 'test' events.\n", " # VERY IMPORTANT! Replce with the following line of code\n", " # to respond to only real 'observation' events.\n", " # if root.attrib['role'] != 'observation': return\n", " if root.attrib['role'] != 'test': return\n", "\n", " # Respond only to 'CBC' events. Change 'CBC' to \"Burst' to respond to only\n", " # unmodeled burst events.\n", " if root.find(\"./What/Param[@name='Group']\").attrib['value'] != 'CBC': return\n", "\n", " skymap, header = get_skymap(root)\n", " prob = prob_observable(skymap, header)\n", " print('Source has a %d%% chance of being observable now' % round(100 * prob))\n", " if prob > 0.5:\n", " pass # FIXME: perform some action\n", "\n", "\n", "# Listen for GCNs until the program is interrupted\n", "# (killed or interrupted with control-C).\n", "gcn.listen(port=8096, handler=process_gcn)" ] } ], "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.3" } }, "nbformat": 4, "nbformat_minor": 2 }