# Basic Observability Calculations¶

Now we are going to teach our GCN handler how to determine whether a gravitational-wave event is observable. We are going to use the astropy.coordinates module. (See also the Astropy example on observation planning in Python.) First, we will need to import a few extra Python modules:

import astropy.coordinates
import astropy.time
import astropy.units as u


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 Astropy to transform those coordinates to a horizontal (altitude–azimuth) 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.

Note

You 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. For more advanced functionality, we recommend the astroplan package.

def prob_observable(m, header):
"""
Determine the integrated probability contained in a gravitational-wave
sky map that is observable from a particular ground-based site at a
particular time.

Bonus: make a plot of probability versus UTC time!
"""

# Determine resolution of sky map
npix = len(m)
nside = hp.npix2nside(npix)

# Get time now
time = astropy.time.Time.now()
# Or at the time of the gravitational-wave event...
# Or at a particular time...
# time = astropy.time.Time('2015-03-01 13:55:27')

# Geodetic coordinates of observatory (example here: Mount Wilson)
observatory = astropy.coordinates.EarthLocation(
lat=34.2247*u.deg, lon=-118.0572*u.deg, height=1742*u.m)

# Alt/az reference frame at observatory, now
frame = astropy.coordinates.AltAz(obstime=time, location=observatory)

# Look up (celestial) spherical polar coordinates of HEALPix grid.
theta, phi = hp.pix2ang(nside, np.arange(npix))
# Convert to RA, Dec.

# Transform grid to alt/az coordinates at observatory, now

# Where is the sun, now?
sun_altaz = astropy.coordinates.get_sun(time).transform_to(altaz)

# How likely is it that the (true, unknown) location of the source
# is within the area that is visible, now? Demand that sun is at
# least 18 degrees below the horizon and that the airmass
# (secant of zenith angle approximation) is at most 2.5.
prob = m[(sun_altaz.alt <= -18*u.deg) & (altaz.secz <= 2.5)].sum()

# Done!
return prob


Finally, we need to update our GCN handler to call this function:

@gcn.handlers.include_notice_types(
gcn.notice_types.LVC_PRELIMINARY,
gcn.notice_types.LVC_INITIAL,
gcn.notice_types.LVC_UPDATE)
# Respond only to 'test' events.
# VERY IMPORTANT! Replce with the following line of code
# to respond to only real 'observation' events.
# if root.attrib['role'] != 'observation':
#    return
if root.attrib['role'] != 'test':
return

# Respond only to 'CBC' events. Change 'CBC' to "Burst'
# to respond to only unmodeled burst events.
if root.find(".//Param[@name='Group']").attrib['value'] != 'CBC':
return

skymap_url = root.find(".//Param[@name='skymap_fits']").attrib['value']

print('Source has a {:d}% chance of being observable now'.format(
int(round(100 * prob))))
if prob > 0.5:
pass # FIXME: perform some action


Let’s run the new GCN handler now…

# Listen for GCNs until the program is interrupted
# (killed or interrupted with control-C).
gcn.listen(handler=process_gcn)


When you run this script, each time you receive a sample LIGO/Virgo GCN Notice, it will print something like the following (note that probability will change as a function of time):

Source has a 76% chance of being observable now