phenomxpy.precession.nnlo

Next-to-next to leading order Euler angles

class phenomxpy.precession.nnlo.NNLO

Bases: object

Class with extra methods for the NNLO angles to be added to pPrec.

compute_nnlo_angles(omega)

Compute analytical NNLO 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

set_NNLO_angles_coefficients()

Set PN coefficients for alpha and gamma.

See Eqs. G9 [6].

phenomxpy.precession.msa

Multi-Scale-Analyis Euler angles

class phenomxpy.precession.msa.MSA

Bases: object

Class 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)

Get c constants from Appendix B (B6, B7, B8) [7] (arXiv:1703.03967).

Constants_d_MSA(LNorm, JNorm, Spl, Spl2, Smi2)

Get d constants from Appendix B (B9, B10, B11) [7] (arXiv:1703.03967).

Get_PN_beta(a, b)

Compute PN spin-orbit couplings.

Get_PN_sigma(a, b)

Compute PN spin-spin couplings.

Get_PN_tau(a, b)

Compute PN spin-spin couplings.

Initialize_MSA(ExpansionOrder=5)

Initialize all the core variables for the MSA system. This will be called first.

JNorm_MSA(L_norm)

Get norm of J using Eq 41 [7] (arXiv:1703.03967).

Lnorm_3PN_of_v(v, v2, L_norm)

Norm of orbital angular momentum L to 3PN order. Eq. 4.7 arXiv:1212.5520.

MSA_Corrections(v, v2, LNorm, JNorm, Spl, Spl2, Smi2, S32)

Get MSA corrections to \(\zeta\) and \(\phi_z\) using Eq. F19 in Appendix F and Eq. 67 [7] (arXiv:1703.03967) respectively.

Psi_MSA(v, v2)

Get \(\psi\) using Eq 51 [7] (arXiv:1703.03967). - 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)

Get \(\dot{\psi}\) using Eq 24 [7] (arXiv:1703.03967). \(\frac{d \psi}{d t} = \frac{A}{2} \sqrt{S_+^2 - S_3^2}\).

Roots_MSA(LNorm, JNorm)

Solve for the roots of Eq 21 [7] (arXiv:1703.03967).

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

SNorm_MSA(v, v2)

Get norm of S, see [7] (arXiv:1703.03967).

Spin_Evolution_Coefficients_MSA(LNorm, JNorm)

Get coefficients for Eq 21 [7] (arXiv:1703.03967).

compute_msa_angles(omega)

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

costhetaLJ(L_norm, J_norm, S_norm)

Calculate (L dot J).

phiz_MSA(v, JNorm)

Get \($\phi_z\) using Eq 66 [7] (arXiv:1703.03967). The coefficients are given in Appendix D (D15 - D26).

psiofv(v, v2, psi0, psi1, psi2)

Equation 51 [7] (arXiv:1703.03967).

zeta_MSA(v)

Eq. F5 [7] (arXiv:1703.03967). \(\zeta_{z,-1} = \eta v^{-3} \sum^5_{n=0} <\Omega_{\zeta}>^n v^n + \zeta_{-1}^0\). Note that the \(<\Omega_{\eta}>^n\) are given by self.Omegazeta(n)_coeff’s as in Eqs. F6-F11.

phenomxpy.precession.numerical

Spin-Taylor numerical Euler angles.

class phenomxpy.precession.numerical.Numerical

Bases: object

Class with extra methods for evolving spin equations and compute numerical angles added to pPrec.

backward_integration(times=None)

Evolve orbit from tref to tmin.

compute_numerical_angles(LNhat, times)

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

delete_evolved_quantities(boolean)

Remove LNhatev, S1/2ev, E1ev, v_interpolant from pPrec object. If the Python garbage collector does act often enough, these quantities could make the memory usage grow to much.

evolve_orbit(tinit, tend, t_eval=None)

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)

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

forward_integration(times=None)

Evolve orbit from tref to t=0.

gamma_integration(times, alpha, cosbeta)

