Sort My Fields [public]
Which of my telescope fields contain the greatest probability of containing the counterpart? The code below reads in a list of fields and a skymap, samples the skymap for the probability at each field center, and sorts into a list with decreasing probability.
# Rank a list of (right ascension, declination) points according to descending
# probability density. See attached file for example input.
import os
import healpy as hp
import numpy as np
# Download list of RA, Dec points.
os.system('curl -O https://gw-astronomy.org/wiki/pub/LV_EM/SortMyFields/CRTS_fields2.txt')
# Download example sky map.
os.system('http://www.ligo.org/scientists/first2years/2015/compare/18951/bayestar.fits.gz')
# Read RA, Dec list.
fields = np.recfromtxt('CRTS_fields2.txt', names='name,ra,dec')
# Read sky map.
skymap = hp.read_map('bayestar.fits.gz', verbose=False)
# Determine HEALPix resolution.
nside = hp.npix2nside(len(skymap))
deg2perpix = hp.nside2pixarea(nside, degrees=True)
# Convert from RA, Dec in degrees to spherical polar coordinates.
theta = 0.5 * np.pi - np.deg2rad(fields['dec'])
phi = np.deg2rad(fields['ra'])
# Interpolate sky map at each point.
# Note that the sky maps are deliberately over-sampled, so alternatively you
# could just evaluate the probability at the nearest neighbor of each point,
# like this:
#
# prob = skymap[hp.ang2pix(theta, phi)]
prob = hp.get_interp_val(skymap, theta, phi)
# Rank points by descending posterior probability density.
prob_fields = sorted(zip(prob, fields), reverse=True)
# Print top 20 fields.
print('Field RA Dec Prob. / deg^2')
for p, field in prob_fields[:20]:
print('{0:s} {1:8.4f} {2:8.4f} {3:.4f}'.format(
field['name'], field['ra'], field['dec'], p / deg2perpix))
--
RoyWilliams - 21 Apr 2015