2Various fitting formulas provided by numerical relativity
4Archisman Ghosh, Nathan K. Johnson-McDaniel, P. Ajith, 2015-04-09
8import scipy.optimize
as so
12 print(
'Cannot import lal SWIG bindings')
18 raise ValueError(
"m1 has to be positive")
20 raise ValueError(
"m2 has to be positive")
22def _check_chi(chi1, chi2):
24 raise ValueError(
"chi1 has to be <= 1")
26 raise ValueError(
"chi2 has to be <= 1")
28 raise ValueError(
"chi1 has to be nonnegative")
30 raise ValueError(
"chi2 has to be nonnegative")
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")
38def _check_mchi(m1,m2,chi1,chi2):
42def _check_mchi_aligned(m1,m2,chi1,chi2):
44 _check_chi_aligned(chi1,chi2)
55 uib2016v1 =
"UIB2016v1"
63 trunc_silent =
"TRUNCATESILENT"
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)))
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))
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))
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).
105 m1, m2 : component masses
111 m1 = np.array(m1,dtype=np.float64)
112 m2 = np.array(m2,dtype=np.float64)
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)))
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).
124 m1, m2 : component masses
128 final dimensionless spin, chif
130 m1 = np.array(m1,dtype=np.float64)
131 m2 = np.array(m2,dtype=np.float64)
133 eta = m1*m2/(m1+m2)**2.
134 return np.sqrt(12.)*eta - 3.871*(eta**2.) + 4.028*(eta**3)
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)
149 a = np.minimum(np.array(a),1.)
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)
155 return 3+z2 - np.sqrt((3.-z1)*(3.+z1+2.*z2))*a_sign
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 """
161 _check_mchi_aligned(m1, m2, chi1, chi2)
165 delta_m = (m1 - m2)/m
167 S = (m1*m1*chi1 + m2*m2*chi2)/m**2
168 Delta = (m2*chi2 - m1*chi1)/m
170 return eta, delta_m, S, Delta
172def _RIT_symm_express(eta, delta_m, S, Delta, coeffs):
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
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).
181 dm2 = delta_m*delta_m
183 Deltadm = Delta*delta_m
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)
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 """
200 J_isco = (3*np.sqrt(r_isco)-2*a_f)*2./np.sqrt(3*r_isco)
203 if version ==
"2014":
223 elif version ==
"2016":
244 raise ValueError(
'Unknown version--should be either "2014" or "2016".')
246 dm2 = delta_m*delta_m
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
253 return abs(a_f-a_f_new)
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")
263 m1, m2 : component masses
264 chi1, chi2 : dimensionless spins of two BHs
268 dimensionless final spin, chif
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
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")
286 m1, m2 : component masses
287 chi1, chi2 : dimensionless spins of two BHs
288 chif: final spin (optional),
if already calculated
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)
299 eta, delta_m, S, Delta = _RIT_setup(m1, m2, chi1, chi2)
304 chif = np.array(chif)
312 if version ==
"2014":
332 elif version ==
"2016":
353 raise ValueError(
'Unknown version--should be either "2014" or "2016".')
356 E_isco = (1. - 2./r_isco + chif/r_isco**1.5)/np.sqrt(1. - 3./r_isco + 2.*chif/r_isco**1.5)
358 dm2 = delta_m*delta_m
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
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
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
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)
413 chif = np.array(chif,dtype=np.float64)
415 _check_mchi(m1,m2,chi1,chi2)
417 chi1proj = chi1*np.cos(tilt1)
418 chi2proj = chi2*np.cos(tilt2)
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.")
424 print(
"Note: Precomputed chif not used by this fit.")
426 elif fitname==bbh_final_state_fits.hlz2014:
428 elif fitname==bbh_final_state_fits.hl2016:
430 elif fitname==bbh_final_state_fits.phenD:
432 print(
"Note: Precomputed chif not used by this fit.")
434 elif fitname==bbh_final_state_fits.uib2016:
436 print(
"Note: Precomputed chif not used by this fit.")
438 elif fitname==bbh_final_state_fits.uib2016v1:
440 print(
"Note: Precomputed chif not used by this fit.")
443 raise ValueError(
"Unrecognized fit name.")
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
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
462 (signed) dimensionless final spin parameter, chif
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)
472 _check_mchi(m1,m2,chi1,chi2)
474 chi1proj = chi1*np.cos(tilt1)
475 chi2proj = chi2*np.cos(tilt2)
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:
483 elif fitname==bbh_final_state_fits.hl2016:
485 elif fitname==bbh_final_state_fits.phenD:
487 elif fitname==bbh_final_state_fits.uib2016:
489 elif fitname==bbh_final_state_fits.uib2016v1:
491 elif fitname==bbh_final_state_fits.hbr2016:
494 raise ValueError(
"Unrecognized fit name.")
496 chif = _truncate_at_Kerr_limit(chif, truncate, fitname)
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.
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)
520 magnitude of the dimensionless final spin parameter, chif
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)
531 _check_mchi(m1,m2,chi1,chi2)
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:
535 elif fitname==bbh_final_state_fits.hbr2016:
539 raise ValueError(
"Unrecognized fit name.")
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)
551 chif = (chifaligned*chifaligned + Sperpmag2/(m1+m2)**4.)**0.5
553 chif = _truncate_at_Kerr_limit(chif, truncate,
"augmented "+fitname)
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
565 m1, m2: component masses
566 chi1, chi2: dimensionless spins of two BHs
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)
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]")
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))
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
607 m1, m2: component masses
608 chi1, chi2: dimensionless spins of two BHs
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)
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]")
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
647 common setup function for UIB final-state
and luminosity fit functions
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)
657 raise ValueError(
"m1 must not be negative")
659 raise ValueError(
"m2 must not be negative")
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]")
669 raise ValueError(
"m1+m2 must be positive")
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)
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)
690 Shat = (chi1*m1sq+chi2*m2sq)/(m1sq+m2sq)
698 chidiff = chi1 - chi2
700 chidiff = np.sign(m1-m2)*chidiff
701 chidiff2 = chidiff*chidiff
706 sqrt1m4eta = (1. - 4.*eta)**0.5
708 return m, eta, eta2, eta3, eta4, Stot, Shat, Shat2, Shat3, Shat4, chidiff, chidiff2, sqrt2, sqrt3, sqrt1m4eta
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
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.
723 m, eta, eta2, eta3, eta4, Stot, Shat, Shat2, Shat3, Shat4, chidiff, chidiff2, sqrt2, sqrt3, sqrt1m4eta = bbh_UIBfits_setup(m1, m2, chi1, chi2)
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
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
756 elif version ==
"v2":
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
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
787 raise ValueError(
'Unknown version -- should be either "v1" or "v2".')
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
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
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.
810 m, eta, eta2, eta3, eta4, Stot, Shat, Shat2, Shat3, Shat4, chidiff, chidiff2, sqrt2, sqrt3, sqrt1m4eta = bbh_UIBfits_setup(m1, m2, chi1, chi2)
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
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
847 elif version ==
"v2":
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
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
882 raise ValueError(
'Unknown version -- should be either "v1" or "v2".')
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
892def _bbh_HBR2016_setup(m1, m2, chi1, chi2):
894 Setup function for the Hofmann, Barausse,
and Rezzolla final spin fits to vectorize the masses
and spins
and calculate the mass ratio.
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)
903 return m1, m2, chi1, chi2, m2/m1
905def _bbh_HBR2016_ell(m1, m2, chi1z, chi2z, version):
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.
910 version can thus be
"M1J2",
"M3J3",
or "M3J4"
915 if version ==
"M1J2":
927 elif version ==
"M3J3":
949 elif version ==
"M3J4":
976 raise ValueError(
'Unknown version--should be either "M1J2", "M3J3", or "M3J4".')
984 atot = (chi1z + q*q*chi2z)/((1.+q)*(1.+q))
985 aeff = atot + xi*eta*(chi1z + chi2z)
989 e_isco = (1. - 2./3./r_isco)**0.5
990 l_isco = 2./(3.*3.**0.5)*(1. + 2.*(3.*r_isco - 2.)**0.5)
995 ksum = k00 + k01*aeff + k02*aeff2 + eta*(k10 + k11*aeff + k12*aeff2)
1000 ksum = ksum + (k03 + eta*k13)*aeff3 + eta2*(k20 + k21*aeff + k22*aeff2 + k23*aeff3) + eta3*(k30 + k31*aeff + k32*aeff2 + k33*aeff3)
1002 ksum = ksum + (k04 + eta*k14 + eta2*k24 + eta3*k34)*aeff3*aeff
1005 ell = abs((l_isco - 2.*atot*(e_isco - 1.)) + eta*ksum)
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.
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.
1016 version can thus be
"M1J2",
"M3J3",
or "M3J4"
1018 m1, m2: component masses
1019 chi1z, chi2z: components of the dimensionless spins of the two BHs along the orbital angular momentum
1024 m1, m2, chi1z, chi2z, q = _bbh_HBR2016_setup(m1, m2, chi1z, chi2z)
1028 atot = (chi1z + chi2z*q*q)/((1.+q)*(1.+q))
1030 ell = _bbh_HBR2016_ell(m1, m2, chi1z, chi2z, version)
1032 return atot + ell/(1./q + 2. + q)
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.
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.
1041 version can thus be
"M1J2",
"M3J3",
or "M3J4"
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
1051 m1, m2, chi1, chi2, q = _bbh_HBR2016_setup(m1, m2, chi1, chi2)
1054 tilt1 = np.array(tilt1,dtype=np.float64)
1055 tilt2 = np.array(tilt2,dtype=np.float64)
1056 phi12 = np.array(phi12,dtype=np.float64)
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
1073 ell = _bbh_HBR2016_ell(m1, m2, chi1*cos_betas, chi2*cos_gammas, version)
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)
1082 return sqrt_arg**0.5/((1.+q)*(1.+q))
1090 t1600018 =
"T1600018"
1094def _rescale_Lpeak(Lpeak):
1096 convert from geometric units (
"Planck luminosity" of c^5/G) to multiples of 10^56 erg/s = 10^49 J/s
1098 LumPl_ergs_per_sec = lal.LUMPL_SI*1e-49
1099 return LumPl_ergs_per_sec*Lpeak
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)
1110 q = np.array(q,dtype=np.float64)
1114 raise ValueError(
"q has to be > 0.")
1116 raise ValueError(
"q has to be <= 1.")
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
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.
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)
1143 m, eta, eta2, eta3, eta4, Stot, Shat, Shat2, Shat3, Shat4, chidiff, chidiff2, sqrt2, sqrt3, sqrt1m4eta =
bbh_UIBfits_setup(m1, m2, chi1, chi2)
1147 chi_eff = (m1*chi1+m2*chi2)/m
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
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
1157 return _rescale_Lpeak(Lpeak)
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
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.
1172 m, eta, eta2, eta3, eta4, Stot, Shat, Shat2, Shat3, Shat4, chidiff, chidiff2, sqrt2, sqrt3, sqrt1m4eta = bbh_UIBfits_setup(m1, m2, chi1, chi2)
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
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)
1207 L0 = 0.016379197203103536
1210 Lpeak = Lpeak*eta2*L0
1213 return _rescale_Lpeak(Lpeak)
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)
1221 m1, m2 : component masses
1222 chi1z, chi2z : components of the dimensionless spins of the two BHs along their orbital angular momentum
1226 peak luminosity
in units of 10^56 ergs/s
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)
1235 eta, delta_m, S, Delta = _RIT_setup(m1, m2, chi1z, chi2z)
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]
1241 Lpeak = _RIT_symm_express(eta, delta_m, S, Delta, coeffs)
1243 return _rescale_Lpeak(Lpeak)
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
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
1260 peak luminosity, Lpeak,
in units of 10^56 ergs/s
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)
1270 _check_mchi(m1,m2,chi1,chi2)
1272 chi1proj = chi1*np.cos(tilt1)
1273 chi2proj = chi2*np.cos(tilt2)
1275 if fitname==bbh_Lpeak_fits.t1600018:
1277 elif fitname==bbh_Lpeak_fits.uib2016:
1279 elif fitname==bbh_Lpeak_fits.hl2016:
1282 raise ValueError(
"Unrecognized fit name.")
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.
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
1306 The shape of m1
is used to determine the shape of the output.
1310 Average of the results
for the given fits
for the chosen quantity
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.")
1316 if quantity
not in [
"Mf",
"af",
"afz",
"Lpeak"]:
1317 raise ValueError(
"Unknown quantity: %s"%quantity)
1321 def _return_quantity(m1, m2, chi1, chi2, tilt1, tilt2, phi12, fitname):
1322 if quantity ==
"Mf":
1324 elif quantity ==
"af":
1326 elif quantity ==
"afz":
1328 elif quantity ==
"Lpeak":
1333 if hasattr(x,
'__len__'):
1340 fits = np.atleast_1d(fits)
1343 raise ValueError(
"The list of fits passed cannot be an empty array.")
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))
1348 data = np.zeros([num_fits, num_data])
1352 if hasattr(m1,
'shape'):
1353 data_shape = m1.shape
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)
1364 data_avg = np.mean(data,axis=0)
1366 return data_avg.reshape(data_shape)
Final mass and spin fits.
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...
def calc_isco_radius(a)
Calculate the ISCO radius of a Kerr BH as a function of the Kerr parameter using eqns.
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...
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,...
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 ...
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...
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...
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...
def bbh_aligned_Lpeak_6mode_SHXJDK(q, chi1, chi2)
wrapper to bbh_peak_luminosity_non_precessing_T1600018() for backwards compatibility
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,...
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...
def bbh_final_spin_projected_spin_Healyetal(m1, m2, chi1, chi2, tilt1, tilt2)
wrapper to bbh_final_spin_projected_spins() for backwards compatibility
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...
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...
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...
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...
def bbh_final_spin_precessing_Healyetal_extension_Minit(m1, m2, chi1, chi2, tilt1, tilt2, phi12)
wrapper to bbh_final_spin_precessing() for backwards compatibility
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...
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 ...
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...
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
def bbh_average_fits_precessing(m1, m2, chi1, chi2, tilt1, tilt2, phi12, quantity, fits)
Convenience functions.
def bbh_UIBfits_setup(m1, m2, chi1, chi2)
common setup function for UIB final-state and luminosity fit functions
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,...