phenomxpy.precession.nnlo
Next-to-next to leasing order Euler angles
- class phenomxpy.precession.nnlo.NNLO[source]
Bases:
objectClass with extra methods for the NNLO angles to be added to
pPrec.
phenomxpy.precession.msa
Multi-Scale-Analyis Euler angles
- class phenomxpy.precession.msa.MSA[source]
Bases:
objectClass with extra methods for the MSA angles to be added to
pPrec.Only
PrecVersion= 223 version is supported.The methods are adapted from PhenomX lalsuite.
- Constants_c_MSA(v, v2, JNorm, Spl2, Smi2)[source]
Get c constants from Appendix B (B6, B7, B8) [8].
- Constants_d_MSA(LNorm, JNorm, Spl, Spl2, Smi2)[source]
Get d constants from Appendix B (B9, B10, B11) [8].
- Initialize_MSA(ExpansionOrder=5)[source]
Initialize all the core variables for the MSA system. This will be called first.
- MSA_Corrections(v, v2, LNorm, JNorm, Spl, Spl2, Smi2, S32)[source]
Get MSA corrections to \(\zeta\) and \(\phi_z\) using Eq. F19 in Appendix F and Eq. 67 [8]. respectively.
- Psi_MSA(v, v2)[source]
Get \(\psi\) using Eq 51 [8]. - Here \(\psi\) is the phase of S as in Eq 23. - Note that the coefficients are defined in Appendix C (C1 and C2).
- Psi_dot_MSA(v, v2, Spl2, S32)[source]
Get \(\dot{\psi}\) using Eq 24 [8]. \(\frac{d \psi}{d t} = \frac{A}{2} \sqrt{S_+^2 - S_3^2}\).
- Roots_MSA(LNorm, JNorm)[source]
Solve for the roots of Eq 21 [8].
Roots for \(\frac{d S^2}{d t^2} = -A^2 (S^2 -S_+^2)(S^2 - S_-^2)(S^2 - S_3^2)\).
- Returns:
\(S_+^2\), \(S_-^2\) and \(S_3^2\).
- Return type:
Tuple with 3 floats or 3 1D ndarrays
- compute_msa_angles(omega, **kwargs)[source]
Compute MSA angles alpha, cosbeta, gamma.
- Parameters:
omega (float, 1D ndarray) – frequency evolution taken from the aligned-spin model.
- Returns:
\(\alpha\), \(\cos \beta\), \(\gamma\)
- Return type:
Tuple with 3 1D ndarrays
phenomxpy.precession.numerical
Spin-Taylor numerical Euler angles.
- class phenomxpy.precession.numerical.Numerical[source]
Bases:
objectClass with extra methods for evolving spin equations and compute numerical angles added to
pPrec.- compute_numerical_angles(LNhat, times)[source]
Compute numerical angles alpha, beta, gamma for a given LN evolution plus ringdown attachment.
The input LNhat is the evolved, normalized one. It is transformed to the J-frame with the angles phiJ_Sf, thetaJ_Sf, kappa.
Then the Euler angles from J to L are:
\[\alpha = \operatorname{atan2}(L_{N,y}, L_{N,x}) \cos(\beta) = L_{N,z}.\]\(\gamma\) is numerically integrated according to the minimal rotation condition.
- if
self.interpolate!= “euler_angles”: times corresponds to the points where LNhat is evaluated plus the points in the ringdown.
- if
self.interpolate== “euler_angles”: LNhat is evaluated in
self.t_eval, times can be any array and the Euler angles are evaluated with CubicSplines.
Sets the angle offsets from J to L0: alphaJtoL0, cosbetaJtoL0, gammaJtoL0
- Returns:
alpha, cosbeta, gamma
- Return type:
Tuple with 3 1D ndarrays
- if
- evolve_orbit(tinit, tend, t_eval=None)[source]
Solve PN spin equations from tinit to tend in a custom time array
t_eval.LNh is the Z axis of the rotating frame, while E1 is the X axis.
The evolution of E1 is not needed for computing the Euler angles. If removed, the solution is slightly different since the solver chooses different adaptive points and interpolation methods. E1 is evolved satisfying the minimal rotation condition. This is used to obtain the quaternions without computing the Euler angles nor integrating gamma. It also prevents the loss of precission due to building CubicSplines, derive them and then numerically integrate them.
- Returns:
LNh, E1, S1, S2, and the evaluated times t_eval.
- Return type:
Tuple with 4 (3, N)-ndarrays and one 1D ndarray
- evolved_final_spin(LNhat, S1, S2)[source]
Compute the final spin from evolved quantities. Eq. 25 arxiv:2105.05872
The aligned final spin contribution is also computed from the evolved quantities and its sign is copied to the precessing one
Update self._pWF.afinal_prec and self.EulerRDslope
- gamma_integration(times, alpha, cosbeta)[source]
Numerical integration of \(\dot{\gamma}\) over the full time array to fulfill minimal rotation condition.
\(\dot{\gamma} = - \dot{\alpha} \cos(\beta)\). Eq. 18 [9].
Different methods selected with option
pPrec.gamma_integration_method:‘piecewise_integral’: analytical integral in each time interval of the 5th order “spline” resulting from the multiplication of two cubic splines: CS.derivative() * CS.
‘boole’: use Boole’s rule, assumes equispaced grid. This is the only one supported for cuda=True and in this case it uses linear interpolation instead of cubic spline.
‘antiderivative’: use .antiderivative() method of CubicSpline. This spline is built after evaluating the alpha.derivative() and cosbeta() CubicSplines. It loses accuracy since higher order from the multiplication of alpha_der * cosbeta are discarded.
- Returns:
Integrated \(\gamma\) satisfying minimal rotation condition.
- Return type:
1D ndarray
- class phenomxpy.precession.numerical.SpinTaylor(eta, m1M, m2M)[source]
Bases:
objectClass to setup SpinTaylorT4 quantities and compute spin derivatives.
- energy_spin_derivative_setup(quadparam1=0, quadparam2=0)[source]
Setup energy and spin derivative coefficients.
Adapted from XLALSimSpinTaylorEnergySpinDerivativeSetup.
- phenomxpy.precession.numerical._spin_derivatives_avg(v, LNh, E1, S1, S2, eta, S1dot3, S2dot3, S1dot4S2Avg, S1dot4S2OAvg, S1dot5, S2dot5, S1dot7S2, S2dot7S1)[source]
Compute right-hand side of ODE.
We call it spin derivatives but it also includes the derivatives of LNh and E1.
Special case of the original function in lalsuite.
- Parameters:
v (float) – \(\sqrt[3]{\omega / 2}\)
LNh (ndarray (3,)) – The angular momentum unit vector.
E1 (ndarray (3,)) – The X axis vector.
S1 (ndarray (3,)) – The spin vector of the first object.
S2 (ndarray (3,)) – The spin vector of the second object.
S1/2dot... (float) – Spin coefficients.
- Returns:
A 1D array with the derivatives of LNh, E1, S1, S2. 12 elements since each vector has 3 components.
- Return type:
1D ndarray
- phenomxpy.precession.numerical.eval_cubic_spline(x, x_breaks, coefs)[source]
Evaluate a cubic spline in a single point with numba.
- phenomxpy.precession.numerical.eval_cubic_spline_array(x_array, x_breaks, coefs)[source]
Evaluate a cubic spline in an array with numba.
- phenomxpy.precession.numerical.integrate_product_cubic_splines_vectorized(cs1, cs2, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/phenomxpy/envs/stable/lib/python3.13/site-packages/numpy/__init__.py'>)[source]
Vectorized integration of the product of two splines.
- Parameters:
cs1 (CubicSpline.derivative) – Derivative of the first cubic spline with CubicSpline.derivative()
cs2 (CubicSpline) – The second cubic spline
- Returns:
Indefinite integral of the product.
- Return type:
1D ndarray
- phenomxpy.precession.numerical.spin_derivatives_avg_numba(v, LNh, E1, S1, S2, eta, S1dot3, S2dot3, S1dot4S2Avg, S1dot4S2OAvg, S1dot5, S2dot5, S1dot7S2, S2dot7S1)
Compute right-hand side of ODE.
We call it spin derivatives but it also includes the derivatives of LNh and E1.
Special case of the original function in lalsuite.
- Parameters:
v (float) – \(\sqrt[3]{\omega / 2}\)
LNh (ndarray (3,)) – The angular momentum unit vector.
E1 (ndarray (3,)) – The X axis vector.
S1 (ndarray (3,)) – The spin vector of the first object.
S2 (ndarray (3,)) – The spin vector of the second object.
S1/2dot... (float) – Spin coefficients.
- Returns:
A 1D array with the derivatives of LNh, E1, S1, S2. 12 elements since each vector has 3 components.
- Return type:
1D ndarray
- phenomxpy.precession.numerical.tphm_spin_derivatives(t, y, v_interpolant, spin_taylor_t4)[source]
Spin derivatives using a v_interpolant without numba.
- phenomxpy.precession.numerical.tphm_spin_derivatives_imr_omega(t, y, pPhase22, spin_taylor_t4)[source]
Spin derivatives using imr_omega without numba.
- phenomxpy.precession.numerical.tphm_spin_derivatives_imr_omega_numba(t, y, pPhase22, eta, S1dot3, S2dot3, S1dot4S2Avg, S1dot4S2OAvg, S1dot5, S2dot5, S1dot7S2, S2dot7S1)[source]
Spin derivatives using imr_omega with numba.
- phenomxpy.precession.numerical.tphm_spin_derivatives_numba(t, y, t_breaks, coeffs, eta, S1dot3, S2dot3, S1dot4S2Avg, S1dot4S2OAvg, S1dot5, S2dot5, S1dot7S2, S2dot7S1)[source]
Spin derivatives using a v_interpolant with numba.
phenomxpy.precession.precession
Define pPrec class and useful function for precession
- phenomxpy.precession.precession.numba_global_rotation(lmax, L0modes, Jmodes, wignerdL2, wignerdL3, wignerdL4, wignerdL5, expAlphaPowers, expGammaPowers)[source]
- phenomxpy.precession.precession.numba_rotation(lmax, Jmodes, CPmodes, wignerdL2, wignerdL3, wignerdL4, wignerdL5, expAlphaPowers, expGammaPowers, mode_array)[source]
- class phenomxpy.precession.precession.pPrec(pwf, prec_version='numerical', final_spin_version=None, afinal_prec=None, add_20_mode=False, numba_rotation='default_modes', store_arrays=False, **kwargs)[source]
Bases:
objectClass for precessing utilities.
- Parameters:
pwf (pWF) – pWF object with waveform arguments
prec_version ({'numerical', 'msa', 'nnlo'}) – Choose the precessing prescription.
final_spin_version (int) –
Options for computing the final spin (0, 1, 2, 3, 4).
When
prec_versionis ‘numerical’,final_spin_versionuses 4. ‘msa’ uses 3. ‘nnlo’ uses 0.0: \(\bar{\chi}_p = \chi_p\) (Eq. 4.7-4.9 [7]).
1: \(\bar{\chi}_p = \chi_{1x}\)
2: \(\bar{\chi}_p = \chi_\perp\) (related to in-plane spin components, Eq. 4.24, 4.28 [7]).
For the options above, \(\bar{\chi}_p\) is then plugged into Eq. 4.25 [7] to obtain the final spin.
3: For MSA angles, uses spin averaged quantities (Eq. 4.23, 4.29 [7]).
4: Numerically evolved (this is done outside the pPrec class).
afinal_prec (float) – Numerical value for the final precessing spin. Skips its internal calculation if not
None.add_20_mode (bool) – Include 20 mode in co-precessing modes
numba_rotation ({'default_modes', 'custom_modes'} or bool) –
Use numba functions to perform the twisting-up rotation of modes.
’default_modes’ or
True: use numba for the default_mode_array defined in phenomt.py’custom_modes’: use custom mode_array
False: do not use numba.
store_arrays (bool) – Store internal arrays like evolved LN, S1/2 for numerical angles or Euler angles.
fall_back_msa_to_nnlo (bool) – Fall back to NNLO angles if MSA initialization fails.
use_final_spin_sign (bool) – Use sign of final spin from MSA or assume positive.
rtol (float) – Relative and absolute tolerances for the ODE solver.
atol (float) – Relative and absolute tolerances for the ODE solver.
use_closest_time_to_tref (bool) – Use closest point in the time array to tref to define there the initial condition for the ODE solver. It only supports equispaced time arrays.
v_function ({'interpolant', 'imr_omega'}) –
Function for v used in the ODE solver.
’interpolant’: compute an interpolant of v.
’imr_omega’: use the pPhase.imr_omega so it is independent of the time array as the interpolant and more accurate.
cubic_interpolation_for_ode (bool) – The ODE solution is interpolated using CubicSpline and evaluated on the cpu, then transferred to the gpu if needed. If False, the interpolation is linear and is done directly on the gpu if requested.
numba_derivatives (bool) – Use numba functions for the right-hand-side of the ODE solver.
interpolate ({‘ode_solution’, ‘euler_angles’} or
None.) –‘ode_solution’: the ODE solution is interpolated from the internal solver points to the requested time array
’euler_angles’: the Euler angles are interpolated from the internal solver points to the requested time array
None: the ODE solution is evaluated in the requested time array
quaternions_from ({'frame', 'euler_angles'}) –
Strategy to compute the quaternions transforming from co-precessing to inertial frame.
’frame’: use the evolved frame Z=LNhatev, X=E1, Y=ZxX to get the quaternions already satisfying the minimal rotation condition. The evolution of E1 is done so it satisfies minimal rotation condition. Avoid computational cost of Euler angles and loss of accuracy due to the use of cubic splines for gamma integration.
’euler_angles’: from the evolution of LNhatev, computes Euler angles alpha, cosbeta, numerically integrate gamma according to minimal rotation condition and then transform to quaternions.
gamma_integration_method ({'piecewise_integral', 'boole', 'antiderivative'}) –
Strategy for the integration of gamma following minimal rotation condition.
’piecewise_integral’: analytical integral in each time interval of the 5th order “spline” resulting from the multiplication of two cubic splines: CS.derivative() * CS.
’boole’: use Boole’s rule, assumes equispaced grid. This is the only one supported for cuda=True and in this case it uses linear interpolation instead of cubic spline.
’antiderivative’: use .antiderivative() method of CubicSpline. This spline is built after evaluating the alpha.derivative() and cosbeta() CubicSplines. It loses accuracy since higher order from the multiplication of alpha_der * cosbeta are discarded.
analytical_RD_gamma (bool) – The ringdown attachment for gamma angles has an analytical form in terms of alpha and beta, which is used for NNLO and MSA but that was not used for the numerical angles in the lalsuite implementation.
- WignerdMatrix(cosBeta, sign=1, lmax=2, global_rotation=False)[source]
Compute Wigner-d matrices coefficients from cosbeta.
Appendix A [7].
Used only for L0/Jmodes, not the default CPmodes.
Coefficients for each l are stored in self.wignerdLi with i = {2, 3, 4, 5}.
global_rotation=True refers to the constant in time rotation from J to L0 frame. This rotation includes all the l=l modes
- _load_from(source_class)[source]
Add specific methods from an input class (NNLO, MSA, Numerical) to the pPrec class
- compute_angular_momentum(v)[source]
Compute orbital angular momentum up to 4PN approximantion.
Eq. 4.13 [7].
- set_angular_momentum_coefficients()[source]
Set coefficients for the PN angular momenum L.
Eqs. G10 [7].
Angular momemtum at fref. Eq. 4.13 arxiv:2004.06503
- set_exp_angle_powers(alpha, gamma)[source]
Set complex exponential values needed for Wigner-D matrices.
- set_final_spin(final_spin_version)[source]
Set final spin according to the final_spin_version
0: \(\bar{\chi}_p = \chi_p\) (Eq. 4.7-4.9 [7]).
1: \(\bar{\chi}_p = \chi_{1x}\)
2: \(\bar{\chi}_p = \chi_\perp\) (related to in-plane spin components, Eq. 4.24, 4.28 [7]).
For the options above, \(\bar{\chi}_p\) is then plugged into Eq. 4.25 [7] to obtain the final spin.
3: For MSA angles, uses spin averaged quantities (Eq. 4.23, 4.29 [7]).
4: Numerically evolved (this is done outside the pPrec class).
phenomxpy.precession.quaternion
Utility functions involving quaternions.
- phenomxpy.precession.quaternion.as_euler_angles(q, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/phenomxpy/envs/stable/lib/python3.13/site-packages/numpy/__init__.py'>)[source]
Get Euler angles from normalized quaternions.
Adapted from M. Boyle to return cosbeta
Returns alpha, cosbeta, gamma
- phenomxpy.precession.quaternion.compute_ylms(qTot, mode_array, use_wigner_from_quaternions, numba_rotation, xp, lmax)[source]
Compute Ylms time series from quaternions for a list of modes.
Used for
computing polarizations_from= ‘CPmodes’.- Parameters:
qTot ((N, 4) ndarray) – Array of quaternions (as float array).
mode_array (list) – List of modes in format [[l,m],[],…] to compute the Ylms for.
use_wigner_from_quaternions (bool) – Compute WignerD matrices (Ylms) directly from quaternions instead of from Euler angles.
numba_rotation (string or bool) –
Use numba functions to perform the twisting-up rotation of modes.
’default_modes’ or
True: use numba for the default_mode_array defined in phenomt.py.’custom_modes’: use numba for a custom mode_array. It is slower than default_modes because of all the if statements.
False: do not use numba.
xp (np/cp) – Module to use for computations on the cpu/gpu.
lmax (int) – Maximum l value for the Ylms to compute (only used for
use_wigner_from_quaternions=False).
- Returns:
Ylm: multidimensional array of Ylms time series. Each row corresponds to a mode in mode_array (same order). gammaTot: float. Global phase factor to correct the strain if
use_wigner_from_quaternions=False.- Return type:
Tuple with ndarray and float
- phenomxpy.precession.quaternion.from_euler_angles(alpha, cosbeta, gamma, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/phenomxpy/envs/stable/lib/python3.13/site-packages/numpy/__init__.py'>, minus_beta=False)[source]
Adapted from quaternion package of Michael Boyle
We use cosbeta instead of beta to speed-up calculations.
Assumes the Euler angles correspond to the quaternion R via
R = exp(alpha*z/2) * exp(beta*y/2) * exp(gamma*z/2)
The angles naturally must be in radians for this to make any sense.
We assume beta [0, pi] => sin(beta/2) >= 0, since this is mostly the case. A minus sign is introduced in q1, q2 if minus_beta=True
- Parameters:
alpha (float or array of floats)
cosbeta (float, or array of floats)
gamma (float, or array of floats)
- Returns:
R
- Return type:
quaternion array
- phenomxpy.precession.quaternion.from_frame(X, Y, Z, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/phenomxpy/envs/stable/lib/python3.13/site-packages/numpy/__init__.py'>)[source]
Compute quaternions from frame.
- phenomxpy.precession.quaternion.from_rotation_matrix(rot, nonorthogonal=True, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/phenomxpy/envs/stable/lib/python3.13/site-packages/numpy/__init__.py'>)[source]
Adapted from M. Boyle quaternion.
Convert input 3x3 rotation matrix to unit quaternion.
For any orthogonal matrix rot, this function returns a quaternion q such that, for every pure-vector quaternion v, we have
q * v * q.conjugate() == rot @ v.vec
Here, @ is the standard python matrix multiplication operator and v.vec is the 3-vector part of the quaternion v. If rot is not orthogonal the “closest” orthogonal matrix is used; see Notes below.
- Parameters:
rot ((..., N, 3, 3) float array) – Each 3x3 matrix represents a rotation by multiplying (from the left) a column vector to produce a rotated column vector. Note that this input may actually have ndims>3; it is just assumed that the last two dimensions have size 3, representing the matrix.
nonorthogonal (bool, optional) – If scipy.linalg is available, use the more robust algorithm of Bar-Itzhack. Default value is True.
- Returns:
q – Unit quaternions resulting in rotations corresponding to input rotations. Output shape is rot.shape[:-2].
- Return type:
array of quaternions
- Raises:
LinAlgError – If any of the eigenvalue solutions does not converge
Notes
By default, if scipy.linalg is available, this function uses Bar-Itzhack’s algorithm to allow for non-orthogonal matrices. [J. Guidance, Vol. 23, No. 6, p. 1085 <http://dx.doi.org/10.2514/2.4654>] This will almost certainly be quite a bit slower than simpler versions, though it will be more robust to numerical errors in the rotation matrix. Also note that Bar-Itzhack uses some pretty weird conventions. The last component of the quaternion appears to represent the scalar, and the quaternion itself is conjugated relative to the convention used throughout this module.
If scipy.linalg is not available or if the optional nonorthogonal parameter is set to False, this function falls back to the possibly faster, but less robust, algorithm of Markley [J. Guidance, Vol. 31, No. 2, p. 440 <http://dx.doi.org/10.2514/1.31730>].
- phenomxpy.precession.quaternion.numba_wigner_from_quaternions_custom_modes_i(wD, R, mode_array)[source]
Compute WignerD matrices from quaternions for a custom mode_array in format [[2,2], [2,-1], …].
Derivation Mathematica notebook.
Modifies input 1D array wD where each element corresponds to one mode, same order as the mode_array
- phenomxpy.precession.quaternion.numba_wigner_from_quaternions_default_modes_i(wD, R)[source]
Compute WignerD matrices from quaternions for default modes: 22, 21, 33, 44, 55 and negatives
Derivation Mathematica notebook.
Modifies the 1D array wD where each element corresponds to one mode, same order as the default mode array
- phenomxpy.precession.quaternion.wigner_from_quaternions(R, mode_array, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/phenomxpy/envs/stable/lib/python3.13/site-packages/numpy/__init__.py'>)[source]
Compute WignerD matrices from quaternions for a custom mode_array with format [[2,2], [2,-1], …].
Derivation Matematica notebook.
Return array where each row corresponds to one mode, same order as the mode_array.