Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInference 4.1.8.1-bb0d041
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
nrutils.py
Go to the documentation of this file.
1"""
2Various fitting formulas provided by numerical relativity
3
4Archisman Ghosh, Nathan K. Johnson-McDaniel, P. Ajith, 2015-04-09
5"""
6
7import numpy as np
8import scipy.optimize as so
9try:
10 import lal
11except ImportError:
12 print('Cannot import lal SWIG bindings')
13
14# Functions to check that input values are physical
15
16def _check_m(m1,m2):
17 if np.any(m1<=0):
18 raise ValueError("m1 has to be positive")
19 if np.any(m2<=0):
20 raise ValueError("m2 has to be positive")
21
22def _check_chi(chi1, chi2):
23 if np.any(chi1>1):
24 raise ValueError("chi1 has to be <= 1")
25 if np.any(chi2>1):
26 raise ValueError("chi2 has to be <= 1")
27 if np.any(chi1<0):
28 raise ValueError("chi1 has to be nonnegative")
29 if np.any(chi2<0):
30 raise ValueError("chi2 has to be nonnegative")
31
32def _check_chi_aligned(chi1, chi2):
33 if np.any(abs(chi1)>1):
34 raise ValueError("chi1 has to be <= 1")
35 if np.any(abs(chi2)>1):
36 raise ValueError("chi2 has to be <= 1")
37
38def _check_mchi(m1,m2,chi1,chi2):
39 _check_m(m1,m2)
40 _check_chi(chi1,chi2)
41
42def _check_mchi_aligned(m1,m2,chi1,chi2):
43 _check_m(m1,m2)
44 _check_chi_aligned(chi1,chi2)
45
46###########################
47# Final mass and spin fits
48###########################
49
50# list of final mass and spin fit names
52 pan2011 = "Pan2011" # Pan et al. [Phys Rev D 84, 124052 (2011)] (mass+spin, non-spinning)
53 hlz2014 = "HLZ2014" # Healy et al. [Phys Rev D 90, 104004 (2014)] (mass+spin, aligned)
54 phenD = "PhenomD" # Husa et al. [Phys Rev D 93, 044006 (2016)] (mass+spin, aligned)
55 uib2016v1 = "UIB2016v1" # Jimenez-Forteza et al. [arXiv:1611.00332v1] (mass+spin, aligned)
56 uib2016 = "UIB2016" # Jimenez-Forteza et al. [LIGO-P1600270-v4 / arXiv:1611.00332v2] (mass+spin, aligned)
57 hbr2016 = "HBR2016" # Hofmann et al. [ApJL 825:L19 (2016)] (spin only, precessing)
58 hl2016 = "HL2016" # Healy and Lousto [arXiv:1610.09713] (mass+spin, aligned)
59
60# list of Kerr truncation behaviours
62 trunc = "TRUNCATE" # truncate fit value to 1.0, with an info message
63 trunc_silent = "TRUNCATESILENT" # truncate fit value to 1.0, without info message
64 keep = "KEEP" # keep fit value, even if superextremal, but still give the info message
65 ign = "IGNORE" # keep fit value, even if superextremal, without info message
66 err = "ERROR" # abort with an error if fit value is superextremal
67
68# Function to truncate final spins at Kerr limit (with various options)
69
70def _truncate_at_Kerr_limit(chif, behavior, fitname="this"):
71 if behavior not in list(bbh_Kerr_trunc_opts.__dict__.values()):
72 raise ValueError("Unknown user option '%s' for Kerr truncation behavior." % behavior)
73 if behavior!=bbh_Kerr_trunc_opts.ign and np.any(abs(chif)>1.0):
74 if isinstance(chif,(list,np.ndarray)):
75 idx_over = np.where(abs(chif)>1.0)
76 if behavior==bbh_Kerr_trunc_opts.trunc or behavior==bbh_Kerr_trunc_opts.trunc_silent:
77 chif_trunc = np.sign(chif[idx_over])*1.0
78 if behavior==bbh_Kerr_trunc_opts.trunc:
79 print("Truncating %d excessive chif values from %s fit to Kerr limit of +-1.0" % (np.size(idx_over), fitname))
80 chif[idx_over] = chif_trunc
81 elif behavior==bbh_Kerr_trunc_opts.keep:
82 print("Note: %s fit predicts %d chif values in excess of the Kerr limit." % (fitname, np.size(idx_over)))
83 elif behavior==bbh_Kerr_trunc_opts.err:
84 raise ValueError("%s fit predicts %d chif values in excess of the Kerr limit." % (fitname, np.size(idx_over)))
85 else:
86 if behavior==bbh_Kerr_trunc_opts.trunc or behavior==bbh_Kerr_trunc_opts.trunc_silent:
87 chif_trunc = np.sign(chif)*1.0
88 if behavior==bbh_Kerr_trunc_opts.trunc:
89 print("Truncating excessive chif of %f from %s fit to Kerr limit of %f" % (chif, fitname, chif_trunc))
90 chif = chif_trunc
91 elif behavior==bbh_Kerr_trunc_opts.keep:
92 print("Note: %s fit predicts chif of %f in excess of the Kerr limit." % (fitname, chif))
93 elif behavior==bbh_Kerr_trunc_opts.err:
94 raise ValueError("%s fit predicts chif of %f in excess of the Kerr limit." % (fitname, chif))
95 return chif
96
97# Final mass and spin functions
98
100 """
101 Calculate the mass of the final BH resulting from the merger of two non-spinning black holes using eqn. 29a of from Pan et al, Phys Rev D 84, 124052 (2011).
102
103 Parameters
104 ----------
105 m1, m2 : component masses
106
107 Returns
108 ------
109 final mass, mf
110 """
111 m1 = np.array(m1,dtype=np.float64)
112 m2 = np.array(m2,dtype=np.float64)
113
114 m = m1 + m2
115 eta = m1*m2/(m1+m2)**2.
116 return m*(1. + (np.sqrt(8./9.)-1.)*eta - 0.4333*(eta**2.) - (0.4392*(eta**3)))
117
119 """
120 Calculate the spin of the final BH resulting from the merger of two non-spinning black holes using eqn. 29b of Pan et al, Phys Rev D 84, 124052 (2011).
121
122 Parameters
123 ----------
124 m1, m2 : component masses
125
126 Returns
127 ------
128 final dimensionless spin, chif
129 """
130 m1 = np.array(m1,dtype=np.float64)
131 m2 = np.array(m2,dtype=np.float64)
132
133 eta = m1*m2/(m1+m2)**2.
134 return np.sqrt(12.)*eta - 3.871*(eta**2.) + 4.028*(eta**3)
135
136def calc_isco_radius(a):
137 """
138 Calculate the ISCO radius of a Kerr BH as a function of the Kerr parameter using eqns. 2.5 and 2.8 from Ori and Thorne, Phys Rev D 62, 24022 (2000)
139
140 Parameters
141 ----------
142 a : Kerr parameter
143
144 Returns
145 -------
146 ISCO radius
147 """
148
149 a = np.minimum(np.array(a),1.) # Only consider a <=1, to avoid numerical problems
150
151 # Ref. Eq. (2.5) of Ori, Thorne Phys Rev D 62 124022 (2000)
152 z1 = 1.+(1.-a**2.)**(1./3)*((1.+a)**(1./3) + (1.-a)**(1./3))
153 z2 = np.sqrt(3.*a**2 + z1**2)
154 a_sign = np.sign(a)
155 return 3+z2 - np.sqrt((3.-z1)*(3.+z1+2.*z2))*a_sign
156
157
158def _RIT_setup(m1, m2, chi1, chi2):
159 """ Internal function to compute the combinations of masses and spins that are used in the RIT fits """
160
161 _check_mchi_aligned(m1, m2, chi1, chi2)
162
163 m = m1+m2
164 eta = m1*m2/m/m
165 delta_m = (m1 - m2)/m
166
167 S = (m1*m1*chi1 + m2*m2*chi2)/m**2 # symmetric spin (dimensionless -- called \tilde{S} in the paper)
168 Delta = (m2*chi2 - m1*chi1)/m # antisymmetric spin (dimensionless -- called tilde{Delta} in the paper
169
170 return eta, delta_m, S, Delta
171
172def _RIT_symm_express(eta, delta_m, S, Delta, coeffs):
173 """
174 Internal function to return the basic even-in-interchange-of-particles expression used in the RIT fits for the final mass, final spin, and peak luminosity
175
176 The coefficients in the expression are passed as a vector as [0 1 2a 2b 2c 2d 3a 3b 3c 3d 4a 4b 4c 4d 4e 4f 4g 4h 4i] (giving the subscripts).
177 """
178
179 S2 = S*S
180
181 dm2 = delta_m*delta_m
182
183 Deltadm = Delta*delta_m
184
185 Delta2 = Delta*Delta
186
187 return 16.*eta*eta*(coeffs[0] + coeffs[1]*S + coeffs[2]*Deltadm + coeffs[3]*S2 + coeffs[4]*Delta2 \
188 + coeffs[5]*dm2 + coeffs[6]*Deltadm*S + coeffs[7]*S*Delta2 + coeffs[8]*S2*S \
189 + coeffs[9]*S*dm2 + coeffs[10]*Deltadm*S2 + coeffs[11]*Delta2*Deltadm \
190 + coeffs[12]*Delta2*Delta2 + coeffs[13]*S2*S2 + coeffs[14]*Delta2*S2 + coeffs[15]*dm2*dm2 + coeffs[16]*Deltadm*dm2 \
191 + coeffs[17]*Delta2*dm2 + coeffs[18]*S2*dm2)
192
193def _final_spin_diff_Healyetal(a_f, eta, delta_m, S, Delta, version):
194 """ Internal function: the final spin with the Healy et al. fits is determined by minimizing this function """
195
196 # calculate ISCO radius
197 r_isco = calc_isco_radius(a_f)
198
199 # angular momentum at ISCO -- Eq.(2.8) of Ori, Thorne Phys Rev D 62 124022 (2000)
200 J_isco = (3*np.sqrt(r_isco)-2*a_f)*2./np.sqrt(3*r_isco)
201
202 # fitting coefficients
203 if version == "2014": # From Table XI of Healy et al Phys Rev D 90, 104004 (2014) [fourth order fits]
204 L0 = 0.686710
205 L1 = 0.613247
206 L2a = -0.145427
207 L2b = -0.115689
208 L2c = -0.005254
209 L2d = 0.801838
210 L3a = -0.073839
211 L3b = 0.004759
212 L3c = -0.078377
213 L3d = 1.585809
214 L4a = -0.003050
215 L4b = -0.002968
216 L4c = 0.004364
217 L4d = -0.047204
218 L4e = -0.053099
219 L4f = 0.953458
220 L4g = -0.067998
221 L4h = 0.001629
222 L4i = -0.066693
223 elif version == "2016": # From Table III of Healy and Lousto arXiv:1610.09713 (values taken from Matlab implementation and thus slightly more precise than the ones in the table)
224 L0 = 0.686732132
225 L1 = 0.613284976
226 L2a = -0.148530075
227 L2b = -0.113826318
228 L2c = -0.00323995784
229 L2d = 0.798011319
230 L3a = -0.0687823713
231 L3b = 0.00129103641
232 L3c = -0.0780143929
233 L3d = 1.55728564
234 L4a = -0.00571010557
235 L4b = 0.005919799
236 L4c = -0.00170575554
237 L4d = -0.0588819084
238 L4e = -0.0101866693
239 L4f = 0.964444768
240 L4g = -0.11088507
241 L4h = -0.00682082169
242 L4i = -0.0816482139
243 else:
244 raise ValueError('Unknown version--should be either "2014" or "2016".')
245
246 dm2 = delta_m*delta_m
247 dm4 = dm2*dm2
248 dm6 = dm4*dm2
249
250 a_f_new = _RIT_symm_express(eta, delta_m, S, Delta, [L0, L1, L2a, L2b, L2c, L2d, L3a, L3b, L3c, L3d, L4a, L4b, L4c, L4d, L4e, L4f, L4g, L4h, L4i]) \
251 + S*(1. + 8.*eta)*dm4 + eta*J_isco*dm6
252
253 return abs(a_f-a_f_new)
254
255def bbh_final_spin_non_precessing_Healyetal(m1, m2, chi1, chi2, version="2014"):
256 """
257 Calculate the spin of the final BH resulting from the merger of two black holes with
258 non-precessing spins using fit from Healy et al Phys Rev D 90, 104004 (2014) (version == "2014")
259 or the small update from Healy and Lousto arXiv:1610.09713 (version == "2016")
260
261 Parameters
262 ----------
263 m1, m2 : component masses
264 chi1, chi2 : dimensionless spins of two BHs
265
266 Returns
267 -------
268 dimensionless final spin, chif
269 """
270 m1 = np.array(m1,dtype=np.float64)
271 m2 = np.array(m2,dtype=np.float64)
272 chi1 = np.array(chi1,dtype=np.float64)
273 chi2 = np.array(chi2,dtype=np.float64)
274 eta, delta_m, S, Delta = _RIT_setup(m1, m2, chi1, chi2)
275 def func(eta, delta_m, S, Delta):
276 return so.leastsq(_final_spin_diff_Healyetal, 0., args=(eta, delta_m, S, Delta, version))[0][0]
277 chif = np.vectorize(func)(eta, delta_m, S, Delta)
278 return chif.item() if np.ndim(chif) == 0 else chif
279
280def bbh_final_mass_non_precessing_Healyetal(m1, m2, chi1, chi2, version="2014", chif=None):
281 """
282 Calculate the mass of the final BH resulting from the merger of two black holes with non-precessing spins using fit from Healy et al Phys Rev D 90, 104004 (2014) (version == "2014") or the small update from Healy and Lousto arXiv:1610.09713 (version == "2016")
283
284 Parameters
285 ----------
286 m1, m2 : component masses
287 chi1, chi2 : dimensionless spins of two BHs
288 chif: final spin (optional), if already calculated
289
290 Returns
291 -------
292 final mass, mf
293 """
294 m1 = np.array(m1,dtype=np.float64)
295 m2 = np.array(m2,dtype=np.float64)
296 chi1 = np.array(chi1,dtype=np.float64)
297 chi2 = np.array(chi2,dtype=np.float64)
298
299 eta, delta_m, S, Delta = _RIT_setup(m1, m2, chi1, chi2)
300
301 if chif is None:
302 chif = bbh_final_spin_non_precessing_Healyetal(m1, m2, chi1, chi2, version)
303 else:
304 chif = np.array(chif)
305
306 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
307 # now compute the final mass
308 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
309 r_isco = calc_isco_radius(chif)
310
311 # fitting coefficients
312 if version == "2014": # From Table XI of Healy et al Phys Rev D 90, 104004 (2014) [fourth order fits]
313 M0 = 0.951507
314 K1 = -0.051379
315 K2a = -0.004804
316 K2b = -0.054522
317 K2c = -0.000022
318 K2d = 1.995246
319 K3a = 0.007064
320 K3b = -0.017599
321 K3c = -0.119175
322 K3d = 0.025000
323 K4a = -0.068981
324 K4b = -0.011383
325 K4c = -0.002284
326 K4d = -0.165658
327 K4e = 0.019403
328 K4f = 2.980990
329 K4g = 0.020250
330 K4h = -0.004091
331 K4i = 0.078441
332 elif version == "2016": # From Table III of Healy and Lousto arXiv:1610.09713 (values taken from Matlab implementation and thus slightly more precise than the ones in the table)
333 M0 = 0.951659087
334 K1 = -0.0511301363
335 K2a = -0.00569897591
336 K2b = -0.0580644933
337 K2c = -0.00186732281
338 K2d = 1.99570464
339 K3a = 0.00499137602
340 K3b = -0.00923776244
341 K3c = -0.120577082
342 K3d = 0.0164168385
343 K4a = -0.0607207285
344 K4b = -0.00179828653
345 K4c = 0.000654388173
346 K4d = -0.156625642
347 K4e = 0.0103033606
348 K4f = 2.97872857
349 K4g = 0.00790433045
350 K4h = 0.000631241195
351 K4i = 0.0844776942
352 else:
353 raise ValueError('Unknown version--should be either "2014" or "2016".')
354
355 # binding energy at ISCO -- Eq.(2.7) of Ori, Thorne Phys Rev D 62 124022 (2000)
356 E_isco = (1. - 2./r_isco + chif/r_isco**1.5)/np.sqrt(1. - 3./r_isco + 2.*chif/r_isco**1.5)
357
358 dm2 = delta_m*delta_m
359 dm6 = dm2*dm2*dm2
360
361 m = m1 + m2
362
363 # final mass -- Eq. (14) of Healy et al Phys Rev D 90, 104004 (2014)
364 mf = _RIT_symm_express(eta, delta_m, S, Delta, [M0, K1, K2a, K2b, K2c, K2d, K3a, K3b, K3c, K3d, K4a, K4b, K4c, K4d, K4e, K4f, K4g, K4h, K4i]) \
365 + (1+eta*(E_isco+11.))*dm6
366
367 return mf*m
368
369def bbh_final_spin_projected_spin_Healyetal(m1, m2, chi1, chi2, tilt1, tilt2):
370 """
371 wrapper to bbh_final_spin_projected_spins() for backwards compatibility
372 """
373 return bbh_final_spin_projected_spins(m1, m2, chi1, chi2, tilt1, tilt2, "HLZ2014")
374
375def bbh_final_mass_projected_spin_Healyetal(m1, m2, chi1, chi2, tilt1, tilt2, chif=None):
376 """
377 wrapper to bbh_final_mass_projected_spins() for backwards compatibility
378 """
379 return bbh_final_mass_projected_spins(m1, m2, chi1, chi2, tilt1, tilt2, "HLZ2014", chif)
380
381def bbh_final_spin_precessing_Healyetal_extension_Minit(m1, m2, chi1, chi2, tilt1, tilt2, phi12):
382 """
383 wrapper to bbh_final_spin_precessing() for backwards compatibility
384 """
385 return bbh_final_spin_precessing(m1, m2, chi1, chi2, tilt1, tilt2, phi12, "HLZ2014")
386
387def bbh_final_mass_projected_spins(m1, m2, chi1, chi2, tilt1, tilt2, fitname, chif=None):
388 """
389 Calculate the mass of the final BH resulting from the merger of two black holes,
390 only using the projected spins along the total angular momentum
391 and some aligned-spin fit from the literature
392
393 Parameters
394 ----------
395 m1, m2 : component masses
396 chi1, chi2 : dimensionless spins of two BHs
397 tilt1, tilt2 : tilts (in radians) in the new spin convention
398 fitname: fit selection currently supports Pan2011 (non-spinning), HLZ2014, PhenomD, UIB2016, HL2016
399 chif: final spin (optional, only used for HLZ2014 and HL2016), if already calculated
400
401 Returns
402 -------
403 final mass, mf
404 """
405
406 m1 = np.array(m1,dtype=np.float64)
407 m2 = np.array(m2,dtype=np.float64)
408 chi1 = np.array(chi1,dtype=np.float64)
409 chi2 = np.array(chi2,dtype=np.float64)
410 tilt1 = np.array(tilt1,dtype=np.float64)
411 tilt2 = np.array(tilt2,dtype=np.float64)
412 if chif is not None:
413 chif = np.array(chif,dtype=np.float64)
414
415 _check_mchi(m1,m2,chi1,chi2) # Check that inputs are physical
416
417 chi1proj = chi1*np.cos(tilt1)
418 chi2proj = chi2*np.cos(tilt2)
419
420 if fitname==bbh_final_state_fits.pan2011:
421 if np.any(chi1!=0) or np.any(chi2!=0) or np.any(tilt1!=0) or np.any(tilt2!=0):
422 print("Note: Pan2011 fit does not use spins.")
423 if chif is not None:
424 print("Note: Precomputed chif not used by this fit.")
426 elif fitname==bbh_final_state_fits.hlz2014:
427 mf = bbh_final_mass_non_precessing_Healyetal(m1, m2, chi1proj, chi2proj, version="2014", chif=chif)
428 elif fitname==bbh_final_state_fits.hl2016:
429 mf = bbh_final_mass_non_precessing_Healyetal(m1, m2, chi1proj, chi2proj, version="2016", chif=chif)
430 elif fitname==bbh_final_state_fits.phenD:
431 if chif is not None:
432 print("Note: Precomputed chif not used by this fit.")
433 mf = bbh_final_mass_non_precessing_Husaetal(m1, m2, chi1proj, chi2proj)
434 elif fitname==bbh_final_state_fits.uib2016:
435 if chif is not None:
436 print("Note: Precomputed chif not used by this fit.")
437 mf = bbh_final_mass_non_precessing_UIB2016(m1, m2, chi1proj, chi2proj, version="v2")
438 elif fitname==bbh_final_state_fits.uib2016v1:
439 if chif is not None:
440 print("Note: Precomputed chif not used by this fit.")
441 mf = bbh_final_mass_non_precessing_UIB2016(m1, m2, chi1proj, chi2proj, version="v1")
442 else:
443 raise ValueError("Unrecognized fit name.")
444
445 return mf
446
447def bbh_final_spin_projected_spins(m1, m2, chi1, chi2, tilt1, tilt2, fitname, truncate=bbh_Kerr_trunc_opts.trunc):
448 """
449 Calculate the (signed) dimensionless spin parameter of the final BH resulting from the merger of two black holes,
450 only using the projected spins along the total angular momentum
451 and some aligned-spin fit from the literature
452
453 Parameters
454 ----------
455 m1, m2 : component masses
456 chi1, chi2 : dimensionless spins of two BHs
457 tilt1, tilt2 : tilts (in radians) in the new spin convention
458 fitname: fit selection currently supports Pan2011 (non-spinning), HLZ2014, PhenomD, UIB2016, HBR2016, HL2016
459
460 Returns
461 -------
462 (signed) dimensionless final spin parameter, chif
463 """
464
465 m1 = np.array(m1,dtype=np.float64)
466 m2 = np.array(m2,dtype=np.float64)
467 chi1 = np.array(chi1,dtype=np.float64)
468 chi2 = np.array(chi2,dtype=np.float64)
469 tilt1 = np.array(tilt1,dtype=np.float64)
470 tilt2 = np.array(tilt2,dtype=np.float64)
471
472 _check_mchi(m1,m2,chi1,chi2) # Check that inputs are physical
473
474 chi1proj = chi1*np.cos(tilt1)
475 chi2proj = chi2*np.cos(tilt2)
476
477 if fitname==bbh_final_state_fits.pan2011:
478 if np.any(chi1!=0) or np.any(chi2!=0) or np.any(tilt1!=0) or np.any(tilt2!=0):
479 print("Note: Pan2011 fit does not use spins.")
481 elif fitname==bbh_final_state_fits.hlz2014:
482 chif = bbh_final_spin_non_precessing_Healyetal(m1, m2, chi1proj, chi2proj, version="2014")
483 elif fitname==bbh_final_state_fits.hl2016:
484 chif = bbh_final_spin_non_precessing_Healyetal(m1, m2, chi1proj, chi2proj, version="2016")
485 elif fitname==bbh_final_state_fits.phenD:
486 chif = bbh_final_spin_non_precessing_Husaetal(m1, m2, chi1proj, chi2proj)
487 elif fitname==bbh_final_state_fits.uib2016:
488 chif = bbh_final_spin_non_precessing_UIB2016(m1, m2, chi1proj, chi2proj, version="v2")
489 elif fitname==bbh_final_state_fits.uib2016v1:
490 chif = bbh_final_spin_non_precessing_UIB2016(m1, m2, chi1proj, chi2proj, version="v1")
491 elif fitname==bbh_final_state_fits.hbr2016:
492 chif = bbh_final_spin_non_precessing_HBR2016(m1, m2, chi1proj, chi2proj)
493 else:
494 raise ValueError("Unrecognized fit name.")
495
496 chif = _truncate_at_Kerr_limit(chif, truncate, fitname) # optionally truncate and/or warn at Kerr limit, depending on user option
497
498 return chif
499
500def bbh_final_spin_precessing(m1, m2, chi1, chi2, tilt1, tilt2, phi12, fitname, truncate=bbh_Kerr_trunc_opts.trunc):
501 """
502 Calculate the magnitude of the dimensionless spin parameter of the final BH resulting from the merger of two black holes,
503 including the in-plane spin components;
504 by either using a precessing fit from the literature;
505 or by first projecting the spins along the angular momentum and using an aligned-spin fit from the literature,
506 then adding in quadrature the in-plane dimensionful spin scaled by the initial mass squared.
507
508 Parameters
509 ----------
510 m1, m2 : component masses
511 chi1, chi2 : dimensionless spins of two BHs
512 tilt1, tilt2 : tilts (in radians) in the new spin convention
513 phi12: angle (in radians) between in-plane spin components
514 fitname: fit selection currently supports Pan2011 (non-spinning), HLZ2014 (aligned+augmentation),
515 PhenomD (aligned+augmentation), UIB2016 (aligned+augmentation),
516 HL2016 (aligned+augmentation), HBR2016 (precessing)
517
518 Returns
519 -------
520 magnitude of the dimensionless final spin parameter, chif
521 """
522
523 m1 = np.array(m1,dtype=np.float64)
524 m2 = np.array(m2,dtype=np.float64)
525 chi1 = np.array(chi1,dtype=np.float64)
526 chi2 = np.array(chi2,dtype=np.float64)
527 tilt1 = np.array(tilt1,dtype=np.float64)
528 tilt2 = np.array(tilt2,dtype=np.float64)
529 phi12 = np.array(phi12,dtype=np.float64)
530
531 _check_mchi(m1,m2,chi1,chi2) # Check that inputs are physical
532
533 if fitname==bbh_final_state_fits.pan2011 or fitname==bbh_final_state_fits.hlz2014 or fitname==bbh_final_state_fits.phenD or fitname==bbh_final_state_fits.uib2016 or fitname==bbh_final_state_fits.uib2016v1 or fitname==bbh_final_state_fits.hl2016:
534 precfit = False
535 elif fitname==bbh_final_state_fits.hbr2016:
536 precfit = True
537 chif = bbh_final_spin_precessing_HBR2016(m1, m2, chi1, chi2, tilt1, tilt2, phi12)
538 else:
539 raise ValueError("Unrecognized fit name.")
540
541 if not precfit:
542 # First compute the final mass and parallel component of the final spin using the aligned components of the initial spins
543 chifaligned = bbh_final_spin_projected_spins(m1, m2, chi1, chi2, tilt1, tilt2, fitname, truncate)
544
545 # Now compute the squared magnitude of the in-plane dimensionful spin, first computing the magnitudes of the initial in-plane spins
546 S1perpmag = m1*m1*chi1*np.sin(tilt1)
547 S2perpmag = m2*m2*chi2*np.sin(tilt2)
548 Sperpmag2 = S1perpmag*S1perpmag + S2perpmag*S2perpmag + 2.*S1perpmag*S2perpmag*np.cos(phi12)
549
550 # Combine together
551 chif = (chifaligned*chifaligned + Sperpmag2/(m1+m2)**4.)**0.5
552
553 chif = _truncate_at_Kerr_limit(chif, truncate, "augmented "+fitname) # optionally truncate and/or warn at Kerr limit, depending on user option
554
555 return chif
556
557def bbh_final_mass_non_precessing_Husaetal(m1, m2, chi1, chi2):
558 """
559 Calculate the mass of the final BH resulting from the
560 merger of two black holes with non-precessing spins using the fits
561 used by IMRPhenomD, given in Eqs. (3.6) and (3.8) of Husa et al.
562 arXiv:1508.07250. Note that Eq. (3.8) gives the radiated energy, not
563 the final mass directly
564
565 m1, m2: component masses
566 chi1, chi2: dimensionless spins of two BHs
567 """
568 m1 = np.array(m1,dtype=np.float64)
569 m2 = np.array(m2,dtype=np.float64)
570 chi1 = np.array(chi1,dtype=np.float64)
571 chi2 = np.array(chi2,dtype=np.float64)
572
573 if np.any(abs(chi1)>1):
574 raise ValueError("chi1 has to be in [-1, 1]")
575 if np.any(abs(chi2)>1):
576 raise ValueError("chi2 has to be in [-1, 1]")
577
578 # binary parameters
579 m = m1+m2
580 msq = m*m
581
582 eta = m1*m2/msq
583 eta2 = eta*eta
584 eta3 = eta2*eta
585 eta4 = eta3*eta
586
587 S1 = chi1*m1**2/msq # spin angular momentum 1 (in m = 1 units)
588 S2 = chi2*m2**2/msq # spin angular momentum 2 (in m = 1 units)
589 S = S1+S2 # total spin
590 Sh = S/(1. - 2.*eta) # rescaled total spin
591
592 # Expressions copied from LALSimIMRPhenomD_internals.c (except with two notation differences: S is capitalized in chif and s -> Sh in Mf, in addition to the "m*(1. - ...)" to obtain the final mass from the radiated mass in m = 1 units which is calculated in the LAL code)
593
594 Mf = m*(1. - ((0.055974469826360077*eta + 0.5809510763115132*eta2 - 0.9606726679372312*eta3 + 3.352411249771192*eta4)*
595 (1. + (-0.0030302335878845507 - 2.0066110851351073*eta + 7.7050567802399215*eta2)*Sh))/(1. + (-0.6714403054720589 - 1.4756929437702908*eta + 7.304676214885011*eta2)*Sh))
596
597 return Mf
598
599def bbh_final_spin_non_precessing_Husaetal(m1, m2, chi1, chi2):
600 """
601 Calculate the spin of the final BH resulting from the
602 merger of two black holes with non-precessing spins using the fits
603 used by IMRPhenomD, given in Eqs. (3.6) and (3.8) of Husa et al.
604 arXiv:1508.07250. Note that Eq. (3.8) gives the radiated energy, not
605 the final mass directly
606
607 m1, m2: component masses
608 chi1, chi2: dimensionless spins of two BHs
609 """
610 # Vectorize the function if arrays are provided as input
611 m1 = np.array(m1,dtype=np.float64)
612 m2 = np.array(m2,dtype=np.float64)
613 chi1 = np.array(chi1,dtype=np.float64)
614 chi2 = np.array(chi2,dtype=np.float64)
615
616 if np.any(abs(chi1)>1):
617 raise ValueError("chi1 has to be in [-1, 1]")
618 if np.any(abs(chi2)>1):
619 raise ValueError("chi2 has to be in [-1, 1]")
620
621 # binary parameters
622 m = m1+m2
623 msq = m*m
624
625 eta = m1*m2/msq
626 eta2 = eta*eta
627 eta3 = eta2*eta
628 eta4 = eta3*eta
629
630 S1 = chi1*m1**2/msq # spin angular momentum 1 (in m = 1 units)
631 S2 = chi2*m2**2/msq # spin angular momentum 2 (in m = 1 units)
632 S = S1+S2 # total spin
633 Sh = S/(1. - 2.*eta) # rescaled total spin
634
635 Ssq = S*S
636 Scu = Ssq*S
637 Squ = Scu*S
638
639 # Expressions copied from LALSimIMRPhenomD_internals.c (except with two notation differences: S is capitalized in chif and s -> Sh in Mf, in addition to the "m*(1. - ...)" to obtain the final mass from the radiated mass in m = 1 units which is calculated in the LAL code)
640
641 chif = 3.4641016151377544*eta - 4.399247300629289*eta2 + 9.397292189321194*eta3 - 13.180949901606242*eta4 + (1 - 0.0850917821418767*eta - 5.837029316602263*eta2)*S + (0.1014665242971878*eta - 2.0967746996832157*eta2)*Ssq + (-1.3546806617824356*eta + 4.108962025369336*eta2)*Scu + (-0.8676969352555539*eta + 2.064046835273906*eta2)*Squ
642
643 return chif
644
645def bbh_UIBfits_setup(m1, m2, chi1, chi2):
646 """
647 common setup function for UIB final-state and luminosity fit functions
648 """
649
650 # Vectorize the function if arrays are provided as input
651 m1 = np.array(m1,dtype=np.float64)
652 m2 = np.array(m2,dtype=np.float64)
653 chi1 = np.array(chi1,dtype=np.float64)
654 chi2 = np.array(chi2,dtype=np.float64)
655
656 if np.any(m1<0):
657 raise ValueError("m1 must not be negative")
658 if np.any(m2<0):
659 raise ValueError("m2 must not be negative")
660
661 if np.any(abs(chi1)>1):
662 raise ValueError("chi1 has to be in [-1, 1]")
663 if np.any(abs(chi2)>1):
664 raise ValueError("chi2 has to be in [-1, 1]")
665
666 # binary masses
667 m = m1+m2
668 if np.any(m<=0):
669 raise ValueError("m1+m2 must be positive")
670 msq = m*m
671 m1sq = m1*m1
672 m2sq = m2*m2
673
674 # symmetric mass ratio
675 eta = m1*m2/msq
676 if np.any(eta>0.25):
677 print("Truncating eta from above to 0.25. This should only be necessary in some rounding corner cases, but better check your m1 and m2 inputs...")
678 eta = np.minimum(eta,0.25)
679 if np.any(eta<0.0):
680 print("Truncating negative eta to 0.0. This should only be necessary in some rounding corner cases, but better check your m1 and m2 inputs...")
681 eta = np.maximum(eta,0.0)
682 eta2 = eta*eta
683 eta3 = eta2*eta
684 eta4 = eta2*eta2
685
686 # spin variables (in m = 1 units)
687 S1 = chi1*m1sq/msq # spin angular momentum 1
688 S2 = chi2*m2sq/msq # spin angular momentum 2
689 Stot = S1+S2 # total spin
690 Shat = (chi1*m1sq+chi2*m2sq)/(m1sq+m2sq) # effective spin, = msq*Stot/(m1sq+m2sq)
691 Shat2 = Shat*Shat
692 Shat3 = Shat2*Shat
693 Shat4 = Shat2*Shat2
694
695 # asymmetric spin combination (spin difference), where the paper assumes m1>m2
696 # to make our implementation symmetric under simultaneous exchange of m1<->m2 and chi1<->chi2,
697 # we flip the sign here when m2>m1
698 chidiff = chi1 - chi2
699 if np.any(m2>m1):
700 chidiff = np.sign(m1-m2)*chidiff
701 chidiff2 = chidiff*chidiff
702
703 # typical squareroots and functions of eta
704 sqrt2 = 2.**0.5
705 sqrt3 = 3.**0.5
706 sqrt1m4eta = (1. - 4.*eta)**0.5
707
708 return m, eta, eta2, eta3, eta4, Stot, Shat, Shat2, Shat3, Shat4, chidiff, chidiff2, sqrt2, sqrt3, sqrt1m4eta
709
710def bbh_final_mass_non_precessing_UIB2016(m1, m2, chi1, chi2, version="v2"):
711 """
712 Calculate the final mass with the aligned-spin NR fit
713 by Xisco Jimenez Forteza, David Keitel, Sascha Husa et al.
714 [LIGO-P1600270] [https://arxiv.org/abs/1611.00332]
715 versions v1 and v2 use the same ansatz,
716 with v2 calibrated to additional SXS and RIT data
717
718 m1, m2: component masses
719 chi1, chi2: dimensionless spins of two BHs
720 Results are symmetric under simultaneous exchange of m1<->m2 and chi1<->chi2.
721 """
722
723 m, eta, eta2, eta3, eta4, Stot, Shat, Shat2, Shat3, Shat4, chidiff, chidiff2, sqrt2, sqrt3, sqrt1m4eta = bbh_UIBfits_setup(m1, m2, chi1, chi2)
724
725 if version == "v1":
726 # rational-function Pade coefficients (exact) from Eq. (22) of 1611.00332v1
727 b10 = 0.487
728 b20 = 0.295
729 b30 = 0.17
730 b50 = -0.0717
731 # fit coefficients from Tables VII-X of 1611.00332v1
732 # values at increased numerical precision copied from
733 # https://gravity.astro.cf.ac.uk/cgit/cardiff_uib_share/tree/Luminosity_and_Radiated_Energy/UIBfits/LALInference/EradUIB2016_pyform_coeffs.txt
734 # git commit 636e5a71462ecc448060926890aa7811948d5a53
735 a2 = 0.5635376058169299
736 a3 = -0.8661680065959881
737 a4 = 3.181941595301782
738 b1 = -0.15800074104558132
739 b2 = -0.15815904609933157
740 b3 = -0.14299315232521553
741 b5 = 8.908772171776285
742 f20 = 3.8071100104582234
743 f30 = 25.99956516423936
744 f50 = 1.552929335555098
745 f10 = 1.7004558922558886
746 f21 = 0.
747 d10 = -0.12282040108157262
748 d11 = -3.499874245551208
749 d20 = 0.014200035799803777
750 d30 = -0.01873720734635449
751 d31 = -5.1830734185518725
752 f11 = 14.39323998088354
753 f31 = -232.25752840151296
754 f51 = -0.8427987782523847
755
756 elif version == "v2":
757 # rational-function Pade coefficients (exact) from Eq. (22) of LIGO-P1600270-v4
758 b10 = 0.346
759 b20 = 0.211
760 b30 = 0.128
761 b50 = -0.212
762 # fit coefficients from Tables VII-X of LIGO-P1600270-v4
763 # values at increased numerical precision copied from
764 # https://dcc.ligo.org/DocDB/0128/P1600270/004/FinalStateUIB2016_suppl_Erad_coeffs.txt
765 a2 = 0.5609904135313374
766 a3 = -0.84667563764404
767 a4 = 3.145145224278187
768 b1 = -0.2091189048177395
769 b2 = -0.19709136361080587
770 b3 = -0.1588185739358418
771 b5 = 2.9852925538232014
772 f20 = 4.271313308472851
773 f30 = 31.08987570280556
774 f50 = 1.5673498395263061
775 f10 = 1.8083565298668276
776 f21 = 0.
777 d10 = -0.09803730445895877
778 d11 = -3.2283713377939134
779 d20 = 0.01118530335431078
780 d30 = -0.01978238971523653
781 d31 = -4.91667749015812
782 f11 = 15.738082204419655
783 f31 = -243.6299258830685
784 f51 = -0.5808669012986468
785
786 else:
787 raise ValueError('Unknown version -- should be either "v1" or "v2".')
788
789 # Calculate the radiated-energy fit from Eq. (28) of LIGO-P1600270-v4
790 Erad = (((1. + -2.0/3.0*sqrt2)*eta + a2*eta2 + a3*eta3 + a4*eta4)*(1. + b10*b1*Shat*(f10 + f11*eta + (16. - 16.*f10 - 4.*f11)*eta2) + b20*b2*Shat2*(f20 + f21*eta + (16. - 16.*f20 - 4.*f21)*eta2) + b30*b3*Shat3*(f30 + f31*eta + (16. - 16.*f30 - 4.*f31)*eta2)))/(1. + b50*b5*Shat*(f50 + f51*eta + (16. - 16.*f50 - 4.*f51)*eta2)) + d10*sqrt1m4eta*eta2*(1. + d11*eta)*chidiff + d30*Shat*sqrt1m4eta*eta*(1. + d31*eta)*chidiff + d20*eta3*chidiff2
791
792 # Convert to actual final mass
793 Mf = m*(1.-Erad)
794
795 return Mf
796
797def bbh_final_spin_non_precessing_UIB2016(m1, m2, chi1, chi2, version="v2"):
798 """
799 Calculate the final spin with the aligned-spin NR fit
800 by Xisco Jimenez Forteza, David Keitel, Sascha Husa et al.
801 [LIGO-P1600270] [https://arxiv.org/abs/1611.00332]
802 versions v1 and v2 use the same ansatz,
803 with v2 calibrated to additional SXS and RIT data
804
805 m1, m2: component masses
806 chi1, chi2: dimensionless spins of two BHs
807 Results are symmetric under simultaneous exchange of m1<->m2 and chi1<->chi2.
808 """
809
810 m, eta, eta2, eta3, eta4, Stot, Shat, Shat2, Shat3, Shat4, chidiff, chidiff2, sqrt2, sqrt3, sqrt1m4eta = bbh_UIBfits_setup(m1, m2, chi1, chi2)
811
812 if version == "v1":
813 # rational-function Pade coefficients (exact) from Eqs. (7) and (8) of 1611.00332v1
814 a20 = 5.28
815 a30 = 1.27
816 a50 = 2.89
817 b10 = -0.194
818 b20 = 0.075
819 b30 = 0.00782
820 b50 = -0.527
821 # fit coefficients from Tables I-IV of 1611.00332v1
822 # values at increased numerical precision copied from
823 # https://gravity.astro.cf.ac.uk/cgit/cardiff_uib_share/tree/Luminosity_and_Radiated_Energy/UIBfits/LALInference/FinalSpinUIB2016_pyform_coeffs.txt
824 # git commit 636e5a71462ecc448060926890aa7811948d5a53
825 a2 = 3.772362507208651
826 a3 = -9.627812453422376
827 a5 = 2.487406038123681
828 b1 = 1.0005294518146604
829 b2 = 0.8823439288807416
830 b3 = 0.7612809461506448
831 b5 = 0.9139185906568779
832 f21 = 8.887933111404559
833 f31 = 23.927104476660883
834 f50 = 1.8981657997557002
835 f11 = 4.411041530972546
836 f52 = 0.
837 d10 = 0.2762804043166152
838 d11 = 11.56198469592321
839 d20 = -0.05975750218477118
840 d30 = 2.7296903488918436
841 d31 = -3.388285154747212
842 f12 = 0.3642180211450878
843 f22 = -40.35359764942015
844 f32 = -178.7813942566548
845 f51 = -5.556957394513334
846
847 elif version == "v2":
848 # rational-function Pade coefficients (exact) from Eqs. (7) and (8) of LIGO-P1600270-v4
849 a20 = 5.24
850 a30 = 1.3
851 a50 = 2.88
852 b10 = -0.194
853 b20 = 0.0851
854 b30 = 0.00954
855 b50 = -0.579
856 # fit coefficients from Tables I-IV of LIGO-P1600270-v4
857 # values at increased numerical precision copied from
858 # https://dcc.ligo.org/DocDB/0128/P1600270/004/FinalStateUIB2016_suppl_spin_coeffs.txt
859 a2 = 3.8326341618708577
860 a3 = -9.487364155598392
861 a5 = 2.5134875145648374
862 b1 = 1.0009563702914628
863 b2 = 0.7877509372255369
864 b3 = 0.6540138407185817
865 b5 = 0.8396665722805308
866 f21 = 8.77367320110712
867 f31 = 22.830033250479833
868 f50 = 1.8804718791591157
869 f11 = 4.409160174224525
870 f52 = 0.
871 d10 = 0.3223660562764661
872 d11 = 9.332575956437443
873 d20 = -0.059808322561702126
874 d30 = 2.3170397514509933
875 d31 = -3.2624649875884852
876 f12 = 0.5118334706832706
877 f22 = -32.060648277652994
878 f32 = -153.83722669033995
879 f51 = -4.770246856212403
880
881 else:
882 raise ValueError('Unknown version -- should be either "v1" or "v2".')
883
884 # Calculate the fit for the Lorb' quantity from Eq. (16) of LIGO-P1600270-v4
885 Lorb = (2.*sqrt3*eta + a20*a2*eta2 + a30*a3*eta3)/(1. + a50*a5*eta) + (b10*b1*Shat*(f11*eta + f12*eta2 + (64. - 16.*f11 - 4.*f12)*eta3) + b20*b2*Shat2*(f21*eta + f22*eta2 + (64. - 16.*f21 - 4.*f22)*eta3) + b30*b3*Shat3*(f31*eta + f32*eta2 + (64. - 16.*f31 - 4.*f32)*eta3))/(1. + b50*b5*Shat*(f50 + f51*eta + f52*eta2 + (64. - 64.*f50 - 16.*f51 - 4.*f52)*eta3)) + d10*sqrt1m4eta*eta2*(1. + d11*eta)*chidiff + d30*Shat*sqrt1m4eta*eta3*(1. + d31*eta)*chidiff + d20*eta3*chidiff2
886
887 # Convert to actual final spin
888 chif = Lorb + Stot
889
890 return chif
891
892def _bbh_HBR2016_setup(m1, m2, chi1, chi2):
893 """
894 Setup function for the Hofmann, Barausse, and Rezzolla final spin fits to vectorize the masses and spins and calculate the mass ratio.
895 """
896
897 # Vectorize if arrays are provided as input
898 m1 = np.array(m1,dtype=np.float64)
899 m2 = np.array(m2,dtype=np.float64)
900 chi1 = np.array(chi1,dtype=np.float64)
901 chi2 = np.array(chi2,dtype=np.float64)
902
903 return m1, m2, chi1, chi2, m2/m1
904
905def _bbh_HBR2016_ell(m1, m2, chi1z, chi2z, version):
906 """
907 Compute the orbital angular momentum ell used by the Hofmann, Barausse, and Rezzolla final spin fit [from ApJL 825, L19 (2016)], henceforth HBR.
908 The three versions available correspond to the three choices of fit coefficients given in Table 1 of that paper, with 6, 16, and 20 coefficients, respectively.
909
910 version can thus be "M1J2", "M3J3", or "M3J4"
911 """
912
913 # Set upper bounds for sums and coefficients from Table 1 in HBR; k00 is calculated from Eq. (11) in the paper
914
915 if version == "M1J2":
916 nM = 1
917 nJ = 2
918
919 k00 = -3.82116
920 k01 = -1.2019
921 k02 = -1.20764
922 k10 = 3.79245
923 k11 = 1.18385
924 k12 = 4.90494
925 xi = 0.41616
926
927 elif version == "M3J3":
928 nM = 3
929 nJ = 3
930
931 k00 = -5.83164
932 k01 = 2.87025
933 k02 = -1.53315
934 k03 = -3.78893
935 k10 = 32.9127
936 k11 = -62.9901
937 k12 = 10.0068
938 k13 = 56.1926
939 k20 = -136.832
940 k21 = 329.32
941 k22 = -13.2034
942 k23 = -252.27
943 k30 = 210.075
944 k31 = -545.35
945 k32 = -3.97509
946 k33 = 368.405
947 xi = 0.463926
948
949 elif version == "M3J4":
950 nM = 3
951 nJ = 4
952
953 k00 = -5.97723
954 k01 = 3.39221
955 k02 = 4.48865
956 k03 = -5.77101
957 k04 = -13.0459
958 k10 = 35.1278
959 k11 = -72.9336
960 k12 = -86.0036
961 k13 = 93.7371
962 k14 = 200.975
963 k20 = -146.822
964 k21 = 387.184
965 k22 = 447.009
966 k23 = -467.383
967 k24 = -884.339
968 k30 = 223.911
969 k31 = -648.502
970 k32 = -697.177
971 k33 = 753.738
972 k34 = 1166.89
973 xi = 0.474046
974
975 else:
976 raise ValueError('Unknown version--should be either "M1J2", "M3J3", or "M3J4".')
977
978 # Calculate eta, atot, and aeff; note that HBR call the symmetric mass ratio nu instead of eta
979
980 m = m1 + m2
981 q = m2/m1
982 eta = m1*m2/(m*m)
983
984 atot = (chi1z + q*q*chi2z)/((1.+q)*(1.+q)) # Eq. (12) in HBR
985 aeff = atot + xi*eta*(chi1z + chi2z) # Inline equation below Eq. (7); see also Eq. (15) for the precessing analogue
986
987 # Calculate ISCO energy and angular momentum
988 r_isco = calc_isco_radius(aeff)
989 e_isco = (1. - 2./3./r_isco)**0.5 # Eq. (2) in HBR
990 l_isco = 2./(3.*3.**0.5)*(1. + 2.*(3.*r_isco - 2.)**0.5) # Eq. (3) in HBR
991
992 # The following expressions give the double sum in Eq. (13) in HBR (without the overall factor of eta that is put in below) specialized to the nJ values for the three versions, i.e., nJ = 2 => M1J2, nJ = 3 => M3J3, nJ = 4 => M3J4
993 if nJ >= 2:
994 aeff2 = aeff*aeff
995 ksum = k00 + k01*aeff + k02*aeff2 + eta*(k10 + k11*aeff + k12*aeff2)
996 if nJ >= 3:
997 eta2 = eta*eta
998 eta3 = eta2*eta
999 aeff3 = aeff2*aeff
1000 ksum = ksum + (k03 + eta*k13)*aeff3 + eta2*(k20 + k21*aeff + k22*aeff2 + k23*aeff3) + eta3*(k30 + k31*aeff + k32*aeff2 + k33*aeff3)
1001 if nJ >= 4:
1002 ksum = ksum + (k04 + eta*k14 + eta2*k24 + eta3*k34)*aeff3*aeff
1003
1004 # Calculate the absolute value of ell
1005 ell = abs((l_isco - 2.*atot*(e_isco - 1.)) + eta*ksum) # Eq. (13) in HBR
1006
1007 return ell
1008
1009def bbh_final_spin_non_precessing_HBR2016(m1, m2, chi1z, chi2z, version="M3J3"):
1010 """
1011 Calculate the (signed) dimensionless spin of the final BH resulting from the
1012 merger of two black holes with aligned spins using the fit from Hofmann, Barausse, and Rezzolla ApJL 825, L19 (2016), henceforth HBR.
1013
1014 The three versions available correspond to the three choices of fit coefficients given in Table 1 of that paper, with 6, 16, and 20 coefficients, respectively.
1015
1016 version can thus be "M1J2", "M3J3", or "M3J4"
1017
1018 m1, m2: component masses
1019 chi1z, chi2z: components of the dimensionless spins of the two BHs along the orbital angular momentum
1020 """
1021
1022 # Calculate q and vectorize the masses and spins if arrays are provided as input
1023
1024 m1, m2, chi1z, chi2z, q = _bbh_HBR2016_setup(m1, m2, chi1z, chi2z)
1025
1026 # Calculate the final spin
1027
1028 atot = (chi1z + chi2z*q*q)/((1.+q)*(1.+q)) # Eq. (12) in HBR
1029
1030 ell = _bbh_HBR2016_ell(m1, m2, chi1z, chi2z, version)
1031
1032 return atot + ell/(1./q + 2. + q) # Eq. (12) in HBR, writing the symmetric mass ratio in terms of q
1033
1034def bbh_final_spin_precessing_HBR2016(m1, m2, chi1, chi2, tilt1, tilt2, phi12, version="M3J3"):
1035 """
1036 Calculate the dimensionless spin of the final BH resulting from the
1037 merger of two black holes with precessing spins using the fit from Hofmann, Barausse, and Rezzolla ApJL 825, L19 (2016), henceforth HBR.
1038
1039 The three versions available correspond to the three choices of fit coefficients given in Table 1 of that paper, with 6, 16, and 20 coefficients, respectively.
1040
1041 version can thus be "M1J2", "M3J3", or "M3J4"
1042
1043 m1, m2: component masses
1044 chi1, chi2: dimensionless spins of two BHs
1045 tilt1, tilt2: tilt angles of the spins from the orbital angular momentum
1046 phi12: angle between in-plane spin components
1047 """
1048
1049 # Calculate q and vectorize the masses and spins if arrays are provided as input
1050
1051 m1, m2, chi1, chi2, q = _bbh_HBR2016_setup(m1, m2, chi1, chi2)
1052
1053 # Vectorize the spin angles if arrays are provided as input
1054 tilt1 = np.array(tilt1,dtype=np.float64)
1055 tilt2 = np.array(tilt2,dtype=np.float64)
1056 phi12 = np.array(phi12,dtype=np.float64)
1057
1058 # Set eps (\epsilon_\beta or \epsilon_\gamma) to the value given below Eq. (18) in HBR
1059
1060 eps = 0.024
1061
1062 # Computing angles defined in Eq. (17) of HBR. The betas and gammas expressions are for the starred quantities computed using the second (approximate) equality in Eq. (18) in HBR
1063 cos_beta = np.cos(tilt1)
1064 cos_betas = np.cos(tilt1 + eps*np.sin(tilt1))
1065 cos_gamma = np.cos(tilt2)
1066 cos_gammas = np.cos(tilt2 + eps*np.sin(tilt2))
1067 cos_alpha = ((1 - cos_beta*cos_beta)*(1 - cos_gamma*cos_gamma))**0.5*np.cos(phi12) + cos_beta*cos_gamma # This rewrites the inner product definition of cos_alpha in terms of cos_beta, cos_gamma, and phi12
1068
1069 # Define a shorthand and compute the final spin
1070
1071 q2 = q*q
1072
1073 ell = _bbh_HBR2016_ell(m1, m2, chi1*cos_betas, chi2*cos_gammas, version)
1074
1075 # Compute the final spin value [Eq. (16) in HBR], truncating the argument of the square root at zero if it becomes negative
1076 sqrt_arg = chi1*chi1 + chi2*chi2*q2*q2 + 2.*chi1*chi2*q2*cos_alpha + 2.*(chi1*cos_betas + chi2*q2*cos_gammas)*ell*q + ell*ell*q2
1077 if np.any(sqrt_arg < 0.):
1078 print(f"bbh_final_spin_precessing_HBR2016(): The argument of the square root is negative for indexes {np.where(sqrt_arg < 0)[0]}; truncating it to zero.")
1079 sqrt_arg = np.clip(sqrt_arg,0.,None)
1080
1081 # Return the final spin value [Eq. (16) in HBR]
1082 return sqrt_arg**0.5/((1.+q)*(1.+q))
1083
1084#######################
1085# Peak luminosity fits
1086#######################
1087
1088# list of peak luminosity fit names
1090 t1600018 = "T1600018" # Jimenez Forteza et al. [LIGO-T1600018 (2016)] (aligned)
1091 uib2016 = "UIB2016" # Keitel et al. [arXiv:1612.09566] (aligned)
1092 hl2016 = "HL2016" # Healy and Lousto [arXiv:1610.09713] (aligned)
1093
1094def _rescale_Lpeak(Lpeak):
1095 """
1096 convert from geometric units ("Planck luminosity" of c^5/G) to multiples of 10^56 erg/s = 10^49 J/s
1097 """
1098 LumPl_ergs_per_sec = lal.LUMPL_SI*1e-49 # Approximate value = 3628.505
1099 return LumPl_ergs_per_sec*Lpeak
1100
1101def bbh_aligned_Lpeak_6mode_SHXJDK(q, chi1, chi2):
1102 """
1103 wrapper to bbh_peak_luminosity_non_precessing_T1600018() for backwards compatibility
1104
1105 q: mass ratio (here m2/m1, where m1>m2)
1106 chi1: the component of the dimensionless spin of m1 along the angular momentum (z)
1107 chi2: the component of the dimensionless spin of m2 along the angular momentum (z)
1108 """
1109 # Vectorize the function if arrays are provided as input
1110 q = np.array(q,dtype=np.float64)
1111
1112 # from bayespputils.py convention is q = m2/m1, where m1>m2.
1113 if np.any(q<=0.):
1114 raise ValueError("q has to be > 0.")
1115 if np.any(q>1.):
1116 raise ValueError("q has to be <= 1.")
1117
1118 # T1600018 fit expects m1>m2;
1119 # implementation here should be symmetric,
1120 # but we follow that convention to be sure
1121 m1 = 1./(1.+q)
1122 m2 = q/(1.+q)
1123
1124 return bbh_peak_luminosity_non_precessing_T1600018(m1, m2, chi1, chi2)
1125
1126def bbh_peak_luminosity_non_precessing_T1600018(m1, m2, chi1, chi2):
1127 """
1128 Calculate the peak luminosity (using modes 22, 21, 33, 32, 44, and 43) of a binary black hole with aligned spins using the fit made by Sascha Husa, Xisco Jimenez Forteza, David Keitel [LIGO-T1500598] using 5th order in chieff and return results in units of 10^56 ergs/s
1129
1130 m1, m2: component masses
1131 chi1: the component of the dimensionless spin of m1 along the angular momentum (z)
1132 chi2: the component of the dimensionless spin of m2 along the angular momentum (z)
1133 Results are symmetric under simultaneous exchange of m1<->m2 and chi1<->chi2.
1134 """
1135
1136 # Vectorize the function if arrays are provided as input
1137 m1 = np.array(m1,dtype=np.float64)
1138 m2 = np.array(m2,dtype=np.float64)
1139 chi1 = np.array(chi1,dtype=np.float64)
1140 chi2 = np.array(chi2,dtype=np.float64)
1141
1142 # Calculate powers of eta and the effective spin S (not used in this fit)
1143 m, eta, eta2, eta3, eta4, Stot, Shat, Shat2, Shat3, Shat4, chidiff, chidiff2, sqrt2, sqrt3, sqrt1m4eta = bbh_UIBfits_setup(m1, m2, chi1, chi2)
1144
1145 dm2 = 1. - 4.*eta # equivalent to sqrt1m4eta**2
1146
1147 chi_eff = (m1*chi1+m2*chi2)/m # equivalent to (q_inv*chi1 + chi2)/(1. + q_inv)
1148 chi_eff2 = chi_eff*chi_eff
1149 chi_eff3 = chi_eff2*chi_eff
1150 chi_eff4 = chi_eff3*chi_eff
1151 chi_eff5 = chi_eff4*chi_eff
1152
1153 # Calculate best fit (from [https://dcc.ligo.org/T1500598-v4])
1154 Lpeak = (0.012851338846828302 + 0.007822265919928252*chi_eff + 0.010221856361035788*chi_eff2 + 0.015805535732661396*chi_eff3 + 0.0011356206806770043*chi_eff4 - 0.009868152529667197*chi_eff5)*eta2 + (0.05681786589129071 - 0.0017473702709303457*chi_eff - 0.10150706091341818*chi_eff2 - 0.2349153289253309*chi_eff3 + 0.015657737820040145*chi_eff4 + 0.19556893194885075*chi_eff5)*eta4 + 0.026161288241420833*dm2**0.541825641769908*eta**3.1629576945611757*chidiff + 0.0007771032100485481*dm2**0.4499151697918658*eta**1.7800346166040835*chidiff2
1155
1156 # Convert to 10^56 ergs/s units
1157 return _rescale_Lpeak(Lpeak)
1158
1159def bbh_peak_luminosity_non_precessing_UIB2016(m1, m2, chi1, chi2):
1160 """
1161 Calculate the peak luminosity with the aligned-spin NR fit
1162 by David Keitel, Xisco Jimenez Forteza, Sascha Husa, Lionel London et al.
1163 [LIGO-P1600279-v5] [https://arxiv.org/abs/1612.09566v1]
1164 using modes up to lmax=6,
1165 and return results in units of 10^56 ergs/s
1166
1167 m1, m2: component masses
1168 chi1, chi2: dimensionless spins of two BHs
1169 Results are symmetric under simultaneous exchange of m1<->m2 and chi1<->chi2.
1170 """
1171
1172 m, eta, eta2, eta3, eta4, Stot, Shat, Shat2, Shat3, Shat4, chidiff, chidiff2, sqrt2, sqrt3, sqrt1m4eta = bbh_UIBfits_setup(m1, m2, chi1, chi2)
1173 eta5 = eta3*eta2
1174
1175 # fit coefficients corresponding to Table I, II, IV,
1176 # exact values corresponding to https://arxiv.org/src/1612.09566v1/anc/LpeakUIB2016_suppl_coeffs.txt
1177 # fi2 coefficients are replaced by functions of fi0, fi1 as in Eq. (10)
1178 a0 = 0.8742169580717333
1179 a1 = -2.111792574893241
1180 a2 = 35.214103272783646
1181 a3 = -244.94930678226913
1182 a4 = 877.1061892200927
1183 a5 = -1172.549896493467
1184 b1 = 0.9800204548606681
1185 b2 = -0.1779843936224084
1186 b4 = 1.7859209418791981
1187 d10 = 3.789116271213293
1188 d20 = 0.40214125006660567
1189 d30 = 4.273116678713487
1190 f10 = 1.6281049269810424
1191 f11 = -3.6322940180721037
1192 f20 = 31.710537408279116
1193 f21 = -273.84758785648336
1194 f30 = -0.23470852321351202
1195 f31 = 6.961626779884965
1196 f40 = 0.21139341988062182
1197 f41 = 1.5255885529750841
1198 f60 = 3.0901740789623453
1199 f61 = -16.66465705511997
1200 f70 = 0.8362061463375388
1201 f71 = 0.
1202
1203 # calculate the Lpeak/(eta2*L0) fit from Eq. (14), using the constraints for fi2 from Eq. (10)
1204 Lpeak = a0 + a1*eta + a2*eta2 + a3*eta3 + a4*eta4 + a5*eta5 + (0.465*b1*Shat*(f10 + f11*eta + (16. - 16.*f10 - 4.*f11)*eta2) + 0.107*b2*Shat2*(f20 + f21*eta + (16. - 16.*f20 - 4.*f21)*eta2) + Shat3*(f30 + f31*eta + (-16.*f30 - 4.*f31)*eta2) + Shat4*(f40 + f41*eta + (-16.*f40 - 4.*f41)*eta2))/(1. - 0.328*b4*Shat*(f60 + f61*eta + (16. - 16.*f60 - 4.*f61)*eta2) + Shat2*(f70 + f71*eta + (-16.*f70 - 4.*f71)*eta2)) + eta3*((d10+d30*Shat)*sqrt1m4eta*chidiff + d20*chidiff2)
1205
1206 # Lpeak(eta=0.25,chi1=chi2=0)/0.25^2
1207 L0 = 0.016379197203103536
1208
1209 # Convert to actual peak luminosity
1210 Lpeak = Lpeak*eta2*L0
1211
1212 # Convert to 10^56 ergs/s units
1213 return _rescale_Lpeak(Lpeak)
1214
1215def bbh_peak_luminosity_non_precessing_Healyetal(m1, m2, chi1z, chi2z):
1216 """
1217 Calculate the peak luminosity of an aligned-spin binary black hole coalescence using the fit from Healy and Lousto arXiv:1610.09713 (henceforth HL)
1218
1219 Parameters
1220 ----------
1221 m1, m2 : component masses
1222 chi1z, chi2z : components of the dimensionless spins of the two BHs along their orbital angular momentum
1223
1224 Returns
1225 -------
1226 peak luminosity in units of 10^56 ergs/s
1227
1228 """
1229
1230 m1 = np.array(m1,dtype=np.float64)
1231 m2 = np.array(m2,dtype=np.float64)
1232 chi1z = np.array(chi1z,dtype=np.float64)
1233 chi2z = np.array(chi2z,dtype=np.float64)
1234
1235 eta, delta_m, S, Delta = _RIT_setup(m1, m2, chi1z, chi2z)
1236
1237 # fitting coefficients (from the RIT group's private Matlab implementation; versions to fewer digits are given in Table IV of HL)
1238 coeffs = [0.00102101737,0.000897428902,-9.77467189e-05,0.000920883841,1.86982704e-05,-0.000391316975,-0.000120214144,0.000148102239,0.00137901473,-0.000493730555,0.000884792724,3.29254224e-07,1.70170729e-05,0.00151484008,-0.000148456828,0.000136659593,0.000160343115,-6.18530577e-05,-0.00103602521]
1239
1240 # Return the peak luminosity in units of 10^56 ergs/s [cf. Eq. (3) in HL]
1241 Lpeak = _RIT_symm_express(eta, delta_m, S, Delta, coeffs)
1242
1243 return _rescale_Lpeak(Lpeak)
1244
1245def bbh_peak_luminosity_projected_spins(m1, m2, chi1, chi2, tilt1, tilt2, fitname):
1246 """
1247 Calculate the peak luminosity of the merger of two black holes,
1248 only using the projected spins along the total angular momentum
1249 and some aligned-spin fit from the literature
1250
1251 Parameters
1252 ----------
1253 m1, m2 : component masses
1254 chi1, chi2 : dimensionless spins of two BHs
1255 tilt1, tilt2 : tilts (in radians) in the new spin convention
1256 fitname: fit selection currently supports T1600018, UIB2016, HL2016
1257
1258 Returns
1259 -------
1260 peak luminosity, Lpeak, in units of 10^56 ergs/s
1261 """
1262
1263 m1 = np.array(m1,dtype=np.float64)
1264 m2 = np.array(m2,dtype=np.float64)
1265 chi1 = np.array(chi1,dtype=np.float64)
1266 chi2 = np.array(chi2,dtype=np.float64)
1267 tilt1 = np.array(tilt1,dtype=np.float64)
1268 tilt2 = np.array(tilt2,dtype=np.float64)
1269
1270 _check_mchi(m1,m2,chi1,chi2) # Check that inputs are physical
1271
1272 chi1proj = chi1*np.cos(tilt1)
1273 chi2proj = chi2*np.cos(tilt2)
1274
1275 if fitname==bbh_Lpeak_fits.t1600018:
1276 Lpeak = bbh_peak_luminosity_non_precessing_T1600018(m1, m2, chi1proj, chi2proj)
1277 elif fitname==bbh_Lpeak_fits.uib2016:
1278 Lpeak = bbh_peak_luminosity_non_precessing_UIB2016(m1, m2, chi1proj, chi2proj)
1279 elif fitname==bbh_Lpeak_fits.hl2016:
1280 Lpeak = bbh_peak_luminosity_non_precessing_Healyetal(m1, m2, chi1proj, chi2proj)
1281 else:
1282 raise ValueError("Unrecognized fit name.")
1283
1284 return Lpeak
1285
1286#########################
1287# Convenience functions
1288#########################
1289
1290def bbh_average_fits_precessing(m1, m2, chi1, chi2, tilt1, tilt2, phi12, quantity, fits):
1291 """
1292 Calculate the average of the predictions of various fits, currently just for the final mass, final spin, or peak luminosity of a quasicircular
1293 binary black hole coalescence. The final spin calculation returns the magnitude of the final spin, including the contribution
1294 from in-plane spins; the other two cases just apply aligned-spin fits to the components of the spins along the orbital angular momentum.
1295
1296 Parameters
1297 ----------
1298 m1, m2 : component masses
1299 chi1, chi2 : dimensionless spins of two BHs
1300 tilt1, tilt2 : tilts (in radians) in the new spin convention
1301 phi12: angle (in radians) between in-plane spin components (only used for the final spin)
1302 quantity: "Mf", "af", "afz", or "Lpeak"
1303 fits: An array of fit names to be used. The possible fit names are those known by bbh_final_mass_projected_spins, bbh_final_spin_precessing,
1304 and bbh_peak_luminosity_projected_spins
1305
1306 The shape of m1 is used to determine the shape of the output.
1307
1308 Returns
1309 -------
1310 Average of the results for the given fits for the chosen quantity
1311 """
1312
1313 if quantity != "af" and max(abs(np.atleast_1d(phi12))) != 0:
1314 print("Note: phi12 is only used for the full final spin calculation.")
1315
1316 if quantity not in ["Mf", "af", "afz", "Lpeak"]:
1317 raise ValueError("Unknown quantity: %s"%quantity)
1318
1319 # Define function to return the appropriate quantity
1320
1321 def _return_quantity(m1, m2, chi1, chi2, tilt1, tilt2, phi12, fitname):
1322 if quantity == "Mf":
1323 return bbh_final_mass_projected_spins(m1, m2, chi1, chi2, tilt1, tilt2, fitname)
1324 elif quantity == "af":
1325 return bbh_final_spin_precessing(m1, m2, chi1, chi2, tilt1, tilt2, phi12, fitname)
1326 elif quantity == "afz":
1327 return bbh_final_spin_projected_spins(m1, m2, chi1, chi2, tilt1, tilt2, fitname)
1328 elif quantity == "Lpeak":
1329 return bbh_peak_luminosity_projected_spins(m1, m2, chi1, chi2, tilt1, tilt2, fitname)
1330
1331 # Define a function that returns the length of an array when passed an array and 1 when not passed an array
1332 def _len_smart(x):
1333 if hasattr(x, '__len__'):
1334 return len(x)
1335 else:
1336 return 1
1337
1338 # Initialize
1339
1340 fits = np.atleast_1d(fits)
1341
1342 if not fits.size: # Check to make sure that fits is not an empty array
1343 raise ValueError("The list of fits passed cannot be an empty array.")
1344
1345 num_fits = len(fits)
1346 num_data = max(_len_smart(m1), _len_smart(m2), _len_smart(chi1), _len_smart(chi2), _len_smart(tilt1), _len_smart(tilt2), _len_smart(phi12))
1347
1348 data = np.zeros([num_fits, num_data])
1349
1350 # Convert the input from bayespputils into a format that is compatible with the following loop
1351
1352 if hasattr(m1, 'shape'):
1353 data_shape = m1.shape
1354 else:
1355 data_shape = -1
1356
1357 # Loop over the fits
1358
1359 for k, fit in enumerate(fits):
1360 data_portion = _return_quantity(m1, m2, chi1, chi2, tilt1, tilt2, phi12, fit)
1361 data[k] = data_portion.reshape(-1)
1362
1363 # Calculate average:
1364 data_avg = np.mean(data,axis=0)
1365
1366 return data_avg.reshape(data_shape)
#define max(a, b)
def bbh_final_mass_non_precessing_Healyetal(m1, m2, chi1, chi2, version="2014", chif=None)
Calculate the mass of the final BH resulting from the merger of two black holes with non-precessing s...
Definition: nrutils.py:293
def calc_isco_radius(a)
Calculate the ISCO radius of a Kerr BH as a function of the Kerr parameter using eqns.
Definition: nrutils.py:147
def bbh_peak_luminosity_non_precessing_T1600018(m1, m2, chi1, chi2)
Calculate the peak luminosity (using modes 22, 21, 33, 32, 44, and 43) of a binary black hole with al...
Definition: nrutils.py:1134
def bbh_final_mass_non_precessing_UIB2016(m1, m2, chi1, chi2, version="v2")
Calculate the final mass with the aligned-spin NR fit by Xisco Jimenez Forteza, David Keitel,...
Definition: nrutils.py:721
def bbh_final_spin_projected_spins(m1, m2, chi1, chi2, tilt1, tilt2, fitname, truncate=bbh_Kerr_trunc_opts.trunc)
Calculate the (signed) dimensionless spin parameter of the final BH resulting from the merger of two ...
Definition: nrutils.py:463
def bbh_final_spin_precessing(m1, m2, chi1, chi2, tilt1, tilt2, phi12, fitname, truncate=bbh_Kerr_trunc_opts.trunc)
Calculate the magnitude of the dimensionless spin parameter of the final BH resulting from the merger...
Definition: nrutils.py:521
def bbh_final_mass_projected_spins(m1, m2, chi1, chi2, tilt1, tilt2, fitname, chif=None)
Calculate the mass of the final BH resulting from the merger of two black holes, only using the proje...
Definition: nrutils.py:404
def bbh_final_mass_non_spinning_Panetal(m1, m2)
Calculate the mass of the final BH resulting from the merger of two non-spinning black holes using eq...
Definition: nrutils.py:110
def bbh_aligned_Lpeak_6mode_SHXJDK(q, chi1, chi2)
wrapper to bbh_peak_luminosity_non_precessing_T1600018() for backwards compatibility
Definition: nrutils.py:1108
def bbh_peak_luminosity_non_precessing_UIB2016(m1, m2, chi1, chi2)
Calculate the peak luminosity with the aligned-spin NR fit by David Keitel, Xisco Jimenez Forteza,...
Definition: nrutils.py:1170
def bbh_final_mass_non_precessing_Husaetal(m1, m2, chi1, chi2)
Calculate the mass of the final BH resulting from the merger of two black holes with non-precessing s...
Definition: nrutils.py:567
def bbh_final_spin_projected_spin_Healyetal(m1, m2, chi1, chi2, tilt1, tilt2)
wrapper to bbh_final_spin_projected_spins() for backwards compatibility
Definition: nrutils.py:372
def bbh_final_spin_non_spinning_Panetal(m1, m2)
Calculate the spin of the final BH resulting from the merger of two non-spinning black holes using eq...
Definition: nrutils.py:129
def bbh_final_spin_non_precessing_Husaetal(m1, m2, chi1, chi2)
Calculate the spin of the final BH resulting from the merger of two black holes with non-precessing s...
Definition: nrutils.py:609
def bbh_final_spin_precessing_HBR2016(m1, m2, chi1, chi2, tilt1, tilt2, phi12, version="M3J3")
Calculate the dimensionless spin of the final BH resulting from the merger of two black holes with pr...
Definition: nrutils.py:1047
def bbh_final_spin_non_precessing_HBR2016(m1, m2, chi1z, chi2z, version="M3J3")
Calculate the (signed) dimensionless spin of the final BH resulting from the merger of two black hole...
Definition: nrutils.py:1020
def bbh_final_spin_precessing_Healyetal_extension_Minit(m1, m2, chi1, chi2, tilt1, tilt2, phi12)
wrapper to bbh_final_spin_precessing() for backwards compatibility
Definition: nrutils.py:384
def bbh_peak_luminosity_non_precessing_Healyetal(m1, m2, chi1z, chi2z)
Calculate the peak luminosity of an aligned-spin binary black hole coalescence using the fit from Hea...
Definition: nrutils.py:1228
def bbh_peak_luminosity_projected_spins(m1, m2, chi1, chi2, tilt1, tilt2, fitname)
Calculate the peak luminosity of the merger of two black holes, only using the projected spins along ...
Definition: nrutils.py:1261
def bbh_final_spin_non_precessing_Healyetal(m1, m2, chi1, chi2, version="2014")
Calculate the spin of the final BH resulting from the merger of two black holes with non-precessing s...
Definition: nrutils.py:269
def bbh_final_mass_projected_spin_Healyetal(m1, m2, chi1, chi2, tilt1, tilt2, chif=None)
wrapper to bbh_final_mass_projected_spins() for backwards compatibility
Definition: nrutils.py:378
def bbh_average_fits_precessing(m1, m2, chi1, chi2, tilt1, tilt2, phi12, quantity, fits)
Convenience functions.
Definition: nrutils.py:1311
def bbh_UIBfits_setup(m1, m2, chi1, chi2)
common setup function for UIB final-state and luminosity fit functions
Definition: nrutils.py:648
def bbh_final_spin_non_precessing_UIB2016(m1, m2, chi1, chi2, version="v2")
Calculate the final spin with the aligned-spin NR fit by Xisco Jimenez Forteza, David Keitel,...
Definition: nrutils.py:808