LALSimulation  5.4.0.1-fe68b98
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 
38 typedef 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 }
53 BBHPhenomParams;
54 
55 /*
56  * private function prototypes; all internal functions use solar masses.
57  *
58  */
59 
60 static BBHPhenomParams *ComputeIMRPhenomAParams(const REAL8 m1, const REAL8 m2);
61 static BBHPhenomParams *ComputeIMRPhenomBParams(const REAL8 m1, const REAL8 m2, const REAL8 chi);
62 
63 static REAL8 EstimateSafeFMinForTD(const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 deltaT);
64 static REAL8 EstimateSafeFMaxForTD(const REAL8 f_max, const REAL8 dt);
65 static REAL8 ComputeTau0(const REAL8 m1, const REAL8 m2, const REAL8 f_min);
66 static size_t EstimateIMRLength(const REAL8 m1, const REAL8 m2, const REAL8 f_min, const REAL8 deltaT);
67 static size_t NextPow2(const size_t n);
68 
69 static REAL8 LorentzianFn(const REAL8 freq, const REAL8 fRing, const REAL8 sigma);
70 
71 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);
72 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);
73 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);
74 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);
75 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);
76 static size_t find_instant_freq(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 target, const size_t start);
77 static size_t find_peak_amp(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc);
78 static int apply_phase_shift(const REAL8TimeSeries *hp, const REAL8TimeSeries *hc, const REAL8 shift);
79 static 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 /* *******************************************************************/
90 static 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 /*********************************************************************/
167 static 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  */
248 static 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  */
258 static 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 
262 static 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  */
271 static 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  */
285 static 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 
293 static 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  */
302 static 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  */
378 static 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  */
474 static 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  */
502 static 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  */
531 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) {
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 */
580 static 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 */
599 static 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 
616 static 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 */
637 static 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  */
662 typedef 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 }
732 IMRDerivCoeffs;
733 
734 static 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 
742 static 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 
749 static 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  */
757 static 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 
954 static 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  */
1140 static 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  */
1286 int 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  */
1341 int 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  */
1426 double 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  */
1447 double 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  */
1459 double 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  */
1479 int 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  */
1574 int 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  */
1627 int 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  */
1667 int 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)
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)
XLALDestroyREAL8TimeSeries
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
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)
p
list p
status
string status
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