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/KAGRA 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...
    # time = astropy.time.Time(header['MJD-OBS'], format='mjd')
    # 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.
    radecs = astropy.coordinates.SkyCoord(
        ra=phi*u.rad, dec=(0.5*np.pi - theta)*u.rad)

    # Transform grid to alt/az coordinates at observatory, now
    altaz = radecs.transform_to(frame)

    # 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)
def process_gcn(payload, root):
    # 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']

    skymap, header = hp.read_map(skymap_url, h=True)
    prob = prob_observable(skymap, header)
    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/KAGRA 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