Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
calc_tilts_prec_avg_regularized.py
Go to the documentation of this file.
1"""
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.
5
6N. K. Johnson-McDaniel, 2021
7"""
8
9import numpy as np
10
11from numpy.linalg import norm
12
13from numpy.polynomial import Polynomial as P
14
15from scipy.special import ellipe, ellipk
16
17from scipy.integrate import ode
18
19from .tilts_at_infinity_utils import *
20
21from warnings import warn
22
23try:
24 from scipy.integrate import solve_ivp
25except ImportError:
26 solve_ivp_avail = False
27else:
28 solve_ivp_avail = True
29
30try:
31 import mpmath as mp
32except ImportError:
33 mpmath_avail = False
34else:
35 mpmath_avail = True
36
37# Set available PN orders for the orbital angular momentum expression in terms of frequency
38
39LPN_orders = [0, 1, 1.5, 2, 2.5]
40
41LPN_order_string = ', '.join(map(str, LPN_orders))
42
43# Define functions
44
45
46def L_vec_from_masses_spins_and_f(m1, m2, f, PNorder=2, S1_vec=np.zeros(3), S2_vec=np.zeros(3),
47 no_mass_spin_input_check=False):
48 """
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>.
51
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.
54
55 Inputs:
56
57 - Required:
58
59 m1, m2: Detector frame masses in kg
60 f: m = 2 GW frequency in Hz
61
62 - Optional:
63
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
70
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
73 """
74
75 # Check inputs
76
77 if not no_mass_spin_input_check:
78 check_masses(m1, m2)
79
80 if f <= 0.:
81 raise ValueError(
82 "The frequency must be positive, while it is f = %f" % f)
83
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))
87
88 M = m1 + m2
89 Msq = M*M
90
91 S1 = norm(S1_vec)
92 S2 = norm(S2_vec)
93
94 if not no_mass_spin_input_check and (S1 > m1*m1/Msq or S2 > m2*m2/Msq):
95 raise ValueError(
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))
99
100 # Compute useful quantities
101
102 X1 = m1/M
103 X2 = m2/M
104
105 X1sq = X1*X1
106 X2sq = X2*X2
107
108 eta = X1*X2
109
110 zhat = np.array([0., 0., 1.])
111
112 S1z_vec = S1_vec[2]*zhat
113 S2z_vec = S2_vec[2]*zhat
114
115 v = (np.pi*M*kg_to_s*f)**(1./3.)
116
117 v2 = v*v
118 v3 = v2*v
119
120 # Put together L to the desired PN order, starting by giving the magnitude and direction of the Newtonian expression
121
122 LN = eta/v
123
124 L = 1.*zhat
125
126 if PNorder >= 1:
127 L += (1.5+eta/6.)*v2*zhat
128 if PNorder >= 1.5:
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.
131 )*v3
132 if PNorder >= 2:
133 L += (3.375 - 2.375*eta + eta*eta/24.)*v2*v2*zhat
134 if PNorder >= 2.5:
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
138 )*v3*v2/48.
139
140 return LN*L
141
142
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):
144 """
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>.
147
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.
150
152
153 Inputs:
154
155 - Required:
156
157 m1, m2: Detector frame masses in kg
158 f: m = 2 GW frequency in Hz
159
160 - Optional:
161
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
168
169 Output: Magnitude of the PN orbital angular momentum in total mass = 1 units
170 """
171
172 return norm(L_vec_from_masses_spins_and_f(m1, m2, f, PNorder=PNorder, S1_vec=S1_vec, S2_vec=S2_vec,
173 no_mass_spin_input_check=no_mass_spin_input_check))
174
175
176def L_from_masses_a_and_ecc(m1, m2, a, e):
177 """
178 Compute the magnitude of the Newtonian orbital angular momentum from the binary's masses, semimajor axis, and
179 eccentricity
180
181 Inputs:
182
183 m1, m2: Source frame masses in kg
184 a: Semimajor axis in m
185 e: eccentricity
186
187 Output: Magnitude of orbital angular momentum in total mass = 1 units
188 """
189
190 # Check inputs
191
192 check_masses(m1, m2)
193
194 if a <= 0.:
195 raise ValueError(
196 "The semimajor axis must be positive, while it is a = %e" % a)
197
198 if e < 0. or e > 1.:
199 raise ValueError(
200 "The eccentricity must be between 0 and 1, while it is e = %f" % e)
201
202 # Compute
203
204 M = m1 + m2
205
206 eta = m1*m2/M/M
207
208 return eta*(a*(1. - e*e)/(M*kg_to_m))**0.5
209
210
211def E_ratio_term(m):
212 """
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)]
214 """
215
216 omm = 1. - m
217 fpm = 4. + m
218
219 return (3.*(1. + 3.*(4. - m)/(omm*fpm))/(16.*omm**1.5) - 0.5)*m/fpm
220
221
222def eq_coeffs(u, kappaxiq, q, S1, S2, xi, S0sq):
223 """
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]
225
226 Inputs:
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)
233
234 Output: The coefficients of the equation (vareps, Btransf, Ctransf, Dtransf)
235 """
236
237 # Define useful quantities
238
239 q2 = q*q
240
241 eps = 1. - q
242
243 opq = 1. + q
244
245 u2 = u*u
246
247 S1sq = S1*S1
248 S2sq = S2*S2
249
250 kappaxiq2 = kappaxiq*kappaxiq
251
252 # Compute the coefficients of the transformed equation q (1 - q^2) u^2 bar{S}^6 + bar{B} bar{S}^4 +
253 # bar{C} bar{S}^2 + bar{D} = 0 [Eqs. (20-22) of the paper <https://dcc.ligo.org/P2100029>, arXiv:2107.11902].
254
255 # Here we refer to the barred coefficients with the suffix "transf" and the coefficient of bar{S}^6 as
256 # vareps, as in Appendix C of the paper <https://dcc.ligo.org/P2100029>, arXiv:2107.11902
257
258 Sigma = 0.5*(S0sq - S1sq - S2sq)
259 Upsilon = opq*(2.*q*Sigma + q2*S1sq + S2sq)
260 zeta = q*(Sigma + q*S1sq)*xi
261
262 vareps = q*opq*eps*u2
263
264 Btransf = 0.25*opq*eps*eps + q*eps * \
265 (xi - 2.*opq*kappaxiq)*u + Upsilon*u2
266
267 Ctransf = eps*(opq*(Sigma + q*kappaxiq2) - q*kappaxiq *
268 xi) + 2.*(zeta - Upsilon*kappaxiq)*u
269
270 Dtransf = opq*(Sigma*Sigma - S1sq*S2sq) + q2*S1sq*xi * \
271 xi/opq - 2.*zeta*kappaxiq + Upsilon*kappaxiq2
272
273 return vareps, Btransf, Ctransf, Dtransf
274
275
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):
277 """
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
281
282 Inputs:
283
284 - Required:
285
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
292
293 - Optional:
294
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
300
301 Output: The roots of the equation ordered from smallest to largest (Ssq3, Ssqm, Ssqp)
302 """
303
304 # Extract coefficients and solve
305
306 vareps, Btransf, Ctransf, Dtransf = coeffs
307
308 if use_mpmath:
309 Btransf_out, Ctransf_out, Dtransf_out = Btransf, Ctransf, Dtransf
310
311 roots, err = mp.polyroots([vareps, Btransf, Ctransf, Dtransf], error=True, maxsteps=polyroots_maxsteps, extraprec=polyroots_extraprec)
312
313 SsqI, SsqII, SsqIII = roots
314
315 # Check that the error on the roots isn't too large
316 if err > root_tol:
317 warn("polyroots error of %e greater than root_tol = %e"%(err, root_tol), RuntimeWarning)
318 else:
319 # The [0] indices here and elsewhere are necessary to avoid a warning
320 Btransf_out, Ctransf_out, Dtransf_out = Btransf[0], Ctransf[0], Dtransf[0]
321
322 cubic = P([Dtransf_out, Ctransf_out, Btransf_out, vareps])
323
324 SsqI, SsqII, SsqIII = cubic.roots()
325
326 # Check that all roots are real
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))
334
335 # Order the roots
336
337 Ssq3, Ssqm, Ssqp = np.sort([SsqI.real, SsqII.real, SsqIII.real])
338
339 # Check for complex conjugates that passed the previous check, but allow for roots that are real and equal,
340 # which occurs to numerical accuracy for close-to-aligned-spin cases
341 if Ssq3.real == Ssqm.real:
342 if Ssq3.imag != 0. or Ssqm.imag != 0.:
343 warn(
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") # This system is nonprecessing to numerical accuracy
347
348 if Ssqm.real == Ssqp.real:
349 if Ssqm.imag != 0. or Ssqp.imag != 0.:
350 warn(
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") # This system is nonprecessing to numerical accuracy
354
355 return Ssq3, Ssqm, Ssqp
356
357
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):
360 """
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).
365
366 Inputs:
367
368 - Required:
369
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
380 together)
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
386
387 - Optional:
388
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
393
394 Output: (<S^2>_prec - S0^2)/(1 - q) [= dkappa_{xi q}/du]
396 """
397
398 # Compute the desired quantity for u > 0, otherwise set to 0
399
400 if u > 0.:
401 # Compute coefficients of equation
402
403 vareps, Btransf, Ctransf, Dtransf = eq_coeffs(
404 u, kappaxiq, q, S1, S2, xi, S0sq)
405
406 # Compute the roots of the equation
407
408 # First, check whether vareps is small enough that we can linearize <S^2>_pr in mm and compute
409 # S_+^2 + S_-^2 to good accuracy by solving the quadratic one obtains by setting eps = 0 in the cubic
410 # (so S_+^2 + S_-^2 = -Ctransf/Btransf).
411
412 # Check if vareps is small enough that we can just solve the quadratic for S_+^2 + S_-^2
413 # This is currently just implemented for the case where Btransf, Ctransf, Dtransf > 0
414 # This is the case where there are issues in practice and also is the one that is easy to analyze
415
416 only_quadratic = False
417
418 # The second check is that vareps is small enough that the square root in z_eps_abs makes sense
419 if min(Btransf, Ctransf, Dtransf) > 0 and 4.*Ctransf*vareps < Btransf*Btransf:
420 z_eps_abs = abs(
421 (-Btransf + (Btransf*Btransf - 4.*Ctransf*vareps)**0.5)/(2.*Ctransf))
422
423 # Error on -Ctransf/Btransf expression for S_+^2 + S_-^2, E_\varepsilon^\text{sum} [Eq. (C5) in the paper <https://dcc.ligo.org/P2100029>, arXiv:2107.11902]
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))
427
428 # Check that the additional inequality is satisfied
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))
432
433 # This bounds S_+^2 - S_-^2 by abs(S_+^2 + S_-^2) in the first part, which is not sharp
434 # We check that m_bound < 1 since E_ratio_term is only real for m_bound < 1
435 # It should be possible to restrict computing E_ratio_term only for smaller m_bound values,
436 # but it's not clear if this is necessary
437 if m_bound < 1.:
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
441
442 if only_quadratic:
443 Ssq_pr = -0.5*Ctransf/Btransf
444
445 # Check that Ssq_pr [= (bar{S}_+^2 + bar{S}_-^2)/2] gives an unbarred value that is nonnegative
446
447 Ssq_pr_unbarred = (1. - q)*Ssq_pr + S0sq
448
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))
452 else:
453 Ssq3, Ssqm, Ssqp = solve_cubic(
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)
456
457 # Check that Ssqp and Ssqm give unbarred values that are nonnegative
458
459 Ssqp_unbarred = (1. - q)*Ssqp + S0sq
460 Ssqm_unbarred = (1. - q)*Ssqm + S0sq
461
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))
465
466 # Compute <S^2>_pr
467
468 if Ssqp == Ssq3:
469 raise RuntimeError(
470 "Degenerate nonprecessing case with all roots equal. Cannot continue with computation.")
471
472 # When abs(Ssq3) is large, mm is small, and one can write Ssq_pr to good accuracy in a form that does not
473 # involve Ssq3 (which has a large residual when put back into the equation when it is large, likely due to
474 # catastrophic cancellation), but the residual is still smaller than Ssq3, so we take the value of Ssq3
475 # output by the solver to give us an accurate enough value to use in checking the size of mm
476
477 mm = (Ssqp - Ssqm)/(Ssqp - Ssq3) # Eq. (25) in Chatziioannou, Klein, Yunes, and Cornish PRD 95, 104004 (2017)
478
479 # Check that the equation is satisfied to good accuracy (by checking that the residual changes sign over an
480 # interval given by root_tol) and compute Ssq_pr
481
482 def cubic(Ssq):
483 if use_mpmath:
484 return mp.polyval([vareps, Btransf, Ctransf, Dtransf], Ssq)
485 else:
486 return ((vareps*Ssq + Btransf)*Ssq + Ctransf)*Ssq + Dtransf
487
488 def cubic_deriv(Ssq):
489 if use_mpmath:
490 return mp.polyval([3*vareps, 2*Btransf, Ctransf], Ssq)
491 else:
492 return (3.*vareps*Ssq + 2.*Btransf)*Ssq + Ctransf
493
494 # If the difference between Ssqp and Ssqm or Ssqm and Ssq3 is smaller than the difference set by root_tol,
495 # shrink root_tol appropriately
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.)
498 else:
499 root_tol_pm = root_tol
500
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.)
503 else:
504 root_tol_m3 = root_tol
505
506 root_tol_message_extra_bit_flag = False
507
508 root_tol_orig = root_tol
509
510 if root_tol_pm < root_tol and root_tol_pm <= root_tol_m3:
511 root_tol = root_tol_pm
512
513 root_tol_message_extra_bit_flag = True
514
515 root_tol_message_extra_bit_insert = "S_+^2 - S_-^2"
516 elif root_tol_m3 < root_tol:
517 root_tol = root_tol_m3
518
519 root_tol_message_extra_bit_flag = True
520
521 root_tol_message_extra_bit_insert = "S_-^2 - S_3^2"
522
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
526 else:
527 root_tol_message_extra_bit = ""
528
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)
533
534 # Check for special cases and if we can linearize
535
536 if use_mpmath:
537 Btransf_out, Ctransf_out, Dtransf_out = Btransf, Ctransf, Dtransf
538 else:
539 Btransf_out, Ctransf_out, Dtransf_out = Btransf[0], Ctransf[0], Dtransf[0]
540
541 if Ssqp == Ssqm:
542 # Since here we have a double root, we will check that the derivative of the cubic changes sign
543 # If the difference between Ssqm and Ssq3 is smaller than the difference set by root_tol,
544 # shrink root_tol appropriately
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.)
547
548 root_tol_message_extra_bit = " the internally modified (since S_-^2 - S_3^2 is small) "\
549 "version of"
550 else:
551 root_tol_m0 = root_tol_orig
552
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)
556
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)
566
567 if root_tol_raise_err:
568 raise RuntimeError(root_tol_message)
569 else:
570 warn(root_tol_message, RuntimeWarning)
571
572 Ssq_pr = Ssqp # Here m = 0
573 elif Ssqm == Ssq3:
574 # Since here we have a double root, we will check that the derivative of the cubic changes sign
575 # If the difference between Ssqp and Ssqm is smaller than the difference set by root_tol,
576 # shrink root_tol appropriately
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.)
579
580 root_tol_message_extra_bit = " the internally modified (since S_+^2 - S_-^2 is small) "\
581 "version of"
582 else:
583 root_tol_m1 = root_tol_orig
584
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)
588
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)
598
599 if root_tol_raise_err:
600 raise RuntimeError(root_tol_message)
601 else:
602 warn(root_tol_message, RuntimeWarning)
603
604 Ssq_pr = Ssq3 # Here m = 1
605 elif E_ratio_term(mm) <= lin_tol*max(1., 1./(Ssqp - Ssqm)): # Linearization check
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)
616
617 if root_tol_raise_err:
618 raise RuntimeError(root_tol_message)
619 else:
620 warn(root_tol_message, RuntimeWarning)
621
622 # Eq. (42) in Chatziioannou, Klein, Yunes, and Cornish PRD 95, 104004 (2017), with the elliptic functions linearized in m
623 Ssq_pr = 0.5*(Ssqp + Ssqm)
624 else:
625 Ssq3_resp, Ssq3_resm = cubic(
626 Ssq3 + max(abs(Ssq3), 1.)*root_tol), cubic(Ssq3 - max(abs(Ssq3), 1.)*root_tol)
627
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,
640 Dtransf_out)
641
642 if root_tol_raise_err:
643 raise RuntimeError(root_tol_message)
644 else:
645 warn(root_tol_message, RuntimeWarning)
646
647 if use_mpmath:
648 ellip_ratio = mp.ellipe(mm)/mp.ellipk(mm)
649 else:
650 ellip_ratio = ellipe(mm)/ellipk(mm)
651
652 # Eq. (17) in the paper <https://dcc.ligo.org/P2100029>, arXiv:2107.11902
653 Ssq_pr = Ssqp + (ellip_ratio - 1.)*(Ssqp - Ssq3)
654 else:
655 Ssq_pr = 0.
656
657 return Ssq_pr
658
659
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):
663 """
664 Compute kappaxiq when the binary has u = uf starting from u = u0
665
666 Inputs:
667
668 - Required:
669
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
676 uf: Final 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
684 together)
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
689 if it is used
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
693
694 - Optional:
695
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.
700
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
708 """
709
710 # Compute initial value of kappaxiq
711
712 kappaxiq0 = S1L
713
714 if use_solve_ivp:
715 # Use scipy.integrate.solve_ivp
716
717 # Set up kwargs for solver (necessary since one gets a warning if one passes min_step with an integrator other than LSODA)
718 solve_ivp_kwargs = {'atol': ode_atol, 'rtol': ode_rtol}
719
720 if method == "LSODA":
721 solve_ivp_kwargs['min_step'] = solve_ivp_lsoda_min_step
722
723 # Solve
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)
726
727 # Return kappaxiq at u = uf
728
729 return soln.y[0]
730 elif use_mpmath:
731 # Use mpmath.odefun
732
733 # Solve (we have to write it in terms of ub := -u so that the final value of ub is greater than the initial value, as
734 # required by mpmath.odefun)
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)
740
741 # Return kappaxiq at u = uf
742
743 return soln(-uf)
744 else:
745 # Use scipy.integrate.ode
746
747 # initialize the ODE solver
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)
753
754 # Solve
755
756 end = False
757
758 while solver.successful() and not end:
759 new_t = solver.t - du
760
761 if new_t <= uf:
762 new_t = uf
763 end = True
764
765 solver.integrate(new_t, step=1)
766
767 if abs(solver.t - uf) > 0.75*du:
768 raise RuntimeError(
769 "Final u of %e is more than 0.75*du different from uf of %e" % (solver.t, uf))
770
771 # Return kappaxiq at u = uf
772
773 return solver.y
774
775
776def 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,
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'):
784 """
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).
792
793 Inputs:
794
795 - Required:
796
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
801 coordinate system
802 fref: Reference frequency, in Hz; ignored when L0_vec is not None
803
804 - Optional initial and final L settings:
805
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,
809 default: None
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
812
813 - Optional settings for the evolution:
814
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),
838 default: "lsoda"
839 use_mpmath_fallback: Use the slow but accurate mpmath integration if the first fallback integration fails (True) or not (False),
840 default: True
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
853 not None, default: 2
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"
873
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)
877 """
878
879 # Set scipy integrators for which this has been tested
880
881 if use_solve_ivp:
882 tested_scipy_integrators = ["RK45", "RK23", "DOP853", "LSODA"]
883 else:
884 tested_scipy_integrators = ["vode", "lsoda", "dopri5", "dop853"]
885
886 tested_scipy_integrator_string = ', '.join(tested_scipy_integrators)
887
888 # Set known failure modes
889
890 known_failure_modes = ['Error', 'NAN', 'None']
891
892 # Check version
893
894 if version != 'v1':
895 raise ValueError("Only version v1 is available at the moment, while version = %s"%version)
896
897 # Make "lsoda" and "dop853" uppercase when passed with use_solve_ivp == True, for compatibility
898
899 methods_to_uppercase = ["lsoda", "dop853"]
900
901 if use_solve_ivp:
902 if method in methods_to_uppercase:
903 method = method.upper()
904
905 if method_fallback in methods_to_uppercase:
906 method_fallback = method_fallback.upper()
907
908 # Check that inputs are physical and otherwise make sense
909
910 check_masses(m1, m2)
911
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:
914 raise ValueError(
915 "chi1_vec and chi2_vec must both be numpy arrays of length 3, while they are", chi1_vec, chi2_vec)
916
917 if L0_vec is not None:
918 if len(L0_vec) != 3:
919 raise ValueError(
920 "L0_vec must be either None or a numpy array of length 3, while it is", L0_vec)
921
922 chi1 = norm(chi1_vec)
923 chi2 = norm(chi2_vec)
924
925 if chi1 == 0. or chi1 > 1. or chi2 == 0. or chi2 > 1.:
926 raise ValueError(
927 "The magnitudes of the spins must both be positive and at most 1, while they are chi1, chi2 = %f, %f" %
928 (chi1, chi2))
929
930 if fref is not None:
931 check_fref(fref, m1, m2, "precession-averaged")
932 elif L0_vec is None:
933 raise ValueError("One of fref and L0_vec must not be None")
934
935 if du <= 0. or du_fallback <= 0.:
936 raise ValueError(
937 "The step in u and its fallback value must both be positive, while they are du, du_fallback = %e, %e" %
938 (du, du_fallback))
939
940 if use_solve_ivp and not solve_ivp_avail:
941 import scipy
942
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__)
945
946 if use_solve_ivp:
947 interface_str = "solve_ivp"
948 else:
949 interface_str = "ode"
950
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))
956
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)
960
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)
964
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)
967
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))
973
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))
977
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))
982
983 # Print a warning if both fref and L0_vec are not None
984
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)
988
989 # Set the failure output and the string corresponding to it
990 # These default to None, since they are not used if failure_mode == 'Error'
991
992 if failure_mode == 'NAN':
993 failure_output = np.nan
994 failure_output_string = 'np.nan'
995 else:
996 failure_output = None
997 failure_output_string = 'None'
998
999 # Set failure message for mpmath evolution
1000
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)
1007
1008 # Check that m1 > m2, and swap the objects if m1 < m2; raise an error if m1 == m2
1009 # This does not rotate the spins, since we do not need to track the orientation of the binary
1010
1011 swap = False
1012
1013 if m1 < m2:
1014 swap = True
1015 m1, m2 = m2, m1
1016 chi1, chi2 = chi2, chi1
1017 chi1_vec, chi2_vec = chi2_vec, chi1_vec
1018
1019 eq_mass_check(m1, m2, Lf)
1020
1021 # Save the masses in kg for later use
1022
1023 m1_kg = m1
1024 m2_kg = m2
1025
1026 # Compute derived quantities (scaling the masses so that their sum is 1)
1027
1028 q = m2/m1
1029
1030 M = m1 + m2
1031
1032 m1 /= M
1033 m2 /= M
1034
1035 S1 = m1*m1*chi1
1036 S2 = m2*m2*chi2
1037
1038 S1_vec = m1*m1*chi1_vec
1039 S2_vec = m2*m2*chi2_vec
1040
1041 # Compute initial orbital angular momentum (in total mass = 1 units)
1042 # We need to pass the masses in kg here to convert the frequency when computing from fref
1043
1044 if L0_vec is not None:
1045 L0 = norm(L0_vec)
1046 else:
1047 if LPNspins:
1048 S1_vec_for_L = S1_vec
1049 S2_vec_for_L = S2_vec
1050 else:
1051 S1_vec_for_L = np.zeros(3)
1052 S2_vec_for_L = np.zeros(3)
1053
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)
1056
1057 L0 = norm(L0_vec)
1058
1059 if L0 == 0.:
1060 raise ValueError("L0_vec must not be the zero vector")
1061
1062 if Lf is not None and Lf <= L0:
1063 raise ValueError(
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))
1066
1067 # Setting ode_atol for the case when ode_atol_inplane_spin_scaling == True
1068
1069 L0_hat = L0_vec/L0
1070
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)
1073
1074 if ode_atol_inplane_spin_scaling:
1075 ode_atol = max(min(chi1_inplane, chi2_inplane)
1076 * ode_atol_base, ode_atol_floor)
1077 else:
1078 ode_atol = ode_atol_base
1079
1080 # Take care of the -1 cases for ode_rtol and lin_tol and their fallback versions
1081
1082 if ode_rtol == -1:
1083 ode_rtol = 0.1*ode_atol
1084
1085 if ode_rtol_fallback == -1:
1086 ode_rtol_fallback = 0.1*ode_atol_fallback
1087
1088 if lin_tol == -1:
1089 lin_tol = ode_atol
1090
1091 if lin_tol_fallback == -1:
1092 lin_tol_fallback = ode_atol_fallback
1093
1094 # Check the tolerance settings
1095
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.:
1098 raise ValueError(
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))
1103
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))
1108
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))
1112
1113 # Compute effective spin [Eq. (12) of Gerosa et al.]
1114
1115 xi = np.dot((1. + q)*S1_vec + (1. + 1./q)*S2_vec, L0_hat)
1116
1117 # Compute projection of S1_vec onto L0_vec
1118
1119 S1L = np.dot(S1_vec, L0_hat)
1120 S2L = np.dot(S2_vec, L0_hat)
1121
1122 # Return the initial tilts when spins are aligned and not in an unstable configuration
1123
1124 cos_tilt1 = np.clip(S1L/S1, -1., 1.)
1125 cos_tilt2 = np.clip(S2L/S2, -1., 1.)
1126
1127 tilt1 = np.arccos(cos_tilt1)
1128 tilt2 = np.arccos(cos_tilt2)
1129
1130 if (tilt1 == np.pi and tilt2 in [0., np.pi]) or (tilt1 == 0. and tilt2 == 0.):
1131 return package_tilts(tilt1, tilt2, Lf, swap=swap)
1132
1133 # Compute initial total dimensionful spin and angular momentum in total mass = 1 units
1134
1135 S0sq = norm(S1_vec + S2_vec)**2
1136
1137 # Compute kappa at infinity or Lf
1138
1139 u0 = 0.5/L0
1140
1141 if Lf is not None:
1142 uf = 0.5/Lf
1143 else:
1144 uf = 0.
1145
1146 mpmath_tilts = False
1147
1148 if method == "odefun":
1149 try:
1150 mpmath_tilts = True
1151
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)
1157
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)
1160
1161 except NonprecessingError as ne:
1162 return package_tilts(tilt1, tilt2, Lf, swap)
1163 else:
1164 try:
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)
1168
1169 except RuntimeError as re:
1170 if use_fallback:
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)
1175
1176 try:
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)
1180
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)
1187
1188 try:
1189 mpmath_tilts = True
1190
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)
1196
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)
1199
1200 except NonprecessingError as ne:
1201 return package_tilts(tilt1, tilt2, Lf, swap)
1202 else:
1203 if not use_mpmath_fallback:
1204 mpmath_bit = "it is disabled"
1205 else:
1206 mpmath_bit = "mpmath is not available"
1207
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)
1211
1212 return evolution_error_handling(failure_mode, failure_message, failure_output, failure_output_string, Lf, swap)
1213
1214 except NonprecessingError as ne:
1215 return package_tilts(tilt1, tilt2, Lf, swap)
1216
1217 else:
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)
1221
1222 return evolution_error_handling(failure_mode, failure_message, failure_output, failure_output_string, Lf, swap)
1223
1224 except NonprecessingError as ne:
1225 return package_tilts(tilt1, tilt2, Lf, swap)
1226
1227 # Compute the desired output. We have already accounted for the case where at least one spin is zero,
1228 # so there are no concerns about division by zero here.
1229
1230 if Lf is None:
1231 # Compute cosines of tilt angles at infinity [Eqs. (15) in the paper <https://dcc.ligo.org/P2100029>, arXiv:2107.11902]
1232
1233 cos_tilt1inf = k_evol/S1
1234 cos_tilt2inf = q*(xi/(1. + q) - k_evol)/S2
1235
1236 # Clip to avoid issues with taking the arccos if there are rounding issues
1237
1238 cos_tilt1inf = np.clip(cos_tilt1inf, -1., 1.)
1239 cos_tilt2inf = np.clip(cos_tilt2inf, -1., 1.)
1240
1241 # Return tilt angles at infinity
1242
1243 if mpmath_tilts:
1244 tilt1inf = float(mp.acos(cos_tilt1inf))
1245 tilt2inf = float(mp.acos(cos_tilt2inf))
1246 else:
1247 tilt1inf = np.arccos(cos_tilt1inf[0])
1248 tilt2inf = np.arccos(cos_tilt2inf[0])
1249
1250 return package_tilts(tilt1inf, tilt2inf, None, swap)
1251 else:
1252 vareps, Btransf, Ctransf, Dtransf = eq_coeffs(
1253 uf, k_evol, q, S1, S2, xi, S0sq)
1254
1255 _, Ssqm, Ssqp = solve_cubic(
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)
1258
1259 if mpmath_tilts:
1260 lin_tol_use = ode_tol_mpmath
1261 else:
1262 lin_tol_use = lin_tol
1263
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)
1267
1268 # Compute cosines of tilt angles at separation [Eqs. (23) in paper <https://dcc.ligo.org/P2100029>, arXiv:2107.11902]
1269
1270 cos_tilt1sep_Sm = (k_evol - uf*Ssqm)/S1
1271 cos_tilt2sep_Sm = q*(xi/(1. + q) - S1*cos_tilt1sep_Sm)/S2
1272
1273 cos_tilt1sep_Sp = (k_evol - uf*Ssqp)/S1
1274 cos_tilt2sep_Sp = q*(xi/(1. + q) - S1*cos_tilt1sep_Sp)/S2
1275
1276 cos_tilt1sep_Savg = (k_evol - uf*Ssqavg)/S1
1277 cos_tilt2sep_Savg = q*(xi/(1. + q) - S1*cos_tilt1sep_Savg)/S2
1278
1279 # Clip to avoid issues with taking the arccos if there are rounding issues
1280
1281 cos_tilt1sep_Sm = np.clip(cos_tilt1sep_Sm, -1., 1.)
1282 cos_tilt2sep_Sp = np.clip(cos_tilt2sep_Sp, -1., 1.)
1283
1284 cos_tilt1sep_Sm = np.clip(cos_tilt1sep_Sm, -1., 1.)
1285 cos_tilt2sep_Sp = np.clip(cos_tilt2sep_Sp, -1., 1.)
1286
1287 cos_tilt1sep_Savg = np.clip(cos_tilt1sep_Savg, -1., 1.)
1288 cos_tilt2sep_Savg = np.clip(cos_tilt2sep_Savg, -1., 1.)
1289
1290 # Return tilt angles at separation
1291
1292 if swap:
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
1296
1297 if mpmath_tilts:
1298 tilt1sepmin = float(mp.acos(cos_tilt1sep_Sm))
1299 tilt1sepmax = float(mp.acos(cos_tilt1sep_Sp))
1300 tilt1sepavg = float(mp.acos(cos_tilt1sep_Savg))
1301
1302 tilt2sepmin = float(mp.acos(cos_tilt2sep_Sp))
1303 tilt2sepmax = float(mp.acos(cos_tilt2sep_Sm))
1304 tilt2sepavg = float(mp.acos(cos_tilt2sep_Savg))
1305 else:
1306 tilt1sepmin = np.arccos(cos_tilt1sep_Sm[0])
1307 tilt1sepmax = np.arccos(cos_tilt1sep_Sp[0])
1308 tilt1sepavg = np.arccos(cos_tilt1sep_Savg[0])
1309
1310 tilt2sepmin = np.arccos(cos_tilt2sep_Sp[0])
1311 tilt2sepmax = np.arccos(cos_tilt2sep_Sm[0])
1312 tilt2sepavg = np.arccos(cos_tilt2sep_Savg[0])
1313
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}
1316
1317
1318def prec_avg_tilt_comp(m1, m2, chi1, chi2, tilt1, tilt2, phi12, fref, Lf=None, **kwargs):
1319 """
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.
1325
1326 Inputs:
1327
1328 - Required:
1329
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
1335
1336 - Optional settings [passed to prec_avg_tilt_comp_vec_inputs(); all but Lf passed through kwargs]:
1337
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
1340
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),
1364 default: "lsoda"
1365 use_mpmath_fallback: Use the slow but accurate mpmath integration if the first fallback integration fails (True) or not (False),
1366 default: True
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"
1398
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)
1402 """
1403
1404 # Return initial tilts for single-spin or nonspinning cases
1405
1406 if chi1 == 0. or chi2 == 0.:
1407 return package_tilts(tilt1, tilt2, Lf, False)
1408
1409 # Check that spins and tilt angles are physical/make sense [other checks carried out in prec_avg_tilt_comp_vec_inputs()]
1410
1411 check_spin_mags(chi1, chi2)
1412
1413 check_tilts(tilt1, tilt2)
1414
1415 # Convert spin magnitudes and angles to vectors
1416
1417 chi1_vec = chi1*np.array([np.sin(tilt1), 0., np.cos(tilt1)])
1418
1419 chi2_vec = chi2*np.array([np.sin(tilt2)*np.cos(phi12),
1420 np.sin(tilt2)*np.sin(phi12), np.cos(tilt2)])
1421
1422 # Pass to prec_avg_tilt_comp_vec_inputs() to do calculation
1423
1424 return prec_avg_tilt_comp_vec_inputs(m1, m2, chi1_vec, chi2_vec, fref, Lf=Lf, **kwargs)
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 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 ...