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
LALSimInspiralTestingGRCorrections.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2017 Walter Del Pozzo
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#include <gsl/gsl_spline.h>
20#include <stdlib.h>
21#include <math.h>
22#include <lal/Date.h>
23#include <lal/FrequencySeries.h>
24#include <lal/LALConstants.h>
25#include <lal/LALDatatypes.h>
26#include <lal/LALSimInspiral.h>
27#include <lal/LALSimInspiralTestingGRCorrections.h>
28#include <lal/Sequence.h>
30#include <lal/LALSimIMR.h>
31
32/*
33 * Peak (angular) frequency of 2,2 mode of GW predicted by fitting NR results. Taken from
34 * Eq. A8 of Bohe et al (2017) (arxiv.org/abs/1611.03703). Used in building SEOBNRv4.
35 */
36static inline REAL8 Get22PeakAngFreq(REAL8 UNUSED eta,
37 REAL8 a) /** Combined spin chi defined in Eq. (2.8) of Bohe et al. (2017) */
38{
39 REAL8 chi = a;
40 REAL8 res;
41 res = 0.5626787200433265 + (-0.08706198756945482 +
42 0.0017434519312586804 * chi) *
43 log (10.26207326082448 -
44 chi * (7.629921628648589 -
45 72.75949266353584 * (-0.25 + eta)) -
46 62.353217004599784 * (-0.25 + eta));
47 return res;
48}
49
50int XLALSimInspiralTestingGRCorrections(COMPLEX16FrequencySeries *htilde, /**< input htilde, will be modified in place */
51 const UINT4 l, /**< degree l of the GW mode */
52 const UINT4 m, /**< order m of the GW mode */
53 const REAL8 m1_SI, /**< mass of companion 1 (kg) */
54 const REAL8 m2_SI, /**< mass of companion 2 (kg) */
55 const REAL8 chi1z, /**< z-component of the dimensionless spin of object 1 */
56 const REAL8 chi2z, /**< z-component of the dimensionless spin of object 2 */
57 const REAL8 f_low, /**< lower GW frequency bound (Hz) */
58 const REAL8 f_ref, /**< reference GW frequency (Hz) */
59 const REAL8 f_window_div_f_Peak, /**< Frequency at which to attach non-GR and GR waveforms, inputted as a fraction of f_Peak (should be between 0 and 1) */
60 const REAL8 NCyclesStep, /**< Number of GW cycles over which to taper the non-GR phase correction */
61 LALDict *LALpars /**< dictionary of additional parameters */
62 )
63{
64 /* check if we have a NULL pnCorrections pointer. If yes, just return */
65 if (LALpars==NULL) return 0;
66 /* external: SI; internal: solar masses */
67 const REAL8 f0 = htilde->f0;
68 const REAL8 deltaF = htilde->deltaF;
69 const REAL8 m1 = m1_SI / LAL_MSUN_SI;
70 const REAL8 m2 = m2_SI / LAL_MSUN_SI;
71 const REAL8 mt = m1 + m2;
72 const REAL8 m_sec = mt * LAL_MTSUN_SI; /* total mass in seconds */
73 const REAL8 eta = m1 * m2 / (mt * mt);
74/* Check that fRef is sane */
75 if( f_ref < f_low )
76 {
77 XLALPrintError("XLAL Error - %s: fRef is smaller than the starting frequency of the waveform, f_low. Please pass in the starting GW frequency instead.\n", __func__);
79 }
80
83
84 REAL8 f22Peak;
85
86 /* Compute the frequency where the amplitude of 2,2 mode peaks differently for BBH, BNS, and NSBH:
87 For BBH, use fit to NR used in construction of SEOBNRv4 (arxiv:1611.03703)
88 For BNS, use fit to NR used in construction of NRTidal models (documented in LALSimNRTunedTides.c, arXiv:1706.02969)
89 For NSBH, use same fit as BBH
90 */
91
92 if(lambda1 == 0.0 && lambda2 == 0.0){
93 f22Peak = Get22PeakAngFreq(eta, 0.5*(chi1z + chi2z) + 0.5*(chi1z - chi2z)*(m1 - m2)/(m1 + m2)/(1. - 2.*eta)) / (2.*LAL_PI * m_sec);
94 //Spin combination defined in Eq. (2.8) of Bohe et al. (2017) (SEOBNRv4 paper)
95 }
96 else if(lambda1 != 0.0 && lambda2 != 0.0)
97 {
98 const REAL8 q = fmax(m1 / m2, m2 / m1);
99 const double kappa2T = XLALSimNRTunedTidesComputeKappa2T(m1_SI, m2_SI, lambda1, lambda2);
100 f22Peak = XLALSimNRTunedTidesMergerFrequency(mt, kappa2T, q) + l*0;
101 }
102 else{
103 f22Peak = Get22PeakAngFreq(eta, 0.5*(chi1z + chi2z) + 0.5*(chi1z - chi2z)*(m1 - m2)/(m1 + m2)/(1. - 2.*eta)) / (2.*LAL_PI * m_sec);
104 }
105
106 INT4 i;
107 INT4 n = (INT4) htilde->data->length;
108
109 /* Indices of f0, f_low, frequency at which non-GR modifications end, and fPeak */
110 INT4 iStart, iStartFinal, iRef, iPeak;
111 iStart = (UINT4) ceil((f_low/2.-f0) / deltaF);
112 iStartFinal = (UINT4) ceil((f_low-f0) / deltaF);
113 iRef = (UINT4) ceil((f_ref*m/2.-f0) / deltaF);
114 iPeak = (UINT4) fmin(ceil((f22Peak - f0) / deltaF),n-1);
115
116 /* Sequence of frequencies where corrections to the model need to be evaluated
117 * Fill with non-zero vals from f0 to fEnd
118 */
119 REAL8Sequence *freqs =NULL;
120 freqs = XLALCreateREAL8Sequence(n);
121
122 for (i = 0; i < n; i++)
123 {
124 freqs->data[i] = f0 + i * deltaF;
125 }
126
127 PNPhasingSeries pfa;
128 const REAL8 qm_def1 = 1.;
129 const REAL8 qm_def2 = 1.;
130
131 XLALSimInspiralPNCorrections(&pfa, m1, m2, chi1z, chi2z, chi1z*chi1z, chi2z*chi2z, chi1z*chi2z, qm_def1, qm_def2, LALpars);
132 XLALSimInspiralPhaseCorrectionsPhasing(htilde,freqs,m,iStart,iRef,iPeak,pfa,m_sec,eta,f_window_div_f_Peak,iStartFinal,NCyclesStep);
134 return 0;
135}
136
137/* Computes the PN coefficients for the non-GR phase correction and stores in pfa.
138 * The TestGRParam values represent fractional deviation from the corresponing PN
139 * coefficients in TaylorF2 expression for the phase, except where the
140 * coefficeint vanishes in GR, in which case the TestGRParam values indicates the
141 * total value of that PN coefficient. The code closely follows XLALSimInspiralPNPhasing_F2()
142 * in LALSimInspiralPNCoefficients.c, but returns only the PN coefficients for the
143 * correction instead of the TaylorF2 value + correction.
144 */
145
146void XLALSimInspiralPNCorrections(PNPhasingSeries *pfa, /**< PN phasing coefficients, to be modified in place */
147 const REAL8 m1, /**< Mass of body 1, in Msol */
148 const REAL8 m2, /**< Mass of body 2, in Msol */
149 const REAL8 chi1L, /**< Component of dimensionless spin 1 along Lhat */
150 const REAL8 chi2L, /**< Component of dimensionless spin 2 along Lhat */
151 const REAL8 chi1sq,/**< Magnitude of dimensionless spin 1 */
152 const REAL8 chi2sq, /**< Magnitude of dimensionless spin 2 */
153 const REAL8 chi1dotchi2, /**< Dot product of dimensionles spin 1 and spin 2 */
154 const REAL8 qm_def1, /**< Quadrupole deformation parameter of body 1 (dimensionless) */
155 const REAL8 qm_def2, /**< Quadrupole deformation parameter of body 2 (dimensionless) */
156 LALDict *LALpars /**< dictionary of additional parameters */
157 )
158{
159 const REAL8 mtot = m1 + m2;
160 const REAL8 eta = m1*m2/mtot/mtot;
161 const REAL8 pfaN = 3.L/(128.L * eta);
162 const REAL8 d = (m1 - m2) / (m1 + m2);
163 const REAL8 m1M = m1/mtot;
164 const REAL8 m2M = m2/mtot;
165 const REAL8 SL = m1M*m1M*chi1L + m2M*m2M*chi2L;
166 const REAL8 dSigmaL = d*(m2M*chi2L - m1M*chi1L);
167 REAL8 pn_sigma = eta * (721.L/48.L*chi1L*chi2L - 247.L/48.L*chi1dotchi2);
168 pn_sigma += (720.L*qm_def1 - 1.L)/96.0L * m1M * m1M * chi1L * chi1L;
169 pn_sigma += (720.L*qm_def2 - 1.L)/96.0L * m2M * m2M * chi2L * chi2L;
170 pn_sigma -= (240.L*qm_def1 - 7.L)/96.0L * m1M * m1M * chi1sq;
171 pn_sigma -= (240.L*qm_def2 - 7.L)/96.0L * m2M * m2M * chi2sq;
172 REAL8 pn_ss3 = (326.75L/1.12L + 557.5L/1.8L*eta)*eta*chi1L*chi2L;
173 pn_ss3 += ((4703.5L/8.4L+2935.L/6.L*m1M-120.L*m1M*m1M)*qm_def1 + (-4108.25L/6.72L-108.5L/1.2L*m1M+125.5L/3.6L*m1M*m1M)) *m1M*m1M * chi1sq;
174 pn_ss3 += ((4703.5L/8.4L+2935.L/6.L*m2M-120.L*m2M*m2M)*qm_def2 + (-4108.25L/6.72L-108.5L/1.2L*m2M+125.5L/3.6L*m2M*m2M)) *m2M*m2M * chi2sq;
175 const REAL8 pn_gamma = (554345.L/1134.L + 110.L*eta/9.L)*SL + (13915.L/84.L - 10.L*eta/3.L)*dSigmaL;
176 /* Introducing here the SIQM test PN cofficients..*/
177 /* See eq. (0.5a) of https://arxiv.org/pdf/1701.06318.pdf */
178 const REAL8 chiS = 0.5L * (chi1L + chi2L);
179 const REAL8 chiA = 0.5L * (chi1L - chi2L);
180 const REAL8 siqm_2pn_kappaS = 50.L * (2.L*eta - 1.L) * (chiS*chiS + chiA*chiA) -100.L * d * chiA * chiS;
181 const REAL8 siqm_2pn_kappaA = -50.L * d * (chiS*chiS + chiA*chiA) + 100.L * (2.L*eta - 1.L) * chiA * chiS;
182 const REAL8 siqm_3pn_kappaS = (chiS*chiS + chiA*chiA) * (26015.L/28.L - 44255.L/21.L * eta - 240.L * eta*eta) + chiA * chiS * d * (26015.L/14.L - 1495.L/3. * eta);
183 const REAL8 siqm_3pn_kappaA = (26015.L/28.L - 1495.L/6.L * eta) * d * (chiS*chiS + chiA*chiA) + (26015.L/14.L - 88510.L/21.L *eta - 480.L *eta*eta) * chiS * chiA;
184
185 /* initialise the PN correction coefficients to 0 identically */
186 memset(pfa, 0, sizeof(PNPhasingSeries));
187
188 /* Non-spinning corrections to phasing terms - see arXiv:0907.0700, Eq. 3.18 */
189 if (XLALDictContains(LALpars,"dchiMinus2")) pfa->vneg[2] = XLALSimInspiralWaveformParamsLookupNonGRDChiMinus2(LALpars);
190 if (XLALDictContains(LALpars,"dchiMinus1")) pfa->vneg[1] = XLALSimInspiralWaveformParamsLookupNonGRDChiMinus1(LALpars);
191
192 if (XLALDictContains(LALpars,"dchi0")) pfa->v[0] = XLALSimInspiralWaveformParamsLookupNonGRDChi0(LALpars);
193 if (XLALDictContains(LALpars,"dchi1")) pfa->v[1] = XLALSimInspiralWaveformParamsLookupNonGRDChi1(LALpars);
194
195 if (XLALDictContains(LALpars,"dchi2"))
196 {
197 pfa->v[2] = 5.L*(743.L/84.L + 11.L * eta)/9.L;
199 }
200
201 if (XLALDictContains(LALpars,"dchi3") || XLALDictContains(LALpars,"dchi3NS") || XLALDictContains(LALpars,"dchi3S"))
202 {
203 pfa->v[3] = -16.L*LAL_PI;
204 pfa->v[3] += 188.L*SL/3.L + 25.L*dSigmaL;
206
207 /* Splitting the deviation parameter into spin and no-spin parts and then adding them here */
208 pfa->v[3] += -16.L*LAL_PI * XLALSimInspiralWaveformParamsLookupNonGRDChi3NS(LALpars);
209 pfa->v[3] += (188.L*SL/3.L + 25.L*dSigmaL) * XLALSimInspiralWaveformParamsLookupNonGRDChi3S(LALpars);
210 }
211
212 if (XLALDictContains(LALpars,"dchi4") || XLALDictContains(LALpars,"dchi4NS") || XLALDictContains(LALpars,"dchi4S") || XLALDictContains(LALpars,"dchikappaS") || XLALDictContains(LALpars,"dchikappaA"))
213 {
214 pfa->v[4] = 5.L*(3058.673L/7.056L + 5429.L/7.L * eta + 617.L * eta*eta)/72.L;
215 pfa->v[4] += -10.L * pn_sigma;
217
218 /* Splitting the deviation parameter into spin and no-spin parts and then adding them here */
219 pfa->v[4] += (5.L*(3058.673L/7.056L + 5429.L/7.L * eta + 617.L * eta*eta)/72.L) * XLALSimInspiralWaveformParamsLookupNonGRDChi4NS(LALpars);
220 pfa->v[4] += -10.L*pn_sigma * XLALSimInspiralWaveformParamsLookupNonGRDChi4S(LALpars);
221
222 /* Adding 2PN SIQM test here..*/
223 pfa->v[4] += siqm_2pn_kappaS * XLALSimInspiralWaveformParamsLookupNonGRDChikappaS(LALpars);
224 pfa->v[4] += siqm_2pn_kappaA * XLALSimInspiralWaveformParamsLookupNonGRDChikappaA(LALpars);
225 }
226
227 if (XLALDictContains(LALpars,"dchi5") || XLALDictContains(LALpars,"dchi5NS") || XLALDictContains(LALpars,"dchi5S"))
228 {
229 pfa->v[5] = 5.L/9.L * (7729.L/84.L - 13.L * eta) * LAL_PI;
230 pfa->v[5] += -1.L * pn_gamma;
232
233 /* Splitting the deviation parameter into spin and no-spin parts and then adding them here */
234 pfa->v[5] += (5.L/9.L * (7729.L/84.L - 13.L * eta) * LAL_PI) * XLALSimInspiralWaveformParamsLookupNonGRDChi5NS(LALpars);
235 pfa->v[5] += -1.L*pn_gamma * XLALSimInspiralWaveformParamsLookupNonGRDChi5S(LALpars);
236 }
237
238 if (XLALDictContains(LALpars,"dchi5l") || XLALDictContains(LALpars,"dchi5lNS") || XLALDictContains(LALpars,"dchi5lS"))
239 {
240 pfa->vlogv[5] = 5.L/3.L * (7729.L/84.L - 13.L * eta) * LAL_PI;
241 pfa->vlogv[5] += -3.L * pn_gamma;
243
244 /* Splitting the deviation parameter into spin and no-spin parts and then adding them here */
245 pfa->vlogv[5] += (5.L/3.L * (7729.L/84.L - 13.L * eta) * LAL_PI) * XLALSimInspiralWaveformParamsLookupNonGRDChi5LNS(LALpars);
246 pfa->vlogv[5] += -3.L*pn_gamma * XLALSimInspiralWaveformParamsLookupNonGRDChi5LS(LALpars);
247 }
248
249 if (XLALDictContains(LALpars,"dchi6") || XLALDictContains(LALpars,"dchi6NS") || XLALDictContains(LALpars,"dchi6S") || XLALDictContains(LALpars,"dchikappaS") || XLALDictContains(LALpars,"dchikappaA"))
250 {
251 pfa->v[6] = (11583.231236531L/4.694215680L - 640.L/3.L * LAL_PI * LAL_PI - 6848.L/21.L*LAL_GAMMA) + eta * (-15737.765635L/3.048192L + 2255./12. * LAL_PI * LAL_PI) + eta*eta * 76055.L/1728.L - eta*eta*eta * 127825.L/1296.L;
252 pfa->v[6] += (-6848.L/21.L)*log(4.);
253 pfa->v[6] += LAL_PI * (3760.L*SL + 1490.L*dSigmaL)/3.L + pn_ss3;
255
256 /* Splitting the deviation parameter into spin and no-spin parts and then adding them here */
257 pfa->v[6] += ((11583.231236531L/4.694215680L - 640.L/3.L * LAL_PI * LAL_PI - 6848.L/21.L*LAL_GAMMA) + eta * (-15737.765635L/3.048192L + 2255./12. * LAL_PI * LAL_PI) + eta*eta * 76055.L/1728.L - eta*eta*eta * 127825.L/1296.L + (-6848.L/21.L)*log(4.)) * XLALSimInspiralWaveformParamsLookupNonGRDChi6NS(LALpars);
258 pfa->v[6] += (LAL_PI * (3760.L*SL + 1490.L*dSigmaL)/3.L + pn_ss3) * XLALSimInspiralWaveformParamsLookupNonGRDChi6S(LALpars);
259
260 /* Adding 3PN SIQM test here..*/
261 pfa->v[6] += siqm_3pn_kappaS * XLALSimInspiralWaveformParamsLookupNonGRDChikappaS(LALpars);
262 pfa->v[6] += siqm_3pn_kappaA * XLALSimInspiralWaveformParamsLookupNonGRDChikappaA(LALpars);
263 }
264
265 if (XLALDictContains(LALpars,"dchi6l"))
266 {
267 pfa->vlogv[6] = -6848.L/21.L;
269 }
270
271 if (XLALDictContains(LALpars,"dchi7")|| XLALDictContains(LALpars,"dchi7NS") || XLALDictContains(LALpars,"dchi7S"))
272 {
273 pfa->v[7] = LAL_PI * ( 77096675.L/254016.L + 378515.L/1512.L * eta - 74045.L/756.L * eta*eta);
274 pfa->v[7] += (-8980424995.L/762048.L + 6586595.L*eta/756.L - 305.L*eta*eta/36.L)*SL - (170978035.L/48384.L - 2876425.L*eta/672.L - 4735.L*eta*eta/144.L) * dSigmaL;
276
277 /* Splitting the deviation parameter into spin and no-spin parts and then adding them here */
278 pfa->v[7] += (LAL_PI * ( 77096675.L/254016.L + 378515.L/1512.L * eta - 74045.L/756.L * eta*eta)) * XLALSimInspiralWaveformParamsLookupNonGRDChi7NS(LALpars);
279 pfa->v[7] += ((-8980424995.L/762048.L + 6586595.L*eta/756.L - 305.L*eta*eta/36.L)*SL - (170978035.L/48384.L - 2876425.L*eta/672.L - 4735.L*eta*eta/144.L) * dSigmaL) * XLALSimInspiralWaveformParamsLookupNonGRDChi7S(LALpars);
280 }
281
282 /* At the very end, multiply everything in the series by the newtonian term */
283 for(int ii = 0; ii <= PN_PHASING_SERIES_MAX_ORDER; ii++)
284 {
285 pfa->v[ii] *= pfaN;
286 pfa->vlogv[ii] *= pfaN;
287 pfa->vlogvsq[ii] *= pfaN;
288 pfa->vneg[ii] *= pfaN;
289 }
290}
291
292/* Compute phase correction that tapers to zero at freqs->data[iEnd] and add to
293 * the phase of of the waveform htilde.
294 */
295
296int XLALSimInspiralPhaseCorrectionsPhasing(COMPLEX16FrequencySeries *htilde, /**< input htilde, will be modified in place */
297 const REAL8Sequence *freqs, /**< frequency series */
298 const UINT4 m, /**< order m of GW mode */
299 const UINT4 iStart, /**< index of starting frequency */
300 const UINT4 iRef, /**< index of reference frequency */
301 const UINT4 iPeak, /**< index of peak frequency */
302 PNPhasingSeries pfa, /**< PN phasing coefficients */
303 const REAL8 mtot, /**< total mass */
304 const REAL8 eta, /**< dimensionless symmetric mass ratio */
305 const REAL8 f_window_div_f_Peak, /**< Frequency at which to attach non-GR and GR waveforms, inputted as a fraction of f_Peak (should be between 0 and 1) */
306 const REAL8 iStartFinal, /**< index of starting frequency for TGR modification */
307 const REAL8 NCyclesStep /**< Choose number of GW cycles over which to taper the non-GR phase correction */
308 )
309{
310
311 UINT4 i;
312 const REAL8 pfaN0 = 3.L/(128.L * eta);
313 const REAL8 piM = LAL_PI * mtot;
314 /*const REAL8 vWindow = cbrt(2*piM * freqs->data[iEnd]/m); //Center of the tapering function in v-space*/
315 const REAL8 vWindow = cbrt(2*piM * f_window_div_f_Peak * freqs->data[iPeak]/2); //Center of the tapering function in v-space
316 const REAL8 width = (NCyclesStep * LAL_PI * vWindow * vWindow * vWindow * vWindow * vWindow * vWindow)/(50. * pfaN0); //Width of tapering step function
317
318
319 /* Compute the second derivative of the phase correction and taper with step function
320 * Phase correction only generated from f_low to f_max
321 */
322 REAL8Sequence *d2phasenonGRdf2Tapered = NULL;
323 d2phasenonGRdf2Tapered = XLALCreateREAL8Sequence( (UINT4) freqs->length );
324 REAL8 f;
325 for ( i = iStart; i < freqs->length; i++)
326 {
327 f = freqs->data[i];
328 if (f>0) d2phasenonGRdf2Tapered->data[i] = PNPhaseSecondDerivative(f, m, pfa, mtot)/ (1. + exp( (cbrt(2*piM * freqs->data[i]/m) - vWindow) / width ) );
329 }
330
331 /* Interpolate and integrate (twice) the tapered second derivative of the phase correction to compute the phase correction
332 * Set the value of the correction and its first derivative to zero at f_ref
333 */
334 REAL8Sequence *dphasenonGRdfTapered = NULL, *phasenonGRTapered = NULL;
335 dphasenonGRdfTapered = XLALCreateREAL8Sequence( (UINT4) freqs->length );
336 phasenonGRTapered = XLALCreateREAL8Sequence( (UINT4) freqs->length );
337
338 gsl_spline *splineTemp = NULL;
339 gsl_interp_accel *acc = NULL;
340 splineTemp = gsl_spline_alloc (gsl_interp_linear, freqs->length);
341 gsl_spline_init(splineTemp, freqs->data, d2phasenonGRdf2Tapered->data, freqs->length);
342
343 dphasenonGRdfTapered->data[iRef] = 0.0;
344 for ( i = iRef; i-- > iStart; )
345 dphasenonGRdfTapered->data[i] = dphasenonGRdfTapered->data[i+1] - gsl_spline_eval_integ(splineTemp, freqs->data[i], freqs->data[i+1], acc);
346 for ( i = iRef+1; i < freqs->length; i++ )
347 dphasenonGRdfTapered->data[i] = dphasenonGRdfTapered->data[i-1] + gsl_spline_eval_integ(splineTemp, freqs->data[i-1], freqs->data[i], acc);
348
349 gsl_spline_init(splineTemp, freqs->data, dphasenonGRdfTapered->data, freqs->length);
350 phasenonGRTapered->data[iRef] = 0.0;
351 for ( i = iRef; i-- > iStart; )
352 phasenonGRTapered->data[i] = phasenonGRTapered->data[i+1] - gsl_spline_eval_integ(splineTemp, freqs->data[i], freqs->data[i+1], acc);
353 for ( i = iRef+1; i < freqs->length; i++ )
354 phasenonGRTapered->data[i] = phasenonGRTapered->data[i-1] + gsl_spline_eval_integ(splineTemp, freqs->data[i-1], freqs->data[i], acc);
355
356
357 /* Compute the first derivative of the phase correction at fPeak and subtract from
358 * the phase correction to ensure the derivative vanishes at fPeak to accomplish alignment
359 * of time-domain waveform.
360 *
361 * Then add phase correction to input waveform phase.
362 */
363 REAL8 PNPhaseRefDerivative = dphasenonGRdfTapered->data[iPeak];
364 /*printf("phase at fref=%f\n", phasenonGRTapered->data[iRef]);*/
365
366 for ( i = iStartFinal; i < freqs->length; i++ ) {
367 REAL8 phasing = phasenonGRTapered->data[i] - PNPhaseRefDerivative * (freqs->data[i]-freqs->data[iRef]) ;
368 htilde->data->data[i] *= cexp(phasing * 1.j);
369 }
370
371 gsl_spline_free(splineTemp);
372 gsl_interp_accel_free(acc);
373 XLALDestroyREAL8Sequence(d2phasenonGRdf2Tapered);
374 XLALDestroyREAL8Sequence(dphasenonGRdfTapered);
375 XLALDestroyREAL8Sequence(phasenonGRTapered);
376 return 0;
377}
378
379REAL8 PNPhase(REAL8 f, /* frequency in Hz */
380 UINT4 m,
381 PNPhasingSeries pfa,
382 const REAL8 mtot) /* total mass in seconds */
383{
384
385 const REAL8 piM = LAL_PI * mtot;
386 const REAL8 pfa7 = pfa.v[7];
387 const REAL8 pfa6 = pfa.v[6];
388 const REAL8 pfl6 = pfa.vlogv[6];
389 const REAL8 pfa5 = pfa.v[5];
390 const REAL8 pfl5 = pfa.vlogv[5];
391 const REAL8 pfa4 = pfa.v[4];
392 const REAL8 pfa3 = pfa.v[3];
393 const REAL8 pfa2 = pfa.v[2];
394 const REAL8 pfa1 = pfa.v[1];
395 const REAL8 pfaN = pfa.v[0];
396 const REAL8 pfaMinus1 = pfa.vneg[1];
397 const REAL8 pfaMinus2 = pfa.vneg[2];
398
399 const REAL8 v = cbrt(2*piM*f/m);
400 const REAL8 logv = log(v);
401 const REAL8 v2 = v * v;
402 const REAL8 v3 = v * v2;
403 const REAL8 v4 = v * v3;
404 const REAL8 v5 = v * v4;
405 const REAL8 v6 = v * v5;
406 const REAL8 v7 = v * v6;
407 REAL8 phasing = 0.0;
408
409 phasing += pfa7 * v7;
410 phasing += (pfa6 + pfl6 * logv) * v6;
411 phasing += (pfa5 + pfl5 * logv) * v5;
412 phasing += pfa4 * v4;
413 phasing += pfa3 * v3;
414 phasing += pfa2 * v2;
415 phasing += pfa1 * v;
416 phasing += pfaN;
417 phasing += pfaMinus1 / v;
418 phasing += pfaMinus2 /v2;
419 phasing /= v5;
420 phasing *=m/2.;
421
422 return -phasing;
423}
424
425/* Derivative of phase with respect to f:
426 * d[ v^(n - 5) ] = (pi * M) / 3 * (n - 5) v^(n - 8)
427 * d[ v^(n - 5) * log(v)] = (pi * M) / 3 * (1 + (n - 5) * log(v)) * v^(n - 8)
428 *
429 * where v = (pi * M * f)^(1 / 3)
430 */
431
432REAL8 PNPhaseDerivative(REAL8 f, /* frequency in Hz */
433 UINT4 m,
434 PNPhasingSeries pfa,
435 const REAL8 mtot) /* total mass in seconds */
436{
437
438 const REAL8 piM = LAL_PI * mtot;
439 const REAL8 pfa7 = pfa.v[7];
440 const REAL8 pfa6 = pfa.v[6];
441 const REAL8 pfl6 = pfa.vlogv[6];
442 //const REAL8 pfa5 = pfa.v[5];
443 const REAL8 pfl5 = pfa.vlogv[5];
444 const REAL8 pfa4 = pfa.v[4];
445 const REAL8 pfa3 = pfa.v[3];
446 const REAL8 pfa2 = pfa.v[2];
447 const REAL8 pfa1 = pfa.v[1];
448 const REAL8 pfaN = pfa.v[0];
449 const REAL8 pfaMinus1 = pfa.vneg[1];
450 const REAL8 pfaMinus2 = pfa.vneg[2];
451
452 const REAL8 v = cbrt(2*piM*f/m);
453 const REAL8 logv = log(v);
454 const REAL8 v2 = v * v;
455 const REAL8 v3 = v * v2;
456 const REAL8 v4 = v * v3;
457 const REAL8 v5 = v * v4;
458// const REAL8 v6 = v * v5;
459// const REAL8 v7 = v * v6;
460 REAL8 phasing = 0.0;
461
462 phasing += 2. * pfa7 * v4;
463 phasing += (1. * pfa6 + 1. * pfl6 + 1. * pfl6 * logv) * v3;
464 phasing += 1. * pfl5 * v2;
465 phasing += -1. * pfa4 * v;
466 phasing += -2. * pfa3;
467 phasing += -3. * pfa2 / v; /* there was a mistake here in the factor. it should be -3 instead of -1.*/
468 phasing += -4. * pfa1 / v2;
469 phasing += -5. * pfaN / v3;
470 phasing += -6. * pfaMinus1 / v4;
471 phasing += -7. * pfaMinus2 /v5;
472 phasing /= v5;
473 phasing *= piM/3.;
474
475 return -phasing;
476}
477
478/* Derivative of phase with respect to f
479 * d2[v^(n - 5)] = (pi * M)^2 / 9 * (n - 5) * (n - 8) * v^(n - 11)
480 * d2[v^(n - 5) * log(v)] = (pi * M)^2 / 9 * ( (2 * n - 13) + (n - 5) * (n - 8) * log(v) ) * v^(n - 11)
481 *
482 * where v = (pi * M * f)^(1 / 3)
483 */
484
485REAL8 PNPhaseSecondDerivative(REAL8 f, /* frequency in Hz */
486 UINT4 m,
487 PNPhasingSeries pfa,
488 const REAL8 mtot) /* total mass in seconds */
489{
490
491 const REAL8 piM = LAL_PI * mtot;
492
493 const REAL8 pfa7 = pfa.v[7];
494 const REAL8 pfa6 = pfa.v[6];
495 const REAL8 pfl6 = pfa.vlogv[6];
496 //const REAL8 pfa5 = pfa.v[5];
497 const REAL8 pfl5 = pfa.vlogv[5];
498 const REAL8 pfa4 = pfa.v[4];
499 const REAL8 pfa3 = pfa.v[3];
500 const REAL8 pfa2 = pfa.v[2];
501 const REAL8 pfa1 = pfa.v[1];
502 const REAL8 pfaN = pfa.v[0];
503 const REAL8 pfaMinus1 = pfa.vneg[1];
504 const REAL8 pfaMinus2 = pfa.vneg[2];
505
506 const REAL8 v = cbrt(2*piM*f/m);
507 const REAL8 logv = log(v);
508 const REAL8 v2 = v * v;
509 const REAL8 v3 = v * v2;
510 const REAL8 v4 = v * v3;
511 const REAL8 v5 = v * v4;
512 const REAL8 v6 = v * v5;
513 const REAL8 v7 = v * v6;
514 const REAL8 v8 = v * v7;
515 REAL8 phasing = 0.0;
516
517 phasing += -2. * pfa7 * v;
518 phasing += (-2. * pfa6 - 1. * pfl6 - 2. * pfl6 * logv);
519 phasing += -3. * pfl5 / v;
520 phasing += 4. * pfa4 / v2 ;
521 phasing += 10. * pfa3 / v3;
522 phasing += 18. * pfa2 / v4;
523 phasing += 28. * pfa1 / v5;
524 phasing += 40. * pfaN / v6;
525 phasing += 54. * pfaMinus1 / v7;
526 phasing += 70. * pfaMinus2 /v8;
527 phasing /= v5;
528 phasing *= piM*piM/9.;
529 phasing *=2./m;
530 return -phasing;
531}
int XLALDictContains(const LALDict *dict, const char *key)
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 vector d(const double L_norm, const double J_norm, const vector roots)
Internal function that returns the coefficients "d_0", "d_2" and "d_4" from 1703.03967 corresponding ...
static REAL8 Get22PeakAngFreq(REAL8 UNUSED eta, REAL8 a)
Combined spin chi defined in Eq.
int XLALSimInspiralPhaseCorrectionsPhasing(COMPLEX16FrequencySeries *htilde, const REAL8Sequence *freqs, const UINT4 m, const UINT4 iStart, const UINT4 iRef, const UINT4 iPeak, PNPhasingSeries pfa, const REAL8 mtot, const REAL8 eta, const REAL8 f_window_div_f_Peak, const REAL8 iStartFinal, const REAL8 NCyclesStep)
REAL8 PNPhase(REAL8 f, UINT4 m, PNPhasingSeries pfa, const REAL8 mtot)
REAL8 PNPhaseDerivative(REAL8 f, UINT4 m, PNPhasingSeries pfa, const REAL8 mtot)
REAL8 PNPhaseSecondDerivative(REAL8 f, UINT4 m, PNPhasingSeries pfa, const REAL8 mtot)
int XLALSimInspiralTestingGRCorrections(COMPLEX16FrequencySeries *htilde, const UINT4 l, const UINT4 m, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi1z, const REAL8 chi2z, const REAL8 f_low, const REAL8 f_ref, const REAL8 f_window_div_f_Peak, const REAL8 NCyclesStep, LALDict *LALpars)
void XLALSimInspiralPNCorrections(PNPhasingSeries *pfa, const REAL8 m1, const REAL8 m2, const REAL8 chi1L, const REAL8 chi2L, const REAL8 chi1sq, const REAL8 chi2sq, const REAL8 chi1dotchi2, const REAL8 qm_def1, const REAL8 qm_def2, LALDict *LALpars)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi4S(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi5(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChiMinus2(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi6(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi3NS(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi6NS(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi5S(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi6L(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi3S(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi0(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupTidalLambda2(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChiMinus1(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi5NS(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi5LNS(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi6S(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi1(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupTidalLambda1(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi4(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi3(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi7S(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi7(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi5L(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChikappaA(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi5LS(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi7NS(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChikappaS(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi2(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupNonGRDChi4NS(LALDict *params)
int l
Definition: bh_qnmode.c:135
double i
Definition: bh_ringdown.c:118
#define LAL_MSUN_SI
#define LAL_PI
#define LAL_MTSUN_SI
#define LAL_GAMMA
double REAL8
uint32_t UINT4
int32_t INT4
#define PN_PHASING_SERIES_MAX_ORDER
Structure for passing around PN phasing coefficients.
static const INT4 m
static const INT4 q
static const INT4 a
void XLALDestroyREAL8Sequence(REAL8Sequence *sequence)
REAL8Sequence * XLALCreateREAL8Sequence(size_t length)
#define XLAL_ERROR(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_EINVAL
COMPLEX16Sequence * data
COMPLEX16 * data
REAL8 vneg[PN_PHASING_SERIES_MAX_ORDER+1]
REAL8 vlogv[PN_PHASING_SERIES_MAX_ORDER+1]
REAL8 vlogvsq[PN_PHASING_SERIES_MAX_ORDER+1]
REAL8 v[PN_PHASING_SERIES_MAX_ORDER+1]
REAL8 * data