-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathIMBH.hs
More file actions
216 lines (170 loc) · 7.31 KB
/
IMBH.hs
File metadata and controls
216 lines (170 loc) · 7.31 KB
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
module IMBH where
{-
// Description
calcEta
Calculation of the symmetric mass ratio
input:
mass1 : mass of a star [unit of solar mass]
mass2 : mass of a star [unit of solar mass]
output:
symmetric mass ratio
calcTotalmass
Calculation of the total mass of a binary system.
input:
mass1 : mass of a star [unit of solar mass]
mass2 : mass of a star [unit of solar mass]
ooutput:
total mass of the binary sistem
rhoImbh
calculation of a signal-to-noise ratio of IMBH with a given binary mass, distance, detector spectrum data.
The lower cut off frequency is fixed to 5 Hz.
input :
mass1 : mass of a star [unit of solar mass]
mass2 : mass of a star [unit of solar mass]
distance : distance from the earth [unit of Mega parsec]
spectrumData : [(frequency, powerspectrum)]
output:
signal-to-noise raio
rhodistImbh
calculation of a SNR*distance of IMBH with a given binary mass, detector spectrum data.
The lower cut off frequency is fixed to 5 Hz.
input :
mass1 : mass of a star [unit of solar mass]
mass2 : mass of a star [unit of solar mass]
spectrumData : [(frequency, powerspectrum)]
output:
SNR*distance
where
SNR : signal-to-noise raio
distance : distance from the earth [unit of Mega parsec]
distImbh
calculation of a distance of IMBH with a given binary mass, SNR=8, detector spectrum data.
The lower cut off frequency is fixed to 5 Hz. The SNR is fixed 8.
input :
mass1 : mass of a star [unit of solar mass]
mass2 : mass of a star [unit of solar mass]
spectrumData : [(frequency, powerspectrum)]
output:
distance [unit of Mpc]
// Test code for IMBHMon
import IMBH
main :: IO Double
main = do
contents <- readFile "prebKAGRA.dat"
let spectrumData = map ((\[x, y] -> (read x::Double, read y**2::Double)) . words) (lines contents)
return $ distImbh 100 100 spectrumData
// reference:
"Template bank for gravitational waveforms from coalescing binary black holes: Nonspinning binaries"
P. Ajith++
Physical Review D, vol. 77, Issue 10, id. 104017
[0] http://adsabs.harvard.edu/abs/2008PhRvD..77j4017A
[1] http://arxiv.org/abs/0710.2335
-}
-- physical cpnstant
msolar_time :: Double
msolar_time = 4.92579497077314E-6
parsec_sec :: Double
parsec_sec = 1.0292712503E8
mparsec_sec :: Double
mparsec_sec = parsec_sec*1E6
data ParamType = F_MERGE | F_RING| SIGMA | F_CUT deriving (Eq)
calcEta :: Double -> Double -> Double
calcEta mass1 mass2 = mass1*mass2/(calcTotalmass mass1 mass2)**2
calcTotalmass :: Double -> Double -> Double
calcTotalmass mass1 mass2 = mass1 + mass2
getWavParam :: (Double, Double, Double) -> Double -> Double -> Double
getWavParam (a, b, c) var_eta var_tmass = (a*var_eta**2 + b*var_eta + c)/(pi*var_tmass)
--from Eq.(4.18) in [1]
setParam :: ParamType -> (Double, Double, Double)
setParam ptype
| ptype == F_MERGE = (2.6740 * 1E-1, 4.4810 * 1E-2, 9.5560 * 1E-2)
| ptype == F_RING = (5.9411 * 1E-1, 8.9794 * 1E-2, 1.9111 * 1E-1)
| ptype == SIGMA = (5.0801 * 1E-1, 7.7515 * 1E-2, 2.2369 * 1E-2)
| ptype == F_CUT = (8.4845 * 1E-1, 1.2848 * 1E-1, 2.7299 * 1E-1)
-- from TABLE I in [1]
getF_merge :: Double -> Double -> Double
getF_merge = getWavParam (setParam F_MERGE)
getF_ring :: Double -> Double -> Double
getF_ring = getWavParam (setParam F_RING)
getSigma :: Double -> Double -> Double
getSigma = getWavParam (setParam SIGMA)
getF_cut :: Double -> Double -> Double
getF_cut = getWavParam (setParam F_CUT)
calcRhoCoef :: Double -> Double -> Double-> Double-> Double
calcRhoCoef var_tmass var_fmerg var_eta dist = var_tmass**(5/6)*var_fmerg**(-7/6)/dist/pi**(2/3)*(5*var_eta/6)**(1/2)
-- from Eq.(B11) in [1]
calcRhoDistCoef :: Double -> Double -> Double -> Double
calcRhoDistCoef var_tmass var_fmerg var_eta = var_tmass**(5/6)*var_fmerg**(-7/6)/pi**(2/3)*(5*var_eta/6)**(1/2)
calcRhoInsp :: Double -> Double -> [(Double, Double)] -> Double
calcRhoInsp var_flow var_fmerg spectrumData = sum $ zipWith (*) shList df
where shList = [(f/var_fmerg)**(-7/3)/sh | (f, sh)<-spectrumData, f>var_flow, f<var_fmerg]
fList = [f | (f, _)<-spectrumData, f>var_flow, f<var_fmerg]
df = zipWith (-) (tail fList) (init fList)
-- from Eq.(B11) in [1]
calcRhoMerg :: Double -> Double -> [(Double, Double)] -> Double
calcRhoMerg var_fmerg var_fring spectrumData = sum $ zipWith (*) shList df
where shList = [(f/var_fmerg)**(-4/3)/sh | (f, sh)<-spectrumData, f>var_fmerg, f<var_fring]
fList = [f | (f, _)<-spectrumData, f>var_fmerg, f<var_fring]
df = zipWith (-) (tail fList) (init fList)
-- from Eq.(B11) in [1]
calcRhoRing :: Double -> Double -> Double -> Double -> [(Double, Double)] -> Double
calcRhoRing var_fmerg var_fring var_fcut sigma spectrumData = w**2 * sum (zipWith (*) shList df)
where w = pi*sigma/2*(var_fring/var_fmerg)**(-2/3)
fList = [f | (f, _)<-spectrumData, f>var_fring, f<var_fcut]
shList = [lorentz f**2/sh | (f, sh)<-spectrumData, f>var_fmerg, f<var_fring]
where lorentz f = 1/(2*pi)*sigma/((f-var_fring)**2+sigma**2/4)
df = zipWith (-) (tail fList) (init fList)
-- from Eq.(B11) in [1]
rhoImbhCore :: Double -> Double -> Double -> Double -> [(Double, Double)]-> Double
rhoImbhCore flow mass1 mass2 distMPC spectrumData = do
let mass1ST = mass1*msolar_time
mass2ST = mass2*msolar_time
eta = calcEta mass1ST mass2ST
tmass = calcTotalmass mass1ST mass2ST
fmerg = getF_merge eta tmass
fring = getF_ring eta tmass
fcut = getF_cut eta tmass
sigma = getSigma eta tmass
dist = distMPC*mparsec_sec
rhoCoef = calcRhoCoef tmass fmerg eta dist
rhoInsp = calcRhoInsp flow fmerg spectrumData
rhoMerg = calcRhoMerg fmerg fring spectrumData
rhoRing = calcRhoRing fmerg fring fcut sigma spectrumData
rhoCoef * (rhoCoef + rhoInsp + rhoMerg + rhoRing)**(1/2)
-- from Eq.(B11) in [1]
rhoImbh :: Double -> Double -> Double -> [(Double, Double)]-> Double
rhoImbh = rhoImbhCore 5
rhodistImbhCore :: Double -> Double -> Double -> [(Double, Double)]-> Double
rhodistImbhCore flow mass1 mass2 spectrumData = do
let mass1ST = mass1*msolar_time
mass2ST = mass2*msolar_time
eta = calcEta mass1ST mass2ST
tmass = calcTotalmass mass1ST mass2ST
fmerg = getF_merge eta tmass
fring = getF_ring eta tmass
fcut = getF_cut eta tmass
sigma = getSigma eta tmass
rhodistCoef = calcRhoDistCoef tmass fmerg eta
rhoInsp = calcRhoInsp flow fmerg spectrumData
rhoMerg = calcRhoMerg fmerg fring spectrumData
rhoRing = calcRhoRing fmerg fring fcut sigma spectrumData
rhodistCoef * (rhoInsp + rhoMerg + rhoRing)**(1/2)
rhodistImbh :: Double -> Double -> [(Double, Double)]-> Double
rhodistImbh = rhodistImbhCore 5
distImbhCore :: Double -> Double -> Double -> [(Double, Double)]-> Double
distImbhCore flow mass1 mass2 spectrumData = do
let mass1ST = mass1*msolar_time
mass2ST = mass2*msolar_time
eta = calcEta mass1ST mass2ST
tmass = calcTotalmass mass1ST mass2ST
fmerg = getF_merge eta tmass
fring = getF_ring eta tmass
fcut = getF_cut eta tmass
sigma = getSigma eta tmass
rhodistCoef = calcRhoDistCoef tmass fmerg eta
rhoInsp = calcRhoInsp flow fmerg spectrumData
rhoMerg = calcRhoMerg fmerg fring spectrumData
rhoRing = calcRhoRing fmerg fring fcut sigma spectrumData
rhodistCoef * (rhoInsp + rhoMerg + rhoRing)**(1/2) /8.0 /mparsec_sec
distImbh :: Double -> Double -> [(Double, Double)]-> Double
distImbh = distImbhCore 5