Numerical integration of \(\dot{\gamma}\) over the full time array to fulfill minimal rotation condition.

\(\dot{\gamma} = - \dot{\alpha} \cos(\beta)\). Eq. 18 [8].

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

set_v_interpolant(times)

Build v_interpolant for the rhs of the ODE solver.

class phenomxpy.precession.numerical.SpinTaylor(eta, m1M, m2M)

Bases: object

Class to setup SpinTaylorT4 quantities and compute spin derivatives.

energy_spin_derivative_setup(quadparam1=0, quadparam2=0)

Setup energy and spin derivative coefficients.

Adapted from XLALSimSpinTaylorEnergySpinDerivativeSetup.

spin_derivatives_avg(v, LNh, E1, S1, S2)

No numba wrapper for the spin derivatives function.

phenomxpy.precession.numerical._spin_derivatives_avg(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.eval_cubic_spline(x, x_breaks, coefs)

Evaluate a cubic spline in a single point with numba.

phenomxpy.precession.numerical.eval_cubic_spline_array(x_array, x_breaks, coefs)

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/latest/lib/python3.13/site-packages/numpy/__init__.py'>)

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.manual_cross(a, b)

Cross product of two 3D vectors.

phenomxpy.precession.numerical.pn_energy_3pn_so(mByM)

Eq. (4.6) of arXiv:1212.5520

phenomxpy.precession.numerical.pn_energy_4PN_qm_s1os1o(mByM)

Eq. (6) of arXiv:astro-ph/0504538v2

phenomxpy.precession.numerical.pn_energy_4PN_qm_s1s1(mByM)

Eq. (6) of arXiv:astro-ph/0504538v2

phenomxpy.precession.numerical.pn_energy_4pn_s1os2o(eta)

Eq. (6) of arXiv:astro-ph/0504538v2

phenomxpy.precession.numerical.pn_energy_4pn_s1s2(eta)

Eq. (6) of arXiv:astro-ph/0504538v2

phenomxpy.precession.numerical.pn_energy_5pn_so(mByM)

Eq. 4.6 of arXiv:1212.5520

phenomxpy.precession.numerical.pn_energy_7pn_so(mByM)

Eq. (4.6) of arXiv:1212.5520 Symbol definitions right above eq. (3.1)

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.spin_dot_3pn(mByM)

dS1, 1.5PN: eq. (8) of gr-qc/0405090

phenomxpy.precession.numerical.spin_dot_4pn_qm_so(mByM)

S1S1 contribution see eq. A.2 of arXiv:1501.01529

phenomxpy.precession.numerical.spin_dot_4pn_s2()

S1S2 contribution see. eq. A.2 of arXiv:1501.01529

phenomxpy.precession.numerical.spin_dot_4pn_s2o()

S1S2 contribution see. eq. A.2 of arXiv:1501.01529

phenomxpy.precession.numerical.spin_dot_5pn(mByM)

dS1, 2.5PN eq. 7.8 of Blanchet et al. gr-qc/0605140

phenomxpy.precession.numerical.spin_dot_7pn(mByM)

dS1, 3.5PN Eq. 3.4 of Bohe’ et al. arXiv:1212.5520

phenomxpy.precession.numerical.tphm_spin_derivatives(t, y, v_interpolant, spin_taylor_t4)

Spin derivatives using a v_interpolant without numba.

phenomxpy.precession.numerical.tphm_spin_derivatives_imr_omega(t, y, pPhase22, spin_taylor_t4)

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)

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)

Spin derivatives using a v_interpolant with numba.

phenomxpy.precession.numerical.unpack(y)

Unpack the state vector y into its components LNh, E1, S1, S2.

phenomxpy.precession.numerical.v_of_omega(omega)

Compute v from a single value of omega.

phenomxpy.precession.precession

Define pPrec class and useful function for precession

class phenomxpy.precession.precession.pPrec(pwf, params: pPrecParams | dict | None = None, **kwargs)

Bases: object

