/
hosts.py
253 lines (191 loc) · 7.17 KB
/
hosts.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
""" Module related to host galaxies of FRBs
Warning: Might get chopped up into pieces sommeday
"""
import numpy as np
import pdb
import os
import pandas
from astropy import units
from astropy.io import fits
from astropy.coordinates import SkyCoord, match_coordinates_sky
from scipy import interpolate
from scipy.integrate import quad
from pkg_resources import resource_filename
def chance_coincidence(rmag, r_i):
"""
Calculate the chance probability of a galaxy to an FRB
Taken from Bloom et al. 2002
https://ui.adsabs.harvard.edu/abs/2002AJ....123.1111B/abstract
..todo.. Expand to allow for other filters
Args:
rmag (float): r-band magnitude
r_i (Angle or Quantity):
Effective radius, angular
Should be the max[2Rhalf, 3 sigma_r0, (R_0^2 + 4 Rhalf^2)^1/2]
See Bloom et al. 2002
Returns:
float: Probability of a chance association
"""
# WHERE DOES THIS EQUATION COME FROM?
sigma = 1. / (3600. ** 2 * 0.334 * np.log(10)) * 10 ** (0.334 * (rmag - 22.963) + 4.320)
# Do it
eta = np.pi * r_i.to('arcsec').value ** 2 * sigma
Pch = 1. - np.exp(-eta)
# Return
return Pch
def chance_dx(rmag):
"""
Returns the angular separation for a secure association (1%)
as in https://ui.adsabs.harvard.edu/abs/2014MNRAS.437.1495T/abstract
Args:
rmag (float):
r-band magnitude
Returns:
Quantity: Angular offset in arcsec
"""
dx = 1.48 * 10**13 * rmag**(-9.53) * units.arcsec
#
return dx
def random_separation(catalog, wcs, npix, trim=1*units.arcmin, ntrial=100):
"""
Find random offsets to
Args:
catalog (astropy.table.Table):
wcs (astropy.WCS.WCS):
npix (int):
trim (astropy.units.Quantity, optional):
ntrial (int, optional):
Returns:
astropy.units.Quantity: angular distances
"""
# Catalog
cat_coord = SkyCoord(ra=catalog['ra'], dec=catalog['dec'], unit='deg')
# Trim
bottom_corner = wcs.pixel_to_world(0, 0)
bottom_offset = bottom_corner.directional_offset_by(-45.*units.deg, trim*np.sqrt(2))
x0,y0 = [float(i) for i in wcs.world_to_pixel(bottom_offset)]
top_corner = wcs.pixel_to_world(npix-1, npix-1)
top_offset = top_corner.directional_offset_by(135.*units.deg, trim*np.sqrt(2))
x1,y1 = [float(i) for i in wcs.world_to_pixel(top_offset)]
# Generate a uniform grid
ndim = int(np.sqrt(ntrial))
xval = np.outer(np.linspace(x0, x1, ndim), np.ones(ndim))
yval = np.outer(np.ones(ndim), np.linspace(y0, y1, ndim))
# Coordinates now
grid_coord = wcs.pixel_to_world(xval.flatten(), yval.flatten())
# Match
idx, d2d, d3d = match_coordinates_sky(grid_coord, cat_coord, nthneighbor=1)
return d2d
def get_R(R_frb, R_0=0.2, R_h=0.25):
"""
Calculates Radius of localisation region in arcsecond
Based on Bloom et al 2002 and Eftekhari et al 2017
Args:
R_frb (float): The 1 sigma localization radius of the FRB
R_0 (float): Radial angular separation between the FRB position and a presumed host
R_h (float): Galaxy half light radius
Returns:
float: radius (in arcseconds)
"""
return np.max([2*R_frb, np.sqrt(R_0**2 + 4*R_h**2)])
def read_r_mags(data_table_path):
"""
Reads data used in Driver et al (2016).
https://iopscience.iop.org/article/10.3847/0004-637X/827/2/108
Args:
data_table_path (string): Path to the fits file with data
Returns:
array: r band magnitudes
array: magnitude bin
array: cosmic variance
"""
table = fits.open(data_table_path)
data = table[1].data
r_mask = np.concatenate((np.where(data['Filtername'] == 'r')[0],
np.where(data['Filtername'] == 'F606W')[0]))
_magbin = data['MagBinCentre'][r_mask]
indexes = _magbin.argsort()
magbin = _magbin[indexes]
r_band_data = data['N(m)'][r_mask][indexes]
cv = data['CosmicVariance'][r_mask][indexes]
mag_uniq = np.unique(magbin)
cvs = []
r_dat = []
for mag in mag_uniq:
loc = np.where(magbin == mag)[0]
if len(loc) > 1:
cv_s = cv[loc]
min_cv_loc = np.where(cv_s == cv_s.min())
r_dat.append(r_band_data[loc[min_cv_loc]][0])
cvs.append(cv_s[min_cv_loc][0])
else:
cvs.append(cv[loc][0])
r_dat.append(r_band_data[loc][0])
cvs = np.array(cvs)
r_dat = np.array(r_dat)
return r_dat, mag_uniq, cvs
def prob_eb17(R_frb, m, R_0=0.2, R_h=0.25, ret_numgal=False):
"""
Calculates chance association probability of a galaxy to an FRB
Taken from:
Section 2 of https://ui.adsabs.harvard.edu/abs/2017ApJ...849..162E/abstract
Args:
R_frb (float): The 1 sigma localization radius of the FRB in arcsec
m (float): r band magnitude of the galaxy
R_0 (float): Radial angular separation between the FRB position and a presumed host
R_h (float): Galaxy half light radius
ret_numgal (bool): to return the number of galaxies along with the chance coincidence probability
Returns:
float: Probability of chance coincidence
"""
r_dat, mag_uniq, cvs = read_r_mags(resource_filename('frb',os.path.join('data','Galaxies','driver2016_t3data.fits')))
spl = interpolate.UnivariateSpline(x=mag_uniq,
y=np.log10(r_dat),
bbox=[-100, 100],
k=3)
def n_gal(m_r):
return 10 ** spl(m_r)
num_dens_gal = quad(n_gal, 0, m)[0]
R = get_R(R_frb, R_0, R_h)
deg2arcsec = 60 * 60
num_gals = np.pi * (R / deg2arcsec) ** 2 * num_dens_gal
if ret_numgal:
return 1 - np.exp(-1 * num_gals), num_gals
else:
return 1 - np.exp(-1 * num_gals)
def load_host_tbl(hosts_file:str=None, host_tbl:pandas.DataFrame=None):
""" Generate a simple host table from a CSV, usually
the public file
Args:
hosts_file (str, optional): [description]. Defaults to None.
host_tbl ([type], optional): [description]. Defaults to None.
Returns:
pandas.DataFrame: [description]
"""
galaxy_path = os.path.join(resource_filename('frb', 'data'),
'Galaxies')
if host_tbl is None:
if hosts_file is None:
hosts_file = os.path.join(galaxy_path, 'public_hosts.csv')
host_tbl = pandas.read_csv(hosts_file)
# Reformat a few columns
sfrbs = [str(ifrb) for ifrb in host_tbl.FRB.values]
host_tbl['FRB'] = sfrbs
# Return
return host_tbl
def load_Mr_pdf(pdf_file:str=None):
""" Load the PDF for Mr for Host galaxies
Generated by M. Bhardwaj from ~20 localized FRBs
Args:
pdf_file (str, optional):
Filename for the PDF. Defaults to None.
Returns:
tuple: Mr, PDF
"""
if pdf_file is None:
pdf_file = os.path.join(resource_filename('frb', 'data'),
'Galaxies', 'PDF_Mr.csv')
#host galaxy M_r
df = pandas.read_csv(pdf_file)
#
return df.Mr.values, df.PDF.values #Host Galaxy PDF using 33 localized FRBs