2Implement the computation of the tilt angles at infinity or bounds and average values of tilts at a finite separation
3from tilt angles at some frequency, taken to be sufficiently small so that precession-averaged evolution is accurate.
4This implementation is described in the paper, <https://dcc.ligo.org/P2100029>, arXiv:2107.11902.
6N. K. Johnson-McDaniel, 2021
11from numpy.linalg
import norm
13from numpy.polynomial
import Polynomial
as P
15from scipy.special
import ellipe, ellipk
17from scipy.integrate
import ode
19from .tilts_at_infinity_utils
import *
21from warnings
import warn
24 from scipy.integrate
import solve_ivp
26 solve_ivp_avail =
False
28 solve_ivp_avail =
True
39LPN_orders = [0, 1, 1.5, 2, 2.5]
41LPN_order_string =
', '.join(map(str, LPN_orders))
47 no_mass_spin_input_check=
False):
49 Compute the PN orbital angular momentum from the binary
's masses and spins and m = 2 GW frequency
50 using the expression in Eq. (4) of <https://dcc.ligo.org/public/0122/T1500554/023/dLdS.pdf>.
52 We orbit-average the contributions
from the spins
as in Eq. (17) of that reference, but the spins are
53 set to zero by default.
59 m1, m2: Detector frame masses
in kg
60 f: m = 2 GW frequency
in Hz
64 PNorder: Order of PN corrections to include--choices are 0, 1, 1.5, 2,
or 2.5, where the half-integer orders are
65 the spin-orbit contributions, default: 2
66 S1_vec, S2_vec: 3d numpy arrays giving the components of the dimensionful spins (
with total mass = 1)
in a coordinate
67 system where the Newtonian orbital angular momentum
is in the z-direction, default: [0., 0., 0.]
68 no_mass_spin_input_check: Turn off mass
and spin input checking
for when this
is called by other functions that have
69 already performed this check, default:
False
71 Output: 3d array giving the PN orbital angular momentum
in total mass = 1 units
in a coordinate system where the
72 Newtonian orbital angular momentum
is in the z-direction
77 if not no_mass_spin_input_check:
82 "The frequency must be positive, while it is f = %f" % f)
84 if PNorder
not in LPN_orders:
85 raise ValueError(
"The PN order to use to compute the orbital angular momentum must be one of %s, while it is "
86 "PNorder = %f" % (LPN_order_string, PNorder))
94 if not no_mass_spin_input_check
and (S1 > m1*m1/Msq
or S2 > m2*m2/Msq):
96 "The magnitudes of the spins must be at most (m1/M)^2 = %f and (m2/M)^2 = %f, respectively, "
97 "where M is the total mass, while they are |S1|, |S2| = %f, %f" %
98 (m1*m1/Msq, m2*m2/Msq, S1, S2))
110 zhat = np.array([0., 0., 1.])
112 S1z_vec = S1_vec[2]*zhat
113 S2z_vec = S2_vec[2]*zhat
115 v = (np.pi*M*kg_to_s*f)**(1./3.)
127 L += (1.5+eta/6.)*v2*zhat
129 L += (-0.25*((1. + 3./X1)*S1_vec + (1. + 3./X2)*S2_vec)
130 - ((1. + 27./X1)*S1z_vec + (1. + 27./X2)*S2z_vec)/12.
133 L += (3.375 - 2.375*eta + eta*eta/24.)*v2*v2*zhat
135 L += ((-71. + 21./X1 - 13.*X1 - 9.*X1sq)*S1_vec + (-71. + 21./X2 - 13.*X2 - 9.*X2sq)*S2_vec
136 + (1. - 147./X1 - 367.*X1/3. + 13.*X1sq/3.)*S1z_vec
137 + (1. - 147./X2 - 367.*X2/3. + 13.*X2sq/3.)*S2z_vec
143def L_from_masses_spins_and_f(m1, m2, f, PNorder=2, S1_vec=np.zeros(3), S2_vec=np.zeros(3), no_mass_spin_input_check=
False):
145 Compute the magnitude of the PN orbital angular momentum from the binary
's masses and spins and m = 2 GW frequency
146 using the expression in Eq. (4) of <https://dcc.ligo.org/public/0122/T1500554/023/dLdS.pdf>.
148 We orbit-average the contributions
from the spins
as in Eq. (17) of that reference, but the spins are
149 set to zero by default.
157 m1, m2: Detector frame masses
in kg
158 f: m = 2 GW frequency
in Hz
162 PNorder: Order of PN corrections to include--choices are 0, 1, 1.5, 2,
or 2.5, where the half-integer orders are
163 the spin-orbit contributions, default: 2
164 S1_vec, S2_vec: 3d numpy arrays giving the components of the dimensionful spins (
with total mass = 1)
in a coordinate
165 system where the Newtonian orbital angular momentum
is in the z-direction, default: [0., 0., 0.]
166 no_mass_spin_input_check: Turn off mass
and spin input checking
for when this
is called by other functions that have
167 already performed this check, default:
False
169 Output: Magnitude of the PN orbital angular momentum
in total mass = 1 units
173 no_mass_spin_input_check=no_mass_spin_input_check))
178 Compute the magnitude of the Newtonian orbital angular momentum from the binary
's masses, semimajor axis, and
183 m1, m2: Source frame masses in kg
184 a: Semimajor axis
in m
187 Output: Magnitude of orbital angular momentum
in total mass = 1 units
196 "The semimajor axis must be positive, while it is a = %e" % a)
200 "The eccentricity must be between 0 and 1, while it is e = %f" % e)
208 return eta*(a*(1. - e*e)/(M*kg_to_m))**0.5
213 Compute m*E_ratio(m) [cf. Eq. (B13) in the paper <https://dcc.ligo.org/P2100029>, arXiv:2107.11902, which gives m^2*E_ratio(m)]
219 return (3.*(1. + 3.*(4. - m)/(omm*fpm))/(16.*omm**1.5) - 0.5)*m/fpm
222def eq_coeffs(u, kappaxiq, q, S1, S2, xi, S0sq):
224 Compute coefficients for the cubic equation
in bar{S}^2 to be solved [Eq. (20)
in the paper <https://dcc.ligo.org/P2100029>, arXiv:2107.11902]
227 u: 1/(2L), where L
is the orbital angular momentum (i.e., the integration variable)
228 kappaxiq: The rescaled variable defined
in Eq. (14) of the paper
229 q: Mass ratio of the binary,
with the < 1 convention
230 S1, S2: Magnitudes of the dimensionful spins of the binary (
with total mass = 1)
231 xi: Conserved effective spin, defined
in Eq. (2) of the paper
232 S0sq: Square of magnitude of the initial total spin (at u = u0,
and with total mass = 1)
234 Output: The coefficients of the equation (vareps, Btransf, Ctransf, Dtransf)
250 kappaxiq2 = kappaxiq*kappaxiq
258 Sigma = 0.5*(S0sq - S1sq - S2sq)
259 Upsilon = opq*(2.*q*Sigma + q2*S1sq + S2sq)
260 zeta = q*(Sigma + q*S1sq)*xi
262 vareps = q*opq*eps*u2
264 Btransf = 0.25*opq*eps*eps + q*eps * \
265 (xi - 2.*opq*kappaxiq)*u + Upsilon*u2
267 Ctransf = eps*(opq*(Sigma + q*kappaxiq2) - q*kappaxiq *
268 xi) + 2.*(zeta - Upsilon*kappaxiq)*u
270 Dtransf = opq*(Sigma*Sigma - S1sq*S2sq) + q2*S1sq*xi * \
271 xi/opq - 2.*zeta*kappaxiq + Upsilon*kappaxiq2
273 return vareps, Btransf, Ctransf, Dtransf
276def solve_cubic(coeffs, u, kappaxiq, q, S1, S2, imag_tol, use_mpmath=False, root_tol=1.e-6, polyroots_maxsteps=50, polyroots_extraprec=50):
278 Solve the cubic equation, check for complex roots,
and order the roots.
279 The inputs besides the coefficients are just needed
in order to be printed
280 if one of the checks
is triggered
286 coeffs: Array of coefficients of cubic [vareps, Btransf, Ctransf, Dtransf]
287 u: 1/(2L), where L
is the orbital angular momentum (i.e., the integration variable)
288 kappaxiq: The rescaled variable defined
in Eq. (14) of the paper <https://dcc.ligo.org/P2100029>, arXiv:2107.11902
289 q: Mass ratio of the binary,
with the < 1 convention
290 S1, S2: Magnitudes of the dimensionful spins of the binary (
with total mass = 1)
291 imag_tol: Tolerance
for the imaginary part of roots of the cubic
295 use_mpmath: If
True, use mpmath.polyroots to solve the cubic instead of numpy
's Polynomial roots method
296 (slower but more accurate), default: False
297 root_tol: Tolerance
for the error estimate
from mpmath.polyroots,
if use_mpmath ==
True, default: 1.e-6
298 polyroots_maxsteps: The maximum number of steps allowed
in mpmath.polyroots
if use_mpmath ==
True, default: 50
299 polyroots_extraprec: The number of extra digits to use
in mpmath.polyroots
if use_mpmath ==
True, default: 50
301 Output: The roots of the equation ordered
from smallest to largest (Ssq3, Ssqm, Ssqp)
306 vareps, Btransf, Ctransf, Dtransf = coeffs
309 Btransf_out, Ctransf_out, Dtransf_out = Btransf, Ctransf, Dtransf
311 roots, err = mp.polyroots([vareps, Btransf, Ctransf, Dtransf], error=
True, maxsteps=polyroots_maxsteps, extraprec=polyroots_extraprec)
313 SsqI, SsqII, SsqIII = roots
317 warn(
"polyroots error of %e greater than root_tol = %e"%(err, root_tol), RuntimeWarning)
320 Btransf_out, Ctransf_out, Dtransf_out = Btransf[0], Ctransf[0], Dtransf[0]
322 cubic = P([Dtransf_out, Ctransf_out, Btransf_out, vareps])
324 SsqI, SsqII, SsqIII = cubic.roots()
327 if max(abs(SsqI.imag), abs(SsqII.imag), abs(SsqIII.imag)) > imag_tol:
328 raise RuntimeError(
"At least one of the roots of the cubic in S^2 has an imaginary part larger than the "
329 "tolerance of %.2g set in the code.\nThe roots are (%.10f, %.10f), (%.10f, %.10f), "
330 "(%.10f, %.10f), and the coefficients are %.10e, %.10e, %.10e, %.10e at u = %.10f, "
331 "kappaxiq = %.10f. q = %f, S1, S2 = %f, %f" %
332 (imag_tol, SsqI.real, SsqI.imag, SsqII.real, SsqII.imag, SsqIII.real, SsqIII.imag,
333 vareps, Btransf_out, Ctransf_out, Dtransf_out, u, kappaxiq, q, S1, S2))
337 Ssq3, Ssqm, Ssqp = np.sort([SsqI.real, SsqII.real, SsqIII.real])
341 if Ssq3.real == Ssqm.real:
342 if Ssq3.imag != 0.
or Ssqm.imag != 0.:
344 r"\bar{S}_3^2 and \bar{S}_-^2 have the same real parts, with values of: %e, %e"%(Ssq3, Ssqm), RuntimeWarning)
345 elif Ssq3 == 0.
and Ssqm == 0.:
346 raise NonprecessingError(
r"The system is nonprecessing: \bar{S}^2_3 = \bar(S}^2_- = 0")
348 if Ssqm.real == Ssqp.real:
349 if Ssqm.imag != 0.
or Ssqp.imag != 0.:
351 r"\bar{S}_-^2 and \bar{S}_+^2 have the same real parts, with values of: %e, %e"%(Ssqm, Ssqp), RuntimeWarning)
352 elif Ssqm == 0.
and Ssqp == 0.:
353 raise NonprecessingError(
r"The system is nonprecessing: \bar{S}^2_- = \bar(S}^2_+ = 0")
355 return Ssq3, Ssqm, Ssqp
358def dkappaxiq_du(u, kappaxiq, q, S1, S2, xi, S0sq, imag_tol, root_tol, lin_tol, root_tol_raise_err, use_mpmath=False,
359 polyroots_maxsteps=50, polyroots_extraprec=50):
361 Compute (<S^2>_prec - S0^2)/(1 - q) [= dkappa_{xi q}/du] using Eq. (42) in Chatziioannou, Klein, Yunes,
and Cornish PRD 95, 104004 (2017)
362 [rewritten
as bar{S}_+^2 + [E(m)/K(m) - 1](bar{S}_+^2 - bar{S}_3^2), where bar{S}_*^2 := (S_*^2 - S_0^2)/(1 - q)].
363 Here kappa_{xi q} := (kappa - S_0^2/(2 L) - q xi/(1 + q))/(1 - q) - xi/4,
and S_0
is the magnitude of the
364 initial total spin (i.e., at u = u0).
370 u: 1/(2L), where L
is the orbital angular momentum (i.e., the integration variable)
371 kappaxiq: The rescaled variable defined above
and in Eq. (14) of the paper <https://dcc.ligo.org/P2100029>, arXiv:2107.11902
372 q: Mass ratio of the binary,
with the < 1 convention
373 S1, S2: Magnitudes of the dimensionful spins of the binary (
with total mass = 1)
374 xi: Conserved effective spin, defined
in Eq. (2) of the paper <https://dcc.ligo.org/P2100029>, arXiv:2107.11902
375 S0sq: Square of magnitude of the initial total spin (at u = u0,
and with total mass = 1)
376 imag_tol: Tolerance
for the imaginary part of roots of the cubic
377 root_tol: Tolerance
for the error
in the roots of the cubic (
is an absolute tolerance
for small values of the roots,
as
378 generally occurs during the beginning of the evolution,
and is a fractional tolerance when the roots are large,
379 which generally occurs towards the end of the evolution; also
is reduced automatically when the roots are close
381 lin_tol: Tolerance
for errors on the r.h.s. of dkappa/du due to linearization
in m (
is a fractional tolerance when m
is
382 large,
as generally occurs during the beginning of the evolution,
and is an absolute tolerance when m
is small,
383 as generally occurs towards the end of the evolution)
384 root_tol_raise_err: Whether to
raise an error (
True)
or just
print a warning (
False) when the error
in the
385 roots of the cubic
is larger than root_tol
389 use_mpmath: If
True, use mpmath functions instead of numpy functions to solve the cubic
and evaluate the elliptic integrals
390 (slower but more accurate), default:
False
391 polyroots_maxsteps: The maximum number of steps allowed
in mpmath.polyroots
if use_mpmath ==
True, default: 50
392 polyroots_extraprec: The number of extra digits to use
in mpmath.polyroots
if use_mpmath ==
True, default: 50
394 Output: (<S^2>_prec - S0^2)/(1 - q) [= dkappa_{xi q}/du]
403 vareps, Btransf, Ctransf, Dtransf =
eq_coeffs(
404 u, kappaxiq, q, S1, S2, xi, S0sq)
416 only_quadratic =
False
419 if min(Btransf, Ctransf, Dtransf) > 0
and 4.*Ctransf*vareps < Btransf*Btransf:
421 (-Btransf + (Btransf*Btransf - 4.*Ctransf*vareps)**0.5)/(2.*Ctransf))
424 sum_err_bound = z_eps_abs * \
425 (abs(Ctransf*Ctransf - Btransf*Dtransf) + Ctransf*Dtransf*z_eps_abs) / \
426 (Btransf*(Btransf - Ctransf*z_eps_abs))
429 if z_eps_abs*(Ctransf/Btransf + sum_err_bound) < 1:
430 m_bound = z_eps_abs*(Ctransf/Btransf + sum_err_bound) / \
431 (1. - z_eps_abs*(Ctransf/Btransf + sum_err_bound))
438 if E_ratio_term(m_bound) <= 0.5*lin_tol*max(1., 1./(Ctransf/Btransf + sum_err_bound)) \
439 and sum_err_bound < lin_tol:
440 only_quadratic =
True
443 Ssq_pr = -0.5*Ctransf/Btransf
447 Ssq_pr_unbarred = (1. - q)*Ssq_pr + S0sq
449 if Ssq_pr_unbarred < 0.:
450 raise RuntimeError(
"The unbarred S_+^2 + S_-^2 value (which was all that was necessary to compute to this accuracy) is negative, "
451 "indicating a problem with the evolution, with a value of %e"%(2.*Ssq_pr_unbarred))
454 [vareps, Btransf, Ctransf, Dtransf], u, kappaxiq, q, S1, S2, imag_tol, use_mpmath=use_mpmath, root_tol=root_tol,
455 polyroots_maxsteps=polyroots_maxsteps, polyroots_extraprec=polyroots_extraprec)
459 Ssqp_unbarred = (1. - q)*Ssqp + S0sq
460 Ssqm_unbarred = (1. - q)*Ssqm + S0sq
462 if Ssqp_unbarred < 0.
or Ssqm_unbarred < 0.:
463 raise RuntimeError(
"At least one of the unbarred S_+^2 and S_-^2 values are negative, indicating a problem with the evolution, "
464 "with values of %e, %e"%(Ssqp_unbarred, Ssqm_unbarred))
470 "Degenerate nonprecessing case with all roots equal. Cannot continue with computation.")
477 mm = (Ssqp - Ssqm)/(Ssqp - Ssq3)
484 return mp.polyval([vareps, Btransf, Ctransf, Dtransf], Ssq)
486 return ((vareps*Ssq + Btransf)*Ssq + Ctransf)*Ssq + Dtransf
488 def cubic_deriv(Ssq):
490 return mp.polyval([3*vareps, 2*Btransf, Ctransf], Ssq)
492 return (3.*vareps*Ssq + 2.*Btransf)*Ssq + Ctransf
496 if Ssqp - Ssqm < 1.1*max(0.5*abs(Ssqp + Ssqm), 1.)*root_tol:
497 root_tol_pm = 0.5*(Ssqp - Ssqm)/max(0.5*abs(Ssqp + Ssqm), 1.)
499 root_tol_pm = root_tol
501 if Ssqm - Ssq3 < 1.1*max(0.5*abs(Ssqm + Ssq3), 1.)*root_tol:
502 root_tol_m3 = 0.5*(Ssqm - Ssq3)/max(0.5*abs(Ssqm + Ssq3), 1.)
504 root_tol_m3 = root_tol
506 root_tol_message_extra_bit_flag =
False
508 root_tol_orig = root_tol
510 if root_tol_pm < root_tol
and root_tol_pm <= root_tol_m3:
511 root_tol = root_tol_pm
513 root_tol_message_extra_bit_flag =
True
515 root_tol_message_extra_bit_insert =
"S_+^2 - S_-^2"
516 elif root_tol_m3 < root_tol:
517 root_tol = root_tol_m3
519 root_tol_message_extra_bit_flag =
True
521 root_tol_message_extra_bit_insert =
"S_-^2 - S_3^2"
523 if root_tol_message_extra_bit_flag:
524 root_tol_message_extra_bit =
" the internally modified (since %s is small) version "\
525 "of"%root_tol_message_extra_bit_insert
527 root_tol_message_extra_bit =
""
529 Ssqp_resp, Ssqp_resm = cubic(
530 Ssqp + max(abs(Ssqp), 1.)*root_tol), cubic(Ssqp - max(abs(Ssqp), 1.)*root_tol)
531 Ssqm_resp, Ssqm_resm = cubic(
532 Ssqm + max(abs(Ssqm), 1.)*root_tol), cubic(Ssqm - max(abs(Ssqm), 1.)*root_tol)
537 Btransf_out, Ctransf_out, Dtransf_out = Btransf, Ctransf, Dtransf
539 Btransf_out, Ctransf_out, Dtransf_out = Btransf[0], Ctransf[0], Dtransf[0]
545 if Ssqm - Ssq3 < 1.1*max(0.5*abs(Ssqm + Ssq3), 1.)*root_tol_orig:
546 root_tol_m0 = 0.5*(Ssqm - Ssq3)/max(0.5*abs(Ssqm + Ssq3), 1.)
548 root_tol_message_extra_bit =
" the internally modified (since S_-^2 - S_3^2 is small) "\
551 root_tol_m0 = root_tol_orig
553 root_tol_message_extra_bit =
""
554 Ssqp_resp, Ssqp_resm = cubic_deriv(
555 Ssqp + max(abs(Ssqp), 1.)*root_tol_m0), cubic_deriv(Ssqp - max(abs(Ssqp), 1.)*root_tol_m0)
557 if np.sign(Ssqp_resp) == np.sign(Ssqp_resm):
558 root_tol_message =
"Error in S_+^2 is greater than%s the tolerance root_tol = %e--the "\
559 "residual does not change sign over the interval determined by root_tol "\
560 "(S_+^2 = S_-^2, so m = 0 and we do not need to consider S_3^2)\n"\
561 "Printing the residuals for S^2 = S_+^2 + max(abs(S_+^2),1)*root_tol, S_+^2 "\
562 "- max(abs(S_+^2),1)*root_tol, as well as the values of u, S_+^2, "\
563 "q*(1 - q^2)*u^2, bar{B}, bar{C}, bar{D}\n%e %e %e %f %e %e %e %e" % \
564 (root_tol_message_extra_bit, root_tol_m0, Ssqp_resp, Ssqp_resm, u, Ssqp, vareps,
565 Btransf_out, Ctransf_out, Dtransf_out)
567 if root_tol_raise_err:
568 raise RuntimeError(root_tol_message)
570 warn(root_tol_message, RuntimeWarning)
577 if Ssqp - Ssqm < 1.1*max(0.5*abs(Ssqp + Ssqm), 1.)*root_tol_orig:
578 root_tol_m1 = 0.5*(Ssqp - Ssqm)/max(0.5*abs(Ssqp + Ssqm), 1.)
580 root_tol_message_extra_bit =
" the internally modified (since S_+^2 - S_-^2 is small) "\
583 root_tol_m1 = root_tol_orig
585 root_tol_message_extra_bit =
""
586 Ssqm_resp, Ssqm_resm = cubic_deriv(
587 Ssqm + max(abs(Ssqm), 1.)*root_tol_m1), cubic_deriv(Ssqm - max(abs(Ssqm), 1.)*root_tol_m1)
589 if np.sign(Ssqm_resp) == np.sign(Ssqm_resm):
590 root_tol_message =
"Error in S_-^2 is greater than%s the tolerance root_tol = %e--the "\
591 "residual does not change sign over the interval determined by root_tol "\
592 "(S_-^2 = S_3^2, so m = 1 and we do not need to consider S_+^2)\n"\
593 "Printing the residuals for S^2 = S_-^2 + max(abs(S_-^2),1)*root_tol, S_-^2 "\
594 "- max(abs(S_-^2),1)*root_tol, as well as the values of u, S_-^2, "\
595 "q*(1 - q^2)*u^2, bar{B}, bar{C}, bar{D}\n%e %e %e %f %e %e %e %e" % \
596 (root_tol_message_extra_bit, root_tol_m1, Ssqm_resp, Ssqm_resm, u, Ssqm, vareps,
597 Btransf_out, Ctransf_out, Dtransf_out)
599 if root_tol_raise_err:
600 raise RuntimeError(root_tol_message)
602 warn(root_tol_message, RuntimeWarning)
605 elif E_ratio_term(mm) <= lin_tol*max(1., 1./(Ssqp - Ssqm)):
606 if np.sign(Ssqp_resp) == np.sign(Ssqp_resm)
or np.sign(Ssqm_resp) == np.sign(Ssqm_resm):
607 root_tol_message =
"Error in S_+^2 and/or S_-^2 is greater than%s the tolerance root_tol = %e--the "\
608 "residual does not change sign over the interval determined by root_tol "\
609 "(m is small enough to be linearized in, so we do not consider S_3^2)\n"\
610 "Printing the residuals for S^2 = S_+^2 + max(abs(S_+^2),1)*root_tol, S_+^2 "\
611 "- max(abs(S_+^2),1)*root_tol, S_-^2 + max(abs(S_-^2),1)*root_tol, "\
612 "S_-^2 - max(abs(S_-^2),1)*root_tol, as well as the values of u, S_+^2, S_-^2, "\
613 "q*(1 - q^2)*u^2, bar{B}, bar{C}, bar{D}\n%e %e %e %e %e %f %f %e %e %e %e" % \
614 (root_tol_message_extra_bit, root_tol, Ssqp_resp, Ssqp_resm, Ssqm_resp, Ssqm_resm,
615 u, Ssqp, Ssqm, vareps, Btransf_out, Ctransf_out, Dtransf_out)
617 if root_tol_raise_err:
618 raise RuntimeError(root_tol_message)
620 warn(root_tol_message, RuntimeWarning)
623 Ssq_pr = 0.5*(Ssqp + Ssqm)
625 Ssq3_resp, Ssq3_resm = cubic(
626 Ssq3 + max(abs(Ssq3), 1.)*root_tol), cubic(Ssq3 - max(abs(Ssq3), 1.)*root_tol)
628 if np.sign(Ssqp_resp) == np.sign(Ssqp_resm)
or np.sign(Ssqm_resp) == np.sign(Ssqm_resm) \
629 or np.sign(Ssq3_resp) == np.sign(Ssq3_resm):
630 root_tol_message =
"Error in S_+^2, S_-^2, and/or S_3^2 is greater than%s the tolerance root_tol = "\
631 "%e--the residual does not change sign over the interval determined by "\
632 "root_tol\nPrinting the residuals for S^2 = S_+^2 + max(abs(S_+^2),1)*root_tol, "\
633 "S_+^2 - max(abs(S_+^2),1)*root_tol, S_-^2 + max(abs(S_-^2),1)*root_tol, "\
634 "S_-^2 - max(abs(S_-^2),1)*root_tol, S_3^2 + max(abs(S_3^2),1)*root_tol, "\
635 "S_3^2 - max(abs(S_-^2),1)*root_tol, as well as the values of u, S_+^2, S_-^2, "\
636 "S_3^2, q*(1 - q^2)*u^2, bar{B}, bar{C}, bar{D}\n"\
637 "%e %e %e %e %e %e %e %f %f %f %e %e %e %e" % \
638 (root_tol_message_extra_bit, root_tol, Ssqp_resp, Ssqp_resm, Ssqm_resp, Ssqm_resm,
639 Ssq3_resp, Ssq3_resm, u, Ssqp, Ssqm, Ssq3, vareps, Btransf_out, Ctransf_out,
642 if root_tol_raise_err:
643 raise RuntimeError(root_tol_message)
645 warn(root_tol_message, RuntimeWarning)
648 ellip_ratio = mp.ellipe(mm)/mp.ellipk(mm)
650 ellip_ratio = ellipe(mm)/ellipk(mm)
653 Ssq_pr = Ssqp + (ellip_ratio - 1.)*(Ssqp - Ssq3)
660def kappaxiq_evol(q, S1, S2, S1L, xi, S0sq, u0, uf, du, method, imag_tol, root_tol, lin_tol, ode_atol, ode_rtol,
661 root_tol_raise_err, use_solve_ivp=False, solve_ivp_lsoda_min_step=0., use_mpmath=False,
662 polyroots_maxsteps=50, polyroots_extraprec=50, odefun_degree=None):
664 Compute kappaxiq when the binary has u = uf starting from u = u0
670 q: Mass ratio of the binary,
with the < 1 convention
671 S1, S2: Magnitudes of the dimensionful spins of the binary (
with total mass = 1)
672 S1L: Component of the dimensionful spin of hole 1 along the orbital angular momentum (
with total mass = 1)
673 xi: Conserved effective spin, defined
in Eq. (2) of the paper <https://dcc.ligo.org/P2100029>, arXiv:2107.11902
674 S0sq: Square of initial dimensionful total spin (
in total mass = 1 units)
675 u0: Initial value of u
677 du: Step
in u
for the evolution
678 method: method to use
for scipy integration (
"vode",
"lsoda",
"dopri5",
or "dop853" for the ode interface
and
679 "RK45",
"RK23",
"DOP853",
or "LSODA" for the solve_ivp interface);
not used
if use_mpmath ==
True
680 imag_tol: Tolerance
for the imaginary part of roots of the cubic
681 root_tol: Tolerance
for the error
in the roots of the cubic (
is an absolute tolerance
for small values of the roots,
as
682 generally occurs during the beginning of the evolution,
and is a fractional tolerance when the roots are large,
683 which generally occurs towards the end of the evolution; also
is reduced automatically when the roots are close
685 lin_tol: Tolerance
for errors on the r.h.s. of dkappa/du due to linearization
in m (
is a fractional tolerance when m
is
686 large,
as generally occurs during the beginning of the evolution,
and is an absolute tolerance when m
is small,
687 as generally occurs towards the end of the evolution)
688 ode_atol: Absolute tolerance
for the default scipy ode solver, also the single tolerance setting
for the mpmath ode solver
690 ode_rtol: Relative tolerance
for the default scipy ode solver
691 root_tol_raise_err: Whether to
raise an error (
True)
or just
print a warning (
False) when the error
in the roots of
692 the cubic
is larger than root_tol
696 use_solve_ivp: If
True, use scipy
's solve_ivp function for the ODE integration, which gives access to a different set of
697 integrators and allows tighter tolerances
for LSODA than the default ode function, but can hang
for certain
698 parameters, default:
False
699 solve_ivp_lsoda_min_step: Minimum step to use
in the solve_ivp LSODA integrator, default: 0.
701 use_mpmath: If
True, use mpmath functions
and ode integration instead of numpy/scipy (slower but more accurate), default:
False
702 polyroots_maxsteps: The maximum number of steps allowed
in mpmath.polyroots
if use_mpmath ==
True, default: 50
703 polyroots_extraprec: The number of extra digits to use
in mpmath.polyroots
if use_mpmath ==
True, default: 50
704 odefun_degree: Degree of polynomial used by mpmath.odefun
if use_mpmath ==
True; calculated automatically by
705 mpmath.odefun
if None, default:
None
707 Output: kappaxiq at Lf
718 solve_ivp_kwargs = {
'atol': ode_atol,
'rtol': ode_rtol}
720 if method ==
"LSODA":
721 solve_ivp_kwargs[
'min_step'] = solve_ivp_lsoda_min_step
724 soln = solve_ivp(dkappaxiq_du, [u0, uf], [kappaxiq0], method=method, t_eval=[uf], args=(q, S1, S2, xi, S0sq, imag_tol,
725 root_tol, lin_tol, root_tol_raise_err), **solve_ivp_kwargs)
735 soln = mp.odefun(
lambda ub, kappaxiq: -
dkappaxiq_du(-ub, kappaxiq, q, S1, S2, xi, S0sq, imag_tol, root_tol,
736 lin_tol, root_tol_raise_err, use_mpmath=
True,
737 polyroots_maxsteps=polyroots_maxsteps,
738 polyroots_extraprec=polyroots_extraprec),
739 -u0, kappaxiq0, ode_atol, odefun_degree)
748 solver = ode(dkappaxiq_du)
749 solver.set_integrator(method, atol=ode_atol, rtol=ode_rtol)
750 solver.set_initial_value(kappaxiq0, u0)
751 solver.set_f_params(q, S1, S2, xi, S0sq, imag_tol,
752 root_tol, lin_tol, root_tol_raise_err)
758 while solver.successful()
and not end:
759 new_t = solver.t - du
765 solver.integrate(new_t, step=1)
767 if abs(solver.t - uf) > 0.75*du:
769 "Final u of %e is more than 0.75*du different from uf of %e" % (solver.t, uf))
777 ode_atol_inplane_spin_scaling=True, ode_atol_floor=1.e-13, ode_rtol=-1, method="lsoda",
778 use_solve_ivp=False, solve_ivp_lsoda_min_step=0., use_fallback=True, du_fallback=1.e-5,
779 ode_atol_fallback=1.e-13, ode_rtol_fallback=-1, method_fallback="lsoda",
780 use_mpmath_fallback=True, mpmath_dps=30, ode_tol_mpmath=1.e-15, odefun_degree=None,
781 polyroots_maxsteps=50, polyroots_extraprec=50, LPNorder=2, LPNspins=True, imag_tol=1.e-6,
782 root_tol=1.e-6, lin_tol=-1, lin_tol_fallback=-1, root_tol_raise_err=False,
783 failure_mode='None', version='v1'):
785 Compute tilt angles at infinite separation or bounds
and average value of tilt angles at finite separation
786 from masses, spins, tilt angles,
and in-plane spin angle at some frequency, taken to be small enough that
787 the precession averaged spin evolution
from Gerosa et al. PRD 92, 064016 (2015) [implemented following
788 Chatziioannou, Klein, Yunes,
and Cornish PRD 95, 104004 (2017)]
is accurate, using the regularized expressions
789 from <https://dcc.ligo.org/P2100029>, arXiv:2107.11902. This version takes vectors of spins
as input (
as opposed
790 to spin angles)
and allows one optionally to specify the initial orbital angular momentum instead of computing it
791 from fref (
and the binary
's parameters).
797 m1, m2: Detector frame masses of the binary, in kg
798 chi1_vec, chi2_vec: 3d numpy arrays
with the dimensionless spins of the binary; these must be
in coordinates where
799 the Newtonian orbital angular momentum
is in the z-direction
if L0_vec
is None (the default),
800 otherwise all that
is necessary
is that chi1_vec, chi2_vec,
and L0_vec are all
in the same
802 fref: Reference frequency,
in Hz; ignored when L0_vec
is not None
804 - Optional initial
and final L settings:
806 L0_vec: 3d numpy array
with the initial orbital angular momentum
in units where the binary
's total mass = 1 and in
807 the same coordinates as the spins; this
is computed
from fref
and the binary parameters
if None is passed
808 and in that case given
in coordinates where the Newtonian orbital angular momentum
is in the z-direction,
810 Lf: Final magnitude of orbital angular momentum
in total mass = 1 units
for output at finite separation,
811 None gives the output at infinity, default:
None
813 - Optional settings
for the evolution:
815 du: Step
in u
for the evolution (external to the ODE integrator), default: 1.e-3
816 ode_atol_base: Base value of absolute tolerance
for ode solver, default: 1.e-8
817 ode_atol_inplane_spin_scaling: Scale ode_atol_base by min(chi1*np.sin(tilt1),chi2*np.sin(tilt2)),
with a floor of
818 ode_atol_floor, to give good accuracy
for small
in-plane spin cases, default:
True
819 ode_atol_floor: Floor value
for absolute tolerance
for ode solver used when ode_atol_inplane_spin_scaling ==
True;
820 should be several orders of magnitude less than ode_atol_base, default: 1.e-13
821 ode_rtol: Relative tolerance
for ode solver, -1 sets this to 0.1*ode_atol, default: -1
822 method: method to use
for integration (either
from scipy
's ode interface, where "vode", "lsoda", "dopri5", and
823 "dop853" have been tested, though
"lsoda" is the recommended one;
"RK45",
"RK23",
"DOP853",
or "LSODA" for
824 scipy
's solve_ivp interface, if use_solve_ivp == True ("lsoda" and "dop853" are equivalent to "LSODA" and
825 "DOP853" here,
for compatibility);
or "odefun" to use mpmath.odefun
with its default Taylor method, which
826 is all that
is available
as of v1.2.0; this uses the mpmath fallback settings
and disables the fallback
827 evolution--it
is also quite slow), default:
"lsoda"
828 use_solve_ivp: If
True, use scipy
's solve_ivp interface and its integrators instead of the default ode interface; note
829 that the solve_ivp integrators apparently hang for certain difficult cases, default:
False
830 solve_ivp_lsoda_min_step: Minimum step to use
in the solve_ivp LSODA integrator, default: 0.
831 use_fallback: Whether to use the fallback integration
if the primary integration fails (
True)
or not (
False),
832 not available
if method ==
"odefun", default:
True
833 du_fallback: Step
in u
for the evolution
if the first evolution fails, default: 1.e-5
834 ode_atol_fallback: Absolute tolerance
for ode solver
if the first evolution fails, default: 1.e-13
835 ode_rtol_fallback: Relative tolerance
for ode solver
if the first evolution fails, -1 sets this to
836 0.1*ode_atol_fallback, default: -1
837 method_fallback: method to use
for integration
if the first evolution fails (any of the scipy integrators available
for method),
839 use_mpmath_fallback: Use the slow but accurate mpmath integration
if the first fallback integration fails (
True)
or not (
False),
841 mpmath_dps: The number of decimal places to use
for the mpmath computations,
if the mpmath integration
is triggered
or
842 method ==
"odefun", default: 30
843 ode_tol_mpmath: Tolerance
for the mpmath ode solver, only used
if the mpmath fallback evolution
is triggered
or
844 method ==
"odefun", default: 1.e-15
845 odefun_degree: Degree of polynomial used by mpmath.odefun
if the mpmath fallback evolution
is triggered
or method ==
"odefun";
846 calculated automatically by mpmath.odefun
if None, default:
None
847 polyroots_maxsteps: The maximum number of steps allowed
in mpmath.polyroots
if the mpmath fallback evolution
is triggered
or
848 method ==
"odefun", default: 50
849 polyroots_extraprec: The number of extra digits to use
in mpmath.polyroots
if the mpmath fallback evolution
is triggered
or
850 method ==
"odefun", default: 50
851 LPNorder: PN order to use
in computing the intial orbital angular momentum
from fref--choices are
852 0, 1, 1.5, 2,
or 2.5, where the half-integer orders are the spin contributions; ignored when L0_vec
is
854 LPNspins: Whether to include the orbit-averaged spin contributions
in the computation of the initial orbital
855 angular momentum; ignored when L0_vec
is not None, default:
True
856 imag_tol: Tolerance
for the imaginary part of roots of the cubic, default: 1.e-6
857 root_tol: Tolerance
for the error
in the roots of the cubic (
is an absolute tolerance
for small values of the roots,
as
858 generally occurs during the beginning of the evolution,
and is a fractional tolerance when the roots are large,
859 which generally occurs towards the end of the evolution; also
is reduced automatically when the roots are close
860 together), default: 1.e-6
861 lin_tol: Tolerance
for errors on the r.h.s. of dkappa/du due to linearization
in m (
is a fractional tolerance when m
is
862 large,
as generally occurs during the beginning of the evolution,
and is an absolute tolerance when m
is small,
863 as generally occurs towards the end of the evolution); -1 sets this to be the same
as ode_atol,
except for the
864 mpmath evolution, where it
is always set to the ode_tol_mpmath value, default: -1
865 lin_tol_fallback: Tolerance
for errors on the r.h.s. of dkappa/du due to linearization
in m
if the first evolution
866 fails, -1 sets this to be the same
as ode_atol_fallback, default: -1
867 root_tol_raise_err: Whether to
raise an error (
True)
or just
print a warning (
False) when the error
in the roots
868 of the cubic
is larger than root_tol, default:
False
869 failure_mode: How the code behaves when the evolution fails (after the fallback evolution,
if use_fallback
is True, the
870 default).
'Error' means that the code raises a RuntimeError,
while 'NAN' and 'None' mean that the code
871 raises a RuntimeWarning
and returns np.nan
or None for the output, respectively, default:
'None'
872 version: Version of the calculation to use--currently only the initial version, v1,
is available, default:
"v1"
874 Output: dictionary
with entries
'tilt1_inf',
'tilt2_inf' for evolution to infinity
and entries
'tilt1_sep_min',
875 'tilt1_sep_max',
'tilt1_sep_avg',
'tilt2_sep_min',
'tilt2_sep_max',
'tilt2_sep_avg' for evolution to a
876 finite separation (i.e., a finite orbital angular momentum)
882 tested_scipy_integrators = [
"RK45",
"RK23",
"DOP853",
"LSODA"]
884 tested_scipy_integrators = [
"vode",
"lsoda",
"dopri5",
"dop853"]
886 tested_scipy_integrator_string =
', '.join(tested_scipy_integrators)
890 known_failure_modes = [
'Error',
'NAN',
'None']
895 raise ValueError(
"Only version v1 is available at the moment, while version = %s"%version)
899 methods_to_uppercase = [
"lsoda",
"dop853"]
902 if method
in methods_to_uppercase:
903 method = method.upper()
905 if method_fallback
in methods_to_uppercase:
906 method_fallback = method_fallback.upper()
912 if type(chi1_vec)
is not np.ndarray
or type(chi2_vec)
is not np.ndarray
or len(np.atleast_1d(chi1_vec)) != 3
or \
913 len(np.atleast_1d(chi2_vec)) != 3:
915 "chi1_vec and chi2_vec must both be numpy arrays of length 3, while they are", chi1_vec, chi2_vec)
917 if L0_vec
is not None:
920 "L0_vec must be either None or a numpy array of length 3, while it is", L0_vec)
922 chi1 = norm(chi1_vec)
923 chi2 = norm(chi2_vec)
925 if chi1 == 0.
or chi1 > 1.
or chi2 == 0.
or chi2 > 1.:
927 "The magnitudes of the spins must both be positive and at most 1, while they are chi1, chi2 = %f, %f" %
931 check_fref(fref, m1, m2,
"precession-averaged")
933 raise ValueError(
"One of fref and L0_vec must not be None")
935 if du <= 0.
or du_fallback <= 0.:
937 "The step in u and its fallback value must both be positive, while they are du, du_fallback = %e, %e" %
940 if use_solve_ivp
and not solve_ivp_avail:
943 raise ValueError(
"solve_ivp is not available in the version of scipy being used (v%s), so use_solve_ivp == True "
944 "is not allowed."%scipy.__version__)
947 interface_str =
"solve_ivp"
949 interface_str =
"ode"
951 if method
not in np.append(tested_scipy_integrators,
"odefun")
or method_fallback
not in tested_scipy_integrators:
952 raise ValueError(
"The integration method must be odefun or one of the tested scipy integrators available with "
953 "the %s interface, %s, and the fallback integration method must be one of the tested scipy "
954 "integrators, while they are method, method_fallback = %s, %s" %
955 (interface_str, tested_scipy_integrator_string, method, method_fallback))
957 if method ==
"odefun" and not mpmath_avail:
958 raise ValueError(
"mpmath is not available, so odefun is not a possible choice for the integration method, "
959 "which must in this case be one of the tested scipy integrators: %s"%tested_scipy_integrator_string)
961 if mpmath_dps <= 0
or not isinstance(mpmath_dps, int):
962 raise ValueError(
"The number of decimal places to use for the mpmath computations must be a natural number, "
963 "while it is mpmath_dps = %f"%mpmath_dps)
965 if odefun_degree
is not None and (odefun_degree <= 0
or not isinstance(odefun_degree, int)):
966 raise ValueError(
"The degree for the mpmath.odefun must be a natural number, while it is odefun_degree = %f"%odefun_degree)
968 if polyroots_maxsteps <= 0
or polyroots_extraprec <= 0
or not isinstance(polyroots_maxsteps, int) \
969 or not isinstance(polyroots_extraprec, int):
970 raise ValueError(
"The maximum number of steps and extra precision for mpmath.polyroots must both be natural "
971 "numbers, while they are polyroots_maxsteps, polyroots_extraprec = %f, %f"%(polyroots_maxsteps,
972 polyroots_extraprec))
974 if L0_vec
is None and LPNorder
not in LPN_orders:
975 raise ValueError(
"The PN order to use to compute the initial orbital angular momentum must be one of %s, "
976 "while it is LPNorder = %f" % (LPN_order_string, LPNorder))
978 if failure_mode
not in known_failure_modes:
979 known_failure_modes_string =
', '.join(known_failure_modes)
980 raise ValueError(
"The failure mode must be one of %s, while it is failure_mode = %s" %
981 (known_failure_modes_string, failure_mode))
985 if fref
is not None and L0_vec
is not None:
986 warn(
"fref and L0_vec are both not None. In this case, the L0_vec value is used to initialize the "
987 "evolution and fref is not used.", ValueWarning)
992 if failure_mode ==
'NAN':
993 failure_output = np.nan
994 failure_output_string =
'np.nan'
996 failure_output =
None
997 failure_output_string =
'None'
1001 def failure_message_mpmath(err):
1002 return "Encountered the exception\n\n%s\n\nusing the mpmath evolution and mpmath_dps = %d, " \
1003 "ode_tol_mpmath = %e, odefun_degree = %s, polyroots_maxsteps = %d, polyroots_extraprec " \
1004 "= %d. However, it may be possible to evolve this case successfully using more digits " \
1005 "and/or a tighter tolerance "%(
format_error(err), mpmath_dps, ode_tol_mpmath, str(odefun_degree),
1006 polyroots_maxsteps, polyroots_extraprec)
1016 chi1, chi2 = chi2, chi1
1017 chi1_vec, chi2_vec = chi2_vec, chi1_vec
1038 S1_vec = m1*m1*chi1_vec
1039 S2_vec = m2*m2*chi2_vec
1044 if L0_vec
is not None:
1048 S1_vec_for_L = S1_vec
1049 S2_vec_for_L = S2_vec
1051 S1_vec_for_L = np.zeros(3)
1052 S2_vec_for_L = np.zeros(3)
1055 m1_kg, m2_kg, fref, PNorder=LPNorder, S1_vec=S1_vec_for_L, S2_vec=S2_vec_for_L, no_mass_spin_input_check=
True)
1060 raise ValueError(
"L0_vec must not be the zero vector")
1062 if Lf
is not None and Lf <= L0:
1064 "The magnitude of the final orbital angular momentum must be greater than that of the "
1065 "initial orbital angular momentum, %f, while it is Lf = %f" % (L0, Lf))
1071 chi1_inplane = norm(chi1_vec - np.dot(chi1_vec, L0_hat)*L0_hat)
1072 chi2_inplane = norm(chi2_vec - np.dot(chi2_vec, L0_hat)*L0_hat)
1074 if ode_atol_inplane_spin_scaling:
1075 ode_atol = max(min(chi1_inplane, chi2_inplane)
1076 * ode_atol_base, ode_atol_floor)
1078 ode_atol = ode_atol_base
1083 ode_rtol = 0.1*ode_atol
1085 if ode_rtol_fallback == -1:
1086 ode_rtol_fallback = 0.1*ode_atol_fallback
1091 if lin_tol_fallback == -1:
1092 lin_tol_fallback = ode_atol_fallback
1096 if ode_atol_base <= 0.
or ode_atol_floor <= 0.
or ode_rtol <= 0.
or ode_atol_fallback <= 0. \
1097 or ode_rtol_fallback <= 0.
or ode_tol_mpmath <= 0.:
1099 "The ode solver tolerance values and their fallback values must all be positive, while they are "
1100 "ode_atol_base, ode_atol_floor, ode_rtol, ode_atol_fallback, ode_rtol_fallback, "
1101 "ode_tol_mpmath = %e, %e, %e, %e, %e, %e" %
1102 (ode_atol_base, ode_atol_floor, ode_rtol, ode_atol_fallback, ode_rtol_fallback, ode_tol_mpmath))
1104 if ode_atol_inplane_spin_scaling
and ode_atol_floor >= ode_atol_base:
1105 raise ValueError(
"ode_atol_floor should be less than ode_atol_base, likely by several orders of magnitude, "
1106 "while they are ode_atol_floor, ode_atol_base = %e, %e" %
1107 (ode_atol_floor, ode_atol_base))
1109 if imag_tol < 0.
or root_tol < 0.
or lin_tol < 0.
or lin_tol_fallback < 0.:
1110 raise ValueError(
"The tolerance values must all be nonnegative, while they are imag_tol, root_tol, lin_tol, "
1111 "lin_tol_fallback = %e, %e, %e, %e" % (imag_tol, root_tol, lin_tol, lin_tol_fallback))
1115 xi = np.dot((1. + q)*S1_vec + (1. + 1./q)*S2_vec, L0_hat)
1119 S1L = np.dot(S1_vec, L0_hat)
1120 S2L = np.dot(S2_vec, L0_hat)
1124 cos_tilt1 = np.clip(S1L/S1, -1., 1.)
1125 cos_tilt2 = np.clip(S2L/S2, -1., 1.)
1127 tilt1 = np.arccos(cos_tilt1)
1128 tilt2 = np.arccos(cos_tilt2)
1130 if (tilt1 == np.pi
and tilt2
in [0., np.pi])
or (tilt1 == 0.
and tilt2 == 0.):
1135 S0sq = norm(S1_vec + S2_vec)**2
1146 mpmath_tilts =
False
1148 if method ==
"odefun":
1152 with mp.workdps(mpmath_dps):
1153 k_evol =
kappaxiq_evol(q, S1, S2, S1L, xi, S0sq, u0, uf, du, method, imag_tol, root_tol, ode_tol_mpmath,
1154 ode_tol_mpmath, ode_rtol, root_tol_raise_err, use_mpmath=
True,
1155 polyroots_maxsteps=polyroots_maxsteps, polyroots_extraprec=polyroots_extraprec,
1156 odefun_degree=odefun_degree)
1158 except (RuntimeError, mp.ctx_base.StandardBaseContext.NoConvergence)
as err:
1159 return evolution_error_handling(failure_mode, failure_message_mpmath(err), failure_output, failure_output_string, Lf, swap)
1161 except NonprecessingError
as ne:
1165 k_evol =
kappaxiq_evol(q, S1, S2, S1L, xi, S0sq, u0, uf, du, method, imag_tol, root_tol, lin_tol,
1166 ode_atol, ode_rtol, root_tol_raise_err, use_solve_ivp=use_solve_ivp,
1167 solve_ivp_lsoda_min_step=solve_ivp_lsoda_min_step)
1169 except RuntimeError
as re:
1171 warn(
"Encountered the exception\n\n%s\n\nusing method %s and du = %e, ode_atol = %e, ode_rtol = %e, "
1172 "lin_tol = %e, trying with %s with du = %e, ode_atol = %e, ode_rtol = %e, lin_tol = %e" %
1173 (
format_error(re), method, du, ode_atol, ode_rtol, lin_tol, method_fallback, du_fallback,
1174 ode_atol_fallback, ode_rtol_fallback, lin_tol_fallback), RuntimeWarning)
1177 k_evol =
kappaxiq_evol(q, S1, S2, S1L, xi, S0sq, u0, uf, du_fallback, method_fallback, imag_tol,
1178 root_tol, lin_tol_fallback, ode_atol_fallback, ode_rtol_fallback,
1179 root_tol_raise_err, solve_ivp_lsoda_min_step=solve_ivp_lsoda_min_step)
1181 except RuntimeError
as re2:
1182 if use_mpmath_fallback
and mpmath_avail:
1183 warn(
"Encountered the exception\n\n%s\n\nin the first fallback evolution. Now trying with the mpmath "
1184 "fallback with mpmath_dps = %d, ode_tol_mpmath = %e, odefun_degree = %s, "
1185 "polyroots_maxsteps = %d, polyroots_extraprec = %d"%(
format_error(re2), mpmath_dps, ode_tol_mpmath, str(odefun_degree),
1186 polyroots_maxsteps, polyroots_extraprec), RuntimeWarning)
1191 with mp.workdps(mpmath_dps):
1192 k_evol =
kappaxiq_evol(q, S1, S2, S1L, xi, S0sq, u0, uf, du, method, imag_tol, root_tol, ode_tol_mpmath,
1193 ode_tol_mpmath, ode_rtol, root_tol_raise_err, use_mpmath=
True,
1194 polyroots_maxsteps=polyroots_maxsteps, polyroots_extraprec=polyroots_extraprec,
1195 odefun_degree=odefun_degree)
1197 except (RuntimeError, mp.ctx_base.StandardBaseContext.NoConvergence)
as err:
1198 return evolution_error_handling(failure_mode, failure_message_mpmath(err), failure_output, failure_output_string, Lf, swap)
1200 except NonprecessingError
as ne:
1203 if not use_mpmath_fallback:
1204 mpmath_bit =
"it is disabled"
1206 mpmath_bit =
"mpmath is not available"
1208 failure_message =
"Fallback evolution failed with the exception\n\n%s\n\nThe mpmath fallback evolution might work in this " \
1209 "case, but %s. Alternatively, you can try adjusting the tolerance settings for the evolution, but this " \
1210 "case may not be possible to evolve accurately with the scipy integration." % (
format_error(re2), mpmath_bit)
1214 except NonprecessingError
as ne:
1218 failure_message =
"Evolution failed with the exception\n\n%s\n\nYou can try adjusting the tolerance settings " \
1219 "for the evolution, but this case may not be possible to evolve accurately with the current " \
1220 "version of the precession-averaged evolution."%
format_error(re)
1224 except NonprecessingError
as ne:
1233 cos_tilt1inf = k_evol/S1
1234 cos_tilt2inf = q*(xi/(1. + q) - k_evol)/S2
1238 cos_tilt1inf = np.clip(cos_tilt1inf, -1., 1.)
1239 cos_tilt2inf = np.clip(cos_tilt2inf, -1., 1.)
1244 tilt1inf = float(mp.acos(cos_tilt1inf))
1245 tilt2inf = float(mp.acos(cos_tilt2inf))
1247 tilt1inf = np.arccos(cos_tilt1inf[0])
1248 tilt2inf = np.arccos(cos_tilt2inf[0])
1252 vareps, Btransf, Ctransf, Dtransf =
eq_coeffs(
1253 uf, k_evol, q, S1, S2, xi, S0sq)
1256 [vareps, Btransf, Ctransf, Dtransf], uf, k_evol, q, S1, S2, imag_tol, use_mpmath=mpmath_tilts,
1257 root_tol=root_tol, polyroots_maxsteps=polyroots_maxsteps, polyroots_extraprec=polyroots_extraprec)
1260 lin_tol_use = ode_tol_mpmath
1262 lin_tol_use = lin_tol
1264 Ssqavg =
dkappaxiq_du(uf, k_evol, q, S1, S2, xi, S0sq, imag_tol, root_tol, lin_tol_use, root_tol_raise_err,
1265 use_mpmath=mpmath_tilts, polyroots_maxsteps=polyroots_maxsteps,
1266 polyroots_extraprec=polyroots_extraprec)
1270 cos_tilt1sep_Sm = (k_evol - uf*Ssqm)/S1
1271 cos_tilt2sep_Sm = q*(xi/(1. + q) - S1*cos_tilt1sep_Sm)/S2
1273 cos_tilt1sep_Sp = (k_evol - uf*Ssqp)/S1
1274 cos_tilt2sep_Sp = q*(xi/(1. + q) - S1*cos_tilt1sep_Sp)/S2
1276 cos_tilt1sep_Savg = (k_evol - uf*Ssqavg)/S1
1277 cos_tilt2sep_Savg = q*(xi/(1. + q) - S1*cos_tilt1sep_Savg)/S2
1281 cos_tilt1sep_Sm = np.clip(cos_tilt1sep_Sm, -1., 1.)
1282 cos_tilt2sep_Sp = np.clip(cos_tilt2sep_Sp, -1., 1.)
1284 cos_tilt1sep_Sm = np.clip(cos_tilt1sep_Sm, -1., 1.)
1285 cos_tilt2sep_Sp = np.clip(cos_tilt2sep_Sp, -1., 1.)
1287 cos_tilt1sep_Savg = np.clip(cos_tilt1sep_Savg, -1., 1.)
1288 cos_tilt2sep_Savg = np.clip(cos_tilt2sep_Savg, -1., 1.)
1293 cos_tilt1sep_Sm, cos_tilt2sep_Sm, cos_tilt1sep_Sp, cos_tilt2sep_Sp = cos_tilt2sep_Sp, cos_tilt1sep_Sp, \
1294 cos_tilt2sep_Sm, cos_tilt1sep_Sm
1295 cos_tilt1sep_Savg, cos_tilt2sep_Savg = cos_tilt2sep_Savg, cos_tilt1sep_Savg
1298 tilt1sepmin = float(mp.acos(cos_tilt1sep_Sm))
1299 tilt1sepmax = float(mp.acos(cos_tilt1sep_Sp))
1300 tilt1sepavg = float(mp.acos(cos_tilt1sep_Savg))
1302 tilt2sepmin = float(mp.acos(cos_tilt2sep_Sp))
1303 tilt2sepmax = float(mp.acos(cos_tilt2sep_Sm))
1304 tilt2sepavg = float(mp.acos(cos_tilt2sep_Savg))
1306 tilt1sepmin = np.arccos(cos_tilt1sep_Sm[0])
1307 tilt1sepmax = np.arccos(cos_tilt1sep_Sp[0])
1308 tilt1sepavg = np.arccos(cos_tilt1sep_Savg[0])
1310 tilt2sepmin = np.arccos(cos_tilt2sep_Sp[0])
1311 tilt2sepmax = np.arccos(cos_tilt2sep_Sm[0])
1312 tilt2sepavg = np.arccos(cos_tilt2sep_Savg[0])
1314 return {
'tilt1_sep_min': tilt1sepmin,
'tilt1_sep_max': tilt1sepmax,
'tilt1_sep_avg': tilt1sepavg,
1315 'tilt2_sep_min': tilt2sepmin,
'tilt2_sep_max': tilt2sepmax,
'tilt2_sep_avg': tilt2sepavg}
1318def prec_avg_tilt_comp(m1, m2, chi1, chi2, tilt1, tilt2, phi12, fref, Lf=None, **kwargs):
1320 Compute tilt angles at infinite separation or bounds
and average value of tilt angles at finite separation
from
1321 masses, spins, tilt angles,
and in-plane spin angle at some frequency, taken to be small enough that the precession
1322 averaged spin evolution
from Gerosa et al. PRD 92, 064016 (2015) [implemented following Chatziioannou, Klein, Yunes,
1323 and Cornish PRD 95, 104004 (2017)]
is accurate, using the regularized expressions. This takes
in spin angles
and a
1324 reference frequency to initialize the evolution.
1330 m1, m2: Detector frame masses of the binary,
in kg
1331 chi1, chi2: Dimensionless spin magnitudes of the binary
1332 tilt1, tilt2: Tilt angles of the binary
's spins (w.r.t. the Newtonian orbital angular momentum) at fref
1333 phi12: Angle between the in-plane components of the spins at fref
1334 fref: Reference frequency,
in Hz
1338 Lf: Final magnitude of orbital angular momentum
in total mass = 1 units
for output at finite separation
1339 None gives the output at infinity, default:
None
1341 du: Step
in u
for the evolution (external to the ODE integrator), default: 1.e-3
1342 ode_atol_base: Base value of absolute tolerance
for ode solver, default: 1.e-8
1343 ode_atol_inplane_spin_scaling: Scale ode_atol_base by min(chi1*np.sin(tilt1),chi2*np.sin(tilt2)),
with a floor of
1344 ode_atol_floor, to give good accuracy
for small
in-plane spin cases, default:
True
1345 ode_atol_floor: Floor value
for absolute tolerance
for ode solver used when ode_atol_inplane_spin_scaling ==
True;
1346 should be several orders of magnitude less than ode_atol_base, default: 1.e-13
1347 ode_rtol: Relative tolerance
for ode solver, -1 sets this to 0.1*ode_atol, default: -1
1348 method: method to use
for integration (either
from scipy
's ode interface, where "vode", "lsoda", "dopri5", and
1349 "dop853" have been tested, though
"lsoda" is the recommended one;
"RK45",
"RK23",
"DOP853",
or "LSODA" for
1350 scipy
's solve_ivp interface, if use_solve_ivp == True ("lsoda" and "dop853" are equivalent to "LSODA" and
1351 "DOP853" here,
for compatibility);
or "odefun" to use mpmath.odefun
with its default Taylor method, which
1352 is all that
is available
as of v1.2.0; this uses the mpmath fallback settings
and disables the fallback
1353 evolution--it
is also quite slow), default:
"lsoda"
1354 use_solve_ivp: If
True, use scipy
's solve_ivp interface and its integrators instead of the default ode interface; note
1355 that the solve_ivp integrators apparently hang for certain difficult cases, default:
False
1356 solve_ivp_lsoda_min_step: Minimum step to use
in the solve_ivp LSODA integrator, default: 0.
1357 use_fallback: Whether to use the fallback integration
if the primary integration fails (
True)
or not (
False),
1358 not available
if method ==
"odefun", default:
True
1359 du_fallback: Step
in u
for the evolution
if the first evolution fails, default: 1.e-5
1360 ode_atol_fallback: Absolute tolerance
for ode solver
if the first evolution fails, default: 1.e-13
1361 ode_rtol_fallback: Relative tolerance
for ode solver
if the first evolution fails, -1 sets this to
1362 0.1*ode_atol_fallback, default: -1
1363 method_fallback: method to use
for integration
if the first evolution fails (any of the scipy integrators available
for method),
1365 use_mpmath_fallback: Use the slow but accurate mpmath integration
if the first fallback integration fails (
True)
or not (
False),
1367 mpmath_dps: The number of decimal places to use
for the mpmath computations,
if the mpmath integration
is triggered
or
1368 method ==
"odefun", default: 30
1369 ode_tol_mpmath: Tolerance
for the mpmath ode solver, only used
if the mpmath fallback evolution
is triggered
or
1370 method ==
"odefun", default: 1.e-15
1371 odefun_degree: Degree of polynomial used by mpmath.odefun
if the mpmath fallback evolution
is triggered
or method ==
"odefun";
1372 calculated automatically by mpmath.odefun
if None, default:
None
1373 polyroots_maxsteps: The maximum number of steps allowed
in mpmath.polyroots
if the mpmath fallback evolution
is triggered
or
1374 method ==
"odefun", default: 50
1375 polyroots_extraprec: The number of extra digits to use
in mpmath.polyroots
if the mpmath fallback evolution
is triggered
or
1376 method ==
"odefun", default: 50
1377 LPNorder: PN order to use
in computing the intial orbital angular momentum
from fref--choices are
1378 0, 1, 1.5, 2,
or 2.5, where the half-integer orders are the spin contributions, default: 2
1379 LPNspins: Whether to include the orbit-averaged spin contributions
in the computation of the initial orbital
1380 angular momentum, default:
True
1381 imag_tol: Tolerance
for the imaginary part of roots of the cubic, default: 1.e-6
1382 root_tol: Tolerance
for the error
in the roots of the cubic (
is an absolute tolerance
for small values of the roots,
as
1383 generally occurs during the beginning of the evolution,
and is a fractional tolerance when the roots are large,
1384 which generally occurs towards the end of the evolution; also
is reduced automatically when the roots are close
1385 together), default: 1.e-6
1386 lin_tol: Tolerance
for errors on the r.h.s. of dkappa/du due to linearization
in m (
is a fractional tolerance when m
is
1387 large,
as generally occurs during the beginning of the evolution,
and is an absolute tolerance when m
is small,
1388 as generally occurs towards the end of the evolution); -1 sets this to be the same
as ode_atol,
except for the
1389 mpmath evolution, where it
is always set to the ode_tol_mpmath value, default: -1
1390 lin_tol_fallback: Tolerance
for errors on the r.h.s. of dkappa/du due to linearization
in m
if the first evolution
1391 fails, -1 sets this to be the same
as ode_atol_fallback, default: -1
1392 root_tol_raise_err: Whether to
raise an error (
True)
or just
print a warning (
False) when the error
in the roots
1393 of the cubic
is larger than root_tol, default:
False
1394 failure_mode: How the code behaves when the evolution fails (after the fallback evolution,
if use_fallback
is True, the
1395 default).
'Error' means that the code raises a RuntimeError,
while 'NAN' and 'None' mean that the code
1396 raises a RuntimeWarning
and returns np.nan
or None for the output, respectively, default:
'None'
1397 version: Version of the calculation to use--currently only the initial version, v1,
is available, default:
"v1"
1399 Output: dictionary
with entries
'tilt1_inf',
'tilt2_inf' for evolution to infinity
and entries
'tilt1_sep_min',
1400 'tilt1_sep_max',
'tilt1_sep_avg',
'tilt2_sep_min',
'tilt2_sep_max',
'tilt2_sep_avg' for evolution to a
1401 finite separation (i.e., a finite orbital angular momentum)
1406 if chi1 == 0.
or chi2 == 0.:
1417 chi1_vec = chi1*np.array([np.sin(tilt1), 0., np.cos(tilt1)])
1419 chi2_vec = chi2*np.array([np.sin(tilt2)*np.cos(phi12),
1420 np.sin(tilt2)*np.sin(phi12), np.cos(tilt2)])
Exception raised when the evolution finds that the system is nonprecessing to numerical accuracy.
def prec_avg_tilt_comp(m1, m2, chi1, chi2, tilt1, tilt2, phi12, fref, Lf=None, **kwargs)
Compute tilt angles at infinite separation or bounds and average value of tilt angles at finite separ...
def solve_cubic(coeffs, u, kappaxiq, q, S1, S2, imag_tol, use_mpmath=False, root_tol=1.e-6, polyroots_maxsteps=50, polyroots_extraprec=50)
Solve the cubic equation, check for complex roots, and order the roots.
def L_vec_from_masses_spins_and_f(m1, m2, f, PNorder=2, S1_vec=np.zeros(3), S2_vec=np.zeros(3), no_mass_spin_input_check=False)
Compute the PN orbital angular momentum from the binary's masses and spins and m = 2 GW frequency usi...
def L_from_masses_a_and_ecc(m1, m2, a, e)
Compute the magnitude of the Newtonian orbital angular momentum from the binary's masses,...
def E_ratio_term(m)
Compute m*E_ratio(m) [cf.
def kappaxiq_evol(q, S1, S2, S1L, xi, S0sq, u0, uf, du, method, imag_tol, root_tol, lin_tol, ode_atol, ode_rtol, root_tol_raise_err, use_solve_ivp=False, solve_ivp_lsoda_min_step=0., use_mpmath=False, polyroots_maxsteps=50, polyroots_extraprec=50, odefun_degree=None)
Compute kappaxiq when the binary has u = uf starting from u = u0.
def eq_coeffs(u, kappaxiq, q, S1, S2, xi, S0sq)
Compute coefficients for the cubic equation in bar{S}^2 to be solved [Eq.
def dkappaxiq_du(u, kappaxiq, q, S1, S2, xi, S0sq, imag_tol, root_tol, lin_tol, root_tol_raise_err, use_mpmath=False, polyroots_maxsteps=50, polyroots_extraprec=50)
Compute (<S^2>prec - S0^2)/(1 - q) [= dkappa{xi q}/du] using Eq.
def prec_avg_tilt_comp_vec_inputs(m1, m2, chi1_vec, chi2_vec, fref, L0_vec=None, Lf=None, du=1.e-3, ode_atol_base=1.e-8, ode_atol_inplane_spin_scaling=True, ode_atol_floor=1.e-13, ode_rtol=-1, method="lsoda", use_solve_ivp=False, solve_ivp_lsoda_min_step=0., use_fallback=True, du_fallback=1.e-5, ode_atol_fallback=1.e-13, ode_rtol_fallback=-1, method_fallback="lsoda", use_mpmath_fallback=True, mpmath_dps=30, ode_tol_mpmath=1.e-15, odefun_degree=None, polyroots_maxsteps=50, polyroots_extraprec=50, LPNorder=2, LPNspins=True, imag_tol=1.e-6, root_tol=1.e-6, lin_tol=-1, lin_tol_fallback=-1, root_tol_raise_err=False, failure_mode='None', version='v1')
Compute tilt angles at infinite separation or bounds and average value of tilt angles at finite separ...
def L_from_masses_spins_and_f(m1, m2, f, PNorder=2, S1_vec=np.zeros(3), S2_vec=np.zeros(3), no_mass_spin_input_check=False)
Compute the magnitude of the PN orbital angular momentum from the binary's masses and spins and m = 2...
def evolution_error_handling(failure_mode, failure_message, failure_output, failure_output_string, Lf, swap, hybrid_evol=False)
Take care of the error message or warning and returning something for the tilts when the precession-a...
def package_tilts(tilt1, tilt2, Lf, swap)
Package tilts to be returned by prec_avg_tilt_comp_vec_inputs() or prec_avg_tilt_comp() depending on ...
def check_tilts(tilt1, tilt2)
def check_fref(fref, m1, m2, evol_type)
def eq_mass_check(m1, m2, Lf)
def check_spin_mags(chi1, chi2)