1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 """
18 This module provides a bunch of user-friendly wrappers to the SWIG-bound REAL8TimeSeries and REAL8FrequencySeries objects and their associated functions.
19 """
20
21 from __future__ import division
22
23 import numpy
24 from scipy import integrate
25 import lal
26 from pylal import git_version
27
28 __author__ = "Duncan M. Macleod <duncan.macleod@ligo.org>"
29 __version__ = git_version.id
30 __date__ = git_version.date
31
32
33
34
35
36 -def inspiralrange(f, S, snr=8, m1=1.4, m2=1.4, fmin=10, fmax=None,\
37 horizon=False):
38 """
39 Calculate the sensitivity distance to an inspiral with the given
40 masses, for the given signal-to-noise ratio. See the following
41 reference for information on the integral performed:
42
43 https://dcc.ligo.org/cgi-bin/private/DocDB/ShowDocument?docid=27267
44
45 @param f: frequency array
46 @type f: C{numpy.array}
47 @param S: power spectral density array
48 @type S: C{numpy.array}
49 @param snr: signal-to-noise ratio at which to calculate range
50 @type snr: C{float}
51 @param m1: mass (in solar masses) of first binary component,
52 default: 1.4
53 @type m1: C{float}
54 @param m2: mass (in solar masses) of second binary component,
55 default: 1.4
56 @type m2: C{float}
57 @param fmin: minimum frequency limit of integral, default: 10 Hz
58 @type fmin: C{float}
59 @param fmax: maximum frequency limit of integral, default: ISCO
60 @type fmax: C{float}
61 @param horizon: return horizon distance in stead of angle-averaged,
62 default: False
63 @type horizon: C{bool}
64
65 @return: sensitive distance to an inspiral (in solar masses) for the
66 given PSD and parameters
67 @rtype: C{float}
68 """
69
70 mtot = m1+m2
71 mchirp = (m1*m2)**(0.6) / mtot**(0.4)
72
73
74 mchirp *= lal.LAL_MSUN_SI
75 pre = (5 * lal.LAL_C_SI**(1/3) *\
76 (mchirp * lal.LAL_G_SI / lal.LAL_C_SI**2)**(5/3) * 1.77**2) /\
77 (96 * lal.LAL_PI ** (4/3) * snr**2)
78
79
80
81 fisco = lal.LAL_C_SI**3\
82 / (lal.LAL_G_SI * lal.LAL_MSUN_SI * 6**1.5 * lal.LAL_PI * mtot)
83 if not fmax:
84 fmax = fisco
85 elif fmax > fisco:
86 warnings.warn("Upper frequency bound greater than %s-%s ISCO "\
87 "frequency of %.2g Hz, using ISCO" % (m1,m2,fisco))
88 fmax = fisco
89
90
91 condition = (f >= fmin) & (f < fmax)
92 integrand = f[condition]**(-7/3)/S[condition]
93 result = integrate.trapz(integrand, f[condition])
94
95 d = (pre*result) ** (1/2) / (lal.LAL_PC_SI*1e6)
96 return d
97
98
99
100
101
102 -def squeezing(f, S, dc=None, opgain=None, cavpf=None, fmin=1000, fmax=2000):
103 """
104 Calculate squeezing factor based on observed noise in given band
105 and predicted shot noise
106
107 @param f: frequency array
108 @type f: C{numpy.array}
109 @param S: power spectral density array
110 @type S: C{numpy.array}
111 @param dc:
112 @type dc: C{float}
113 @param opgain:
114 @type opgain: C{float}
115 @param cavpf:
116 @type cavpf: C{float}
117 @param fmin: minimum frequency limit of integral, default: 1000 Hz
118 @type fmin: C{float}
119 @param fmax: maximum frequency limit of integral, default: 2000 Hz
120 @type fmax: C{float}
121
122 @return: squeezing ratio in dB
123 @rtype: C{float}
124 """
125
126
127 condition = (f >= fmin) & (f < fmax)
128
129
130 model = abs(1+1j*f[condition]/cavpf) * (dc)**(1/2) /opgain
131
132
133 d = numpy.median(model /S[condition]**(1/2))
134
135
136 d = 20 * numpy.log10(d)
137
138 return d
139
140
141
142
143
145 """
146 Calculate the sensitive distance to a GW burst with the given intrinsic
147 energy for the given signal-to-noise ratio snr, as a function of f.
148
149 @param f: frequency of interest
150 @type f: C{float} or C{numpy.array}
151 @param S: power spectral density
152 @type S: C{numpy.array}
153 @param snr: signal-to-noise ratio of burst
154 @type snr: C{float}
155 @param E: instrinsic energy of burst, default: grb-like 0.01
156 @type E: C{float}
157
158 @return: sensitive distance in pc at which a GW burst with the given
159 energy would be detected with the given SNR, as a function of
160 it's frequency
161 @rtype: C{float}
162 """
163 A = ((lal.LAL_G_SI * E * lal.LAL_MSUN_SI * 0.4)\
164 / (lal.LAL_PI**2 * lal.LAL_C_SI))**(1/2) / lal.LAL_PC_SI
165 return A / (snr * S**(1/2) * f)
166
167 -def burstrange(f, S, snr=8, E=1e-2, fmin=0, fmax=None, unit="Mpc"):
168 """
169 Calculate the sensitive distance to a GW burst with the given intrinsic
170 energy for the given signal-to-noise ratio snr, integrated over frequency.
171
172 @param f: frequency of interest
173 @type f: C{float} or C{numpy.array}
174 @param S: power spectral density
175 @type S: C{numpy.array}
176 @param snr: signal-to-noise ratio of burst, default: 8
177 @type snr: C{float}
178 @param E: instrinsic energy of burst, default: grb-like 0.01
179 @type E: C{float}
180 @param fmin: minimum frequency limit of integral, default: 10 Hz
181 @type fmin: C{float}
182 @param fmax: maximum frequency limit of integral, default: ISCO
183 @type fmax: C{float}
184
185 @return: sensitive distance at which a GW burst with the given
186 energy would be detected with the given SNR, integrated over
187 frequency.
188 @rtype: C{float}
189 """
190 if not fmin:
191 fmin = f.min()
192 if not fmax:
193 fmax = f.max()
194
195
196 condition = (f >= fmin) & (f < fmax)
197
198
199 integrand = fdependent_burstrange(f[condition], S[condition], snr, E)**3
200 result = integrate.trapz(integrand, f[condition])
201
202 d = (result / (fmax-fmin))**(1/3)
203 if unit == "Mpc":
204 d = d/1e6
205 elif unit == "kpc":
206 d = d/1e3
207 elif unit == "pc":
208 d = d
209 else:
210 raise ValueError("Unrecognized unit: %s" % unit)
211
212 return d
213