Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-5e288d3
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimIMRPhenom.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2011 P. Ajith, Nickolas Fotopoulos
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 <complex.h>
22
23#include <lal/LALStdlib.h>
24#include <lal/LALSimIMR.h>
25#include <lal/LALConstants.h>
26#include <lal/Date.h>
27#include <lal/FrequencySeries.h>
28#include <lal/StringInput.h>
29#include <lal/TimeSeries.h>
30#include <lal/TimeFreqFFT.h>
31#include <lal/Units.h>
32#include <lal/XLALError.h>
33#include <gsl/gsl_blas.h>
34#include <gsl/gsl_linalg.h>
35#include <gsl/gsl_permutation.h>
36#include <gsl/gsl_complex.h>
37
38typedef struct tagBBHPhenomParams{
39 REAL8 fMerger;
40 REAL8 fRing;
41 REAL8 fCut;
42 REAL8 sigma;
43 REAL8 psi0;
44 REAL8 psi1;
45 REAL8 psi2;
46 REAL8 psi3;
47 REAL8 psi4;
48 REAL8 psi5;
49 REAL8 psi6;
50 REAL8 psi7;
51 REAL8 psi8;
52}
53BBHPhenomParams;
54
55/*
56 * private function prototypes; all internal functions use solar masses.
57 *
58 */
59
60static BBHPhenomParams *ComputeIMRPhenomAParams(const REAL8 m1, const REAL8 m2);
61static BBHPhenomParams *ComputeIMRPhenomBParams(const REAL8 m1, const REAL8 m2, const REAL8 chi);
62
63static REAL8 EstimateSafeFMinForTD(const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 deltaT);
64static REAL8 EstimateSafeFMaxForTD(const REAL8 f_max, const REAL8 dt);
65static REAL8 ComputeTau0(const REAL8 m1, const REAL8 m2, const REAL8 f_min);
66static size_t EstimateIMRLength(const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 deltaT);
67static size_t NextPow2(const size_t n);
68
69static REAL8 LorentzianFn(const REAL8 freq, const REAL8 fRing, const REAL8 sigma);
70
71static int IMRPhenomAGenerateFD(COMPLEX16FrequencySeries **htilde, const REAL8 phi0, const REAL8 deltaF, const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, const BBHPhenomParams *params);
72static int IMRPhenomBGenerateFD(COMPLEX16FrequencySeries **htilde, const REAL8 phi0, const REAL8 deltaF, const REAL8 m1, const REAL8 m2, const REAL8 chi, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, const BBHPhenomParams *params);
73static int IMRPhenomAGenerateTD(REAL8TimeSeries **h, const REAL8 phiPeak, const REAL8 deltaT, const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, const BBHPhenomParams *params);
74static int IMRPhenomBGenerateTD(REAL8TimeSeries **h, const REAL8 phiPeak, const REAL8 deltaT, const REAL8 m1, const REAL8 m2, const REAL8 chi, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, const BBHPhenomParams *params);
75static int FDToTD(REAL8TimeSeries **signalTD, const COMPLEX16FrequencySeries *signalFD, const REAL8 totalMass, const REAL8 deltaT, const REAL8 f_min, const REAL8 f_max, const REAL8 f_min_wide, const REAL8 f_max_wide);
76static size_t find_instant_freq(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 target, const size_t start);
77static size_t find_peak_amp(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc);
78static int apply_phase_shift(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 shift);
79static int apply_inclination(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 inclination);
80
81/* INTERNAL ROUTINES */
82
83/* *******************************************************************/
84/* Compute phenomenological parameters for non-spinning binaries */
85/* Ref. Eq.(4.18) of http://arxiv.org/pdf/0710.2335 and */
86/* Table I of http://arxiv.org/pdf/0712.0343 */
87/* */
88/* Takes solar masses. */
89/* *******************************************************************/
90static BBHPhenomParams *ComputeIMRPhenomAParams(const REAL8 m1, const REAL8 m2) {
91 /* calculate the total mass and symmetric mass ratio */
92 const REAL8 totalMass = m1 + m2;
93 const REAL8 eta = m1 * m2 / (totalMass * totalMass);
94 const REAL8 piM = totalMass * LAL_PI * LAL_MTSUN_SI;
95 const REAL8 etap2 = eta*eta;
96
97 const REAL8 fMerg_a = 6.6389e-01;
98 const REAL8 fMerg_b = -1.0321e-01;
99 const REAL8 fMerg_c = 1.0979e-01;
100
101 const REAL8 fRing_a = 1.3278e+00;
102 const REAL8 fRing_b = -2.0642e-01;
103 const REAL8 fRing_c = 2.1957e-01;
104
105 const REAL8 sigma_a = 1.1383e+00;
106 const REAL8 sigma_b = -1.7700e-01;
107 const REAL8 sigma_c = 4.6834e-02;
108
109 const REAL8 fCut_a = 1.7086e+00;
110 const REAL8 fCut_b = -2.6592e-01;
111 const REAL8 fCut_c = 2.8236e-01;
112
113 const REAL8 psi0_a = -1.5829e-01;
114 const REAL8 psi0_b = 8.7016e-02;
115 const REAL8 psi0_c = -3.3382e-02;
116
117 const REAL8 psi2_a = 3.2967e+01;
118 const REAL8 psi2_b = -1.9000e+01;
119 const REAL8 psi2_c = 2.1345e+00;
120
121 const REAL8 psi3_a = -3.0849e+02;
122 const REAL8 psi3_b = 1.8211e+02;
123 const REAL8 psi3_c = -2.1727e+01;
124
125 const REAL8 psi4_a = 1.1525e+03;
126 const REAL8 psi4_b = -7.1477e+02;
127 const REAL8 psi4_c = 9.9692e+01;
128
129 const REAL8 psi6_a = 1.2057e+03;
130 const REAL8 psi6_b = -8.4233e+02;
131 const REAL8 psi6_c = 1.8046e+02;
132
133 const REAL8 psi7_a = 0.;
134 const REAL8 psi7_b = 0.;
135 const REAL8 psi7_c = 0.;
136
137 BBHPhenomParams *phenParams = (BBHPhenomParams *) XLALMalloc(sizeof(BBHPhenomParams));
138 if (!phenParams) XLAL_ERROR_NULL(XLAL_EFUNC);
139 memset(phenParams, 0, sizeof(BBHPhenomParams));
140
141 /* Evaluate the polynomials. See Eq. (4.18) of P. Ajith et al
142 * arXiv:0710.2335 [gr-qc] */
143 phenParams->fCut = (fCut_a*etap2 + fCut_b*eta + fCut_c)/piM;
144 phenParams->fMerger = (fMerg_a*etap2 + fMerg_b*eta + fMerg_c)/piM;
145 phenParams->fRing = (fRing_a*etap2 + fRing_b*eta + fRing_c)/piM;
146 phenParams->sigma = (sigma_a*etap2 + sigma_b*eta + sigma_c)/piM;
147
148 phenParams->psi0 = (psi0_a*etap2 + psi0_b*eta + psi0_c)/(eta*pow(piM, 5./3.));
149 phenParams->psi1 = 0.;
150 phenParams->psi2 = (psi2_a*etap2 + psi2_b*eta + psi2_c)/(eta*pow(piM, 3./3.));
151 phenParams->psi3 = (psi3_a*etap2 + psi3_b*eta + psi3_c)/(eta*pow(piM, 2./3.));
152 phenParams->psi4 = (psi4_a*etap2 + psi4_b*eta + psi4_c)/(eta*cbrt(piM));
153 phenParams->psi5 = 0.;
154 phenParams->psi6 = (psi6_a*etap2 + psi6_b*eta + psi6_c)/(eta/cbrt(piM));
155 phenParams->psi7 = (psi7_a*etap2 + psi7_b*eta + psi7_c)/(eta*pow(piM, -2./3.));
156
157 return phenParams;
158}
159
160/*********************************************************************/
161/* Compute phenomenological parameters for aligned-spin binaries */
162/* Ref. Eq.(2) and Table I of http://arxiv.org/pdf/0909.2867 */
163/* */
164/* Takes solar masses. Populates and returns a new BBHPhenomParams */
165/* structure. */
166/*********************************************************************/
167static BBHPhenomParams *ComputeIMRPhenomBParams(const REAL8 m1, const REAL8 m2, const REAL8 chi) {
168 /* calculate the total mass and symmetric mass ratio */
169 const REAL8 totalMass = m1 + m2;
170 const REAL8 eta = m1 * m2 / (totalMass * totalMass);
171 const REAL8 piM = totalMass * LAL_PI * LAL_MTSUN_SI;
172
173 /* spinning phenomenological waveforms */
174 const REAL8 etap2 = eta*eta;
175 const REAL8 chip2 = chi*chi;
176 const REAL8 etap3 = etap2*eta;
177 const REAL8 etap2chi = etap2*chi;
178 const REAL8 etachip2 = eta*chip2;
179 const REAL8 etachi = eta*chi;
180
181 BBHPhenomParams *phenParams = (BBHPhenomParams *) XLALMalloc(sizeof(BBHPhenomParams));
182 if (!phenParams) XLAL_ERROR_NULL(XLAL_EFUNC);
183 memset(phenParams, 0, sizeof(BBHPhenomParams));
184
185 phenParams->psi0 = 3./(128.*eta);
186
187 phenParams->psi2 = 3715./756. +
188 -9.2091e+02*eta + 4.9213e+02*etachi + 1.3503e+02*etachip2 +
189 6.7419e+03*etap2 + -1.0534e+03*etap2chi +
190 -1.3397e+04*etap3 ;
191
192 phenParams->psi3 = -16.*LAL_PI + 113.*chi/3. +
193 1.7022e+04*eta + -9.5659e+03*etachi + -2.1821e+03*etachip2 +
194 -1.2137e+05*etap2 + 2.0752e+04*etap2chi +
195 2.3859e+05*etap3 ;
196
197 phenParams->psi4 = 15293365./508032. - 405.*chip2/8. +
198 -1.2544e+05*eta + 7.5066e+04*etachi + 1.3382e+04*etachip2 +
199 8.7354e+05*etap2 + -1.6573e+05*etap2chi +
200 -1.6936e+06*etap3 ;
201
202 phenParams->psi6 = -8.8977e+05*eta + 6.3102e+05*etachi + 5.0676e+04*etachip2 +
203 5.9808e+06*etap2 + -1.4148e+06*etap2chi +
204 -1.1280e+07*etap3 ;
205
206 phenParams->psi7 = 8.6960e+05*eta + -6.7098e+05*etachi + -3.0082e+04*etachip2 +
207 -5.8379e+06*etap2 + 1.5145e+06*etap2chi +
208 1.0891e+07*etap3 ;
209
210 phenParams->psi8 = -3.6600e+05*eta + 3.0670e+05*etachi + 6.3176e+02*etachip2 +
211 2.4265e+06*etap2 + -7.2180e+05*etap2chi +
212 -4.5524e+06*etap3;
213
214 phenParams->fMerger = 1. - 4.4547*pow(1.-chi,0.217) + 3.521*pow(1.-chi,0.26) +
215 6.4365e-01*eta + 8.2696e-01*etachi + -2.7063e-01*etachip2 +
216 -5.8218e-02*etap2 + -3.9346e+00*etap2chi +
217 -7.0916e+00*etap3 ;
218
219 phenParams->fRing = (1. - 0.63*pow(1.-chi,0.3))/2. +
220 1.4690e-01*eta + -1.2281e-01*etachi + -2.6091e-02*etachip2 +
221 -2.4900e-02*etap2 + 1.7013e-01*etap2chi +
222 2.3252e+00*etap3 ;
223
224 phenParams->sigma = (1. - 0.63*pow(1.-chi,0.3))*pow(1.-chi,0.45)/4. +
225 -4.0979e-01*eta + -3.5226e-02*etachi + 1.0082e-01*etachip2 +
226 1.8286e+00*etap2 + -2.0169e-02*etap2chi +
227 -2.8698e+00*etap3 ;
228
229 phenParams->fCut = 3.2361e-01 + 4.8935e-02*chi + 1.3463e-02*chip2 +
230 -1.3313e-01*eta + -8.1719e-02*etachi + 1.4512e-01*etachip2 +
231 -2.7140e-01*etap2 + 1.2788e-01*etap2chi +
232 4.9220e+00*etap3 ;
233
234 phenParams->fCut /= piM;
235 phenParams->fMerger/= piM;
236 phenParams->fRing /= piM;
237 phenParams->sigma /= piM;
238
239 phenParams->psi1 = 0.;
240 phenParams->psi5 = 0.;
241
242 return phenParams;
243}
244
245/*
246 * Return tau0, the Newtonian chirp length estimate. Ref. Eq.(6) of B. S. Sathyaprakash, http://arxiv.org/abs/gr-qc/9411043v1
247 */
248static REAL8 ComputeTau0(const REAL8 m1, const REAL8 m2, const REAL8 f_min) {
249 const REAL8 totalMass = m1 + m2;
250 const REAL8 eta = m1 * m2 / (totalMass * totalMass);
251 return 5. * totalMass * LAL_MTSUN_SI / (256. * eta * pow(LAL_PI * totalMass * LAL_MTSUN_SI * f_min, 8./3.));
252}
253
254/*
255 * Estimate the length of a TD vector that can hold the waveform as the Newtonian
256 * chirp time tau0 plus 1000 M.
257 */
258static size_t EstimateIMRLength(const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 deltaT) {
259 return (size_t) floor((ComputeTau0(m1, m2, f_min) + 1000 * (m1 + m2) * LAL_MTSUN_SI) / deltaT);
260}
261
262static size_t NextPow2(const size_t n) {
263 return 1 << (size_t) ceil(log2(n));
264}
265
266/*
267 * Find a lower value for f_min (using the definition of Newtonian chirp
268 * time) such that the waveform has a minimum length of tau0. This is
269 * necessary to avoid FFT artifacts.
270 */
271static REAL8 EstimateSafeFMinForTD(const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 deltaT) {
272 const REAL8 totalMass = m1 + m2;
273 const REAL8 eta = m1 * m2 / (totalMass * totalMass);
274 REAL8 tau0 = deltaT * NextPow2(1.025 * EstimateIMRLength(m1, m2, f_min, deltaT));
275 REAL8 temp_f_min = pow((tau0 * 256. * eta * pow(totalMass * LAL_MTSUN_SI, 5./3.) / 5.), -3./8.) / LAL_PI;
276 if (temp_f_min > f_min) temp_f_min = f_min;
277 if (temp_f_min < 0.5) temp_f_min = 0.5;
278 return temp_f_min;
279}
280
281/*
282 * Find a higher value of f_max so that we can safely apply a window later.
283 * The safety factor 1.025 is an empirical estimation
284 */
285static REAL8 EstimateSafeFMaxForTD(const REAL8 f_max, const REAL8 deltaT) {
286 REAL8 temp_f_max = 1.025 * f_max;
287
288 /* make sure that these frequencies are not too out of range */
289 if (temp_f_max > 2. / deltaT - 100.) temp_f_max = 2. / deltaT - 100.;
290 return temp_f_max;
291}
292
293static REAL8 LorentzianFn(const REAL8 freq, const REAL8 fRing, const REAL8 sigma) {
294 return sigma / (LAL_TWOPI * ((freq - fRing)*(freq - fRing)
295 + sigma*sigma / 4.0));
296}
297
298
299/*
300 * Private function to generate IMRPhenomA frequency-domain waveforms given coefficients
301 */
302static int IMRPhenomAGenerateFD(
303 COMPLEX16FrequencySeries **htilde, /**< FD waveform */
304 const REAL8 phi0, /**< orbital phase at peak (rad) */
305 const REAL8 deltaF, /**< frequency resolution */
306 const REAL8 m1, /**< mass of companion 1 [solar masses] */
307 const REAL8 m2, /**< mass of companion 2 [solar masses] */
308 const REAL8 f_min, /**< start frequency */
309 const REAL8 f_max, /**< end frequency */
310 const REAL8 distance, /**< distance of source */
311 const BBHPhenomParams *params /**< from ComputeIMRPhenomAParams */
312) {
313 const LIGOTimeGPS ligotimegps_zero = {0, 0};
314 COMPLEX16 *data;
315 size_t i;
316
317 const REAL8 fMerg = params->fMerger;
318 const REAL8 fRing = params->fRing;
319 const REAL8 sigma = params->sigma;
320 const REAL8 totalMass = m1 + m2;
321 const REAL8 eta = m1 * m2 / (totalMass * totalMass);
322
323 /* compute the amplitude pre-factor */
324 const REAL8 amp0 = - pow(LAL_MTSUN_SI*totalMass, 5./6.) * pow(fMerg, -7./6.)
325 / pow(LAL_PI, 2./3.) * sqrt(5. * eta / 24.) / (distance / LAL_C_SI);
326
327 /* allocate htilde */
328 size_t n = NextPow2(f_max / deltaF) + 1;
329 *htilde = XLALCreateCOMPLEX16FrequencySeries("htilde: FD waveform", &ligotimegps_zero, 0.0, deltaF, &lalStrainUnit, n);
330 memset((*htilde)->data->data, 0, n * sizeof(COMPLEX16));
331 XLALUnitMultiply(&((*htilde)->sampleUnits), &((*htilde)->sampleUnits), &lalSecondUnit);
332 if (!(*htilde)) XLAL_ERROR(XLAL_EFUNC);
333
334 /* now generate the waveform at all frequency bins except DC and Nyquist */
335 data = (*htilde)->data->data;
336 for (i=1; i < n - 1; i++) {
337 REAL8 ampEff, psiEff;
338 /* Fourier frequency corresponding to this bin */
339 REAL8 f = i * deltaF;
340 REAL8 cbrt_f = cbrt(f);
341 REAL8 fNorm = f / fMerg;
342
343 /* compute the amplitude */
344 if ((f < f_min) || (f > f_max)) continue;
345 else if (f <= fMerg) ampEff = amp0 * pow(fNorm, -7./6.);
346 else if ((f > fMerg) & (f <= fRing)) ampEff = amp0 * pow(fNorm, -2./3.);
347 else if (f > fRing)
348 ampEff = amp0 * LAL_PI_2 * pow(fRing / fMerg, -2./3.) * sigma
349 * LorentzianFn(f, fRing, sigma);
350 else {
351 XLALDestroyCOMPLEX16FrequencySeries(*htilde);
352 *htilde = NULL;
353 XLAL_ERROR(XLAL_EDOM);
354 }
355
356 /* now compute the phase */
357 psiEff = - 2.*phi0
358 + params->psi0 / (f * f) * cbrt_f /* f^-5/3 */
359 + params->psi1 / (f * cbrt_f) /* f^-4/3 */
360 + params->psi2 / f /* f^-3/3 */
361 + params->psi3 / f * cbrt_f /* f^-2/3 */
362 + params->psi4 / cbrt_f /* f^-1/3 */
363 + params->psi5 /* f^0/3 */
364 + params->psi6 * cbrt_f /* f^1/3 */
365 + params->psi7 * f / cbrt_f; /* f^2/3 */
366
367 /* generate the waveform */
368 data[i] = ampEff * cos(psiEff);
369 data[i] += -I * ampEff * sin(psiEff);
370 }
371
372 return XLAL_SUCCESS;
373}
374
375/*
376 * Private function to generate IMRPhenomB frequency-domain waveforms given coefficients
377 */
378static int IMRPhenomBGenerateFD(
379 COMPLEX16FrequencySeries **htilde, /**< FD waveform */
380 const REAL8 phi0, /**< orbital phase at peak (rad) */
381 const REAL8 deltaF, /**< frequency resolution */
382 const REAL8 m1, /**< mass of companion 1 [solar masses] */
383 const REAL8 m2, /**< mass of companion 2 [solar masses] */
384 const REAL8 chi, /**< mass-weighted aligned-spin parameter */
385 const REAL8 f_min, /**< start frequency */
386 const REAL8 f_max, /**< end frequency */
387 const REAL8 distance, /**< distance of source */
388 const BBHPhenomParams *params /**< from ComputeIMRPhenomBParams */
389) {
390 static LIGOTimeGPS ligotimegps_zero = {0, 0};
391 size_t i;
392
393 const REAL8 fMerg = params->fMerger;
394 const REAL8 fRing = params->fRing;
395 const REAL8 sigma = params->sigma;
396 const REAL8 totalMass = m1 + m2;
397 const REAL8 eta = m1 * m2 / (totalMass * totalMass);
398 const REAL8 piM = LAL_PI * totalMass * LAL_MTSUN_SI;
399
400 /* compute the amplitude pre-factor */
401 REAL8 amp0 = - pow(LAL_MTSUN_SI*totalMass, 5./6.) * pow(fMerg, -7./6.)
402 / pow(LAL_PI, 2./3.) * sqrt(5. * eta / 24.) / (distance / LAL_C_SI);
403
404 /***********************************************************************/
405 /* these are the parameters required for the "new" phenomenological IMR
406 * waveforms*/
407 /***********************************************************************/
408
409 /* PN corrections to the frequency domain amplitude of the (2,2) mode */
410 const REAL8 alpha2 = -323./224. + 451.*eta/168.;
411 const REAL8 alpha3 = (27./8. - 11.*eta/6.)*chi;
412
413 /* leading order power law of the merger amplitude */
414 static REAL8 mergPower = -2./3.;
415
416 /* spin-dependent corrections to the merger amplitude */
417 const REAL8 epsilon_1 = 1.4547*chi - 1.8897;
418 const REAL8 epsilon_2 = -1.8153*chi + 1.6557;
419
420 /* normalisation constant of the inspiral amplitude */
421 const REAL8 vMerg = cbrt(piM * fMerg);
422 const REAL8 vRing = cbrt(piM * fRing);
423
424 REAL8 w1 = (1. + alpha2 * vMerg * vMerg + alpha3 * piM * fMerg) / (1. + epsilon_1 * vMerg + epsilon_2 * vMerg * vMerg);
425 REAL8 w2 = w1 * (LAL_PI * sigma / 2.) * pow(fRing / fMerg, mergPower)
426 * (1. + epsilon_1 * vRing + epsilon_2 * vRing * vRing);
427
428 /* allocate htilde */
429 size_t n = NextPow2(f_max / deltaF) + 1;
430 *htilde = XLALCreateCOMPLEX16FrequencySeries("htilde: FD waveform", &ligotimegps_zero, 0.0, deltaF, &lalStrainUnit, n);
431 memset((*htilde)->data->data, 0, n * sizeof(COMPLEX16));
432 XLALUnitMultiply(&((*htilde)->sampleUnits), &((*htilde)->sampleUnits), &lalSecondUnit);
433 if (!(*htilde)) XLAL_ERROR(XLAL_EFUNC);
434
435 /* now generate the waveform */
436 size_t ind_max = (size_t) (f_max / deltaF);
437 for (i = (size_t) (f_min / deltaF); i < ind_max; i++) {
438 REAL8 ampEff, psiEff;
439 REAL8 v, v2, v3, v4, v5, v6, v7, v8;
440
441 /* Fourier frequency corresponding to this bin */
442 REAL8 f = i * deltaF;
443
444 /* PN expansion parameter */
445 v = cbrt(piM * f);
446 v2 = v*v; v3 = v2*v; v4 = v2*v2; v5 = v4*v; v6 = v3*v3; v7 = v6*v, v8 = v7*v;
447
448 /* compute the amplitude */
449 if (f <= fMerg)
450 ampEff = pow(f / fMerg, -7./6.)*(1. + alpha2 * v2 + alpha3 * v3);
451 else if (f > fRing)
452 ampEff = w2 * LorentzianFn(f, fRing, sigma);
453 else /* fMerg < f <= fRing */
454 ampEff = w1 * pow(f / fMerg, mergPower) * (1. + epsilon_1 * v + epsilon_2 * v2);
455
456 /* now compute the phase */
457 psiEff = -2.*phi0
458 + 3./(128.*eta*v5)*(1 + params->psi2*v2
459 + params->psi3*v3 + params->psi4*v4
460 + params->psi5*v5 + params->psi6*v6
461 + params->psi7*v7 + params->psi8*v8);
462
463 /* generate the waveform */
464 ((*htilde)->data->data)[i] = amp0 * ampEff * cos(psiEff);
465 ((*htilde)->data->data)[i] += -I * amp0 * ampEff * sin(psiEff);
466 }
467
468 return XLAL_SUCCESS;
469}
470
471/*
472 * Private function to generate time-domain waveforms given coefficients
473 */
474static int IMRPhenomAGenerateTD(REAL8TimeSeries **h, const REAL8 phi0, const REAL8 deltaT, const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, const BBHPhenomParams *params) {
475 COMPLEX16FrequencySeries *htilde=NULL;
476 /* We will generate the waveform from a frequency which is lower than the
477 * f_min chosen. Also the cutoff frequency may be higher than the f_max. We
478 * will later apply a window function, and truncate the time-domain waveform
479 * below an instantaneous frequency f_min. */
480 REAL8 f_min_wide = EstimateSafeFMinForTD(m1, m2, f_min, deltaT);
481 const REAL8 f_max_wide = 0.5 / deltaT;
482 REAL8 deltaF;
483
484 if (EstimateSafeFMaxForTD(f_max, deltaT) > f_max_wide)
485 XLALPrintWarning("Warning: sampling rate too low to capture chosen f_max\n");
486 deltaF = 1. / (deltaT * NextPow2(EstimateIMRLength(m1, m2, f_min_wide, deltaT)));
487
488 /* generate in frequency domain */
489 if (IMRPhenomAGenerateFD(&htilde, phi0, deltaF, m1, m2, f_min_wide, f_max_wide, distance, params)) XLAL_ERROR(XLAL_EFUNC);
490
491 /* convert to time domain */
492 FDToTD(h, htilde, m1 + m2, deltaT, f_min, f_max, f_min_wide, f_max_wide);
493 XLALDestroyCOMPLEX16FrequencySeries(htilde);
494 if (!*h) XLAL_ERROR(XLAL_EFUNC);
495
496 return XLAL_SUCCESS;
497}
498
499/*
500 * Private function to generate time-domain waveforms given coefficients
501 */
502static int IMRPhenomBGenerateTD(REAL8TimeSeries **h, const REAL8 phi0, const REAL8 deltaT, const REAL8 m1, const REAL8 m2, const REAL8 chi, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, const BBHPhenomParams *params) {
503 REAL8 deltaF;
504 COMPLEX16FrequencySeries *htilde=NULL;
505 /* We will generate the waveform from a frequency which is lower than the
506 * f_min chosen. Also the cutoff frequency is higher than the f_max. We
507 * will later apply a window function, and truncate the time-domain waveform
508 * below an instantaneous frequency f_min. */
509 REAL8 f_min_wide = EstimateSafeFMinForTD(m1, m2, f_min, deltaT);
510 const REAL8 f_max_wide = 0.5 / deltaT;
511 if (EstimateSafeFMaxForTD(f_max, deltaT) > f_max_wide)
512 XLALPrintWarning("Warning: sampling rate (%" LAL_REAL8_FORMAT " Hz) too low for expected spectral content (%" LAL_REAL8_FORMAT " Hz) \n", deltaT, EstimateSafeFMaxForTD(f_max, deltaT));
513 deltaF = 1. / (deltaT * NextPow2(EstimateIMRLength(m1, m2, f_min_wide, deltaT)));
514
515 /* generate in frequency domain */
516 if (IMRPhenomBGenerateFD(&htilde, phi0, deltaF, m1, m2, chi, f_min_wide, f_max_wide, distance, params)) XLAL_ERROR(XLAL_EFUNC);
517
518 /* convert to time domain */
519 FDToTD(h, htilde, m1 + m2, deltaT, f_min, f_max, f_min_wide, f_max_wide);
520 XLALDestroyCOMPLEX16FrequencySeries(htilde);
521 if (!*h) XLAL_ERROR(XLAL_EFUNC);
522
523 return XLAL_SUCCESS;
524}
525
526/*
527 * Window and IFFT a FD waveform to TD, then window in TD.
528 * Requires that the FD waveform be generated outside of f_min and f_max.
529 * FD waveform is modified.
530 */
531static int FDToTD(REAL8TimeSeries **signalTD, const COMPLEX16FrequencySeries *signalFD, const REAL8 totalMass, const REAL8 deltaT, const REAL8 f_min, const REAL8 f_max, const REAL8 f_min_wide, const REAL8 f_max_wide) {
532 const LIGOTimeGPS gpstime_zero = {0, 0};
533 const size_t nf = signalFD->data->length;
534 const size_t nt = 2 * (nf - 1);
535 const REAL8 windowLength = 20. * totalMass * LAL_MTSUN_SI / deltaT;
536 const REAL8 winFLo = (f_min + f_min_wide) / 2.;
537 REAL8 winFHi = (f_max + f_max_wide) / 2.;
538 COMPLEX16 *FDdata = signalFD->data->data;
539 REAL8FFTPlan *revPlan;
540 REAL8 *TDdata;
541 size_t k;
542
543 /* check inputs */
544 if (f_min_wide >= f_min) XLAL_ERROR(XLAL_EDOM);
545
546 /* apply the softening window function */
547 if (winFHi > 0.5 / deltaT) winFHi = 0.5 / deltaT;
548 for (k = nf;k--;) {
549 const REAL8 f = k / (deltaT * nt);
550 REAL8 softWin = (1. + tanh(f - winFLo))
551 * (1. - tanh(f - winFHi)) / 4.;
552 FDdata[k] *= softWin;
553 }
554
555 /* allocate output */
556 *signalTD = XLALCreateREAL8TimeSeries("h", &gpstime_zero, 0.0, deltaT, &lalStrainUnit, nt);
557
558 /* Inverse Fourier transform */
559 revPlan = XLALCreateReverseREAL8FFTPlan(nt, 1);
560 if (!revPlan) {
561 XLALDestroyREAL8TimeSeries(*signalTD);
562 *signalTD = NULL;
563 XLAL_ERROR(XLAL_EFUNC);
564 }
565 XLALREAL8FreqTimeFFT(*signalTD, signalFD, revPlan);
566 XLALDestroyREAL8FFTPlan(revPlan);
567 if (!(*signalTD)) XLAL_ERROR(XLAL_EFUNC);
568
569 /* apply a linearly decreasing window at the end
570 * of the waveform in order to avoid edge effects. */
571 if (windowLength > (*signalTD)->data->length) XLAL_ERROR(XLAL_ERANGE);
572 TDdata = (*signalTD)->data->data;
573 for (k = windowLength; k--;)
574 TDdata[nt-k-1] *= k / windowLength;
575
576 return XLAL_SUCCESS;
577}
578
579/* return the index before the instantaneous frequency rises past target */
580static size_t find_instant_freq(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 target, const size_t start) {
581 size_t k = start + 1;
582 const size_t n = hp->data->length - 1;
583
584 /* complex_h := hp + i hc := A e^{-i \phi}. Hence \phi = arctan(-hc/hp)
585 and F = d \phi/dt = (hp hc' - hc hp') / (2 \pi A^2)
586 We use first order finite difference to compute the derivatives hp' and hc'*/
587 for (; k < n; k++) {
588 const REAL8 hpDot = (hp->data->data[k+1] - hp->data->data[k-1]) / (2 * hp->deltaT);
589 const REAL8 hcDot = (hc->data->data[k+1] - hc->data->data[k-1]) / (2 * hc->deltaT);
590 REAL8 f = hcDot * hp->data->data[k] - hpDot * hc->data->data[k];
591 f /= LAL_TWOPI;
592 f /= hp->data->data[k] * hp->data->data[k] + hc->data->data[k] * hc->data->data[k];
593 if (f >= target) return k - 1;
594 }
595 XLAL_ERROR(XLAL_EDOM);
596}
597
598/* Return the index of the sample at with the peak amplitude */
599static size_t find_peak_amp(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc) {
600 const REAL8 *hpdata = hp->data->data;
601 const REAL8 *hcdata = hc->data->data;
602 size_t k = hp->data->length;
603 size_t peak_ind = -1;
604 REAL8 peak_amp_sq = 0.;
605
606 for (;k--;) {
607 const REAL8 amp_sq = hpdata[k] * hpdata[k] + hcdata[k] * hcdata[k];
608 if (amp_sq > peak_amp_sq) {
609 peak_ind = k;
610 peak_amp_sq = amp_sq;
611 }
612 }
613 return peak_ind;
614}
615
616static int apply_phase_shift(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 shift) {
617 REAL8 *hpdata = hp->data->data;
618 REAL8 *hcdata = hc->data->data;
619 size_t k = hp->data->length;
620 const double cs = cos(shift);
621 const double ss = sin(shift);
622
623 for (;k--;) {
624 const REAL8 temp_hpdata = hpdata[k] * cs - hcdata[k] * ss;
625 hcdata[k] = hpdata[k] * ss + hcdata[k] * cs;
626 hpdata[k] = temp_hpdata;
627 }
628 return 0;
629}
630
631/* Apply the inclination-angle weighting to the two polarizations
632* Ref. Eq.(63) of B.S. Sathyaprakash and Bernard F. Schutz,
633* “Physics, Astrophysics and Cosmology with Gravitational Waves”,
634* Living Rev. Relativity, 12, (2009), 2. [Online Article]: cited 27 Nov 2013,
635* http://www.livingreviews.org/lrr-2009-2
636*/
637static int apply_inclination(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 inclination) {
638 REAL8 inclFacPlus, inclFacCross;
639 REAL8 *hpdata = hp->data->data;
640 REAL8 *hcdata = hc->data->data;
641 size_t k = hp->data->length;
642
643 inclFacCross = cos(inclination);
644 inclFacPlus = 0.5 * (1. + inclFacCross * inclFacCross);
645 for (;k--;) {
646 hpdata[k] *= inclFacPlus;
647 hcdata[k] *= inclFacCross;
648 }
649
650 return XLAL_SUCCESS;
651}
652
653/* ***********************************************************************************/
654/* END OF THE REVIEWED CODE: Below is the code for generating the IMRPhenomB metric */
655/* ***********************************************************************************/
656
657/*
658 * Structure containing the coefficients of various powers of the Fourier frequency
659 * f in the derivative of the amplitude and phase of IMRPhenomB w.r.t. the parameters
660 * of the binary
661 */
662typedef struct tagIMRDerivCoeffs{
663 REAL8 dA1dM_0; /*Coef of k = 0 term in (2k-7)/6 expansion of inspiral phase of amplitude*/
664 REAL8 dA1dM_1; /*Coef of k = 2 term in (2k-7)/6 expansion of inspiral phase of amplitude*/
665 REAL8 dA1dM_2; /*Coef of k = 3 term in (2k-7)/6 expansion of inspiral phase of amplitude*/
666 REAL8 dA1deta_0;
667 REAL8 dA1deta_1;
668 REAL8 dA1deta_2;
669 REAL8 dA1dchi_0;
670 REAL8 dA1dchi_1;
671 REAL8 dA1dchi_2;
672 REAL8 dA2dM_0; /*Coef of k = 0 term in (k-2)/3 expansion of merger phase of amplitude*/
673 REAL8 dA2dM_1; /*Coef of k = 1 term in (k-2)/3 expansion of merger phase of amplitude*/
674 REAL8 dA2dM_2; /*Coef of k = 2 term in (k-2)/3 expansion of merger phase of amplitude*/
675 REAL8 dA2deta_0;
676 REAL8 dA2deta_1;
677 REAL8 dA2deta_2;
678 REAL8 dA2dchi_0;
679 REAL8 dA2dchi_1;
680 REAL8 dA2dchi_2;
681 REAL8 dA3dMnum_0;
682 REAL8 dA3dMnum_1;
683 REAL8 dA3dMnum_2;
684 REAL8 dA3detanum_0; /*Coef of k = 0 in pow(f,k) expansion of numerator of derivative of ringdown amplitude */
685 REAL8 dA3detanum_1;
686 REAL8 dA3detanum_2;
687 REAL8 dA3dchinum_0;
688 REAL8 dA3dchinum_1;
689 REAL8 dA3dchinum_2;
690 REAL8 dA3denom_0;
691 REAL8 dA3denom_1;
692 REAL8 dA3denom_2;
693 REAL8 dA3denom_3;
694 REAL8 dA3denom_4;
695 REAL8 dPsidM_0; /*Coef of k = 0 in frequency series of Phase */
696 REAL8 dPsidM_1;
697 REAL8 dPsidM_2;
698 REAL8 dPsidM_3;
699 REAL8 dPsidM_4;
700 REAL8 dPsidM_5;
701 REAL8 dPsidM_6;
702 REAL8 dPsidM_7;
703 REAL8 dPsideta_0;
704 REAL8 dPsideta_1;
705 REAL8 dPsideta_2;
706 REAL8 dPsideta_3;
707 REAL8 dPsideta_4;
708 REAL8 dPsideta_5;
709 REAL8 dPsideta_6;
710 REAL8 dPsideta_7;
711 REAL8 dPsidchi_0;
712 REAL8 dPsidchi_1;
713 REAL8 dPsidchi_2;
714 REAL8 dPsidchi_3;
715 REAL8 dPsidchi_4;
716 REAL8 dPsidchi_5;
717 REAL8 dPsidchi_6;
718 REAL8 dPsidchi_7;
719 REAL8 Wm;
720 REAL8 Wr;
721 REAL8 dWmdM;
722 REAL8 dWmdEta;
723 REAL8 dWmdChi;
724 REAL8 dWrdM;
725 REAL8 dWrdEta;
726 REAL8 dWrdChi;
727 REAL8 alpha2;
728 REAL8 alpha3;
729 REAL8 eps1;
730 REAL8 eps2;
731}
732IMRDerivCoeffs;
733
734static gsl_matrix *XLALSimIMRPhenomBFisherMatrix(
735 const REAL8 m1, /**< component mass 1 (M_sun) */
736 const REAL8 m2, /**< component mass 2 (M_sun) */
737 const REAL8 chi, /**< effective spin parameter of IMRPhenomB: chi = (m1 chi1 + m2 chi2)/(m1+m2) */
738 const REAL8 fLow, /**< low-frequency cutoff (Hz) */
739 const REAL8FrequencySeries *Sh /**< PSD in strain per root Hertz */
740);
741
742static IMRDerivCoeffs *XLALComputeIMRPhenomBDerivativeCoeffs(
743 const REAL8 m1, /**< component mass 1 (M_sun) */
744 const REAL8 m2, /**< component mass 2 (M_sun) */
745 const REAL8 chi, /**< effective spin parameter of IMRPhenomB: chi = (m1 chi1 + m2 chi2)/(m1+m2) */
746 BBHPhenomParams *params /**< phenomenological parameters of IMRPhenomB */
747);
748
749static gsl_matrix *XLALSimIMRPhenomBProjectExtrinsicParam(
750 gsl_matrix *g /**< Fisher matrix of IMRPhenomB */
751);
752
753/*
754 * Compute the coefficients of the series expansion (in Fourier frequency) of the derivatives of the
755 * IMRPhenomB waveform with respect to parameters M, eta, chi
756 */
757static IMRDerivCoeffs *XLALComputeIMRPhenomBDerivativeCoeffs(
758 const REAL8 m1, /**< component mass 1 (M_sun) */
759 const REAL8 m2, /**< component mass 2 (M_sun) */
760 const REAL8 chi, /**< effective spin parameter of IMRPhenomB: chi = (m1 chi1 + m2 chi2)/(m1+m2) */
761 BBHPhenomParams *params /**< phenomenological parameters of IMRPhenomB */
762){
763 /* allocate memory for the structure IMRDerivCoeffs */
764 IMRDerivCoeffs *derCoeffs = (IMRDerivCoeffs *) XLALMalloc(sizeof(IMRDerivCoeffs));
765 if (!derCoeffs) XLAL_ERROR_NULL(XLAL_EFUNC);
766 memset(derCoeffs, 0, sizeof(IMRDerivCoeffs));
767
768 /* define mass and symmetric mass ratio */
769 REAL8 M = (m1 + m2)*LAL_MTSUN_SI;
770 REAL8 eta = m1*m2/((m1 + m2)*(m1 + m2));
771
772 /* extract various phenomenological parameters */
773 REAL8 f1 = params->fMerger;
774 REAL8 f2 = params->fRing;
775 REAL8 sigma = params->sigma;
776
777 /*Define various terms*/
778 REAL8 M_p2 = M*M;
779 REAL8 M_p1by6 = sqrt(cbrt(M));
780 REAL8 M_p5by6 = M/M_p1by6;
781 REAL8 M_p7by6 = M*M_p1by6;
782 REAL8 M_p1by3 = cbrt(M);
783 REAL8 M_p2by3 = M_p1by3*M_p1by3;
784 REAL8 sqrt_eta = sqrt(eta);
785 REAL8 eta_p2 = eta*eta;
786 REAL8 chi_x_eta = chi*eta;
787 REAL8 chi_p2 = chi*chi;
788 REAL8 f2_p2 = f2*f2;
789 REAL8 f2_p3 = f2_p2*f2;
790 REAL8 f2_p4 = f2_p3*f2;
791 REAL8 f1_p7by6 = f1*sqrt(cbrt(f1));
792 REAL8 f1_p3by2 = sqrt(f1*f1*f1);
793 REAL8 f1_p3 = f1*f1*f1;
794 REAL8 sigma_p2 = sigma*sigma;
795 REAL8 sigma_p3 = sigma_p2*sigma;
796 REAL8 f1M_p1by3 = cbrt(f1*M);
797 REAL8 f1M_p2by3 = f1M_p1by3*f1M_p1by3;
798 REAL8 f2M_p1by3 = cbrt(f2*M);
799 REAL8 f2M_p2by3 = f2M_p1by3*f2M_p1by3;
800 REAL8 f1byf2_p1by3 = cbrt(f1/f2);
801 REAL8 f1byf2_p2by3 = f1byf2_p1by3*f1byf2_p1by3;
802
803 /* derivatives w.r.t physical parameters of fMerger */
804 REAL8 dF1dM = -f1/M;
805 REAL8 dF1dEta = (0.3183098861837907*(0.64365 + 0.82696*chi - 0.27063*chi_p2 - 0.116436*eta - 7.8692*chi_x_eta - 21.2748*eta_p2))/M;
806 REAL8 dF1dChi = (0.3183098861837907*(0.9666699/pow(1. - chi,0.783) - 0.91546/pow(1. - chi,0.74)))/M + (0.3183098861837907*(0.82696*eta - 0.54126*chi_x_eta - 3.9346*eta_p2))/M;
807
808 /* derivatives w.r.t physical parameters of fRing */
809 REAL8 dF2dM = -f2/M;
810 REAL8 dF2dEta = (0.3183098861837907*(0.1469 - 0.12281*chi - 0.026091*chi_p2 - 0.0498*eta + 0.34026*chi_x_eta + 6.9756*eta_p2))/M;
811 REAL8 dF2dChi = 0.03008028424436822/(pow(1. - chi,0.7)*M) + (0.3183098861837907*(-0.12281*eta - 0.052182*chi_x_eta + 0.17013*eta_p2))/M;
812
813 /* derivatives of w.r.t physical parameters Sigma parameter */
814 REAL8 dSigmadM = -sigma/M;
815 REAL8 dSigmadEta = (0.3183098861837907*(-0.40979 - 0.035226*chi + 0.10082*chi_p2 + 3.6572*eta - 0.040338*chi_x_eta - 8.6094*eta_p2))/M;
816 REAL8 dSigmadChi = (-0.03580986219567645*(1. - 0.63*pow(1. - chi,0.3)))/(pow(1. - chi,0.55)*M) + 0.01504014212218411/(pow(1. - chi,0.24999999999999994)*M) + (0.3183098861837907*(-0.035226*eta + 0.20164*chi_x_eta - 0.020169*eta_p2))/M;
817
818 /* PN corrections to the frequency domain amplitude of the (2,2) mode */
819 derCoeffs->alpha2 = -323./224. + 451.*eta/168.;
820 derCoeffs->alpha3 = (27./8. - 11.*eta/6.)*chi;
821
822 /* spin-dependent corrections to the merger amplitude */
823 derCoeffs->eps1 = 1.4547*chi - 1.8897;
824 derCoeffs->eps2 = -1.8153*chi + 1.6557;
825
826 /* derivatives of PN corrections to the frequency domain amplitude of the (2,2) mode */
827 REAL8 dAlpha2dEta = 2.6845238095238093;
828 REAL8 dAlpha3dEta = -1.8333333333333333*chi;
829 REAL8 dAlpha3dChi = 3.375 - 1.8333333333333333*eta;
830
831 /* derivatives of spin-dependent corrections to the merger amplitude */
832 REAL8 dEps1dChi = 1.4547;
833 REAL8 dEps2dChi = -1.8153;
834
835 /* Recurring expressions in normalization */
836 REAL8 norm_expr_1 = (1. + 1.4645918875615231*derCoeffs->eps1*f2M_p1by3 + 2.1450293971110255*derCoeffs->eps2*f2M_p2by3);
837 REAL8 norm_expr_2 = 0.33424583931521507*sqrt_eta*f1byf2_p2by3*M_p5by6;
838 REAL8 norm_expr_3 = 1. + 1.4645918875615231*derCoeffs->eps1*f1M_p1by3 + 2.1450293971110255*derCoeffs->eps2*f1M_p2by3;
839 REAL8 norm_expr_4 = sqrt_eta*f1byf2_p2by3*M_p5by6*norm_expr_1/f1_p7by6;
840 REAL8 norm_expr_5 = 0.22283055954347672*sqrt_eta*M_p5by6*norm_expr_1*sigma;
841
842 /* normalisation constant of the inspiral/merger amplitude */
843 derCoeffs->Wm = (1. + 3.141592653589793*derCoeffs->alpha3*f1*M + 2.1450293971110255*derCoeffs->alpha2*f1M_p2by3)/(norm_expr_3);
844 derCoeffs->Wr = derCoeffs->Wm*(norm_expr_2*norm_expr_1*sigma)/f1_p7by6;
845
846 /* compute the derivatives of the phenomenological parameters w.r.t the physical parameters.*/
847 REAL8 dWmDen = norm_expr_3*norm_expr_3; /* This is the common denominatior to dWmdM and dWmdeta dWmdchi */
848
849 derCoeffs->dWmdM = 0.;
850 derCoeffs->dWmdEta = ((-0.48819729585384103*dF1dEta*M*(derCoeffs->eps1 + 2.9291837751230463*derCoeffs->eps2*f1M_p1by3)*(1. + 3.141592653589793*derCoeffs->alpha3*f1*M + 2.1450293971110255*derCoeffs->alpha2*f1M_p2by3))/f1M_p2by3 + (3.141592653589793*derCoeffs->alpha3*dF1dEta*M + 3.141592653589793*dAlpha3dEta*f1*M + (1.430019598074017*derCoeffs->alpha2*dF1dEta*M)/f1M_p1by3 + 2.1450293971110255*dAlpha2dEta*f1M_p2by3)*(norm_expr_3))/dWmDen;
851 derCoeffs->dWmdChi = ((3.141592653589793*derCoeffs->alpha3*dF1dChi*M + 3.141592653589793*dAlpha3dChi*f1*M + (1.430019598074017*derCoeffs->alpha2*dF1dChi*M)/f1M_p1by3)*(norm_expr_3) - (0.48819729585384103*M*(1. + 3.141592653589793*derCoeffs->alpha3*f1*M + 2.1450293971110255*derCoeffs->alpha2*f1M_p2by3)*(3.*f1*(dEps1dChi + 1.4645918875615231*dEps2dChi*f1M_p1by3) + dF1dChi*(derCoeffs->eps1 + 2.9291837751230463*derCoeffs->eps2*f1M_p1by3)))/f1M_p2by3)/dWmDen;
852
853 derCoeffs->dWrdM = derCoeffs->Wm*((0.33424583931521507*dSigmadM*norm_expr_4) + (0.27853819942934593*sqrt_eta*f1byf2_p2by3*norm_expr_1*sigma)/(f1_p7by6*M_p1by6) + (norm_expr_5*((-1.*dF2dM*f1)/f2_p2 + dF1dM/f2))/(f1_p7by6*f1byf2_p1by3) - (0.38995347920108425*dF1dM*norm_expr_4*sigma)/(f1) + (norm_expr_2*((0.48819729585384103*derCoeffs->eps1*(f2 + dF2dM*M))/f2M_p2by3 + (1.430019598074017*derCoeffs->eps2*(f2 + dF2dM*M))/f2M_p1by3)*sigma)/f1_p7by6);
854 derCoeffs->dWrdEta = derCoeffs->Wm*((0.33424583931521507*dSigmadEta*norm_expr_4) + (norm_expr_2*((0.48819729585384103*dF2dEta*derCoeffs->eps1*M)/f2M_p2by3 + (1.430019598074017*dF2dEta*derCoeffs->eps2*M)/f2M_p1by3)*sigma)/f1_p7by6 + (norm_expr_5*((-1.*dF2dEta*f1)/f2_p2 + dF1dEta/f2))/(f1_p7by6*f1byf2_p1by3) - (0.38995347920108425*dF1dEta*norm_expr_4*sigma)/(f1) + (0.16712291965760753*f1byf2_p2by3*M_p5by6*norm_expr_1*sigma)/(sqrt_eta*f1_p7by6)) + derCoeffs->Wr*derCoeffs->dWmdEta/derCoeffs->Wm;
855 derCoeffs->dWrdChi = derCoeffs->Wm*((0.33424583931521507*dSigmadChi*norm_expr_4) + (norm_expr_2*((0.48819729585384103*dF2dChi*derCoeffs->eps1*M)/f2M_p2by3 + (1.430019598074017*dF2dChi*derCoeffs->eps2*M)/f2M_p1by3 + 1.4645918875615231*dEps1dChi*f2M_p1by3 + 2.1450293971110255*dEps2dChi*f2M_p2by3)*sigma)/f1_p7by6 + (norm_expr_5*((-1.*dF2dChi*f1)/f2_p2 + dF1dChi/f2))/(f1_p7by6*f1byf2_p1by3) - (0.38995347920108425*dF1dChi*norm_expr_4*sigma)/(f1)) + derCoeffs->Wr*derCoeffs->dWmdChi/derCoeffs->Wm;
856
857 /* ------------------------------------------------------------------- */
858 /* compute the coefficients of the series expansion of the derivatives */
859 /* ------------------------------------------------------------------- */
860
861 /* Coefficients of the derivatives of the inspiral amplitude */
862 derCoeffs->dA1dM_0 = (0.1773229251163862*sqrt_eta)/M_p1by6;
863 derCoeffs->dA1dM_1 = 0.6846531968814576*derCoeffs->alpha2*sqrt(eta*M);
864 derCoeffs->dA1dM_2 = 1.2255680774891218*derCoeffs->alpha3*sqrt_eta*M_p5by6;
865
866 derCoeffs->dA1deta_0 = (0.10639375506983172*M_p5by6)/sqrt_eta;
867 derCoeffs->dA1deta_1 = (0.22821773229381923*(derCoeffs->alpha2 + 2.*dAlpha2dEta*eta)*sqrt(M*M*M))/sqrt_eta;
868 derCoeffs->dA1deta_2 = (0.33424583931521507*(derCoeffs->alpha3 + 2.*dAlpha3dEta*eta)*M_p5by6*M)/sqrt_eta;
869
870 derCoeffs->dA1dchi_0 = 0.;
871 derCoeffs->dA1dchi_1 = 0.;
872 derCoeffs->dA1dchi_2 = 0.6684916786304301*dAlpha3dChi*sqrt_eta*M_p5by6*M;
873
874 /*Coefficients of the derivatives of the merger amplitude */
875 derCoeffs->dA2dM_0 = (0.03546458502327724*sqrt_eta*(-3.*dF1dM*M*derCoeffs->Wm + f1*(6.*derCoeffs->dWmdM*M + 5.*derCoeffs->Wm)))/(f1_p3by2*M_p1by6);
876 derCoeffs->dA2dM_1 = (0.05194114352082774*derCoeffs->eps1*sqrt_eta*M_p1by6*(-3.*dF1dM*M*derCoeffs->Wm + f1*(6.*derCoeffs->dWmdM*M + 7.*derCoeffs->Wm)))/f1_p3by2;
877 derCoeffs->dA2dM_2 = (0.22821773229381923*derCoeffs->eps2*sqrt(eta*M)*(-1.*dF1dM*M*derCoeffs->Wm + f1*(2.*derCoeffs->dWmdM*M + 3.*derCoeffs->Wm)))/f1_p3by2;
878
879 derCoeffs->dA2deta_0 = (0.10639375506983172*M_p5by6*(-1.*dF1dEta*eta*derCoeffs->Wm + f1*(2.*derCoeffs->dWmdEta*eta + derCoeffs->Wm)))/sqrt(eta*f1_p3);
880 derCoeffs->dA2deta_1 = (0.1558234305624832*derCoeffs->eps1*M_p7by6*(-1.*dF1dEta*eta*derCoeffs->Wm + f1*(2.*derCoeffs->dWmdEta*eta + derCoeffs->Wm)))/sqrt(eta*f1_p3);
881 derCoeffs->dA2deta_2 = (0.22821773229381923*derCoeffs->eps2*sqrt(M*M*M)*(-1.*dF1dEta*eta*derCoeffs->Wm + f1*(2.*derCoeffs->dWmdEta*eta + derCoeffs->Wm)))/sqrt(eta*f1_p3);
882
883 derCoeffs->dA2dchi_0 = (0.10639375506983172*sqrt_eta*M_p5by6*(2.*derCoeffs->dWmdChi*f1 - 1.*dF1dChi*derCoeffs->Wm))/f1_p3by2;
884 derCoeffs->dA2dchi_1 = (0.1558234305624832*sqrt_eta*M_p7by6*(-1.*dF1dChi*derCoeffs->eps1*derCoeffs->Wm + 2.*f1*(derCoeffs->dWmdChi*derCoeffs->eps1 + dEps1dChi*derCoeffs->Wm)))/f1_p3by2;
885 derCoeffs->dA2dchi_2 = 0.22821773229381923*sqrt_eta*sqrt(M*M*M/(f1*f1*f1))*(-1.*dF1dChi*derCoeffs->eps2*derCoeffs->Wm + 2.*f1*(derCoeffs->dWmdChi*derCoeffs->eps2 + dEps2dChi*derCoeffs->Wm));
886
887 /* Coefficients of the derivatives of the ringdown amplitude (numerator) */
888 derCoeffs->dA3dMnum_0 = 8.*derCoeffs->dWrdM*f2_p2*sigma + 2.*derCoeffs->dWrdM*sigma_p3 + (8.*dSigmadM*f2_p2 - 16.*dF2dM*f2*sigma - 2.*dSigmadM*sigma_p2)*derCoeffs->Wr;
889 derCoeffs->dA3dMnum_1 = -16.*derCoeffs->dWrdM*f2*sigma - (16.*dSigmadM*f2 - 16.*dF2dM*sigma)*derCoeffs->Wr;
890 derCoeffs->dA3dMnum_2 = 8.*derCoeffs->dWrdM*sigma + 8.*dSigmadM*derCoeffs->Wr;
891
892 derCoeffs->dA3detanum_0 = 8.*derCoeffs->dWrdEta*f2_p2*sigma + 2.*derCoeffs->dWrdEta*sigma_p3 + (8.*dSigmadEta*f2_p2 - 16.*dF2dEta*f2*sigma - 2.*dSigmadEta*sigma_p2)*derCoeffs->Wr;
893 derCoeffs->dA3detanum_1 = -16.*derCoeffs->dWrdEta*f2*sigma - (16.*dSigmadEta*f2 - 16.*dF2dEta*sigma)*derCoeffs->Wr;
894 derCoeffs->dA3detanum_2 = 8.*derCoeffs->dWrdEta*sigma + 8.*dSigmadEta*derCoeffs->Wr;
895
896 derCoeffs->dA3dchinum_0 = 8.*derCoeffs->dWrdChi*f2_p2*sigma + 2.*derCoeffs->dWrdChi*sigma_p3 + (8.*dSigmadChi*f2_p2 - 16.*dF2dChi*f2*sigma - 2.*dSigmadChi*sigma_p2)*derCoeffs->Wr;
897 derCoeffs->dA3dchinum_1 = -16.*derCoeffs->dWrdChi*f2*sigma - (16.*dSigmadChi*f2 - 16.*dF2dChi*sigma)*derCoeffs->Wr;
898 derCoeffs->dA3dchinum_2 = 8.*derCoeffs->dWrdChi*sigma + 8.*dSigmadChi*derCoeffs->Wr;
899
900 /* Coefficients of the derivatives of the ringdown amplitude (denominator) */
901 derCoeffs->dA3denom_0 = 50.26548245743669*f2_p4 + 25.132741228718345*f2_p2*sigma_p2 + 3.141592653589793*sigma_p2*sigma_p2;
902 derCoeffs->dA3denom_1 = -201.06192982974676*f2_p3 - 50.26548245743669*f2*sigma_p2;
903 derCoeffs->dA3denom_2 = 301.59289474462014*f2_p2 + 25.132741228718345*sigma_p2;
904 derCoeffs->dA3denom_3 = -201.06192982974676*f2;
905 derCoeffs->dA3denom_4 = 50.26548245743669;
906
907 /* Coefficients of the derivatives of the phase */
908 REAL8 dPsi2dEta = -920.91 + 492.13*chi + 135.03*chi_p2 + 13483.8*eta - 2106.8*chi_x_eta - 40191.*eta_p2;
909 REAL8 dPsi2dChi = 492.13*eta + 270.06*chi_x_eta - 1053.4*eta_p2;
910 REAL8 dPsi3dEta = 17022. - 9565.9*chi - 2182.1*chi_p2 - 242740.*eta + 41504.*chi_x_eta + 715770.*eta_p2;
911 REAL8 dPsi3dChi = 37.666666666666664 - 9565.9*eta - 4364.2*chi_x_eta + 20752.*eta_p2;
912 REAL8 dPsi4dEta = -125440. + 75066.*chi + 13382.*chi_p2 + 1.74708e6*eta - 331460.*chi_x_eta - 5.0808e6*eta_p2;
913 REAL8 dPsi4dChi = -101.25*chi + 75066.*eta + 26764.*chi_x_eta - 165730.*eta_p2;
914 REAL8 dPsi6dEta = -889770. + 631020.*chi + 50676.*chi_p2 + 1.19616e7*eta - 2.8296e6*chi_x_eta - 3.383999999999999e7*eta_p2;
915 REAL8 dPsi6dChi = 631020.*eta + 101352.*chi_x_eta - 1.4148e6*eta_p2;
916 REAL8 dPsi7dEta = 869600. - 670980.*chi - 30082.*chi_p2 - 1.16758e7*eta + 3.029e6*chi_x_eta + 3.2673e7*eta_p2;
917 REAL8 dPsi7dChi = -670980.*eta - 60164.*chi_x_eta + 1.5145e6*eta_p2;
918
919 derCoeffs->dPsidM_0 = -0.005796647796902313/(eta*M_p2by3*M*M);
920 derCoeffs->dPsidM_1 = 0.;
921 derCoeffs->dPsidM_2 = (-0.007460387957432594*params->psi2)/(eta*M_p2);
922 derCoeffs->dPsidM_3 = (-0.007284282453678307*params->psi3)/(eta*M_p5by6*M_p5by6);
923 derCoeffs->dPsidM_4 = (-0.005334250494181998*params->psi4)/(eta*M_p1by3*M);
924 derCoeffs->dPsidM_5 = 0.;
925 derCoeffs->dPsidM_6 = (0.0114421241215744*params->psi6)/(eta*M_p2by3);
926 derCoeffs->dPsidM_7 = (0.03351608432985977*params->psi7)/(eta*M_p1by3);
927
928 derCoeffs->dPsideta_0 = -0.0034779886781413872/(eta_p2*M_p5by6*M_p5by6);
929 derCoeffs->dPsideta_1 = 0.;
930 derCoeffs->dPsideta_2 = (0.007460387957432594*(dPsi2dEta*eta - 1.*params->psi2))/(eta_p2*M);
931 derCoeffs->dPsideta_3 = (0.010926423680517461*(dPsi3dEta*eta - 1.*params->psi3))/(eta_p2*M_p2by3);
932 derCoeffs->dPsideta_4 = (0.016002751482545992*(dPsi4dEta*eta - 1.*params->psi4))/(eta_p2*M_p1by3);
933 derCoeffs->dPsideta_5 = 0.;
934 derCoeffs->dPsideta_6 = (0.0343263723647232*M_p1by3*(dPsi6dEta*eta - 1.*params->psi6))/eta_p2;
935 derCoeffs->dPsideta_7 = (0.050274126494789656*M_p2by3*(dPsi7dEta*eta - 1.*params->psi7))/eta_p2;
936
937 derCoeffs->dPsidchi_0 = 0.;
938 derCoeffs->dPsidchi_1 = 0.;
939 derCoeffs->dPsidchi_2 = (0.007460387957432594*dPsi2dChi)/(eta*M);
940 derCoeffs->dPsidchi_3 = (0.010926423680517461*dPsi3dChi)/(eta*M_p2by3);
941 derCoeffs->dPsidchi_4 = (0.016002751482545992*dPsi4dChi)/(eta*M_p1by3);
942 derCoeffs->dPsidchi_5 = 0.;
943 derCoeffs->dPsidchi_6 = (0.0343263723647232*dPsi6dChi*M_p1by3)/eta;
944 derCoeffs->dPsidchi_7 = (0.050274126494789656*dPsi7dChi*M_p2by3)/eta;
945
946 return derCoeffs;
947};
948
949/*
950 * Compute the Fisher information matrix of the IMRPhenomB waveform in the
951 * M, eta, chi, t0, phi0 parameter space, for an SNR=1/sqrt(2) signal.
952 */
953
954static gsl_matrix *XLALSimIMRPhenomBFisherMatrix(
955 const REAL8 m1, /**< component mass 1 (M_sun) */
956 const REAL8 m2, /**< component mass 2 (M_sun) */
957 const REAL8 chi, /**< effective spin parameter of IMRPhenomB: chi = (m1 chi1 + m2 chi2)/(m1+m2) */
958 const REAL8 fLow, /**< low-frequency cutoff (Hz) */
959 const REAL8FrequencySeries *Sh /**< PSD in strain per root Hertz */
960){
961
962 IMRDerivCoeffs *coef;
963 BBHPhenomParams *params;
964
965 /* check inputs */
966 if (m1 <= 0) XLAL_ERROR_NULL(XLAL_EDOM);
967 if (m2 <= 0) XLAL_ERROR_NULL(XLAL_EDOM);
968 if (chi <= -1.) XLAL_ERROR_NULL(XLAL_EDOM);
969 if (chi >= 1.) XLAL_ERROR_NULL(XLAL_EDOM);
970 if (fLow <= 0) XLAL_ERROR_NULL(XLAL_EDOM);
971
972 /* define physical parameters and required coefficients */
973 REAL8 M = (m1 + m2)*LAL_MTSUN_SI;
974 REAL8 eta = m1*m2/((m1+m2)*(m1+m2));
975 REAL8 piM = LAL_PI * M;
976 REAL8 AmpCoef = cbrt(1./(LAL_PI*LAL_PI))*sqrt(5.*eta/24.)*pow(M,5.0/6.0);
977
978 /* compute the phenomenological parameters describing the IMRPhenomB waveform */
979 params = ComputeIMRPhenomBParams(m1, m2, chi);
980
981 /* extract various phenomenological parameters */
982 REAL8 f1 = params->fMerger;
983 REAL8 f2 = params->fRing;
984 REAL8 sigma = params->sigma;
985 REAL8 f3 = params->fCut;
986
987 /* create a view of the PSD between flow and fCut */
988 REAL8 df = Sh->deltaF;
989 size_t nBins = (f3 - fLow) / df;
990
991 /* compute the coefficients of the series expansion of the derivatives */
992 coef = XLALComputeIMRPhenomBDerivativeCoeffs(m1, m2, chi, params);
993
994 /* initialise the components of the Fisher matrix */
995 REAL8 gamma_MM = 0.;
996 REAL8 gamma_MEta= 0.;
997 REAL8 gamma_MChi = 0.;
998 REAL8 gamma_MT0 = 0.;
999 REAL8 gamma_MPhi0 = 0.;
1000 REAL8 gamma_EtaEta = 0.;
1001 REAL8 gamma_EtaChi = 0.;
1002 REAL8 gamma_EtaT0 = 0.;
1003 REAL8 gamma_EtaPhi0= 0.;
1004 REAL8 gamma_ChiChi = 0.;
1005 REAL8 gamma_ChiT0 = 0.;
1006 REAL8 gamma_ChiPhi0 = 0.;
1007 REAL8 gamma_T0T0 = 0.;
1008 REAL8 gamma_T0Phi0 = 0.;
1009 REAL8 gamma_Phi0Phi0 = 0.;
1010
1011 REAL8 amp, dAdM, dAdEta, dAdChi, dA3Den;
1012 REAL8 f1_pm7by6 = 1./(f1*cbrt(sqrt(f1)));
1013 REAL8 f1_p2by3 = cbrt(f1*f1);
1014 REAL8 sigma_p2 = sigma*sigma;
1015 REAL8 amp_merg_const = coef->Wm*AmpCoef*f1_pm7by6;
1016 REAL8 amp_ring_const = (coef->Wr*sigma)/LAL_TWOPI;
1017
1018 /* compute derivatives over a freq vector from fLow to fCut with frequency resolution df
1019 * * and use this to compute the Fisher matrix */
1020 REAL8 f = fLow;
1021
1022 for (size_t k=0; k<nBins; k++) {
1023
1024 REAL8 v = cbrt(piM*f);
1025 REAL8 v_p2 = v*v;
1026 REAL8 v_p3 = v_p2*v;
1027 REAL8 f_p2 = f*f;
1028 REAL8 f_p3 = f_p2*f;
1029 REAL8 f_p4 = f_p3*f;
1030 REAL8 f_p1by3 = cbrt(f);
1031 REAL8 f_p2by3 = f_p1by3*f_p1by3;
1032
1033 REAL8 f_pm1 = 1./f;
1034 REAL8 f_pm1by2 = sqrt(f_pm1);
1035 REAL8 f_pm1by3 = 1./f_p1by3;
1036 REAL8 f_pm2by3 = 1./f_p2by3;
1037
1038 REAL8 f_pm1by6 = cbrt(sqrt(f_pm1));
1039 REAL8 f_pm5by3 = f_pm2by3*f_pm1;
1040 REAL8 f_pm7by6 = f_pm1by6*f_pm1;
1041
1042 REAL8 fbyf1_pm2by3 = f1_p2by3*f_pm2by3;
1043
1044 /* compute derivatives of the amplitude w.r.t M, eta and chi */
1045 if (f <= f1) {
1046 amp = AmpCoef*f_pm7by6*(1. + coef->alpha2*v_p2 + coef->alpha3*v_p3);
1047 dAdM = coef->dA1dM_0*f_pm7by6 + coef->dA1dM_1*f_pm1by2 + coef->dA1dM_2*f_pm1by6;
1048 dAdEta = coef->dA1deta_0*f_pm7by6 + coef->dA1deta_1*f_pm1by2 + coef->dA1deta_2*f_pm1by6;
1049 dAdChi = coef->dA1dchi_0*f_pm7by6 + coef->dA1dchi_1*f_pm1by2 + coef->dA1dchi_2*f_pm1by6;
1050 }
1051 else if ((f1<f) && (f<=f2)) {
1052 amp = amp_merg_const*fbyf1_pm2by3*(1. + coef->eps1*v + coef->eps2*v_p2);
1053 dAdM = coef->dA2dM_0*f_pm2by3 + coef->dA2dM_1*f_pm1by3 + coef->dA2dM_2;
1054 dAdEta = coef->dA2deta_0*f_pm2by3 + coef->dA2deta_1*f_pm1by3 + coef->dA2deta_2;
1055 dAdChi = coef->dA2dchi_0*f_pm2by3 + coef->dA2dchi_1*f_pm1by3 + coef->dA2dchi_2;
1056 }
1057 else {
1058 amp = amp_ring_const/((f-f2)*(f-f2) + sigma_p2*0.25);
1059 dA3Den = coef->dA3denom_0 + coef->dA3denom_1*f + coef->dA3denom_2*f_p2 + coef->dA3denom_3*f_p3 + coef->dA3denom_4*f_p4;
1060 dAdM = (coef->dA3dMnum_0 + coef->dA3dMnum_1*f + coef->dA3dMnum_2*f_p2)/dA3Den;
1061 dAdEta = (coef->dA3detanum_0 + coef->dA3detanum_1*f + coef->dA3detanum_2*f_p2)/dA3Den;
1062 dAdChi = (coef->dA3dchinum_0 + coef->dA3dchinum_1*f + coef->dA3dchinum_2*f_p2)/dA3Den;
1063 }
1064
1065 /* compute derivatives of the phase w.r.t M, eta and chi */
1066 REAL8 dPsidM = coef->dPsidM_0*f_pm5by3 + coef->dPsidM_2*f_pm1 + coef->dPsidM_3*f_pm2by3 + coef->dPsidM_4*f_pm1by3 + coef->dPsidM_6*f_p1by3 + coef->dPsidM_7*f_p2by3;
1067 REAL8 dPsidEta = coef->dPsideta_0*f_pm5by3 + coef->dPsideta_2*f_pm1 + coef->dPsideta_3*f_pm2by3 + coef->dPsideta_4*f_pm1by3 + coef->dPsideta_6*f_p1by3 + coef->dPsideta_7*f_p2by3;
1068 REAL8 dPsidChi = coef->dPsidchi_0*f_pm5by3 + coef->dPsidchi_2*f_pm1 + coef->dPsidchi_3*f_pm2by3 + coef->dPsidchi_4*f_pm1by3 + coef->dPsidchi_6*f_p1by3 + coef->dPsidchi_7*f_p2by3;
1069 REAL8 dPsidT0 = LAL_TWOPI * f;
1070
1071 /* compute the elements of the Fisher matrix in M, eta, chi */
1072 REAL8 amp_p2 = amp*amp;
1073 REAL8 ampSqr_x_dPsiM = amp_p2*dPsidM;
1074 REAL8 ampSqr_x_dPsiEta = amp_p2*dPsidEta;
1075 REAL8 ampSqr_x_dPsidChi = amp_p2*dPsidChi;
1076 REAL8 psd_pm1 = 1./Sh->data->data[k];
1077 REAL8 ampSqr_by_Sh = amp_p2*psd_pm1;
1078
1079 gamma_MM += (ampSqr_x_dPsiM*dPsidM + dAdM*dAdM)*psd_pm1;
1080 gamma_MEta += (ampSqr_x_dPsiM*dPsidEta + dAdM*dAdEta)*psd_pm1;
1081 gamma_MChi += (ampSqr_x_dPsiM*dPsidChi + dAdM*dAdChi)*psd_pm1;
1082 gamma_MT0 += ampSqr_x_dPsiM*dPsidT0*psd_pm1;
1083 gamma_MPhi0 += ampSqr_x_dPsiM*psd_pm1;
1084 gamma_EtaEta += (ampSqr_x_dPsiEta*dPsidEta + dAdEta*dAdEta)*psd_pm1;
1085 gamma_EtaChi += (ampSqr_x_dPsiEta*dPsidChi + dAdEta*dAdChi)*psd_pm1;
1086 gamma_EtaT0 += ampSqr_x_dPsiEta*dPsidT0*psd_pm1;
1087 gamma_EtaPhi0 += ampSqr_x_dPsiEta*psd_pm1;
1088 gamma_ChiChi += (ampSqr_x_dPsidChi*dPsidChi + dAdChi*dAdChi)*psd_pm1;
1089 gamma_ChiT0 += ampSqr_x_dPsidChi*dPsidT0*psd_pm1;
1090 gamma_ChiPhi0 += ampSqr_x_dPsidChi*psd_pm1;
1091 gamma_T0T0 += dPsidT0*dPsidT0*ampSqr_by_Sh;
1092 gamma_T0Phi0 += dPsidT0*ampSqr_by_Sh;
1093 gamma_Phi0Phi0 += ampSqr_by_Sh;
1094
1095 f += df;
1096 }
1097
1098 /* this is actually twice the h-square 2*||h||^2 */
1099 REAL8 hSqr = 2.*gamma_Phi0Phi0;
1100
1101 /* fill the gsl matrix containing the Fisher matrix in (M, eta, chi, t0, phi0) */
1102 gsl_matrix * g = gsl_matrix_calloc (5, 5);
1103
1104 gsl_matrix_set (g, 0,0, gamma_MM/hSqr);
1105 gsl_matrix_set (g, 0,1, gamma_MEta/hSqr);
1106 gsl_matrix_set (g, 0,2, gamma_MChi/hSqr);
1107 gsl_matrix_set (g, 0,3, gamma_MT0/hSqr);
1108 gsl_matrix_set (g, 0,4, gamma_MPhi0/hSqr);
1109
1110 gsl_matrix_set (g, 1,0, gsl_matrix_get(g, 0,1));
1111 gsl_matrix_set (g, 1,1, gamma_EtaEta/hSqr);
1112 gsl_matrix_set (g, 1,2, gamma_EtaChi/hSqr);
1113 gsl_matrix_set (g, 1,3, gamma_EtaT0/hSqr);
1114 gsl_matrix_set (g, 1,4, gamma_EtaPhi0/hSqr);
1115
1116 gsl_matrix_set (g, 2,0, gsl_matrix_get(g, 0,2));
1117 gsl_matrix_set (g, 2,1, gsl_matrix_get(g, 1,2));
1118 gsl_matrix_set (g, 2,2, gamma_ChiChi/hSqr);
1119 gsl_matrix_set (g, 2,3, gamma_ChiT0/hSqr);
1120 gsl_matrix_set (g, 2,4, gamma_ChiPhi0/hSqr);
1121
1122 gsl_matrix_set (g, 3,0, gsl_matrix_get(g, 0,3));
1123 gsl_matrix_set (g, 3,1, gsl_matrix_get(g, 1,3));
1124 gsl_matrix_set (g, 3,2, gsl_matrix_get(g, 2,3));
1125 gsl_matrix_set (g, 3,3, gamma_T0T0/hSqr);
1126 gsl_matrix_set (g, 3,4, gamma_T0Phi0/hSqr);
1127
1128 gsl_matrix_set (g, 4,0, gsl_matrix_get(g, 0,4));
1129 gsl_matrix_set (g, 4,1, gsl_matrix_get(g, 1,4));
1130 gsl_matrix_set (g, 4,2, gsl_matrix_get(g, 2,4));
1131 gsl_matrix_set (g, 4,3, gsl_matrix_get(g, 3,4));
1132 gsl_matrix_set (g, 4,4, gamma_Phi0Phi0/hSqr);
1133
1134 return g;
1135};
1136
1137/*
1138 * Project the Fisher matrix orthogonal to the extrnsic parameters t0 and phi0
1139 */
1140static gsl_matrix *XLALSimIMRPhenomBProjectExtrinsicParam(
1141 gsl_matrix *g /**< Fisher matrix of IMRPhenomB */
1142){
1143 int s = 0;
1144
1145 /* Form submatrices g1, g2, g3, g4, defined as:
1146 * * g = [ g1 g2
1147 * * g4 g3 ] */
1148 gsl_matrix_view g1v = gsl_matrix_submatrix (g, 0, 0, 3, 3);
1149 gsl_matrix_view g2v = gsl_matrix_submatrix (g, 0, 3, 3, 2);
1150 gsl_matrix_view g3v = gsl_matrix_submatrix (g, 3, 3, 2, 2);
1151 gsl_matrix_view g4v = gsl_matrix_submatrix (g, 3, 0, 2, 3);
1152
1153 gsl_matrix * g1 = gsl_matrix_calloc (3, 3);
1154 gsl_matrix * g2 = gsl_matrix_calloc (3, 2);
1155 gsl_matrix * g3 = gsl_matrix_calloc (2, 2);
1156 gsl_matrix * g4 = gsl_matrix_calloc (2, 3);
1157 gsl_matrix * g3invg4 = gsl_matrix_calloc (2, 3);
1158 gsl_matrix * g2g3invg4 = gsl_matrix_calloc (3, 3);
1159
1160 /* Project out the t0 and phi0 dimensions: gamma = g1 - g2 g3^{-1} g4 */
1161 gsl_matrix * g3inv = gsl_matrix_calloc (2, 2);
1162 gsl_permutation * p = gsl_permutation_calloc (2);
1163
1164 gsl_matrix_memcpy (g1, &g1v.matrix);
1165 gsl_matrix_memcpy (g2, &g2v.matrix);
1166 gsl_matrix_memcpy (g3, &g3v.matrix);
1167 gsl_matrix_memcpy (g4, &g4v.matrix);
1168 gsl_matrix_free (g);
1169
1170 gsl_linalg_LU_decomp (g3, p, &s);
1171 gsl_linalg_LU_invert (g3, p, g3inv);
1172 gsl_permutation_free (p);
1173 gsl_matrix_free (g3);
1174
1175 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, g3inv, g4, 0.0, g3invg4);
1176 gsl_matrix_free (g4);
1177 gsl_matrix_free (g3inv);
1178 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, g2, g3invg4, 0.0, g2g3invg4);
1179 gsl_matrix_free (g2);
1180 gsl_matrix_free (g3invg4);
1181
1182 gsl_matrix_sub (g1, g2g3invg4);
1183 gsl_matrix_free (g2g3invg4);
1184 return g1;
1185
1186};
1187
1188/* EXTERNAL ROUTINES */
1189
1190/**
1191 * @addtogroup LALSimIMRPhenom_c
1192 * @brief Routines to produce IMRPhenom-family of phenomenological
1193 * inspiral-merger-ringdown waveforms.
1194 *
1195 * These are frequency-domain models for compact binaries at comparable masses,
1196 * tuned to numerical-relativity simulations.
1197 * * IMRPhenomA models non-spinning binaries.
1198 * * IMRPhenomB/C/D model spinning, but non-precessing binaries.
1199 * * Versions A/B/C should not be used anymore
1200 * unless there are very specific reasons.
1201 * * IMRPhenomP models precessing binaries,
1202 * based on IMRPhenomC
1203 * (outdated, should not be used unless for very specific reasons).
1204 * * IMRPhenomPv2 models precessing binaries,
1205 * based on IMRPhenomD.
1206 * * IMRPhenomHM models spinning, non-precessing binaries,
1207 * based on IMRPhenomD that also includes higher order modes.
1208 * * IMRPhenomPv3 models precessing binaries,
1209 * based on IMRPhenomD but with a different precession prescription
1210 * than IMRPhenomPv2.
1211 * * IMRPhenomPv3HM models precessing binaries with higher-order modes,
1212 * based on IMRPhenomHM.
1213 * * The IMRPhenomX family includes updated models for
1214 * aligned-spin binaries with/without higher modes
1215 * (IMRPhenomXAS and IMRPhenomXHM),
1216 * and for precessing binaries with/without higher modes
1217 * (IMRPhenomXP and IMRPhenomXPHM).
1218 * See @ref LALSimIMRPhenomX_c
1219 * for details.
1220 * * IMRPhenomPv2_NRTidal models precessing binaries,
1221 * adds NR-tuned tidal effects
1222 * to IMRPhenomD phasing and twists the waveform
1223 * to generate the corresponding precessing waveform.
1224 * Two flavors of NRTidal models are available: original (_NRTidal, based on https://arxiv.org/pdf/1706.02969.pdf) and an improved version 2 (_NRTidalv2, based on https://arxiv.org/pdf/1905.06011.pdf).
1225 * * IMRPhenomNSBH models single-spin, non-precessing neutron-star-black-hole
1226 * binaries, based on the amplitude of IMRPhenomC and the phase of
1227 * IMRPhenomD_NRTidalv2
1228 *
1229 * @review IMRPhenomB routines reviewed by Frank Ohme, P. Ajith, Alex Nitz
1230 * and Riccardo Sturani. The review concluded with git hash
1231 * 43ce3b0a8753eb266d75a43ba94b6fb6412121d0 (May 2014).
1232 *
1233 * @review IMRPhenomD routines reviewed by Alex Nielsen, Carl Haster,
1234 * Sebastian Khan, Sascha Husa, Frank Ohme, Mark Hannam, Ofek Brinholtz, Lionel London and
1235 * David Keitel.
1236 * The review concluded with git hash
1237 * db16d17013531cd10451c7d0c6906972ce731866 (Oct/Nov 2015).
1238 *
1239 * @review original IMRPhenomP not reviewed, nor going to be.
1240 * IMRPhenomPv2 reviewed by Capano, Pürrer, Bohe et al. Conludeded
1241 * with git hash 1354291cf6a897995a04cd12dce42b7acaca7b34 (May 2016)
1242 *
1243 * @review IMRPhenomPv2_NRTidal reviewed by Ohme, Haney, Khan, Samajdar,
1244 * Riemenschneider, Setyawati, Hinderer. Concluded with git hash
1245 * f15615215a7e70488d32137a827d63192cbe3ef6 (February 2019).
1246 *
1247 * @review IMRPhneomPv2_NRTidalv2 reviewed by Haney, Ossokine, Pürrer. Concluded
1248 * January 2020, for details and final git hash see review wiki:
1249 * https://git.ligo.org/waveforms/reviews/nrtidal_v2/-/wikis/home
1250 *
1251 * @review IMRPhenomHM review wiki page can be found here
1252 * https://git.ligo.org/waveforms/reviews/phenomhm/wikis/home
1253 *
1254 * @review IMRPhenomNSBH review by Frank Ohme, Tim Dietrich, Shrobana Ghosh,
1255 * Andrew Matas, Jonathan Thompson, Edward Fauchon-Jones. The review concluded
1256 * on 3 February 2020. The review documentation, resources, and final git hash
1257 * can be found at https://git.ligo.org/waveforms/reviews/nsbh-models/wikis/home.
1258 *
1259 * @review IMRPhenomPv3 review:
1260 * https://git.ligo.org/waveforms/reviews/phenompv3hm/-/wikis/home
1261 *
1262 * @review IMRPhenomPv3HM review:
1263 * https://git.ligo.org/waveforms/reviews/phenompv3hm/-/wikis/home
1264 *
1265 * @review IMRPhenomX family review:
1266 * https://git.ligo.org/waveforms/reviews/imrphenomx/-/wikis/home
1267 * @{
1268 */
1269
1270/**
1271 * @name Routines for IMR Phenomenological Model "A"
1272 * @{
1273 */
1274
1275/**
1276 * Driver routine to compute the non-spinning, inspiral-merger-ringdown
1277 * phenomenological waveform IMRPhenomA in the frequency domain.
1278 *
1279 * Reference:
1280 * - Waveform: Eq.(4.13) and (4.16) of http://arxiv.org/pdf/0710.2335
1281 * - Coefficients: Eq.(4.18) of http://arxiv.org/pdf/0710.2335 and
1282 * Table I of http://arxiv.org/pdf/0712.0343
1283 *
1284 * All input parameters should be SI units.
1285 */
1286int XLALSimIMRPhenomAGenerateFD(
1287 COMPLEX16FrequencySeries **htilde, /**< FD waveform */
1288 const REAL8 phi0, /**< orbital phase at peak (rad) */
1289 const REAL8 deltaF, /**< frequency resolution (Hz) */
1290 const REAL8 m1_SI, /**< mass of companion 1 (kg) */
1291 const REAL8 m2_SI, /**< mass of companion 2 (kg) */
1292 const REAL8 f_min, /**< starting GW frequency (Hz) */
1293 const REAL8 f_max, /**< end frequency; if 0, set to fCut */
1294 const REAL8 distance /**< distance of source (m) */
1295) {
1296 BBHPhenomParams *params;
1297 REAL8 f_max_prime;
1298
1299 /* external: SI; internal: solar masses */
1300 const REAL8 m1 = m1_SI / LAL_MSUN_SI;
1301 const REAL8 m2 = m2_SI / LAL_MSUN_SI;
1302
1303 /* check inputs for sanity */
1304 if (*htilde) XLAL_ERROR(XLAL_EFAULT);
1305 if (deltaF <= 0) XLAL_ERROR(XLAL_EDOM);
1306 if (m1 < 0) XLAL_ERROR(XLAL_EDOM);
1307 if (m2 < 0) XLAL_ERROR(XLAL_EDOM);
1308 if (f_min <= 0) XLAL_ERROR(XLAL_EDOM);
1309 if (f_max < 0) XLAL_ERROR(XLAL_EDOM);
1310 if (distance <= 0) XLAL_ERROR(XLAL_EDOM);
1311
1312 /* phenomenological parameters*/
1313 params = ComputeIMRPhenomAParams(m1, m2);
1314 if (!params) XLAL_ERROR(XLAL_EFUNC);
1315 if (params->fCut <= f_min) {
1316 XLALPrintError("fCut <= f_min");
1317 XLAL_ERROR(XLAL_EDOM);
1318 }
1319
1320 /* default f_max to params->fCut */
1321 f_max_prime = f_max ? f_max : params->fCut;
1322 if (f_max_prime <= f_min) {
1323 XLALPrintError("f_max <= f_min");
1324 XLAL_ERROR(XLAL_EDOM);
1325 }
1326
1327 return IMRPhenomAGenerateFD(htilde, phi0, deltaF, m1, m2, f_min, f_max_prime, distance, params);
1328}
1329
1330/**
1331 * Driver routine to compute the non-spinning, inspiral-merger-ringdown
1332 * phenomenological waveform IMRPhenomA in the time domain.
1333 *
1334 * Reference:
1335 * - Waveform: Eq.(4.13) and (4.16) of http://arxiv.org/pdf/0710.2335
1336 * - Coefficients: Eq.(4.18) of http://arxiv.org/pdf/0710.2335 and
1337 * Table I of http://arxiv.org/pdf/0712.0343
1338 *
1339 * All input parameters should be in SI units. Angles should be in radians.
1340 */
1341int XLALSimIMRPhenomAGenerateTD(
1342 REAL8TimeSeries **hplus, /**< +-polarization waveform */
1343 REAL8TimeSeries **hcross, /**< x-polarization waveform */
1344 const REAL8 phiPeak, /**< orbital phase at peak (rad) */
1345 const REAL8 deltaT, /**< sampling interval (s) */
1346 const REAL8 m1_SI, /**< mass of companion 1 (kg) */
1347 const REAL8 m2_SI, /**< mass of companion 2 (kg) */
1348 const REAL8 f_min, /**< starting GW frequency (Hz) */
1349 const REAL8 f_max, /**< end GW frequency; 0 defaults to ringdown cutoff freq */
1350 const REAL8 distance, /**< distance of source (m) */
1351 const REAL8 inclination /**< inclination of source (rad) */
1352) {
1353 BBHPhenomParams *params;
1354 size_t cut_ind, peak_ind;
1355 REAL8 peak_phase; /* measured, not intended */
1356 REAL8 f_max_prime;
1357
1358 /* external: SI; internal: solar masses */
1359 const REAL8 m1 = m1_SI / LAL_MSUN_SI;
1360 const REAL8 m2 = m2_SI / LAL_MSUN_SI;
1361
1362 /* check inputs for sanity */
1363 if (*hplus) XLAL_ERROR(XLAL_EFAULT);
1364 if (*hcross) XLAL_ERROR(XLAL_EFAULT);
1365 if (deltaT <= 0) XLAL_ERROR(XLAL_EDOM);
1366 if (m1 < 0) XLAL_ERROR(XLAL_EDOM);
1367 if (m2 < 0) XLAL_ERROR(XLAL_EDOM);
1368 if (f_min <= 0) XLAL_ERROR(XLAL_EDOM);
1369 if (f_max < 0) XLAL_ERROR(XLAL_EDOM);
1370 if (distance <= 0) XLAL_ERROR(XLAL_EDOM);
1371
1372 /* phenomenological parameters*/
1373 params = ComputeIMRPhenomAParams(m1, m2);
1374 if (!params) XLAL_ERROR(XLAL_EFUNC);
1375 if (params->fCut <= f_min) {
1376 XLALPrintError("fCut <= f_min");
1377 XLAL_ERROR(XLAL_EDOM);
1378 }
1379
1380 /* default f_max to params->fCut */
1381 f_max_prime = f_max ? f_max : params->fCut;
1382 if (f_max_prime <= f_min) {
1383 XLALPrintError("f_max <= f_min");
1384 XLAL_ERROR(XLAL_EDOM);
1385 }
1386
1387 /* generate hplus */
1388 IMRPhenomAGenerateTD(hplus, 0, deltaT, m1, m2, f_min, f_max_prime, distance, params);
1389 if (!(*hplus)) {
1390 XLALFree(params);
1391 XLAL_ERROR(XLAL_EFUNC);
1392 }
1393
1394 /* generate hcross, which is hplus w/ GW phase shifted by -pi/2
1395 * <==> orb. phase shifted by -pi/4 */
1396 IMRPhenomAGenerateTD(hcross, -LAL_PI_4, deltaT, m1, m2, f_min, f_max_prime, distance, params);
1397 XLALFree(params);
1398 if (!(*hcross)) {
1399 XLALDestroyREAL8TimeSeries(*hplus);
1400 *hplus = NULL;
1401 XLAL_ERROR(XLAL_EFUNC);
1402 }
1403
1404 /* clip the parts below f_min */
1405 cut_ind = find_instant_freq(*hplus, *hcross, f_min, (*hplus)->data->length - EstimateIMRLength(m1, m2, f_min, deltaT) + EstimateIMRLength(m1, m2, f_max_prime, deltaT));
1406 *hplus = XLALResizeREAL8TimeSeries(*hplus, cut_ind, (*hplus)->data->length - cut_ind);
1407 *hcross = XLALResizeREAL8TimeSeries(*hcross, cut_ind, (*hcross)->data->length - cut_ind);
1408 if (!(*hplus) || !(*hcross))
1409 XLAL_ERROR(XLAL_EFUNC);
1410
1411 /* set phase and time at peak */
1412 peak_ind = find_peak_amp(*hplus, *hcross);
1413 peak_phase = atan2((*hcross)->data->data[peak_ind], (*hplus)->data->data[peak_ind]);
1414 // NB: factor of 2 b/c phiPeak is *orbital* phase, and we're shifting GW phase
1415 apply_phase_shift(*hplus, *hcross, 2.*phiPeak - peak_phase);
1416 XLALGPSSetREAL8(&((*hplus)->epoch), -(peak_ind * deltaT));
1417 XLALGPSSetREAL8(&((*hcross)->epoch), -(peak_ind * deltaT));
1418
1419 /* apply inclination */
1420 return apply_inclination(*hplus, *hcross, inclination);
1421}
1422
1423/**
1424 * Compute the default final frequency
1425 */
1426double XLALSimIMRPhenomAGetFinalFreq(
1427 const REAL8 m1,
1428 const REAL8 m2
1429) {
1430 BBHPhenomParams *phenomParams;
1431 phenomParams = ComputeIMRPhenomAParams(m1, m2);
1432 return phenomParams->fCut;
1433}
1434
1435/** @} */
1436
1437/**
1438 * @name Routines for IMR Phenomenological Model "B"
1439 * @{
1440 */
1441
1442/**
1443 * Compute the dimensionless, spin-aligned parameter chi as used in the
1444 * IMRPhenomB waveform. This is different from chi in SpinTaylorRedSpin!
1445 * Reference: http://arxiv.org/pdf/0909.2867, paragraph 3.
1446 */
1447double XLALSimIMRPhenomBComputeChi(
1448 const REAL8 m1, /**< mass of companion 1 */
1449 const REAL8 m2, /**< mass of companion 2 */
1450 const REAL8 s1z, /**< spin of companion 1 */
1451 const REAL8 s2z /**< spin of companion 2 */
1452) {
1453 return (m1 * s1z + m2 * s2z) / (m1 + m2);
1454}
1455
1456/**
1457 * Compute the default final frequency
1458 */
1459double XLALSimIMRPhenomBGetFinalFreq(
1460 const REAL8 m1,
1461 const REAL8 m2,
1462 const REAL8 chi
1463) {
1464 BBHPhenomParams *phenomParams;
1465 phenomParams = ComputeIMRPhenomBParams(m1, m2, chi);
1466 return phenomParams->fCut;
1467}
1468
1469/**
1470 * Driver routine to compute the spin-aligned, inspiral-merger-ringdown
1471 * phenomenological waveform IMRPhenomB in the time domain.
1472 *
1473 * Reference: http://arxiv.org/pdf/0909.2867
1474 * - Waveform: Eq.(1)
1475 * - Coefficients: Eq.(2) and Table I
1476 *
1477 * All input parameters should be in SI units. Angles should be in radians.
1478 */
1479int XLALSimIMRPhenomBGenerateTD(
1480 REAL8TimeSeries **hplus, /**< +-polarization waveform */
1481 REAL8TimeSeries **hcross, /**< x-polarization waveform */
1482 const REAL8 phiPeak, /**< orbital phase at peak (rad) */
1483 const REAL8 deltaT, /**< sampling interval (s) */
1484 const REAL8 m1_SI, /**< mass of companion 1 (kg) */
1485 const REAL8 m2_SI, /**< mass of companion 2 (kg) */
1486 const REAL8 chi, /**< mass-weighted aligned-spin parameter */
1487 const REAL8 f_min, /**< starting GW frequency (Hz) */
1488 const REAL8 f_max, /**< end GW frequency; 0 defaults to ringdown cutoff freq */
1489 const REAL8 distance, /**< distance of source (m) */
1490 const REAL8 inclination /**< inclination of source (rad) */
1491) {
1492 BBHPhenomParams *params;
1493 size_t cut_ind, peak_ind;
1494 REAL8 peak_phase; /* measured, not intended */
1495 REAL8 f_max_prime;
1496
1497 /* external: SI; internal: solar masses */
1498 const REAL8 m1 = m1_SI / LAL_MSUN_SI;
1499 const REAL8 m2 = m2_SI / LAL_MSUN_SI;
1500
1501 /* check inputs for sanity */
1502 if (*hplus) XLAL_ERROR(XLAL_EFAULT);
1503 if (*hcross) XLAL_ERROR(XLAL_EFAULT);
1504 if (deltaT <= 0) XLAL_ERROR(XLAL_EDOM);
1505 if (m1 < 0) XLAL_ERROR(XLAL_EDOM);
1506 if (m2 < 0) XLAL_ERROR(XLAL_EDOM);
1507 if (fabs(chi) > 1) XLAL_ERROR(XLAL_EDOM);
1508 if (f_min <= 0) XLAL_ERROR(XLAL_EDOM);
1509 if (f_max < 0) XLAL_ERROR(XLAL_EDOM);
1510 if (distance <= 0) XLAL_ERROR(XLAL_EDOM);
1511
1512 /* phenomenological parameters*/
1513 params = ComputeIMRPhenomBParams(m1, m2, chi);
1514 if (!params) XLAL_ERROR(XLAL_EFUNC);
1515 if (params->fCut <= f_min) {
1516 XLALPrintError("fCut <= f_min");
1517 XLAL_ERROR(XLAL_EDOM);
1518 }
1519
1520 /* default f_max to params->fCut */
1521 f_max_prime = f_max ? f_max : params->fCut;
1522 if (f_max_prime <= f_min) {
1523 XLALPrintError("f_max <= f_min");
1524 XLAL_ERROR(XLAL_EDOM);
1525 }
1526
1527 /* generate plus */
1528 IMRPhenomBGenerateTD(hplus, 0., deltaT, m1, m2, chi, f_min, f_max_prime, distance, params);
1529 if (!(*hplus)) {
1530 XLALFree(params);
1531 XLAL_ERROR(XLAL_EFUNC);
1532 }
1533
1534 /* generate hcross, which is hplus w/ GW phase shifted by -pi/2
1535 * <==> orb. phase shifted by -pi/4 */
1536 IMRPhenomBGenerateTD(hcross, -LAL_PI_4, deltaT, m1, m2, chi, f_min, f_max_prime, distance, params);
1537 XLALFree(params);
1538 if (!(*hcross)) {
1539 XLALDestroyREAL8TimeSeries(*hplus);
1540 *hplus = NULL;
1541 XLAL_ERROR(XLAL_EFUNC);
1542 }
1543
1544 /* clip the parts below f_min */
1545 cut_ind = find_instant_freq(*hplus, *hcross, f_min, (*hplus)->data->length - EstimateIMRLength(m1, m2, f_min, deltaT) + EstimateIMRLength(m1, m2, f_max_prime, deltaT));
1546 *hplus = XLALResizeREAL8TimeSeries(*hplus, cut_ind, (*hplus)->data->length - cut_ind);
1547 *hcross = XLALResizeREAL8TimeSeries(*hcross, cut_ind, (*hcross)->data->length - cut_ind);
1548 if (!(*hplus) || !(*hcross))
1549 XLAL_ERROR(XLAL_EFUNC);
1550
1551 /* set phase and time at peak */
1552 peak_ind = find_peak_amp(*hplus, *hcross);
1553 peak_phase = atan2((*hcross)->data->data[peak_ind], (*hplus)->data->data[peak_ind]);
1554 // NB: factor of 2 b/c phiPeak is *orbital* phase, and we're shifting GW phase
1555 apply_phase_shift(*hplus, *hcross, 2.*phiPeak - peak_phase);
1556 XLALGPSSetREAL8(&((*hplus)->epoch), -(peak_ind * deltaT));
1557 XLALGPSSetREAL8(&((*hcross)->epoch), -(peak_ind * deltaT));
1558
1559 /* apply inclination */
1560 return apply_inclination(*hplus, *hcross, inclination);
1561}
1562
1563
1564/**
1565 * Driver routine to compute the spin-aligned, inspiral-merger-ringdown
1566 * phenomenological waveform IMRPhenomB in the frequency domain.
1567 *
1568 * Reference: http://arxiv.org/pdf/0909.2867
1569 * - Waveform: Eq.(1)
1570 * - Coefficients: Eq.(2) and Table I
1571 *
1572 * All input parameters should be in SI units. Angles should be in radians.
1573 */
1574int XLALSimIMRPhenomBGenerateFD(
1575 COMPLEX16FrequencySeries **htilde, /**< FD waveform */
1576 const REAL8 phi0, /**< orbital phase at peak (rad) */
1577 const REAL8 deltaF, /**< sampling interval (Hz) */
1578 const REAL8 m1_SI, /**< mass of companion 1 (kg) */
1579 const REAL8 m2_SI, /**< mass of companion 2 (kg) */
1580 const REAL8 chi, /**< mass-weighted aligned-spin parameter */
1581 const REAL8 f_min, /**< starting GW frequency (Hz) */
1582 const REAL8 f_max, /**< end frequency; 0 defaults to ringdown cutoff freq */
1583 const REAL8 distance /**< distance of source (m) */
1584) {
1585 BBHPhenomParams *params;
1586 int status;
1587 REAL8 f_max_prime;
1588
1589 /* external: SI; internal: solar masses */
1590 const REAL8 m1 = m1_SI / LAL_MSUN_SI;
1591 const REAL8 m2 = m2_SI / LAL_MSUN_SI;
1592
1593 /* check inputs for sanity */
1594 if (*htilde) XLAL_ERROR(XLAL_EFAULT);
1595 if (deltaF <= 0) XLAL_ERROR(XLAL_EDOM);
1596 if (m1 < 0) XLAL_ERROR(XLAL_EDOM);
1597 if (m2 < 0) XLAL_ERROR(XLAL_EDOM);
1598 if (fabs(chi) > 1) XLAL_ERROR(XLAL_EDOM);
1599 if (f_min <= 0) XLAL_ERROR(XLAL_EDOM);
1600 if (f_max < 0) XLAL_ERROR(XLAL_EDOM);
1601 if (distance <= 0) XLAL_ERROR(XLAL_EDOM);
1602
1603 /* phenomenological parameters*/
1604 params = ComputeIMRPhenomBParams(m1, m2, chi);
1605 if (!params) XLAL_ERROR(XLAL_EFUNC);
1606 if (params->fCut <= f_min) {
1607 XLALPrintError("fCut <= f_min");
1608 XLAL_ERROR(XLAL_EDOM);
1609 }
1610
1611 /* default f_max to params->fCut */
1612 f_max_prime = f_max ? f_max : params->fCut;
1613 if (f_max_prime <= f_min) {
1614 XLALPrintError("f_max <= f_min");
1615 XLAL_ERROR(XLAL_EDOM);
1616 }
1617
1618 status = IMRPhenomBGenerateFD(htilde, phi0, deltaF, m1, m2, chi, f_min, f_max_prime, distance, params);
1619 LALFree(params);
1620 return status;
1621}
1622
1623/**
1624 * Compute the template-space metric of the IMRPhenomB waveform in the
1625 * M, eta, chi coordinates
1626 */
1627int XLALSimIMRPhenomBMetricInMEtaChi(
1628 REAL8 *gamma00, /**< template metric coeff. 00 in PN Chirp Time */
1629 REAL8 *gamma01, /**< template metric coeff. 01/10 PN Chirp Time */
1630 REAL8 *gamma02, /**< template metric coeff. 01/10 PN Chirp Time */
1631 REAL8 *gamma11, /**< template metric coeff. 11 in PN Chirp Time */
1632 REAL8 *gamma12, /**< template metric coeff. 01/10 PN Chirp Time */
1633 REAL8 *gamma22, /**< template metric coeff. 01/10 PN Chirp Time */
1634 const REAL8 m1_SI, /**< component mass 1 (kg) */
1635 const REAL8 m2_SI, /**< component mass 2 (kg) */
1636 const REAL8 chi, /**< effective spin parameter of IMRPhenomB: chi = (m1 chi1 + m2 chi2)/(m1+m2) */
1637 const REAL8 fLow, /**< low-frequency cutoff (Hz) */
1638 const REAL8FrequencySeries *Sh /**< PSD in strain per root Hertz */
1639){
1640
1641 /* external: SI; internal: solar masses */
1642 const REAL8 m1 = m1_SI / LAL_MSUN_SI;
1643 const REAL8 m2 = m2_SI / LAL_MSUN_SI;
1644
1645 /* compute the Fisher matrix in (M, eta, chi, t0, phi0) coords */
1646 gsl_matrix * g = XLALSimIMRPhenomBFisherMatrix(m1, m2, chi, fLow, Sh);
1647
1648 /* Project out t0 and phi0 */
1649 gsl_matrix * g1= XLALSimIMRPhenomBProjectExtrinsicParam(g);
1650
1651 *gamma00 = gsl_matrix_get(g1, 0, 0);
1652 *gamma01 = gsl_matrix_get(g1, 0, 1);
1653 *gamma02 = gsl_matrix_get(g1, 0, 2);
1654 *gamma11 = gsl_matrix_get(g1, 1, 1);
1655 *gamma12 = gsl_matrix_get(g1, 1, 2);
1656 *gamma22 = gsl_matrix_get(g1, 2, 2);
1657
1658 gsl_matrix_free (g1);
1659 return XLAL_SUCCESS;
1660};
1661
1662
1663/**
1664 * Compute the template-space metric of the IMRPhenomB waveform in the
1665 * theta0, theta3, theta3S coordinates
1666 */
1667int XLALSimIMRPhenomBMetricInTheta0Theta3Theta3S(
1668 REAL8 *gamma00, /**< template metric coeff. 00 in PN Chirp Time */
1669 REAL8 *gamma01, /**< template metric coeff. 01/10 PN Chirp Time */
1670 REAL8 *gamma02, /**< template metric coeff. 01/10 PN Chirp Time */
1671 REAL8 *gamma11, /**< template metric coeff. 11 in PN Chirp Time */
1672 REAL8 *gamma12, /**< template metric coeff. 01/10 PN Chirp Time */
1673 REAL8 *gamma22, /**< template metric coeff. 01/10 PN Chirp Time */
1674 const REAL8 m1_SI, /**< component mass 1 (kg) */
1675 const REAL8 m2_SI, /**< component mass 2 (kg) */
1676 const REAL8 chi, /**< effective spin parameter of IMRPhenomB: chi = (m1 chi1 + m2 chi2)/(m1+m2) */
1677 const REAL8 fLow, /**< low-frequency cutoff (Hz) */
1678 const REAL8FrequencySeries *Sh /**< PSD in strain per root Hertz */
1679){
1680
1681 /* external: SI; internal: solar masses */
1682 const REAL8 m1 = m1_SI / LAL_MSUN_SI;
1683 const REAL8 m2 = m2_SI / LAL_MSUN_SI;
1684
1685 const REAL8 M = (m1+m2)*LAL_MTSUN_SI;
1686 const REAL8 eta = m1*m2/((m1+m2)*(m1+m2));
1687 const REAL8 v0 = cbrt(LAL_PI*M*fLow);
1688
1689 /* compute the Fisher matrix in (M, eta, chi, t0, phi0) coords */
1690 gsl_matrix * gMass = XLALSimIMRPhenomBFisherMatrix(m1, m2, chi, fLow, Sh);
1691
1692 const REAL8 theta0 = 5.0/(128.0*eta*v0*v0*v0*v0*v0);
1693 const REAL8 theta3 = (LAL_PI/(4.0*eta*v0*v0));
1694 const REAL8 theta3S = (LAL_PI/(4.0*v0*v0))*(17022. - 9565.9*chi);
1695 const REAL8 theta3_p2 = theta3*theta3;
1696 const REAL8 theta0_p2 = theta0*theta0;
1697 const REAL8 theta0_p2by3 = cbrt(theta3_p2);
1698
1699 /* Co-ordinate Derivatives for Jacobian matrix */
1700 const REAL8 dMdTheta0 = (-0.015831434944115277*theta3)/(fLow*theta0_p2);
1701 const REAL8 dMdTheta3 = 0.015831434944115277/(fLow*theta0);
1702 const REAL8 dMdTheta3S = 0.;
1703 const REAL8 dEtadTheta0 = 3.8715528021485643/(cbrt(theta0)*theta3*theta0_p2by3);
1704 const REAL8 dEtadTheta3 = (-9.678882005371412*cbrt(theta0_p2))/(theta3_p2*theta0_p2by3);
1705 const REAL8 dEtadTheta3S = 0.;
1706 const REAL8 dChidTheta0 = (0.000012000696668619612*theta0_p2by3*theta3S)/(theta0*cbrt(theta0_p2));
1707 const REAL8 dChidTheta3 = (-0.000012000696668619612*theta3S)/(cbrt(theta0_p2)*cbrt(theta3));
1708 const REAL8 dChidTheta3S = (-0.00001800104500292942*theta0_p2by3)/cbrt(theta0_p2);
1709
1710 /* Define the Jacobian transformation matrix which would give the Fisher Matrix in theta co-ordinates */
1711 gsl_matrix * jaco = gsl_matrix_calloc (5, 5);
1712 gsl_matrix * gPre = gsl_matrix_calloc (5, 5);
1713 gsl_matrix * g = gsl_matrix_calloc (5, 5);
1714
1715 gsl_matrix_set (jaco, 0,0 , dMdTheta0);
1716 gsl_matrix_set (jaco, 0,1 , dEtadTheta0);
1717 gsl_matrix_set (jaco, 0,2 , dChidTheta0);
1718 gsl_matrix_set (jaco, 0,3 , 0.);
1719 gsl_matrix_set (jaco, 0,4 , 0.);
1720
1721 gsl_matrix_set (jaco, 1,0 , dMdTheta3);
1722 gsl_matrix_set (jaco, 1,1 , dEtadTheta3);
1723 gsl_matrix_set (jaco, 1,2 , dChidTheta3);
1724 gsl_matrix_set (jaco, 1,3 , 0.);
1725 gsl_matrix_set (jaco, 1,4 , 0.);
1726
1727 gsl_matrix_set (jaco, 2,0 , dMdTheta3S);
1728 gsl_matrix_set (jaco, 2,1 , dEtadTheta3S);
1729 gsl_matrix_set (jaco, 2,2 , dChidTheta3S);
1730 gsl_matrix_set (jaco, 2,3 , 0.);
1731 gsl_matrix_set (jaco, 2,4 , 0.);
1732
1733 gsl_matrix_set (jaco, 3,0 , 0.);
1734 gsl_matrix_set (jaco, 3,1 , 0.);
1735 gsl_matrix_set (jaco, 3,2 , 0.);
1736 gsl_matrix_set (jaco, 3,3 , 1.);
1737 gsl_matrix_set (jaco, 3,4 , 0.);
1738
1739 gsl_matrix_set (jaco, 4,0 , 0.);
1740 gsl_matrix_set (jaco, 4,1 , 0.);
1741 gsl_matrix_set (jaco, 4,2 , 0.);
1742 gsl_matrix_set (jaco, 4,3 , 0.);
1743 gsl_matrix_set (jaco, 4,4 , 1.);
1744
1745 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, gMass, jaco, 0.0, gPre);
1746 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, jaco, gPre, 0.0, g);
1747 gsl_matrix_free (jaco);
1748 gsl_matrix_free (gPre);
1749
1750 /* Project out t0 and phi0 */
1751 gsl_matrix * g1 = XLALSimIMRPhenomBProjectExtrinsicParam(g);
1752
1753 *gamma00 = gsl_matrix_get(g1, 0, 0);
1754 *gamma01 = gsl_matrix_get(g1, 0, 1);
1755 *gamma02 = gsl_matrix_get(g1, 0, 2);
1756 *gamma11 = gsl_matrix_get(g1, 1, 1);
1757 *gamma12 = gsl_matrix_get(g1, 1, 2);
1758 *gamma22 = gsl_matrix_get(g1, 2, 2);
1759
1760 gsl_matrix_free (g1);
1761 return XLAL_SUCCESS;
1762};
1763
1764/** @} */
1765/** @} */
LALFree
#define LALFree(p)
g3
const double g3
Definition: LALSimIMRLackeyTidal2013.c:57
g2
const double g2
Definition: LALSimIMRLackeyTidal2013.c:56
g1
const double g1
Definition: LALSimIMRLackeyTidal2013.c:55
ComputeIMRPhenomBParams
static BBHPhenomParams * ComputeIMRPhenomBParams(const REAL8 m1, const REAL8 m2, const REAL8 chi)
Definition: LALSimIMRPhenom.c:167
ComputeTau0
static REAL8 ComputeTau0(const REAL8 m1, const REAL8 m2, const REAL8 f_min)
Definition: LALSimIMRPhenom.c:248
EstimateSafeFMinForTD
static REAL8 EstimateSafeFMinForTD(const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 deltaT)
Definition: LALSimIMRPhenom.c:271
find_peak_amp
static size_t find_peak_amp(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc)
Definition: LALSimIMRPhenom.c:599
IMRPhenomAGenerateTD
static int IMRPhenomAGenerateTD(REAL8TimeSeries **h, const REAL8 phiPeak, const REAL8 deltaT, const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, const BBHPhenomParams *params)
Definition: LALSimIMRPhenom.c:474
ComputeIMRPhenomAParams
static BBHPhenomParams * ComputeIMRPhenomAParams(const REAL8 m1, const REAL8 m2)
Definition: LALSimIMRPhenom.c:90
apply_phase_shift
static int apply_phase_shift(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 shift)
Definition: LALSimIMRPhenom.c:616
XLALSimIMRPhenomBProjectExtrinsicParam
static gsl_matrix * XLALSimIMRPhenomBProjectExtrinsicParam(gsl_matrix *g)
Definition: LALSimIMRPhenom.c:1140
IMRPhenomBGenerateTD
static int IMRPhenomBGenerateTD(REAL8TimeSeries **h, const REAL8 phiPeak, const REAL8 deltaT, const REAL8 m1, const REAL8 m2, const REAL8 chi, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, const BBHPhenomParams *params)
Definition: LALSimIMRPhenom.c:502
FDToTD
static int FDToTD(REAL8TimeSeries **signalTD, const COMPLEX16FrequencySeries *signalFD, const REAL8 totalMass, const REAL8 deltaT, const REAL8 f_min, const REAL8 f_max, const REAL8 f_min_wide, const REAL8 f_max_wide)
Definition: LALSimIMRPhenom.c:531
XLALSimIMRPhenomBFisherMatrix
static gsl_matrix * XLALSimIMRPhenomBFisherMatrix(const REAL8 m1, const REAL8 m2, const REAL8 chi, const REAL8 fLow, const REAL8FrequencySeries *Sh)
Definition: LALSimIMRPhenom.c:954
find_instant_freq
static size_t find_instant_freq(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 target, const size_t start)
Definition: LALSimIMRPhenom.c:580
apply_inclination
static int apply_inclination(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 inclination)
Definition: LALSimIMRPhenom.c:637
IMRPhenomBGenerateFD
static int IMRPhenomBGenerateFD(COMPLEX16FrequencySeries **htilde, const REAL8 phi0, const REAL8 deltaF, const REAL8 m1, const REAL8 m2, const REAL8 chi, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, const BBHPhenomParams *params)
Definition: LALSimIMRPhenom.c:378
LorentzianFn
static REAL8 LorentzianFn(const REAL8 freq, const REAL8 fRing, const REAL8 sigma)
Definition: LALSimIMRPhenom.c:293
EstimateIMRLength
static size_t EstimateIMRLength(const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 deltaT)
Definition: LALSimIMRPhenom.c:258
EstimateSafeFMaxForTD
static REAL8 EstimateSafeFMaxForTD(const REAL8 f_max, const REAL8 dt)
Definition: LALSimIMRPhenom.c:285
NextPow2
static size_t NextPow2(const size_t n)
Definition: LALSimIMRPhenom.c:262
IMRPhenomAGenerateFD
static int IMRPhenomAGenerateFD(COMPLEX16FrequencySeries **htilde, const REAL8 phi0, const REAL8 deltaF, const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, const BBHPhenomParams *params)
Definition: LALSimIMRPhenom.c:302
XLALComputeIMRPhenomBDerivativeCoeffs
static IMRDerivCoeffs * XLALComputeIMRPhenomBDerivativeCoeffs(const REAL8 m1, const REAL8 m2, const REAL8 chi, BBHPhenomParams *params)
Definition: LALSimIMRPhenom.c:757
sigma
static double sigma(const double a, const double b, const sysq *system)
Internal function that computes the spin-spin couplings.
Definition: LALSimInspiralFDPrecAngles_internals.c:806
g4
REAL8 g4
Definition: SEOBNRv2_opt_3derivs.h:420
s
int s
Definition: bh_qnmode.c:137
M
REAL8 M
Definition: bh_qnmode.c:133
dt
double dt
Definition: bh_ringdown.c:113
i
double i
Definition: bh_ringdown.c:118
w2
const double w2
Definition: exact_derivatives-HrealHeader.c:55
data
sigmaKerr data[0]
Definition: exact_derivatives-HrealHeader.c:4
XLALCreateCOMPLEX16FrequencySeries
COMPLEX16FrequencySeries * XLALCreateCOMPLEX16FrequencySeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaF, const LALUnit *sampleUnits, size_t length)
XLALDestroyCOMPLEX16FrequencySeries
void XLALDestroyCOMPLEX16FrequencySeries(COMPLEX16FrequencySeries *series)
LAL_PI_2
#define LAL_PI_2
LAL_C_SI
#define LAL_C_SI
LAL_MSUN_SI
#define LAL_MSUN_SI
LAL_PI
#define LAL_PI
LAL_TWOPI
#define LAL_TWOPI
LAL_MTSUN_SI
#define LAL_MTSUN_SI
LAL_PI_4
#define LAL_PI_4
COMPLEX16
double complex COMPLEX16
REAL8
double REAL8
XLALMalloc
void * XLALMalloc(size_t n)
XLALFree
void XLALFree(void *p)
XLALSimIMRPhenomAGetFinalFreq
double XLALSimIMRPhenomAGetFinalFreq(const REAL8 m1, const REAL8 m2)
Compute the default final frequency.
Definition: LALSimIMRPhenom.c:1426
XLALSimIMRPhenomBGenerateFD
int XLALSimIMRPhenomBGenerateFD(COMPLEX16FrequencySeries **htilde, const REAL8 phi0, const REAL8 deltaF, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi, const REAL8 f_min, const REAL8 f_max, const REAL8 distance)
Driver routine to compute the spin-aligned, inspiral-merger-ringdown phenomenological waveform IMRPhe...
Definition: LALSimIMRPhenom.c:1574
XLALSimIMRPhenomBMetricInTheta0Theta3Theta3S
int XLALSimIMRPhenomBMetricInTheta0Theta3Theta3S(REAL8 *gamma00, REAL8 *gamma01, REAL8 *gamma02, REAL8 *gamma11, REAL8 *gamma12, REAL8 *gamma22, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi, const REAL8 fLow, const REAL8FrequencySeries *Sh)
Compute the template-space metric of the IMRPhenomB waveform in the theta0, theta3,...
Definition: LALSimIMRPhenom.c:1667
XLALSimIMRPhenomBComputeChi
double XLALSimIMRPhenomBComputeChi(const REAL8 m1, const REAL8 m2, const REAL8 s1z, const REAL8 s2z)
Compute the dimensionless, spin-aligned parameter chi as used in the IMRPhenomB waveform.
Definition: LALSimIMRPhenom.c:1447
XLALSimIMRPhenomAGenerateFD
int XLALSimIMRPhenomAGenerateFD(COMPLEX16FrequencySeries **htilde, const REAL8 phi0, const REAL8 deltaF, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 f_min, const REAL8 f_max, const REAL8 distance)
Driver routine to compute the non-spinning, inspiral-merger-ringdown phenomenological waveform IMRPhe...
Definition: LALSimIMRPhenom.c:1286
XLALSimIMRPhenomAGenerateTD
int XLALSimIMRPhenomAGenerateTD(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, const REAL8 phiPeak, const REAL8 deltaT, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, const REAL8 inclination)
Driver routine to compute the non-spinning, inspiral-merger-ringdown phenomenological waveform IMRPhe...
Definition: LALSimIMRPhenom.c:1341
XLALSimIMRPhenomBGetFinalFreq
double XLALSimIMRPhenomBGetFinalFreq(const REAL8 m1, const REAL8 m2, const REAL8 chi)
Compute the default final frequency.
Definition: LALSimIMRPhenom.c:1459
XLALSimIMRPhenomBGenerateTD
int XLALSimIMRPhenomBGenerateTD(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, const REAL8 phiPeak, const REAL8 deltaT, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi, const REAL8 f_min, const REAL8 f_max, const REAL8 distance, const REAL8 inclination)
Driver routine to compute the spin-aligned, inspiral-merger-ringdown phenomenological waveform IMRPhe...
Definition: LALSimIMRPhenom.c:1479
XLALSimIMRPhenomBMetricInMEtaChi
int XLALSimIMRPhenomBMetricInMEtaChi(REAL8 *gamma00, REAL8 *gamma01, REAL8 *gamma02, REAL8 *gamma11, REAL8 *gamma12, REAL8 *gamma22, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi, const REAL8 fLow, const REAL8FrequencySeries *Sh)
Compute the template-space metric of the IMRPhenomB waveform in the M, eta, chi coordinates.
Definition: LALSimIMRPhenom.c:1627
LAL_REAL8_FORMAT
#define LAL_REAL8_FORMAT
XLALDestroyREAL8FFTPlan
void XLALDestroyREAL8FFTPlan(REAL8FFTPlan *plan)
XLALCreateReverseREAL8FFTPlan
REAL8FFTPlan * XLALCreateReverseREAL8FFTPlan(UINT4 size, int measurelvl)
XLALREAL8FreqTimeFFT
int XLALREAL8FreqTimeFFT(REAL8TimeSeries *tser, const COMPLEX16FrequencySeries *freq, const REAL8FFTPlan *plan)
XLALDestroyREAL8TimeSeries
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
XLALResizeREAL8TimeSeries
REAL8TimeSeries * XLALResizeREAL8TimeSeries(REAL8TimeSeries *series, int first, size_t length)
XLALCreateREAL8TimeSeries
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
lalStrainUnit
const LALUnit lalStrainUnit
lalSecondUnit
const LALUnit lalSecondUnit
XLALUnitMultiply
LALUnit * XLALUnitMultiply(LALUnit *output, const LALUnit *unit1, const LALUnit *unit2)
XLAL_ERROR_NULL
#define XLAL_ERROR_NULL(...)
XLAL_ERROR
#define XLAL_ERROR(...)
XLALPrintError
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLALPrintWarning
int int XLALPrintWarning(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFAULT
XLAL_ERANGE
XLAL_ERANGE
XLAL_EFUNC
XLAL_EFUNC
XLAL_EDOM
XLAL_EDOM
XLALGPSSetREAL8
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
status
string status
p
p
COMPLEX16FrequencySeries::data
COMPLEX16Sequence * data
COMPLEX16Vector::data
COMPLEX16 * data
REAL8FrequencySeries::data
REAL8Sequence * data
REAL8TimeSeries::data
REAL8Sequence * data
REAL8Vector::data
REAL8 * data
params
Definition: burst.c:245
f_min
double f_min
Definition: unicorn.c:22
deltaT
double deltaT
Definition: unicorn.c:24
f_max
double f_max
Definition: unicorn.c:23