Class 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_version is ‘numerical’, final_spin_version uses 4. ‘msa’ uses 3. ‘nnlo’ uses 0.

  • 0: \(\bar{\chi}_p = \chi_p\) Eq. 4.7-4.9 [6] (arxiv:2004.06503).

  • 1: \(\bar{\chi}_p = \chi_{1x}\)

  • 2: \(\bar{\chi}_p = \chi_\perp\) related to in-plane spin components, Eq. 4.24, 4.28 [6] (arxiv:2004.06503).

  For the options above, \(\bar{\chi}_p\) is then plugged into Eq. 4.25 [6] (arxiv:2004.06503) to obtain the final spin.

  • 3: For MSA angles, uses spin averaged quantities. Eq. 4.23, 4.29 [6] (arxiv:2004.06503).

  • 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 and not much slower.

• 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.

• polarizations_from ({'L0modes', 'Jmodes', 'CPmodes'}) –

  Strategy to compute the polarizations that should return the same output.

  • LOmodes: method employed in LAL involving two rotations: CP->J, J->L0

  • Jmodes: similar to PhenomX, rotate modes from CP->J and compute polarizations using theta_JN instead of inclination

  • CPmodes (default): use efficient method with quaternions from M. Boyle, Appendix B of [9].

• use_wigner_from_quaternions (bool) – If False, then use the method of seobnrv5.

• compute_CPmodes_at_once (bool) – Compute all the co-precessing modes at once, this is simpler but uses more memory. If False, compute one mode and add its contribution to the polarizations in a loop over modes.

• compute_ylms_at_once (bool) – Same as for the CPmodes but for the Ylms time series.

WignerDMatrix(l, mp, m)

Compute Wigner-D matrix.

WignerdMatrix(cosBeta, sign=1, lmax=2, global_rotation=False)

Compute Wigner-d matrices coefficients from cosbeta.

Appendix A [6] (arxiv:2004.06503).

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)

Add specific methods from an input class (NNLO, MSA, Numerical) to the pPrec class

compute_angular_momentum(v)

Compute orbital angular momentum up to 4PN approximantion.

Eq. 4.13 [6] (arxiv:2004.06503).

get_chip()

Compute \(\chi_p\). Originally in Eq.3.3-3.4 [10] (arxiv:1408.1810).

get_euler_slope(af, mf)

Slope for the alpha angle in the ringdown attachment.

rotate_y(angle, old)

Rotate vector around y-axis.

rotate_z(angle, old)

Rotate vector around z-axis.

set_angular_momentum_coefficients()

Set coefficients for the PN angular momenum L.

Coefficients in Eqs. G10 [6] (arxiv:2004.06503).

Reconstruction of L in Eq. 4.13.

set_exp_angle_powers(alpha, gamma)

Set complex exponential values needed for Wigner-D matrices.

set_final_spin(final_spin_version)

Set final spin according to the final_spin_version

• 0: \(\bar{\chi}_p = \chi_p\) Eq. 4.7-4.9 [6] (arxiv:2004.06503).

• 1: \(\bar{\chi}_p = \chi_{1x}\)

• 2: \(\bar{\chi}_p = \chi_\perp\) related to in-plane spin components, Eq. 4.24, 4.28 [6] (arxiv:2004.06503).

For the options above, \(\bar{\chi}_p\) is then plugged into Eq. 4.25 [6] (arxiv:2004.06503) to obtain the final spin.

• 3: For MSA angles, uses spin averaged quantities. Eq. 4.23, 4.29 [6] (arxiv:2004.06503).

• 4: Numerically evolved (this is done outside the pPrec class).

set_initial_quantities()

Angles and vectors at fref

set_sqrt()

Set useful values of square roots used in the Wigner-d matrices.

