From 1d5845de477cbb658af5a4c6b1de8b25ae544d52 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Thu, 20 Aug 2026 13:48:58 -0400 Subject: [PATCH] Propagate tidal parameters through GWSignal --- .../Code/RIFT/physics/GWSignal.py | 95 +++++-------- .../Code/test/waveform/test_gwsignal_tides.py | 127 ++++++++++++++++++ 2 files changed, 162 insertions(+), 60 deletions(-) create mode 100644 MonteCarloMarginalizeCode/Code/test/waveform/test_gwsignal_tides.py diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/GWSignal.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/GWSignal.py index d6998bfa6..5d3e40387 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/GWSignal.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/GWSignal.py @@ -34,6 +34,38 @@ print("GWsignal import failed") +def _waveform_parameter_dict(P, taper, **kwargs): + """Translate RIFT waveform parameters to the GWSignal convention. + + Tidal deformabilities are dimensionless in both RIFT and GWSignal. Keep + them in this common constructor so mode and polarization generation cannot + silently disagree about whether matter effects are enabled. + """ + python_dict = {'mass1' : P.m1/lal.MSUN_SI * u.solMass, + 'mass2' : P.m2/lal.MSUN_SI * u.solMass, + 'spin1x' : P.s1x*u.dimensionless_unscaled, + 'spin1y' : P.s1y*u.dimensionless_unscaled, + 'spin1z' : P.s1z*u.dimensionless_unscaled, + 'spin2x' : P.s2x*u.dimensionless_unscaled, + 'spin2y' : P.s2y*u.dimensionless_unscaled, + 'spin2z' : P.s2z*u.dimensionless_unscaled, + 'lambda1' : P.lambda1*u.dimensionless_unscaled, + 'lambda2' : P.lambda2*u.dimensionless_unscaled, + 'deltaT' : P.deltaT*u.s, + 'f22_start' : P.fmin*u.Hz, + 'f22_ref': P.fref*u.Hz, + 'phi_ref' : P.phiref*u.rad, + 'distance' : P.dist/(1e6*lal.PC_SI)*u.Mpc, + 'inclination' : P.incl*u.rad, + 'eccentricity' : P.eccentricity*u.dimensionless_unscaled, + 'longAscNodes' : P.psi*u.rad, + 'meanPerAno' : P.meanPerAno*u.rad, + 'condition' : taper} + if 'lmax_nyquist' in kwargs: + python_dict['lmax_nyquist'] = kwargs['lmax_nyquist'] + return python_dict + + def std_and_conj_hlmoff(P, Lmax=2,approx_string=None,**kwargs): hlms = hlmoft(P, Lmax,approx_string=approx_string,**kwargs) hlmsF = {} @@ -64,26 +96,7 @@ def hlmoft(P, Lmax=2,approx_string=None,no_trust_align_method=None,internal_phas taper=0 if P.taper != lalsim.SIM_INSPIRAL_TAPER_NONE: taper = 1 - python_dict = {'mass1' : P.m1/lal.MSUN_SI * u.solMass, - 'mass2' : P.m2/lal.MSUN_SI * u.solMass, - 'spin1x' : P.s1x*u.dimensionless_unscaled, - 'spin1y' : P.s1y*u.dimensionless_unscaled, - 'spin1z' : P.s1z*u.dimensionless_unscaled, - 'spin2x' : P.s2x*u.dimensionless_unscaled, - 'spin2y' : P.s2y*u.dimensionless_unscaled, - 'spin2z' : P.s2z*u.dimensionless_unscaled, - 'deltaT' : P.deltaT*u.s, - 'f22_start' : P.fmin*u.Hz, - 'f22_ref': P.fref*u.Hz, - 'phi_ref' : P.phiref*u.rad, - 'distance' : P.dist/(1e6*lal.PC_SI)*u.Mpc, - 'inclination' : P.incl*u.rad, - 'eccentricity' : P.eccentricity*u.dimensionless_unscaled, - 'longAscNodes' : P.psi*u.rad, - 'meanPerAno' : P.meanPerAno*u.rad, - 'condition' : taper } - if 'lmax_nyquist' in kwargs: - python_dict['lmax_nyquist'] = kwargs['lmax_nyquist'] + python_dict = _waveform_parameter_dict(P, taper, **kwargs) # if needed # lal_dict = gws.core.utils.to_lal_dict(python_dict) @@ -192,26 +205,7 @@ def hoft(P, Fp=None, Fc=None,approx_string=None, **kwargs): taper=0 if P.taper != lalsim.SIM_INSPIRAL_TAPER_NONE: taper = 1 - python_dict = {'mass1' : P.m1/lal.MSUN_SI * u.solMass, - 'mass2' : P.m2/lal.MSUN_SI * u.solMass, - 'spin1x' : P.s1x*u.dimensionless_unscaled, - 'spin1y' : P.s1y*u.dimensionless_unscaled, - 'spin1z' : P.s1z*u.dimensionless_unscaled, - 'spin2x' : P.s2x*u.dimensionless_unscaled, - 'spin2y' : P.s2y*u.dimensionless_unscaled, - 'spin2z' : P.s2z*u.dimensionless_unscaled, - 'deltaT' : P.deltaT*u.s, - 'f22_start' : P.fmin*u.Hz, - 'f22_ref': P.fref*u.Hz, - 'phi_ref' : P.phiref*u.rad, - 'distance' : P.dist/(1e6*lal.PC_SI)*u.Mpc, - 'inclination' : P.incl*u.rad, - 'eccentricity' : P.eccentricity*u.dimensionless_unscaled, - 'longAscNodes' : P.psi*u.rad, - 'meanPerAno' : P.meanPerAno*u.rad, - 'condition' : taper} - if 'lmax_nyquist' in kwargs: - python_dict['lmax_nyquist'] = kwargs['lmax_nyquist'] + python_dict = _waveform_parameter_dict(P, taper, **kwargs) # if needed # lal_dict = gws.core.utils.to_lal_dict(python_dict) @@ -293,26 +287,7 @@ def complex_hoft(P, Fp=None, Fc=None,approx_string=None,sgn=-1, **kwargs): taper=0 if P.taper != lalsim.SIM_INSPIRAL_TAPER_NONE: taper = 1 - python_dict = {'mass1' : P.m1/lal.MSUN_SI * u.solMass, - 'mass2' : P.m2/lal.MSUN_SI * u.solMass, - 'spin1x' : P.s1x*u.dimensionless_unscaled, - 'spin1y' : P.s1y*u.dimensionless_unscaled, - 'spin1z' : P.s1z*u.dimensionless_unscaled, - 'spin2x' : P.s2x*u.dimensionless_unscaled, - 'spin2y' : P.s2y*u.dimensionless_unscaled, - 'spin2z' : P.s2z*u.dimensionless_unscaled, - 'deltaT' : P.deltaT*u.s, - 'f22_start' : P.fmin*u.Hz, - 'f22_ref': P.fref*u.Hz, - 'phi_ref' : P.phiref*u.rad, - 'distance' : P.dist/(1e6*lal.PC_SI)*u.Mpc, - 'inclination' : P.incl*u.rad, - 'eccentricity' : P.eccentricity*u.dimensionless_unscaled, - 'longAscNodes' : P.psi*u.rad, - 'meanPerAno' : P.meanPerAno*u.rad, - 'condition' : taper} - if 'lmax_nyquist' in kwargs: - python_dict['lmax_nyquist'] = kwargs['lmax_nyquist'] + python_dict = _waveform_parameter_dict(P, taper, **kwargs) # if needed # lal_dict = gws.core.utils.to_lal_dict(python_dict) diff --git a/MonteCarloMarginalizeCode/Code/test/waveform/test_gwsignal_tides.py b/MonteCarloMarginalizeCode/Code/test/waveform/test_gwsignal_tides.py new file mode 100644 index 000000000..09f0b17ec --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/waveform/test_gwsignal_tides.py @@ -0,0 +1,127 @@ +"""Regression tests for tidal parameters in RIFT's GWSignal interface.""" + +import numpy as np +import pytest + +import lal +import astropy.units as u + +import RIFT.lalsimutils as lalsimutils +import RIFT.physics.GWSignal as rift_gws + + +def _tidal_params(lambda1=500.0, lambda2=800.0, s1z=0.05, s2z=-0.03): + return lalsimutils.ChooseWaveformParams( + m1=1.4 * lal.MSUN_SI, + m2=1.3 * lal.MSUN_SI, + s1z=s1z, + s2z=s2z, + lambda1=lambda1, + lambda2=lambda2, + deltaT=1.0 / 4096.0, + fmin=40.0, + fref=40.0, + dist=100.0e6 * lal.PC_SI, + taper=lalsimutils.lsu_TAPER_NONE, + ) + + +def _normalized_tail_difference(first, second): + first = np.asarray(first) + second = np.asarray(second) + assert np.all(np.isfinite(first)) + assert np.all(np.isfinite(second)) + + n_common = min(first.size, second.size) + first = first[-n_common:] + second = second[-n_common:] + first_norm = np.linalg.norm(first) + second_norm = np.linalg.norm(second) + assert first_norm > 0.0 + assert second_norm > 0.0 + return np.linalg.norm(first / first_norm - second / second_norm) + + +def test_waveform_parameter_dict_propagates_tidal_deformabilities(): + params = rift_gws._waveform_parameter_dict( + _tidal_params(), taper=0, lmax_nyquist=4 + ) + + assert params["lambda1"].unit == u.dimensionless_unscaled + assert params["lambda2"].unit == u.dimensionless_unscaled + assert params["lambda1"].value == pytest.approx(500.0) + assert params["lambda2"].value == pytest.approx(800.0) + assert params["lmax_nyquist"] == 4 + + +def test_public_tidal_waveform_responds_to_rift_lambdas(): + """Run a public, valid tidal waveform through the GWSignal boundary.""" + if not rift_gws.has_gws: + pytest.skip("GWSignal is unavailable") + + tidal = rift_gws.hoft( + _tidal_params(s1z=0.0, s2z=0.0), + Fp=1.0, + Fc=0.0, + approx_string="TaylorT4", + ) + point_mass = rift_gws.hoft( + _tidal_params(lambda1=0.0, lambda2=0.0, s1z=0.0, s2z=0.0), + Fp=1.0, + Fc=0.0, + approx_string="TaylorT4", + ) + + assert ( + _normalized_tail_difference(tidal.data.data, point_mass.data.data) > 1.0e-6 + ) + + +def test_seobnrv5thm_gwsignal_modes_respond_to_tides(): + """Exercise a real mode-generating tidal model when its review build exists. + + The test intentionally skips in ordinary environments until SEOBNRv5THM is + exposed as a GWSignal model. In the review/production environment it guards + against the dangerous failure mode where nonzero RIFT lambdas produce the + same modes as the black-hole limit. + """ + if not rift_gws.has_gws: + pytest.skip("GWSignal is unavailable") + + try: + tidal_generator = rift_gws.gws.models.gwsignal_get_waveform_generator( + "SEOBNRv5THM" + ) + except ValueError as exc: + if "Approximant not implemented in GWSignal" not in str(exc): + raise + pytest.skip("SEOBNRv5THM is not exposed through GWSignal: {}".format(exc)) + + metadata = tidal_generator.metadata + if callable(metadata): + metadata = metadata() + assert metadata["modes"], "SEOBNRv5THM must expose native modes through GWSignal" + + tidal_dict = rift_gws._waveform_parameter_dict(_tidal_params(), taper=0) + point_mass_dict = rift_gws._waveform_parameter_dict( + _tidal_params(lambda1=0.0, lambda2=0.0), taper=0 + ) + # Use distinct generators so the comparison cannot pass or fail because a + # plugin retains waveform state between calls. + point_mass_generator = rift_gws.gws.models.gwsignal_get_waveform_generator( + "SEOBNRv5THM" + ) + tidal_modes = rift_gws.wfm.GenerateTDModes(tidal_dict, tidal_generator) + point_mass_modes = rift_gws.wfm.GenerateTDModes( + point_mass_dict, point_mass_generator + ) + + assert (2, 2) in tidal_modes + assert (2, 2) in point_mass_modes + # End-align, where tidal dephasing is largest, and normalize out distance. + assert ( + _normalized_tail_difference( + tidal_modes[(2, 2)].value, point_mass_modes[(2, 2)].value + ) + > 1.0e-6 + )