Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimIMRPhenomD_NRTidal.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015 Sebastian Khan, Michael Puerrer
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with with program; see the file COPYING. If not, write to the
16 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17 * MA 02110-1301 USA
18 */
19
20#include <math.h>
21#include <lal/LALSimIMR.h>
22#include <lal/FrequencySeries.h>
23#include <lal/Sequence.h>
24#include <lal/Units.h>
25#include <lal/LALConstants.h>
26
27
28#ifndef _OPENMP
29#define omp ignore
30#endif
31
32
33// internal function
35 COMPLEX16FrequencySeries **htilde, /**< Output: Frequency-domain waveform h+ */
36 REAL8 phiRef, /**< Phase at reference time */
37 REAL8 fRef, /**< Reference frequency (Hz); 0 defaults to fLow */
38 REAL8 distance, /**< Distance of source (m) */
39 REAL8 m1_SI, /**< Mass of neutron star 1 (kg) */
40 REAL8 m2_SI, /**< Mass of neutron star 2 (kg) */
41 REAL8 chi1, /**< Dimensionless aligned component spin of NS 1 */
42 REAL8 chi2, /**< Dimensionless aligned component spin of NS 2 */
43 REAL8 lambda1, /**< Dimensionless tidal deformability of NS 1 */
44 REAL8 lambda2, /**< Dimensionless tidal deformability of NS 2 */
45 LALDict *extraParams, /**< linked list containing the extra testing GR parameters */
46 const REAL8Sequence *freqs_in, /**< Frequency points at which to evaluate the waveform (Hz) */
47 REAL8 deltaF, /**< Sampling frequency (Hz) */
48 NRTidal_version_type NRTidal_version /**< Version of NRTides; can be any one of NRTidal_V (arXiv:1706.02969), NRTidalv2_V (arXiv:1905.06011) or NRTidalv2NoAmpCorr_V (arXiv:1905.06011, without amplitude corrections) */
49);
50
51// Implementation //////////////////////////////////////////////////////////////
52
54 COMPLEX16FrequencySeries **htilde, /**< Output: Frequency-domain waveform h+ */
55 REAL8 phiRef, /**< Phase at reference time */
56 REAL8 fRef, /**< Reference frequency (Hz); 0 defaults to fLow */
57 REAL8 distance, /**< Distance of source (m) */
58 REAL8 m1_SI, /**< Mass of neutron star 1 (kg) */
59 REAL8 m2_SI, /**< Mass of neutron star 2 (kg) */
60 REAL8 chi1, /**< Dimensionless aligned component spin of NS 1 */
61 REAL8 chi2, /**< Dimensionless aligned component spin of NS 2 */
62 REAL8 lambda1, /**< Dimensionless tidal deformability of NS 1 */
63 REAL8 lambda2, /**< Dimensionless tidal deformability of NS 2 */
64 LALDict *extraParams, /**< linked list containing the extra testing GR parameters */
65 const REAL8Sequence *freqs_in, /**< Frequency points at which to evaluate the waveform (Hz) */
66 REAL8 deltaF, /**< Sampling frequency (Hz) */
67 NRTidal_version_type NRTidal_version /**< Version of NRTides; can be one of NRTidal or NRTidalv2NoAmpCorr */ )
68{
69 /* Check output arrays */
70 if(!htilde) XLAL_ERROR(XLAL_EFAULT);
71 if(*htilde) {
72 XLALPrintError("(*htilde) is supposed to be NULL, but got %p",(*htilde));
74 }
75
76 if (!freqs_in) XLAL_ERROR(XLAL_EFAULT);
77 double fLow = freqs_in->data[0];
78 double fHigh = freqs_in->data[freqs_in->length - 1];
79 if(fRef == 0.0)
80 fRef = fLow;
81
84
85 /* Internally we need m1 > m2, so change around if this is not the case */
86 if (m1_SI < m2_SI) {
87 // Swap m1 and m2
88 double m1temp = m1_SI;
89 double chi1temp = chi1;
90 double lambda1temp = lambda1;
91 double dquadmon1temp = dquadmon1;
92 m1_SI = m2_SI;
93 chi1 = chi2;
94 m2_SI = m1temp;
95 chi2 = chi1temp;
96
97 if (lambda1 != lambda2){
98 lambda1 = lambda2;
100 lambda2 = lambda1temp;
102 }
103 if (dquadmon1 != dquadmon2) {
104 dquadmon1 = dquadmon2;
105 XLALSimInspiralWaveformParamsInsertdQuadMon1(extraParams, dquadmon1);
106 dquadmon2 = dquadmon1temp;
107 XLALSimInspiralWaveformParamsInsertdQuadMon2(extraParams, dquadmon2);
108 }
109 }
110
111 // Call IMRPhenomD. We call either the FrequencySequence version
112 // or the regular LAL version depending on how we've been called.
113
114 int ret = XLAL_SUCCESS;
115 if (deltaF > 0)
116 {
117 // if using a uniform frequency series then we only need to generate
118 // phenomD upto a bit beyond the BNS merger frequency.
119 // if asked for a frequency beyond NRTIDAL_FMAX then the
120 // returned waveform contains frequencies up to the input fHigh but
121 // only contains zeros beyond NRTIDAL_FMAX
122 double f_max_nr_tidal = fHigh;
123 /**< tidal coupling constant.*/
124 const double kappa2T = XLALSimNRTunedTidesComputeKappa2T(m1_SI, m2_SI, lambda1, lambda2);
125 /* Prepare tapering of amplitude beyond merger frequency */
126 const double fHz_mrg = XLALSimNRTunedTidesMergerFrequency( (m1_SI+m2_SI)/LAL_MSUN_SI, kappa2T, m1_SI/m2_SI);
127 const double NRTIDAL_FMAX = 1.3*fHz_mrg;
128
129 if ( ( fHigh > NRTIDAL_FMAX ) || ( fHigh == 0.0 ) )
130 {
131 // only generate upto NRTIDAL_FMAX
132 f_max_nr_tidal = NRTIDAL_FMAX;
133 }
134
136 htilde,
137 phiRef, fRef, deltaF,
138 m1_SI, m2_SI,
139 chi1, chi2,
140 fLow, f_max_nr_tidal,
141 distance,
142 extraParams, NRTidal_version);
143
144 // if uniform sampling and fHigh > NRTIDAL_FMAX then resize htilde
145 // so that it goes up to the user fHigh but is filled with zeros
146 // beyond NRTIDAL_FMAX
147 if (fHigh > NRTIDAL_FMAX)
148 {
149 // resize
150 // n_full is the next power of 2 +1.
151 size_t n_full = (size_t) pow(2,ceil(log2(fHigh / deltaF))) + 1;
152 *htilde = XLALResizeCOMPLEX16FrequencySeries(*htilde, 0, n_full);
153 XLAL_CHECK ( *htilde, XLAL_ENOMEM, "Failed to resize waveform COMPLEX16FrequencySeries");
154 }
155
156 } else {
158 htilde,
159 freqs_in,
160 phiRef, fRef,
161 m1_SI, m2_SI,
162 chi1, chi2,
163 distance,
164 extraParams, NRTidal_version);
165 }
166
167 XLAL_CHECK(XLAL_SUCCESS == ret, ret, "Failed to generate IMRPhenomD waveform.");
168
169 UINT4 offset;
170 REAL8Sequence *freqs = NULL;
171 REAL8Sequence *phi_tidal = NULL;
172 REAL8Sequence *amp_tidal = NULL;
173 REAL8Sequence *planck_taper = NULL;
174
175 /* Initialising terms for HO spin terms added in PhenomD_NRTidalv2 */
176 REAL8 SS_3p5PN = 0., SSS_3p5PN = 0.;
177 const REAL8 m1 = m1_SI / LAL_MSUN_SI;
178 const REAL8 m2 = m2_SI / LAL_MSUN_SI;
179 const REAL8 M = m1 + m2;
180 const REAL8 m_sec = M * LAL_MTSUN_SI; /* Total mass in seconds */
181 REAL8 eta = m1 * m2 / (M*M); /* Symmetric mass-ratio */
182 const REAL8 piM = LAL_PI * m_sec;
183 REAL8 X_A = m1/M;
184 REAL8 X_B = m2/M;
185 REAL8 pn_fac = 3.*pow(piM,2./3.)/(128.*eta);
186 /* End of initialising new parameters */
187
188 if (deltaF > 0) { // uniform frequencies
189 // Recreate freqs using only the lower and upper bounds
190 UINT4 iStart = (UINT4) (fLow / deltaF);
191 UINT4 iStop = (*htilde)->data->length - 1; // use the length calculated in the ROM function
192 freqs = XLALCreateREAL8Sequence(iStop - iStart);
193 if (!freqs) XLAL_ERROR(XLAL_EFUNC, "Frequency array allocation failed.");
194 for (UINT4 i=iStart; i<iStop; i++)
195 freqs->data[i-iStart] = i*deltaF;
196
197 offset = iStart;
198 }
199 else { // unequally spaced frequency sequence
200 freqs = XLALCreateREAL8Sequence(freqs_in->length);
201 if (!freqs) XLAL_ERROR(XLAL_EFUNC, "Frequency array allocation failed.");
202 for (UINT4 i=0; i<freqs_in->length; i++)
203 freqs->data[i] = freqs_in->data[i]; // just copy input
204 offset = 0;
205 }
206 COMPLEX16 *data=(*htilde)->data->data;
207
208 // Get FD tidal phase correction and amplitude factor from arXiv:1706.02969
209 amp_tidal = XLALCreateREAL8Sequence(freqs->length);
210 phi_tidal = XLALCreateREAL8Sequence(freqs->length);
211 planck_taper = XLALCreateREAL8Sequence(freqs->length);
212 if (NRTidal_version == NRTidalv2_V) {
213 ret = XLALSimNRTunedTidesFDTidalPhaseFrequencySeries(phi_tidal, amp_tidal, planck_taper, freqs, m1_SI, m2_SI, lambda1, lambda2, 0.0, 0.0, NRTidalv2NoAmpCorr_V);
214 XLAL_CHECK(XLAL_SUCCESS == ret, ret, "XLALSimNRTunedTidesFDTidalPhaseFrequencySeries Failed.");
215 XLALSimInspiralGetHOSpinTerms(&SS_3p5PN, &SSS_3p5PN, X_A, X_B, chi1, chi2, dquadmon1+1., dquadmon2+1.);
216 }
217 else {
218 XLALSimNRTunedTidesFDTidalPhaseFrequencySeries(phi_tidal, amp_tidal, planck_taper, freqs, m1_SI, m2_SI, lambda1, lambda2, 0.0, 0.0, NRTidal_version);
219 XLAL_CHECK(XLAL_SUCCESS == ret, ret, "XLALSimNRTunedTidesFDTidalPhaseFrequencySeries Failed.");
220 }
221
222 // Assemble waveform from amplitude and phase
223 for (size_t i=0; i<freqs->length; i++) { // loop over frequency points in sequence
224 int j = i + offset; // shift index for frequency series if needed
225 // Apply tidal phase correction and amplitude taper
226 double f = freqs->data[i];
227 COMPLEX16 Corr = planck_taper->data[i] * cexp(-I*phi_tidal->data[i] - I*pn_fac*(SS_3p5PN + SSS_3p5PN)*pow(f,2./3.));
228 data[j] *= Corr;
229 }
230
232 XLALDestroyREAL8Sequence(phi_tidal);
233 XLALDestroyREAL8Sequence(amp_tidal);
234 XLALDestroyREAL8Sequence(planck_taper);
235
236 return XLAL_SUCCESS;
237}
238
239
240/**
241 * @addtogroup LALSimIMRTIDAL_c
242 *
243 * @{
244 *
245 * @name IMRPhenomD_NRTidal
246 *
247 * @author Sebastian Khan, Michael Puerrer
248 *
249 * @brief C code for IMRPhenomD Husa:2015iqa, Khan:2015jqa with added
250 * tidal phase correction from arXiv:1706.02969.
251 *
252 * This is a frequency domain model that adds tidal modifications of
253 * the phasing to the IMRPhenomD model.
254 *
255 * @note Parameter ranges:
256 * * ? <= eta <= 0.25
257 * * 0 <= Lambda_i <= ?
258 * * -1 <= chi_i <= 1
259 *
260 * Aligned component spin on neutron stars.
261 * Symmetric mass-ratio eta = m1*m2/(m1+m2)^2.
262 * Total mass Mtot.
263 *
264 * @{
265 */
266
267
268/**
269 * Compute waveform in LAL format at specified frequencies for the IMRPhenomD_NRTidal
270 * tidal model based on IMRPhenomD.
271 *
272 * XLALSimIMRIMRPhenomDNRTidal() returns the plus and cross polarizations as a complex
273 * frequency series with equal spacing deltaF and contains zeros from zero frequency
274 * to the starting frequency and zeros beyond the cutoff frequency in the ringdown.
275 *
276 * In contrast, XLALSimIMRIMRPhenomDNRTidalFrequencySequence() returns a
277 * complex frequency series with entries exactly at the frequencies specified in
278 * the sequence freqs (which can be unequally spaced). No zeros are added.
279 *
280 * If XLALSimIMRIMRPhenomDNRTidalFrequencySequence() is called with frequencies that
281 * are beyond the maxium allowed geometric frequency for the ROM, zero strain is returned.
282 * It is not assumed that the frequency sequence is ordered.
283 *
284 * This function is designed as an entry point for reduced order quadratures.
285 */
287 COMPLEX16FrequencySeries **htilde, /**< Output: Frequency-domain waveform h+ */
288 const REAL8Sequence *freqs, /**< Frequency points at which to evaluate the waveform (Hz) */
289 REAL8 phiRef, /**< Phase at reference time */
290 REAL8 fRef, /**< Reference frequency (Hz); 0 defaults to fLow */
291 REAL8 distance, /**< Distance of source (m) */
292 REAL8 m1_SI, /**< Mass of neutron star 1 (kg) */
293 REAL8 m2_SI, /**< Mass of neutron star 2 (kg) */
294 REAL8 chi1, /**< Dimensionless aligned component spin of NS 1 */
295 REAL8 chi2, /**< Dimensionless aligned component spin of NS 2 */
296 REAL8 lambda1, /**< Dimensionless tidal deformability of NS 1 */
297 REAL8 lambda2, /**< Dimensionless tidal deformability of NS 2 */
298 LALDict *extraParams, /**< linked list containing the extra testing GR parameters */
299 NRTidal_version_type NRTidal_version /**< Version of NRTides; can be any one of NRTidal_V (arXiv:1706.02969), NRTidalv2_V (arXiv:1905.06011) or NRTidalv2NoAmpCorr_V (arXiv:1905.06011, without amplitude corrections) */
300) {
301 if (!freqs) XLAL_ERROR(XLAL_EFAULT);
302
303 // Call the internal core function with deltaF = 0 to indicate that freqs is non-uniformly
304 // spaced and we want the strain only at these frequencies
305 int retcode = IMRPhenomD_NRTidal_Core(htilde,
306 phiRef, fRef, distance, m1_SI, m2_SI, chi1, chi2, lambda1, lambda2, extraParams, freqs, 0, NRTidal_version);
307
308 return(retcode);
309}
310
311
312/**
313 * Compute waveform in LAL format for the IMRPhenomD_NRTidal
314 * tidal model based on IMRPhenomD.
315 *
316 * Returns the plus and cross polarizations as a complex frequency series with
317 * equal spacing deltaF and contains zeros from zero frequency to the starting
318 * frequency fLow and zeros beyond the cutoff frequency in the ringdown.
319 */
321 COMPLEX16FrequencySeries **htilde, /**< Output: Frequency-domain waveform h+ */
322 REAL8 phiRef, /**< Phase at reference time */
323 REAL8 deltaF, /**< Sampling frequency (Hz) */
324 REAL8 fLow, /**< Starting GW frequency (Hz) */
325 REAL8 fHigh, /**< End frequency; 0 defaults to Mf=0.14 */
326 REAL8 fRef, /**< Reference frequency (Hz); 0 defaults to fLow */
327 REAL8 distance, /**< Distance of source (m) */
328 REAL8 m1_SI, /**< Mass of neutron star 1 (kg) */
329 REAL8 m2_SI, /**< Mass of neutron star 2 (kg) */
330 REAL8 chi1, /**< Dimensionless aligned component spin of NS 1 */
331 REAL8 chi2, /**< Dimensionless aligned component spin of NS 2 */
332 REAL8 lambda1, /**< Dimensionless tidal deformability of NS 1 */
333 REAL8 lambda2, /**< Dimensionless tidal deformability of NS 2 */
334 LALDict *extraParams, /**< linked list containing the extra testing GR parameters */
335 NRTidal_version_type NRTidal_version /**< Version of NRTides; can be any one of NRTidal_V (arXiv:1706.02969), NRTidalv2_V (arXiv:1905.06011) or NRTidalv2NoAmpCorr_V (arXiv:1905.06011, without amplitude corrections) */
336) {
337 // Use fLow, fHigh, deltaF to compute freqs sequence
338 // Instead of building a full sequence we only transfer the boundaries and let
339 // the internal core function do the rest (and properly take care of corner cases).
341 freqs->data[0] = fLow;
342 freqs->data[1] = fHigh;
343
344 int retcode = IMRPhenomD_NRTidal_Core(htilde,
345 phiRef, fRef, distance, m1_SI, m2_SI, chi1, chi2, lambda1, lambda2, extraParams, freqs, deltaF, NRTidal_version);
346
348
349 return(retcode);
350}
351
352/** @} */
353
354/** @} */
void XLALSimInspiralGetHOSpinTerms(REAL8 *SS_3p5PN, REAL8 *SSS_3p5PN, REAL8 X_A, REAL8 X_B, REAL8 chi1, REAL8 chi2, REAL8 quadparam1, REAL8 quadparam2)
Function to add 3.5PN spin-squared and 3.5PN spin-cubed terms.
int XLALSimNRTunedTidesFDTidalPhaseFrequencySeries(const REAL8Sequence *phi_tidal, const REAL8Sequence *amp_tidal, const REAL8Sequence *planck_taper, const REAL8Sequence *fHz, REAL8 m1_SI, REAL8 m2_SI, REAL8 lambda1, REAL8 lambda2, REAL8 chi1, REAL8 chi2, NRTidal_version_type NRTidal_version)
Function to call the frequency domain tidal correction over an array of input frequencies.
double XLALSimNRTunedTidesMergerFrequency(const REAL8 mtot_MSUN, const REAL8 kappa2T, const REAL8 q)
compute the merger frequency of a BNS system.
double XLALSimNRTunedTidesComputeKappa2T(REAL8 m1_SI, REAL8 m2_SI, REAL8 lambda1, REAL8 lambda2)
convenient function to compute tidal coupling constant.
static int IMRPhenomD_NRTidal_Core(COMPLEX16FrequencySeries **htilde, REAL8 phiRef, REAL8 fRef, REAL8 distance, REAL8 m1_SI, REAL8 m2_SI, REAL8 chi1, REAL8 chi2, REAL8 lambda1, REAL8 lambda2, LALDict *extraParams, const REAL8Sequence *freqs_in, REAL8 deltaF, NRTidal_version_type NRTidal_version)
REAL8 XLALSimInspiralWaveformParamsLookupdQuadMon1(LALDict *params)
int XLALSimInspiralWaveformParamsInsertTidalLambda1(LALDict *params, REAL8 value)
int XLALSimInspiralWaveformParamsInsertTidalLambda2(LALDict *params, REAL8 value)
REAL8 XLALSimInspiralWaveformParamsLookupdQuadMon2(LALDict *params)
int XLALSimInspiralWaveformParamsInsertdQuadMon2(LALDict *params, REAL8 value)
int XLALSimInspiralWaveformParamsInsertdQuadMon1(LALDict *params, REAL8 value)
REAL8 M
Definition: bh_qnmode.c:133
double i
Definition: bh_ringdown.c:118
sigmaKerr data[0]
COMPLEX16FrequencySeries * XLALResizeCOMPLEX16FrequencySeries(COMPLEX16FrequencySeries *series, int first, size_t length)
#define LAL_MSUN_SI
#define LAL_PI
#define LAL_MTSUN_SI
double complex COMPLEX16
double REAL8
uint32_t UINT4
NRTidal_version_type
Definition: LALSimIMR.h:80
@ NRTidalv2NoAmpCorr_V
version NRTidalv2, without amplitude corrections
Definition: LALSimIMR.h:84
@ NRTidalv2_V
version NRTidalv2: https://arxiv.org/abs/1905.06011
Definition: LALSimIMR.h:82
int XLALSimIMRPhenomDGenerateFD(COMPLEX16FrequencySeries **htilde, const REAL8 phi0, const REAL8 fRef, const REAL8 deltaF, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi1, const REAL8 chi2, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, LALDict *extraParams, NRTidal_version_type NRTidal_version)
Driver routine to compute the spin-aligned, inspiral-merger-ringdown phenomenological waveform IMRPhe...
int XLALSimIMRPhenomDFrequencySequence(COMPLEX16FrequencySeries **htilde, const REAL8Sequence *freqs, const REAL8 phi0, const REAL8 fRef_in, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi1, const REAL8 chi2, const REAL8 distance, LALDict *extraParams, NRTidal_version_type NRTidal_version)
Compute waveform in LAL format at specified frequencies for the IMRPhenomD model.
int XLALSimIMRPhenomDNRTidalFrequencySequence(COMPLEX16FrequencySeries **htilde, const REAL8Sequence *freqs, REAL8 phiRef, REAL8 fRef, REAL8 distance, REAL8 m1_SI, REAL8 m2_SI, REAL8 chi1, REAL8 chi2, REAL8 lambda1, REAL8 lambda2, LALDict *extraParams, NRTidal_version_type NRTidal_version)
Compute waveform in LAL format at specified frequencies for the IMRPhenomD_NRTidal tidal model based ...
int XLALSimIMRPhenomDNRTidal(COMPLEX16FrequencySeries **htilde, REAL8 phiRef, REAL8 deltaF, REAL8 fLow, REAL8 fHigh, REAL8 fRef, REAL8 distance, REAL8 m1_SI, REAL8 m2_SI, REAL8 chi1, REAL8 chi2, REAL8 lambda1, REAL8 lambda2, LALDict *extraParams, NRTidal_version_type NRTidal_version)
Compute waveform in LAL format for the IMRPhenomD_NRTidal tidal model based on IMRPhenomD.
void XLALDestroyREAL8Sequence(REAL8Sequence *sequence)
REAL8Sequence * XLALCreateREAL8Sequence(size_t length)
#define XLAL_ERROR(...)
#define XLAL_CHECK(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFUNC
REAL8 * data