pylal.spawaveApp

1 #spawaveApp.py 2 #Created by Satya Mohapatra on 2/26/10. 3 #Copyright (C) 2010 Satya Mohapatra 4 # This program is free software; you can redistribute it and/or modify it 5 # under the terms of the GNU General Public License as published by the 6 # Free Software Foundation; either version 2 of the License, or (at your 7 # option) any later version. 8 # 9 # This program is distributed in the hope that it will be useful, but 10 # WITHOUT ANY WARRANTY; without even the implied warranty of 11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 12 # Public License for more details. 13 # 14 # You should have received a copy of the GNU General Public License along 15 # with this program; if not, write to the Free Software Foundation, Inc., 16 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 17 18 19 # 20 # ============================================================================= 21 # 22 # Preamble 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 #interpolatation = interpolate.splrep(design[:,0],design[:,1],s=0) 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 #N = interpolate.splev(x,interpolatation) 84 SNR = math.sqrt(numpy.sum(numpy.divide(numpy.square(S),numpy.square(N))))/d 85 return SNR

86

87 -def IMRhrss(m1,m2,spin1z,spin2z,d):

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

112 -def IMRpeakAmp(m1,m2,spin1z,spin2z,d):

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 #s = numpy.abs(s) 134 s = numpy.real(s) 135 max = numpy.max(s)/d 136 return max

137

138 -def IMRtargetburstfreq(m1,m2,spin1z,spin2z):

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 #S = numpy.real(s) 159 S = numpy.abs(s) 160 x = scipy.linspace(fFinal, imrfFinal, numpy.size(S)) 161 N = interpolatation(x) 162 #N = interpolate.splev(x,interpolatation) 163 ratio = numpy.divide(numpy.square(S),numpy.square(N)) 164 #ratio = numpy.divide(S,N) 165 maxindex = numpy.argmax(ratio) 166 maxfreq = x[maxindex] 167 return maxfreq

168 169

170 -def IMRfinalspin(m1, m2,spin1z,spin2z):

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

Source Code for Module pylal.spawaveApp