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 This module is an enhancement based on the module spawaveform available in pylal.
29 This module can do the following:
30 1) Find the expected SNR of a spin-aligned IMR waveform (see Ajith et al.) with respect to initial LIGO SRD.
31 submodule:IMRsnr
32 2) Find the hrss of the waveform.
33 submodule:IMRhrss
34 3) Find the peak amplitude of the waveform.
35 submodule:IMRpeakAmp
36 4) Find the target burst frequency of the waveform:
37 submodule:IMRtargetburstfreq
38
39 """
40
41
42
43 from pylal import spawaveform
44 import math
45 from scipy import interpolate
46 import scipy
47 import numpy
48 import time
49
50
51 from pylal import git_version
52 __author__ = "Satya Mohapatra <satya@physics.umass.edu>"
53 __version__ = "git id %s" % git_version.id
54 __date__ = git_version.date
55
56 design = numpy.loadtxt('ligo4k_design.txt')
57 interpolatation = interpolate.interp1d(design[:,0],design[:,1])
58
59
60 -def IMRsnr(m1,m2,spin1z,spin2z,d):
61
62 """
63 IMRsnr finds the SNR of the waveform with Initial LIGO SRD for a given source parameters and the source distance.
64 usage: IMRsnr(m1,m2,spin1z,spin2z,distance)
65 e.g.
66 spawaveApp.IMRsnr(30,40,0.45,0.5,100)
67
68 """
69
70 chi = spawaveform.computechi(m1, m2, spin1z, spin2z)
71 imrfFinal = spawaveform.imrffinal(m1, m2, chi, 'fcut')
72 fLower = 10.0
73 order = 7
74 dur = 2**numpy.ceil(numpy.log2(spawaveform.chirptime(m1,m2,order,fLower)))
75 sr = 2**numpy.ceil(numpy.log2(imrfFinal*2))
76 deltaF = 1.0 / dur
77 deltaT = 1.0 / sr
78 s = numpy.empty(sr * dur, 'complex128')
79 spawaveform.imrwaveform(m1, m2, deltaF, fLower, s, spin1z, spin2z)
80 S = numpy.abs(s)
81 x = scipy.linspace(fLower, imrfFinal, numpy.size(S))
82 N = interpolatation(x)
83
84 SNR = math.sqrt(numpy.sum(numpy.divide(numpy.square(S),numpy.square(N))))/d
85 return SNR
86
88
89 """
90 IMRhrss finds the SNR of the waveform for a given source parameters and the source distance.
91 usage: IMRhrss(m1,m2,spin1z,spin2z,distance)
92 e.g.
93 spawaveApp.IMRhrss(30,40,0.45,0.5,100)
94
95 """
96
97 chi = spawaveform.computechi(m1, m2, spin1z, spin2z)
98 imrfFinal = spawaveform.imrffinal(m1, m2, chi, 'fcut')
99 fLower = 10.0
100 order = 7
101 dur = 2**numpy.ceil(numpy.log2(spawaveform.chirptime(m1,m2,order,fLower)))
102 sr = 2**numpy.ceil(numpy.log2(imrfFinal*2))
103 deltaF = 1.0 / dur
104 deltaT = 1.0 / sr
105 s = numpy.empty(sr * dur, 'complex128')
106 spawaveform.imrwaveform(m1, m2, deltaF, fLower, s, spin1z, spin2z)
107 s = numpy.abs(s)
108 s = numpy.square(s)
109 hrss = math.sqrt(numpy.sum(s))/d
110 return hrss
111
113
114 """
115 IMRpeakAmp finds the peak amplitude of the waveform for a given source parameters and the source distance.
116 usage: IMRpeakAmp(m1,m2,spin1z,spin2z,distance)
117 e.g.
118 spawaveApp.IMRpeakAmp(30,40,0.45,0.5,100)
119
120 """
121
122 chi = spawaveform.computechi(m1, m2, spin1z, spin2z)
123 imrfFinal = spawaveform.imrffinal(m1, m2, chi, 'fcut')
124 fLower = 10.0
125 order = 7
126 dur = 2**numpy.ceil(numpy.log2(spawaveform.chirptime(m1,m2,order,fLower)))
127 sr = 2**numpy.ceil(numpy.log2(imrfFinal*2))
128 deltaF = 1.0 / dur
129 deltaT = 1.0 / sr
130 s = numpy.empty(sr * dur, 'complex128')
131 spawaveform.imrwaveform(m1, m2, deltaF, fLower, s, spin1z, spin2z)
132 s = scipy.ifft(s)
133
134 s = numpy.real(s)
135 max = numpy.max(s)/d
136 return max
137
139
140 """
141 IMRtargetburstfreq finds the peak amplitude of the waveform for a given source parameters and the source distance.
142 usage: IMRtargetburstfreq(m1,m2,spin1z,spin2z)
143 e.g.
144 spawaveApp.IMRtargetburstfreq(30,40,0.45,0.5)
145
146 """
147 chi = spawaveform.computechi(m1, m2, spin1z, spin2z)
148 fFinal = spawaveform.ffinal(m1,m2,'schwarz_isco')
149 imrfFinal = spawaveform.imrffinal(m1, m2, chi, 'fcut')
150 fLower = 40.0
151 order = 7
152 dur = 2**numpy.ceil(numpy.log2(spawaveform.chirptime(m1,m2,order,fFinal)))
153 sr = 2**numpy.ceil(numpy.log2(imrfFinal*2))
154 deltaF = 1.0 / dur
155 deltaT = 1.0 / sr
156 s = numpy.empty(sr * dur, 'complex128')
157 spawaveform.imrwaveform(m1, m2, deltaF, fFinal, s, spin1z, spin2z)
158
159 S = numpy.abs(s)
160 x = scipy.linspace(fFinal, imrfFinal, numpy.size(S))
161 N = interpolatation(x)
162
163 ratio = numpy.divide(numpy.square(S),numpy.square(N))
164
165 maxindex = numpy.argmax(ratio)
166 maxfreq = x[maxindex]
167 return maxfreq
168
169
171
172 """
173 The final spin calculation is based on Rezolla et al.
174 IMRfinalspin final spin of the end product black hole.
175 usage: IMRfinalspin(m1,m2,spin1z,spin2z)
176 e.g.
177 spawaveApp.IMRfinalspin(30.,40.,0.45,0.5)
178
179 """
180
181 s4 = -0.129
182 s5 = -0.384
183 t0 = -2.686
184 t2 = -3.454
185 t3 = 2.353
186 spin1x = 0.
187 spin1y = 0.
188 spin2x = 0.
189 spin2y = 0.
190 M = m1 + m2
191 q = m2/m1
192 eta = m1*m2/(m1+m2)**2
193 a1 = math.sqrt(spin1x**2 + spin1y**2 + spin1z**2)
194 a2 = math.sqrt(spin2x**2 + spin2y**2 + spin2z**2)
195 if (a1 != 0) and (a2 != 0): cosa = (spin1x*spin2x + spin1y*spin2y + spin1z*spin2z)/(a1*a2)
196 else:cosa = 0
197 if a1 != 0: cosb = spin1z/a1
198 else: cosb = 0
199 if a2 != 0: cosc = spin2z/a2
200 else: cosc = 0
201 l = (s4/(1+q**2)**2 * (a1**2 + (a2**2)*(q**4) + 2*a1*a2*(q**2)*cosa) + (s5*eta + t0 + 2)/(1+q**2) * (a1*cosb + a2*(q**2)*cosc) + 2*math.sqrt(3.) + t2*eta + t3*(eta**2))
202 afin = 1/(1+q)**2 * math.sqrt(a1**2 + (a2**2)*(q**4) + 2*a1*a2*(q**2)*cosa + 2*(a1*cosb + a2*(q**2)*cosc)*l*q + (l**2)*(q**2))
203 return afin
204