Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-3a66518
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
test_SEOBNRv4PHM_vs_4HM_ringdown.py
Go to the documentation of this file.
1# -*- coding: utf-8 -*-
2#
3# Copyright (C) 2020 Serguei Ossokine
4#
5# This program is free software: you can redistribute it and/or modify
6# it under the terms of the GNU General Public License as published by
7# the Free Software Foundation, either version 3 of the License, or
8# (at your option) any later version.
9#
10# This program is distributed in the hope that it will be useful,
11# but WITHOUT ANY WARRANTY; without even the implied warranty of
12# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13# GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License
16# along with this program. If not, see <http: //www.gnu.org/licenses/>.
17
18"""Test that the aligned-spin limit of SEOBNRv4PHM is the same as SEOBNRv4HM
19In particular we consider cases where the final spin is *negative* with respect to the inital
20 orbital angular momentum. """
21
22import sys
23import pytest
24from scipy.interpolate import InterpolatedUnivariateSpline
25import lalsimulation as ls
26import lal
27import numpy as np
28
29# -- test functions ---------------------
30
31
32def get_SEOBNRv4HM_modes(q, M, chi1, chi2, f_start, distance, deltaT):
33 """Generate SEOBNRv4HM modes"""
34 m1SI = lal.MSUN_SI * q * M / (1.0 + q)
35 m2SI = lal.MSUN_SI * M / (1.0 + q)
36 nqcCoeffsInput = lal.CreateREAL8Vector(10)
37
38 PAparams = lal.CreateDict()
39 lal.DictInsertUINT4Value(PAparams, "PAFlag", 0)
40
41 TGRparams = lal.CreateDict()
42 lal.DictInsertREAL8Value(TGRparams, "domega220", 0.)
43 lal.DictInsertREAL8Value(TGRparams, "dtau220", 0.)
44 lal.DictInsertREAL8Value(TGRparams, "domega210", 0.)
45 lal.DictInsertREAL8Value(TGRparams, "dtau210", 0.)
46 lal.DictInsertREAL8Value(TGRparams, "domega330", 0.)
47 lal.DictInsertREAL8Value(TGRparams, "dtau330", 0.)
48 lal.DictInsertREAL8Value(TGRparams, "domega440", 0.)
49 lal.DictInsertREAL8Value(TGRparams, "dtau440", 0.)
50 lal.DictInsertREAL8Value(TGRparams, "domega550", 0.)
51 lal.DictInsertREAL8Value(TGRparams, "dtau550", 0.)
52
53 sphtseries, dyn, dynHI = ls.SimIMRSpinAlignedEOBModes(
54 deltaT,
55 m1SI,
56 m2SI,
57 f_start,
58 distance,
59 chi1,
60 chi2,
61 41,
62 0.0,
63 0.0,
64 0.0,
65 0.0,
66 0.0,
67 0.0,
68 0.0,
69 0.0,
70 1.0,
71 1.0,
72 nqcCoeffsInput,
73 0, # nqcFlag
74 PAparams,
75 TGRparams,
76 )
77
78 # The minus sign in front of the modes takes into account the fact that the polarization basis in EOB
79 # conventions is different wrt the one in LAL
80 hI = {}
81 modes = [(2, 2), (2, 1), (3, 3), (4, 4), (5, 5)]
82 for lm in modes:
83 hI[lm] = np.trim_zeros(
84 -1 * ls.SphHarmTimeSeriesGetMode(sphtseries, lm[0], lm[1]).data.data, "b"
85 )
86
87 return hI
88
89
90def get_SEOBNRv4PHM_modes(q, M, chi1, chi2, f_start, distance, deltaT):
91 """Generate SEOBNRv4PHM modes"""
92
93 m1SI = lal.MSUN_SI * q * M / (1.0 + q)
94 m2SI = lal.MSUN_SI * M / (1.0 + q)
95 approx = ls.SEOBNRv4PHM
96 hlm = ls.SimInspiralChooseTDModes(0.,
97 deltaT,
98 m1SI,
99 m2SI,
100 chi1[0],
101 chi1[1],
102 chi1[2],
103 chi2[0],
104 chi2[1],
105 chi2[2],
106 f_start,
107 f_start,
108 distance,
109 None,
110 5,
111 approx,
112 )
113 hI = {}
114 modes = [(2, 2), (2, 1), (3, 3), (4, 4), (5, 5)]
115 for lm in modes:
116 hI[lm] = ls.SphHarmTimeSeriesGetMode(hlm, lm[0], lm[1]).data.data
117
118 return hI
119
120
121@pytest.mark.parametrize("q", [5.0, 10.0, 15.0, 20.0])
123 """ Test that the (2,2) mode frequency in the ringdown is within in 1Hz between v4HM and v4PHM.
124 In particular, consider cases where the final spin is negative with respect to the initial
125 orbital angular momentum. Of course considers only aligned-spin cases"""
126 # Step 0: set parameters
127
128 chi1_x, chi1_y, chi1_z = 0.0, 0.0, -0.95
129 chi2_x, chi2_y, chi2_z = 0.0, 0.0, -0.95
130 m_total = 150 # Feducial total mass in solar masses
131 deltaT = 1.0 / 16384
132 distance = 500 * 1e6 * lal.PC_SI # 500 Mpc in m
133 f_start22 = 20
134
135 # Step 1: Get the SEOBNRv4PHM modes
137 q,
138 m_total,
139 [chi1_x, chi2_y, chi1_z],
140 [chi1_x, chi2_y, chi2_z],
141 f_start22 * 0.5,
142 distance,
143 deltaT,
144 )
145
146 # Step 2: Get the SEOBNRv4HM modes
148 q, m_total, chi1_z, chi2_z, f_start22 * 0.5, distance, deltaT
149 )
150
151 time = deltaT * np.arange(len(hlmA[(2, 2)]))
152 time2 = deltaT * np.arange(len(hlmP[(2, 2)]))
153 # Compute the (2,2) mode frequencies
154 mode = (2, 2)
155
156 phaseA = np.unwrap(np.angle(hlmA[mode]))
157 phaseP = np.unwrap(np.angle(hlmP[mode]))
158
159 intrpA = InterpolatedUnivariateSpline(time, phaseA)
160 intrpP = InterpolatedUnivariateSpline(time2, phaseP)
161
162 omega_A = intrpA.derivative()
163 omega_P = intrpP.derivative()
164
165 # Consider the frequency in the ringdown
166 amp = np.abs(hlmA[mode].data)
167 am = np.argmax(amp)
168
169 omega_A_rd = omega_A(time[am : am + 350])
170 omega_P_rd = omega_P(time[am : am + 350])
171
172 # Ensure the difference is small. Note here atol=1 corresponds to 1 Hz difference on frequncies of hundreds of Hz. These small
173 # diffrences are expected due to some minor code differences.
174 # Before a bug was fixed the difference would be huge, i.e. ~100s Hz
175 assert np.allclose(
176 omega_A_rd, omega_P_rd, atol=1
177 ), "The frequencies of the (2,2) mode don't agree between SEOBNRv4PHM and SEOBNRv4HM!"
178
179
180# -- run the tests ------------------------------
181
182if __name__ == "__main__":
183 args = sys.argv[1:] or ["-v", "-rs", "--junit-xml=junit-SEOBNRv4PHM_vs_4HM_ringdown.xml"]
184 sys.exit(pytest.main(args=[__file__] + args))
def get_SEOBNRv4PHM_modes(q, M, chi1, chi2, f_start, distance, deltaT)
Generate SEOBNRv4PHM modes.
def get_SEOBNRv4HM_modes(q, M, chi1, chi2, f_start, distance, deltaT)
Generate SEOBNRv4HM modes.
def test_aligned_spin_limit_ringdown(q)
Test that the (2,2) mode frequency in the ringdown is within in 1Hz between v4HM and v4PHM.