class phenomxpy.precession.precession.pPrecParams(*, prec_version: Literal['numerical', 'msa', 'nnlo'] = 'numerical', final_spin_version: int | None = None, afinal_prec: float | None = None, add_20_mode: bool = False, numba_rotation: Literal['default_modes', 'custom_modes'] | None = 'default_modes', store_arrays: bool = False, save_omega: bool | None = None, fall_back_msa_to_nnlo: bool = True, use_final_spin_sign: bool = True, use_closest_time_to_tref: bool = False, v_function: Literal['interpolant', 'imr_omega'] = 'imr_omega', cubic_interpolation_for_ode: bool = True, numba_derivatives: bool = True, interpolate: Literal['ode_solution', 'euler_angles'] | None = 'ode_solution', quaternions_from: Literal['frame', 'euler_angles'] = 'frame', gamma_integration_method: Literal['piecewise_integral', 'boole', 'antiderivative'] = 'piecewise_integral', analytical_RD_gamma: bool = True, ode_solver: str = 'RK45', atol_prec: float = 1e-12, rtol_prec: float = 1e-12, polarizations_from: Literal['CPmodes', 'Jmodes', 'L0modes'] = 'CPmodes', compute_CPmodes_at_once: bool = False, compute_ylms_at_once: bool = True, use_wigner_from_quaternions: bool = True)

Bases: BaseModel

Auxiliary pydantic class to handle parameters supported by pPrec class

_abc_impl = <_abc._abc_data object>

add_20_mode: bool

afinal_prec: float | None

analytical_RD_gamma: bool

atol_prec: float

compute_CPmodes_at_once: bool

compute_ylms_at_once: bool

cubic_interpolation_for_ode: bool

fall_back_msa_to_nnlo: bool

final_spin_version: int | None

gamma_integration_method: Literal['piecewise_integral', 'boole', 'antiderivative']

interpolate: Literal['ode_solution', 'euler_angles'] | None

model_config = {}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

classmethod normalize_mode(v)

Lowercase the value for prec_version option

numba_derivatives: bool

numba_rotation: Literal['default_modes', 'custom_modes'] | None

ode_solver: str

polarizations_from: Literal['CPmodes', 'Jmodes', 'L0modes']

prec_version: Literal['numerical', 'msa', 'nnlo']

quaternions_from: Literal['frame', 'euler_angles']

rtol_prec: float

save_omega: bool | None

store_arrays: bool

use_closest_time_to_tref: bool

use_final_spin_sign: bool

use_wigner_from_quaternions: bool

v_function: Literal['interpolant', 'imr_omega']

classmethod warn_unknown_fields(values)

Warning about unknown arguments provided to pPrec

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/latest/lib/python3.13/site-packages/numpy/__init__.py'>)

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)

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.ell_prefactor(l: int)

Prefactor needed for the Wigner coefficients from quaternions for each l.

Derivation Matematica notebook.

phenomxpy.precession.quaternion.from_euler_angles(alpha, cosbeta, gamma, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/phenomxpy/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>, minus_beta=False)

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/latest/lib/python3.13/site-packages/numpy/__init__.py'>)

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/latest/lib/python3.13/site-packages/numpy/__init__.py'>)

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(R, mode_array)

Evaluate wigner coefficients from an array of quaternions. Custom array of co-precessing modes.

phenomxpy.precession.quaternion.numba_wigner_from_quaternions_custom_modes_i(wD, R, mode_array)

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(R)

Evaluate wigner coefficients from an array of quaternions. Only when using the co-precessing default modes.

phenomxpy.precession.quaternion.numba_wigner_from_quaternions_default_modes_i(wD, R)

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.quaternion_invert(q, xp)

Compute the inverse of a quaternion.

Parameters:

q – (…, 4) CuPy array of quaternions [w, x, y, z]

Returns:

(…, 4) CuPy array of inverted quaternions

phenomxpy.precession.quaternion.quaternion_multiply(q, r, xp)

Multiply two quaternions elementwise.

Parameters:

• q – (…, 4) CuPy array representing quaternions [w, x, y, z]

• r – (…, 4) CuPy array representing quaternions [w, x, y, z]

Returns:

(…, 4) CuPy array of quaternion products

phenomxpy.precession.quaternion.wigner_from_quaternions(R, mode_array, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/phenomxpy/envs/latest/lib/python3.13/site-packages/numpy/__init__.py'>)

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.