# Copyright (C) 2023, 2024, 2025 Cecilio García Quirós
"""
IMRPhenomTP
-----------
.. autoclass:: phenomxpy.phenomt.phenomtp.IMRPhenomTP
:members:
:private-members:
:undoc-members:
:show-inheritance:
:noindex:
IMRPhenomTPHM
-------------
.. autoclass:: phenomxpy.phenomt.phenomtp.IMRPhenomTPHM
:members:
:private-members:
:undoc-members:
:show-inheritance:
:noindex:
"""
import numpy as np
try:
import cupy as cp
except ImportError:
cp = None
import math
from scipy.spatial.transform import Rotation
from phenomxpy.phenomt.internals import pWFHM, pPhase, pWFParams
from phenomxpy.phenomt.phenomt import _PhenomT, IMRPhenomTHM, THMParams
from phenomxpy.utils import SpinWeightedSphericalHarmonic, rotate_by_polarization_angle
from phenomxpy.precession.precession import pPrec, pPrecParams
from phenomxpy.precession.quaternion import (
from_euler_angles,
from_frame,
as_euler_angles,
compute_ylms,
quaternion_invert,
quaternion_multiply,
)
import logging
logger = logging.getLogger(__name__)
class TPHMParams(THMParams, pPrecParams):
pass
[docs]
class IMRPhenomTPHM(_PhenomT):
"""
Class for the IMRPhenomTPHM model :cite:`phenomtphm`.
Multipolar quasicircular precessing-spin waveform model.
The mode_array is specified in the co-precessing frame.
"""
def __init__(self, params: TPHMParams | dict | None = None, **kwargs):
# Allow either a Params object, a dictionary or kwargs
if params is not None:
if isinstance(params, dict):
tphm_params = TPHMParams(**params)
else:
tphm_params = params
else:
tphm_params = TPHMParams(**kwargs)
tphm_params = tphm_params.model_dump()
# Split parameters for each class
pwf_params = pWFParams.model_construct(**{k: v for k, v in tphm_params.items() if k in pWFParams.model_fields})
pprec_params = pPrecParams.model_construct(**{k: v for k, v in tphm_params.items() if k in pPrecParams.model_fields})
thm_params = THMParams.model_construct(**{k: v for k, v in tphm_params.items() if k in THMParams.model_fields})
# Set pWF struct
self.pWF = pWFHM(mode=[2, 2], params=pwf_params)
# Set pPrec struct.
# This updates pWF.afinal_prec according to the final_spin_version, except for version 4 which is done next
self.pPrec = pPrec(self.pWF, params=pprec_params)
# Copy module for easier use
self.xp = self.pPrec.xp
# ONLY for NUMERICAL ANGLES
if self.pPrec.prec_version == "numerical":
# Set pPhase22 for computing the v_function used in the ODE solver
self.pPrec._pPhase22 = pPhase(self.pWF)
if self.pPrec.final_spin_version == 4 and self.pPrec.input_afinal_prec is None:
# Solve PN spin equations to obtain the final spin, recycle the evolved quantities later to compute the Euler angles
# If delta_t is not set (=0), there is not a time array where to compute the interpolant for v(t),
# in that case, switch to v_function = "imr_omega" which does not need a pre-defined time array.
if self.pWF.delta_t > 0:
times = self.set_time_array()
else:
times = np.array([self.pWF.tref, 0.0])
logger.info("Switching to v_function='imr_omega' since delta_t is not provided.")
self.pPrec.v_function = "imr_omega"
# Build interpolant for v(t) used in the rhs of the ODE solver.
# The other option is to use "imr_omega" function directly (that is more precise).
if self.pPrec.v_function == "interpolant":
self.pPrec.set_v_interpolant(times)
# Forward Integration: from tref to t=0
self.pPrec.forward_integration(times=times)
# Compute final spin from last point (t=0) in the evolution
self.pPrec.evolved_final_spin(self.pPrec.LNhatev, self.pPrec.S1ev, self.pPrec.S2ev)
# Initialize co-precessing class
self.phenTHM = IMRPhenomTHM(pWF_input=self.pWF, params=thm_params)
# Populate global pWF with values computed in co-precessing
key = "22"
if self.pWF.tmin is None:
self.pWF.tmin = self.phenTHM.phenT_classes[key].pWF.tmin
if self.pWF.tref is None:
self.pWF.tref = self.phenTHM.phenT_classes[key].pWF.tref
if self.pWF.delta_t > 0:
# Set equispaced time array lengths
if hasattr(self.pWF, "length") is False:
self.pWF.length = self.phenTHM.phenT_classes[key].pWF.length
if hasattr(self.pWF, "len_neg") is False:
self.pWF.len_neg = self.phenTHM.phenT_classes[key].pWF.len_neg
if hasattr(self.pWF, "len_pos") is False:
self.pWF.len_pos = self.phenTHM.phenT_classes[key].pWF.len_pos
# Check consistency time array with co-precessing
if math.isclose(self.pWF.tmin, self.phenTHM.phenT_classes[key].pWF.tmin, rel_tol=1e-15) is False:
raise AssertionError(
f"Inconsistent tmin with co-precessing {self.pWF.tmin}!={self.phenTHM.phenT_classes[key].pWF.tmin}"
)
if math.isclose(self.pWF.tref, self.phenTHM.phenT_classes[key].pWF.tref, rel_tol=1e-15) is False:
raise AssertionError(
f"Inconsistent tref with co-precessing {self.pWF.tref}!={self.phenTHM.phenT_classes[key].pWF.tref}"
)
if (
self.pWF.delta_t > 0
and math.isclose(self.pWF.length, self.phenTHM.phenT_classes[key].pWF.length, rel_tol=1e-15) is False
):
raise AssertionError(
f"Inconsistent length with co-precessing {self.pWF.length}!={self.phenTHM.phenT_classes[key].pWF.length}"
)
# Set epoch
self.epoch = self.phenTHM.epoch
[docs]
def compute_euler_angles(self, times=None):
r"""
Compute Euler angles alpha, cosbeta, gamma for a time array.
3 prescriptions supported:
- NNLO
- MSA
- Numerical/SpinTaylor (default)
Parameters
----------
times: 1D ndarray
Time array where to evaluate the angles. If ``None``, use equispaced grid ``self.times``.
Returns
-------
Tuple with 3 1D ndarrays
:math:`\\alpha(t)`, :math:`\cos(\\beta(t))`, :math:`\\gamma(t)`,
"""
# Variable to decide if we need to evolve the orbit in the new time array (numerical angles only)
new_times = True if times is not None else False
# Set time array
# set_time_array was called at the initialization step for the case of uniform grids, here we account for the custom time array case.
times = self.set_time_array(times=times)
if times[0] < self.pWF.tmin:
raise ValueError("Time array starts before tmin")
# Compute omega(t) in the time array
omega = self.phenTHM.phenT_classes["22"].pPhase.imr_omega(times[times <= 0])
##########################
# NNLO - MSA ANGLES #
##########################
if self.pPrec.prec_version in ["nnlo", "msa"]:
# Choose proper function to evaluate the angles: compute_nnlo_angles or compute_msa_angles
euler_function = getattr(self.pPrec, "compute_" + self.pPrec.prec_version + "_angles")
# Compute angles in time array up to mergeer t<=0
alpha, cosbeta, gamma = euler_function(omega)
# Values at t=0 for the ringdown attachment
self.pPrec.omegaPeak = self.phenTHM.phenT_classes["22"].pPhase.omegaPeak
self.pPrec.alphaRD0, self.pPrec.cosbetaRD0, self.pPrec.gammaRD0 = euler_function(self.pPrec.omegaPeak)
# Angles for global rotation from J to L0
self.pPrec.alphaJtoL0, self.pPrec.cosbetaJtoL0, self.pPrec.gammaJtoL0 = euler_function(self.pPrec.omegaRef)
# Ringdown part
alphaRD = self.pPrec.alphaRD0 + self.pPrec.EulerRDslope * times[times > 0]
cosbetaRD = self.xp.full(len(times[times > 0]), self.pPrec.cosbetaRD0)
gammaRD = self.pPrec.gammaRD0 - alphaRD * cosbetaRD + self.pPrec.alphaRD0 * self.pPrec.cosbetaRD0
# Attach ringdown
alpha = self.xp.concatenate((alpha, alphaRD))
cosbeta = self.xp.concatenate((cosbeta, cosbetaRD))
gamma = self.xp.concatenate((gamma, gammaRD))
##########################
# NUMERICAL ANGLES #
##########################
else:
###############################################
# FORWARD INTEGRATION: from tref to t=0 #
###############################################
# The PN spin equations have not been yet evolved in the initialization of the class if final_spin_version = 4,
# or a input_afinal_prec is given.
# If that is not the case, we evolve them now from tref to merger (t=0)
if self.pPrec.final_spin_version != 4 or self.pPrec.input_afinal_prec is not None or new_times:
if self.pPrec.v_function == "interpolant" and hasattr(self.pPrec, "v_interpolant") is False:
self.pPrec.set_v_interpolant(times)
self.pPrec.forward_integration(times)
#################################################
# BACKWARD INTEGRATION: from tref to tmin #
#################################################
if self.pWF.tmin < self.pWF.tref:
LNhatev, _ = self.pPrec.backward_integration(times)
else:
LNhatev = self.pPrec.LNhatev
# Explictly delete evolved quantities to prevent memory growth
self.pPrec.delete_evolved_quantities(not self.pPrec.store_arrays)
# Compute angles in the full time array (PN evolved + Ringdown attachment)
alpha, cosbeta, gamma = self.pPrec.compute_numerical_angles(LNhatev, times)
# Optionally store arrays.
if self.pPrec.store_arrays is True:
self.alpha, self.cosbeta, self.gamma = alpha, cosbeta, gamma
return alpha, cosbeta, gamma
[docs]
def compute_quaternions(self, times=None):
"""
Compute quaternions in time array to transform from L to L0 frame and align with the line-of-sight.
For the NNLO and MSA prescriptions, the quaternions are computed from the Euler angles.
For the numerical prescription, the default option is to compute the quaternions directly from the evolution of the L frame:
the Z axis is LNhatev, the X axis is named E1ev (same name as in lalsuite), and the Y axis is obtained by completing the triad with the cross product Y=ZxX.
There is also the option to compute the quaternions from the numerical Euler angles. However, this option involves more computational steps and potential
loss of accuracy due to extra numerical integrations and spline derivatives.
Parameters
----------
times: 1D ndarray
Time array to compute quaternions for. If ``None``, use internal equispaced grid ``self.times``.
torch_single_precision: bool
If ``True``, use single precision for ``torch`` tensors, slightly faster.
Returns
-------
(N, 4) ndarray
Quaternions evaluated in time array.
"""
torch_single_precision = self.pWF.extra_options.get("torch_single_precision", False)
# Quaternion corresponding to the alignment with the line-of-sight
qL02I0 = from_euler_angles(np.pi / 2 - self.pWF.phi_ref, np.cos(self.pWF.inclination), 0.0, xp=self.xp)
################################
# Quaternions from frame #
################################
# Only for numerical angles.
# Skips euler angles computation, calculating the quaternions directly from L, E1 evolution
if self.pPrec.prec_version == "numerical" and self.pPrec.quaternions_from == "frame":
# Variable to decide if we need to evolve the orbit in the new time array (numerical angles only)
new_times = True if times is not None else False
# Set time array
# set_time_array was called at the initialization step for the case of uniform grids, here we account for the custom time array case.
times = self.set_time_array(times=times)
if times[0] < self.pWF.tmin:
raise ValueError(f"Time array starts before tmin {times[0]} < {self.pWF.tmin}")
#################################################
# FORWARD INTEGRATION: from tref to t=0 #
#################################################
# PN spin equations have not been yet evolved if final_spin_version != 4,
# or if the final spin is provided through input_afinal_prec. Do forward integration
if self.pPrec.final_spin_version != 4 or self.pPrec.input_afinal_prec is not None or new_times:
if self.pPrec.v_function == "interpolant" and hasattr(self.pPrec, "v_interpolant") is False:
self.pPrec.set_v_interpolant(times)
self.pPrec.forward_integration(times)
#################################################
# BACKWARD INTEGRATION: from tref to tmin #
#################################################
if self.pWF.tmin < self.pWF.tref:
LNhatev, E1ev = self.pPrec.backward_integration(times)
else:
LNhatev = self.pPrec.LNhatev
E1ev = self.pPrec.E1ev
# Explictly delete evolved quantities to prevent memory growth
self.pPrec.delete_evolved_quantities(not self.pPrec.store_arrays)
# Set axis of evolved frame
X = E1ev
Z = LNhatev
Y = self.xp.cross(Z, X, axis=0)
# Compute quaternions from evolved frame
qtwist = from_frame(X, Y, Z, xp=self.xp)
# Ringdown attachment
# Need to transform to the J-frame where the Euler angles for the attachment are defined
qL0toJ = from_euler_angles(
-self.pPrec.kappa, np.cos(-self.pPrec.thetaJ_Sf), -self.pPrec.phiJ_Sf, minus_beta=True, xp=self.xp
)
alphaRD0, cosbetaRD0, gammaRD0 = as_euler_angles(quaternion_multiply(qL0toJ, qtwist[-1], xp=self.xp))
# Prevent rounding error for cosbetaRD0
if abs(cosbetaRD0) > 1:
cosbetaRD0 = np.copysign(1, cosbetaRD0)
# Remove t=0 if it was added in forward_integration
if self.pPrec.t0_added_to_array:
qtwist = qtwist[:-1]
RDtimes = times[times > 0]
alphaRD = alphaRD0 + self.pPrec.EulerRDslope * RDtimes
cosbetaRD = self.xp.full(len(alphaRD), cosbetaRD0)
gammaRD = gammaRD0 - self.pPrec.EulerRDslope * cosbetaRD0 * RDtimes
# These values will not coincide with the method from Euler angles, but the final rotation is consistent
self.pPrec.alphaRD0 = alphaRD0
self.pPrec.cosbetaRD0 = cosbetaRD0
self.pPrec.gammaRD0 = gammaRD0
q1 = quaternion_invert(from_euler_angles(alphaRD, cosbetaRD, gammaRD, xp=self.xp), xp=self.xp)
qtwistRD = quaternion_multiply(q1, qL0toJ, xp=self.xp)
# Attach RD part
q1 = quaternion_invert(qtwist, xp=self.xp)
qtwist = self.xp.concatenate((q1, qtwistRD), axis=0)
# Final quaternion
qTot = quaternion_multiply(qtwist, qL02I0, xp=self.xp)
######################################
# Quaternions from Euler angles #
######################################
else:
# Compute Euler angles
alpha, cosbeta, gamma = self.compute_euler_angles(times=times)
# Final quaternion
qP2J = from_euler_angles(alpha, cosbeta, gamma, xp=self.xp)
qJ2L0 = from_euler_angles(self.pPrec.alphaJtoL0, self.pPrec.cosbetaJtoL0, self.pPrec.gammaJtoL0, xp=self.xp)
qP2Jb = quaternion_invert(qP2J, xp=self.xp)
qTot = quaternion_multiply(qP2Jb, quaternion_multiply(qJ2L0, qL02I0, xp=self.xp), xp=self.xp)
if self.pPrec.store_arrays is True:
self.qTot = qTot
return qTot
[docs]
def compute_CPmodes(self, times=None):
"""
Compute hlms(t) in the co-precessing frame for a given time array. Use internal array if ``times`` is ``None``.
Parameters
----------
times: 1D ndarray
Time array where to evaluate the modes
Returns
-------
Tuple with dictionary of modes and time array
CPmodes(t), t
"""
CPmodes, times = self.phenTHM.compute_hlms(times=times)
if self.pPrec.store_arrays is True:
self.CPmodes = CPmodes
return CPmodes, times
[docs]
def compute_Jmodes(self, times=None):
"""
Compute hlms(t) in the J-frame for a given time array. Use internal array if ``times`` is ``None``.
Parameters
----------
times: 1D ndarray
Time array where to evaluate the modes
Returns
-------
Tuple with dictionary of modes and time array
Jmodes(t), t
"""
# Compute Euler angles
alpha, cosbeta, gamma = self.compute_euler_angles(times=times)
# Compute CPmodes
CPmodes, times = self.compute_CPmodes(times=times)
# Initialize Jmodes dictionary to zero
Jmodes = {}
n_inertial_modes = 0
for l in range(2, self.phenTHM.lmax + 1):
for m in range(-l, l + 1):
n_inertial_modes = n_inertial_modes + 1
mode_string = str(l) + str(m)
Jmodes[mode_string] = self.xp.zeros(len(alpha), dtype=self.xp.cdouble)
# Compute Wigner-d matrices for Euler angle rotations
self.pPrec.set_exp_angle_powers(gamma, alpha)
self.pPrec.WignerdMatrix(cosbeta, sign=1, lmax=self.phenTHM.lmax)
# Rotation from CP to J frame
for l in range(2, self.phenTHM.lmax + 1):
for m in range(-l, l + 1):
Jmode = str(l) + str(m)
for mp in range(-l, l + 1):
comode = str(l) + str(mp)
if comode in CPmodes:
WignerD = self.pPrec.WignerDMatrix(l, mp, m)
Jmodes[Jmode] += WignerD * CPmodes[comode]
return Jmodes, times
[docs]
def compute_L0modes(self, times=None):
"""
Compute hlms(t) in the L0-frame for a given time array. Use internal array if ``times`` is ``None``..
Parameters
----------
times: 1D ndarray
Time array where to evaluate the modes
Returns
-------
Tuple with dictionary of modes and time array
L0modes(t), t
"""
# Compute Jmodes
Jmodes, times = self.compute_Jmodes(times=times)
# Compute Wigner-d matrices for Euler angle rotation
self.pPrec.set_exp_angle_powers(-self.pPrec.alphaJtoL0, -self.pPrec.gammaJtoL0)
self.pPrec.WignerdMatrix(self.pPrec.cosbetaJtoL0, sign=-1, global_rotation=True, lmax=self.phenTHM.lmax)
# Rotation from J to L0 frame
# Get length of first value from the dict
N = len(next(iter(Jmodes.values())))
# Initialize L0modes to zero
L0modes = {}
for l in range(2, self.phenTHM.lmax + 1):
for m in range(-l, l + 1):
mode_string = str(l) + str(m)
L0modes[mode_string] = self.xp.zeros(N, dtype=self.xp.cdouble)
# Modes rotation
for l in range(2, self.phenTHM.lmax + 1):
for m in range(-l, l + 1):
L0mode = str(l) + str(m)
for mp in range(-l, l + 1):
Jmode = str(l) + str(mp)
WignerD = self.pPrec.WignerDMatrix(l, mp, m)
L0modes[L0mode] += WignerD * Jmodes[Jmode]
return L0modes, times
[docs]
def compute_polarizations(self, times=None):
"""
Compute polarizations hp(t), hc(t) for a given time array. Use internal array if ``times`` is ``None``.
Parameters
----------
times: 1D ndarray
Time array where to evaluate the polarizations
Returns
-------
Tuple with 3 1D ndarrays
hp(t), hc(t), t
"""
# Read extra options
polarizations_from = self.pPrec.polarizations_from
compute_CPmodes_at_once = self.pPrec.compute_CPmodes_at_once
compute_ylms_at_once = self.pPrec.compute_ylms_at_once
use_wigner_from_quaternions = self.pPrec.use_wigner_from_quaternions
####################################
# polarizations_from L0modes #
####################################
if polarizations_from == "L0modes":
# Compute L0modes
hlms, times = self.compute_L0modes(times=times)
# Angles for the Ylms
polar, azimuthal = self.pWF.inclination, np.pi / 2 - self.pWF.phi_ref
####################################
# polarizations_from Jmodes #
####################################
elif polarizations_from == "Jmodes":
# Compute Jmodes
hlms, times = self.compute_Jmodes(times=times)
# Compute polar angles for Ylms
N = [
np.sin(self.pWF.inclination) * np.cos(np.pi / 2 - self.pWF.phi_ref),
np.sin(self.pWF.inclination) * np.sin(np.pi / 2 - self.pWF.phi_ref),
np.cos(self.pWF.inclination),
]
r = Rotation.from_euler(
"zyz", [self.pPrec.gammaJtoL0, np.arccos(self.pPrec.cosbetaJtoL0), self.pPrec.alphaJtoL0], degrees=False
)
N_Jf = r.apply(N)
polar, azimuthal = np.arccos(N_Jf[2]), np.arctan2(N_Jf[1], N_Jf[0])
self.pPrec.J0x_Sf = self.pPrec._pWF.m1_2 * self.pPrec._pWF.s1[0] + self.pPrec._pWF.m2_2 * self.pPrec._pWF.s2[0]
self.pPrec.J0y_Sf = self.pPrec._pWF.m1_2 * self.pPrec._pWF.s1[1] + self.pPrec._pWF.m2_2 * self.pPrec._pWF.s2[1]
self.pPrec.J0z_Sf = (
self.pPrec._pWF.m1_2 * self.pPrec._pWF.s1[2] + self.pPrec._pWF.m2_2 * self.pPrec._pWF.s2[2] + self.pPrec.LRef
)
J0 = [self.pPrec.J0x_Sf, self.pPrec.J0y_Sf, self.pPrec.J0z_Sf]
J0_norm = np.linalg.norm(J0)
J_L0 = r.inv().apply([0, 0, J0_norm])
yN = np.cross(np.array([0, 0, 1]), N)
yN /= np.linalg.norm(yN)
xN = np.cross(yN, N)
xN /= np.linalg.norm(xN)
yNp = np.cross(J_L0, N)
yNp /= np.linalg.norm(yNp)
xNp = np.cross(yNp, N)
xNp /= np.linalg.norm(xNp)
# Eq. C22 arxiv:2004.06503
self.pPrec.zeta = np.arctan2(np.dot(xN, yNp), np.dot(xN, xNp))
####################################
# polarizations_from CPmodes #
####################################
elif polarizations_from == "CPmodes":
############################
# Compute quaternions #
############################
qTot = self.compute_quaternions(times=times)
############################
# Ylms from quaternions #
############################
if compute_ylms_at_once:
# If mode_array is not the default, then use custom_modes numba_rotation function
if self.phenTHM.mode_array != self.phenTHM.default_mode_array and (
self.pPrec.numba_rotation in ["default_modes", True]
):
self.pPrec.numba_rotation = "custom_modes"
Ylm, gammaTot = compute_ylms(
qTot,
self.phenTHM.mode_array,
use_wigner_from_quaternions,
self.pPrec.numba_rotation,
self.xp,
self.phenTHM.lmax,
)
##################################
# Compute strain from CPmodes #
##################################
# Intialize strain
strain = self.xp.zeros(len(qTot), dtype=complex)
# Compute all CPmodes at once. Simple way, uses more memory
if compute_CPmodes_at_once:
# Compute co-precessing modes
CPmodes, times = self.compute_CPmodes(times=times)
# Loop over modes and add strain contribution
for i, [l, m] in enumerate(self.phenTHM.mode_array):
mode = str(l) + str(m)
# Build strain
strain += Ylm[i] * CPmodes[mode]
# Another option is to generate hlms one by one to save memory (default)
else:
cache = None
modes_already_computed = []
for idx, [l, m] in enumerate(self.phenTHM.mode_array):
if [l, m] not in modes_already_computed:
mode = str(l) + str(m)
# NOTE only valid for equatorial symmetry
# Skip the evaluation of the opposite mode by using symmetry
# marray_ylms is to compute the Ylms one by one to save memory
if m != 0 and [l, -m] in self.phenTHM.mode_array:
add_opposite_mode = True
marray_ylms = [[l, m], [l, -m]]
else:
add_opposite_mode = False
marray_ylms = [[l, m]]
modes_already_computed.append([l, -m])
# Compute hlm
key = str(l) + str(abs(m))
phen = self.phenTHM.phenT_classes[key]
if key == "22":
hlm, times, cache = phen.compute_hlm(times=times, return_cache=True)
else:
hlm, times = phen.compute_hlm(times=times, cache=cache)
# Apply equatorial symmetry for negative modes
if m < 0:
hlm = (-1) ** l * self.xp.conj(hlm)
# Compute Ylm if not computed all at once (this saves memory but it is slower)
if compute_ylms_at_once is False:
Ylm, gammaTot = compute_ylms(
qTot,
marray_ylms,
use_wigner_from_quaternions,
self.pPrec.numba_rotation,
self.xp,
self.phenTHM.lmax,
)
i, i_opposite = 0, 1
# All Ylms are already computed. Get proper indeces
else:
i = idx
# Find index in custom mode array for the opposite mode
if add_opposite_mode:
matches = np.all(np.array(self.phenTHM.mode_array) == [l, -m], axis=1)
i_opposite = np.where(matches)[0][0]
# Add mode contribution to strain
strain += Ylm[i] * hlm
# Add opposite mode. NOTE only valid for equatorial symmetry
if add_opposite_mode:
strain += Ylm[i_opposite] * (-1) ** l * self.xp.conj(hlm)
# Correction needed for the seobnrv5 method. Eq.28 arxiv:2303.18046
if use_wigner_from_quaternions is False:
strain *= self.xp.exp(2j * gammaTot)
# Option not supported for polarizations_from
else:
raise ValueError(
f"polarizations_from={polarizations_from} not supported. Available options are CPmodes (default), L0modes (LAL) or Jmodes (PhenomX)"
)
########################################
# Compute strain for Jmodes, L0modes #
########################################
if polarizations_from != "CPmodes":
keys = list(hlms.keys())
strain = self.xp.zeros(len(hlms[keys[0]]), dtype=self.xp.cdouble)
for key, hlm in hlms.items():
ell = int(key[0])
emm = int(key[1:])
Ylm = SpinWeightedSphericalHarmonic(polar, azimuthal, ell, emm)
strain += hlm * Ylm
# Zeta correction for Jmodes. Eqs. D6-D7 arxiv:2004.06503
if polarizations_from == "Jmodes":
strain *= self.xp.exp(1j * self.pPrec.zeta * 2)
####################################
# Get polarizations from strain #
####################################
# Polarizations from complex strain
hp = self.xp.real(strain)
hc = -self.xp.imag(strain)
# Rotate polarizations by polarization_angle
if self.pWF.polarization_angle != 0:
hp, hc = rotate_by_polarization_angle(hp, hc, self.pWF.polarization_angle)
# Condition polarizations
if self.pWF.condition:
hp, hc, times = self.condition_polarizations(hp, hc)
return hp, hc, times
[docs]
class IMRPhenomTP(IMRPhenomTPHM):
"""
Class for the IMRPhenomTP model :cite:`phenomtphm`.
Precessing model with only the 22 mode in the co-precessing frame.
Wrapper to the ``IMRPhenomTPHM`` class called with only the 22, 2-2 modes.
"""
def __init__(self, **kwargs):
not_allowed_arguments = ["mode_array", "add_20_mode"]
for key in not_allowed_arguments:
if key in kwargs:
raise KeyError(f"{key} argument not allowed for IMRPhenomTP")
super().__init__(**kwargs, mode_array=[[2, 2], [2, -2]])