Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInspiral 5.0.3.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALPhenomWaveform.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2008-2009 P. Ajith, Badri Krishnan, Lucia Santamaria, 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
21#include <lal/Units.h>
22#include <lal/LALInspiral.h>
23#include <lal/FindRoot.h>
24#include <lal/SeqFactories.h>
25#include <lal/NRWaveInject.h>
26#include <lal/RealFFT.h>
27#include <lal/TimeFreqFFT.h>
28#include <lal/TimeSeries.h>
29#include <lal/SeqFactories.h>
30#include <lal/VectorOps.h>
31#include <lal/BBHPhenomCoeffs.h>
32#include <lal/AVFactories.h>
33
34#include <math.h>
35
36#ifdef __GNUC__
37#define UNUSED __attribute__ ((unused))
38#else
39#define UNUSED
40#endif
41
42typedef struct tagBBHPhenomParams{
43 REAL8 fMerger;
44 REAL8 fRing;
45 REAL8 fCut;
46 REAL8 sigma;
47 REAL8 psi0;
48 REAL8 psi1;
49 REAL8 psi2;
50 REAL8 psi3;
51 REAL8 psi4;
52 REAL8 psi5;
53 REAL8 psi6;
54 REAL8 psi7;
55 REAL8 psi8;
56}
58
59
60static void XLALComputePhenomParams(
61 BBHPhenomParams *phenParams,
63
64static void XLALComputePhenomParams2(
65 BBHPhenomParams *phenParams,
67
69 REAL8 freq,
70 REAL8 fRing,
71 REAL8 sigma);
72
73static void XLALBBHPhenWaveFD (
75 InspiralTemplate *insp_template,
76 REAL4Vector *signalvec);
77
78static void XLALBBHPhenWaveFD2 (
80 InspiralTemplate *insp_template,
81 REAL4Vector *signalvec);
82
83static void XLALComputeInstantFreq(
84 REAL4Vector *Freq,
85 REAL4Vector *hp,
86 REAL4Vector *hc,
87 REAL8 dt);
88
90 REAL4Vector *h,
91 REAL4Vector *freq,
92 REAL8 cutFreq,
93 REAL8 deltaT);
94
95/*
96 *
97 * Wrapper functions to be called from LALInspiralWave
98 *
99 */
100
101
102/* A = non-spinning binaries. Ref. http://arxiv.org/pdf/0710.2335 */
104 REAL4Vector *signalvec, /**< output array */
105 InspiralTemplate *params /**< inspiral parameters */
106 ) {
107 BBHPhenomParams phenParams;
108
109 /* check inputs */
110 if (!signalvec || !(signalvec->data) || !params) {
111 XLALPrintError(LALINSPIRALH_MSGENULL);
113 }
114 if ((signalvec->length < 2)
115 || (params->mass1 <= 0) || (params->mass2 <= 0)
116 || (params->spin1[0] != 0) || (params->spin1[1] != 0) || (params->spin1[2] != 0)
117 || (params->spin2[0] != 0) || (params->spin2[1] != 0) || (params->spin2[2] != 0)) {
118 XLALPrintError(LALINSPIRALH_MSGECHOICE);
120 }
121
122 /* compute the phenomenological parameters */
123 XLALComputePhenomParams(&phenParams, params);
124
125 /* generate the phenomenological waveform in frequency domain */
126 XLALBBHPhenWaveFD (&phenParams, params, signalvec);
127 return XLAL_SUCCESS;
128}
129
130/* B = aligned-spin binaries. Ref. http://arxiv.org/pdf/0909.2867 */
132 REAL4Vector *signalvec, /**< output array */
133 InspiralTemplate *params /**< inspiral parameters */
134 ) {
135 BBHPhenomParams phenParams;
136
137 /* check inputs */
138 if (!signalvec || !(signalvec->data) || !params) {
139 XLALPrintError(LALINSPIRALH_MSGENULL);
141 }
142 if ((signalvec->length < 2)
143 || (params->mass1 <= 0) || (params->mass2 <= 0)
144 || (params->spin1[0] != 0) || (params->spin1[1] != 0)
145 || (params->spin2[0] != 0) || (params->spin2[1] != 0)) {
146 XLALPrintError(LALINSPIRALH_MSGECHOICE);
148 }
149
150 /* compute the phenomenological parameters */
151 XLALComputePhenomParams2(&phenParams, params);
152
153 /* generate the phenomenological waveform in frequency domain */
154 XLALBBHPhenWaveFD2 (&phenParams, params, signalvec);
155
156 return XLAL_SUCCESS;
157}
158
159/* A = non-spinning binaries. Ref. http://arxiv.org/pdf/0710.2335 */
161 REAL4Vector *signalvec1,
162 REAL4Vector *signalvec2,
164 if (!signalvec1 || !signalvec2 ||
165 !(signalvec1->data) || !(signalvec2->data)) {
166 XLALPrintError(LALINSPIRALH_MSGENULL);
168 }
169
170 /* Initially the waveforms are empty */
171 memset(signalvec1->data, 0, signalvec1->length * sizeof(REAL4));
172 memset(signalvec2->data, 0, signalvec2->length * sizeof(REAL4));
173
174 /* generate one waveform with startPhase specified by the user */
175 if (XLALBBHPhenWaveAFreqDom(signalvec1, params))
177
178 /* generate another waveform orthogonal to it */
179 params->startPhase += LAL_PI_2;
180 if (XLALBBHPhenWaveAFreqDom(signalvec2, params))
182
183 return 0;
184}
185
186/* B = aligned-spin binaries. Ref. http://arxiv.org/pdf/0909.2867 */
188 REAL4Vector *signalvec1,
189 REAL4Vector *signalvec2,
191 if (!signalvec1 || !signalvec2
192 || !(signalvec1->data) || !(signalvec2->data)
193 || !params) {
194 XLALPrintError(LALINSPIRALH_MSGENULL);
196 }
197
198 /* Initially the waveforms are empty */
199 memset(signalvec1->data, 0, signalvec1->length * sizeof(REAL4));
200 memset(signalvec2->data, 0, signalvec2->length * sizeof(REAL4));
201
202 /* generate one waveform with startPhase specified by the user */
203 if (XLALBBHPhenWaveBFreqDom(signalvec1, params))
205
206 /* generate another waveform orthogonal to it */
207 params->startPhase += LAL_PI_2;
208 if (XLALBBHPhenWaveBFreqDom(signalvec2, params))
210
211 return 0;
212}
213
215 REAL4Vector *signalvec1,
216 InspiralTemplate *insp_template) {
217 UINT4 count;
218
219 /* check inputs */
220 if (!signalvec1 || !(signalvec1->data)) {
221 XLALPrintError(LALINSPIRALH_MSGENULL);
223 }
224 if (signalvec1->length <= 2) {
225 XLALPrintError(LALINSPIRALH_MSGECHOICE);
227 }
228 /* Initially the waveforms are empty */
229 memset(signalvec1->data, 0, signalvec1->length * sizeof(REAL4));
230
231 if (XLALBBHPhenTimeDomEngine(signalvec1, NULL, NULL, NULL, NULL, NULL, &count, insp_template))
233
234 return 0;
235}
236
238 REAL4Vector *signalvec1,
239 REAL4Vector *signalvec2,
240 InspiralTemplate *insp_template) {
241 UINT4 count;
242
243 /* check inputs */
244 if (!signalvec1 || !(signalvec1->data)
245 || !signalvec2 || !(signalvec2->data)
246 || !insp_template) {
247 XLALPrintError(LALINSPIRALH_MSGENULL);
249 }
250 if ((signalvec1->length <= 2) || (signalvec2->length <= 2)) {
251 XLALPrintError(LALINSPIRALH_MSGECHOICE);
253 }
254
255 /* Initially the waveforms are empty */
256 memset(signalvec1->data, 0, signalvec1->length * sizeof(REAL4));
257 memset(signalvec2->data, 0, signalvec2->length * sizeof(REAL4));
258
259 return XLALBBHPhenTimeDomEngine(signalvec1, signalvec2, NULL, NULL, NULL, NULL, &count, insp_template);
260}
261
263 REAL4Vector *signalvec1, /**< optional output waveform with phi_c = 0 */
264 REAL4Vector *signalvec2, /**< optional output waveform with phi_c = pi/2 */
265 REAL4Vector *h, /**< optional output waveforms, alternating h+ and hx components */
266 REAL4Vector *aVec, /**< optional output inst. amplitude, alternating A+ and Ax components, assuming that h+ & hx have equal ...*/
267 REAL4Vector *freqVec, /**< optional output instant. freq */
268 REAL8Vector *phiVec, /**< optional output phase evolution */
269 UINT4 *countback, /**< output number of non-zero samples */
270 InspiralTemplate *params) /**< UNDOCUMENTED */ {
271 REAL8 dt, cosI, fLower, peakAmp, fCut, fRes, f, totalMass, softWin, z1, z2;
272 REAL8 fLowerOrig, eta, tau0, winFLo, sigLo, sigHi, tF0, expectedAmplRatio;
273 REAL8 phaseShift, sig1, sig2, startPhaseOrig, phiC;
274 REAL4 windowLength;
275 UINT4 i, j, k, l, n, peakAmpIdx, sigLength, iLower;
276 REAL4Vector *signalFD1 = NULL, *signalFD2 = NULL, *signalTD1 = NULL, *signalTD2 = NULL, *a=NULL, *fVec=NULL;
277 REAL8Vector *phi=NULL;
278 REAL4FFTPlan *revPlan = NULL;
279 BBHPhenomParams phenParams;
280
281 /* check inputs */
282 if (!params) {
283 XLALPrintError(LALINSPIRALH_MSGENULL);
285 }
286 if ((params->nStartPad < 0) || (params->nEndPad < 0)
287 || (params->fLower <= 0) || (params->tSampling <= 0)
288 || (params->mass1 < 0) || (params->mass2 < 0)
289 || (params->spin1[0] != 0) || (params->spin1[1] != 0)
290 || (params->spin2[0] != 0) || (params->spin2[1] != 0)) {
291 XLALPrintError(LALINSPIRALH_MSGECHOICE);
293 }
294
295 dt = 1./params->tSampling;
296 cosI = cos( params->inclination );
297
298 /* compute the phenomenological parameters */
299 switch (params->approximant) {
300 /* non-spinning binaries. Ref. http://arxiv.org/pdf/0710.2335 */
301 case IMRPhenomA:
302 if ((params->spin1[2] != 0) || (params->spin2[2] != 0)) {
303 XLALPrintError(LALINSPIRALH_MSGECHOICE);
305 }
306 XLALComputePhenomParams(&phenParams, params);
307 break;
308 /* aligned-spin binaries. Ref. http://arxiv.org/abs/0909.2867 */
309 case IMRPhenomB:
310 XLALComputePhenomParams2(&phenParams, params);
311 break;
312 default:
313 XLALPrintError(LALINSPIRALH_MSGESWITCH);
315 }
316
317 totalMass = params->mass1 + params->mass2;
318 eta = params->mass1*params->mass2/pow(totalMass,2.);
319
320 /* we will generate the waveform from a frequency which is lower than the
321 * fLower chosen. Also the cutoff frequency is higher than the fCut. We
322 * will later apply a window function, and truncate the time-domain waveform
323 * below an instantaneous frequency fLower */
324 fLowerOrig = params->fLower; /* this is the low-freq set by the user */
325 startPhaseOrig = params->startPhase; /* this is the coalescence phase set by the user*/
326
327 /* Find an optimum value for fLower (using the definition of Newtonian chirp time)
328 * such that the waveform has a minimum length of tau0. This is necessary to avoid
329 * FFT artifacts */
330 tau0 = 64.;
331 fLower = pow((tau0*256.*eta*pow(totalMass*LAL_MTSUN_SI,5./3.)/5.),-3./8.)/LAL_PI;
332 fCut = (1.025)*phenParams.fCut;
333
334 /* make sure that these frequencies are not too out of range */
335 if (fLower > fLowerOrig) fLower = fLowerOrig;
336 if (fLower < 0.5) fLower = 0.5;
337 if (fCut > params->tSampling/2.-100.) fCut = params->tSampling/2.-100.;
338
339 /* generate waveforms over this frequency range */
340 params->fLower = fLower;
341 phenParams.fCut = params->tSampling/2.;
342
343 /* allocate memory for the freq and domain tempaltes - templates are generated
344 * using a lower value of fLower than prescribed by the user */
345 tau0 = 5.*totalMass*LAL_MTSUN_SI/(256.*eta*pow(LAL_PI*totalMass*LAL_MTSUN_SI*fLower, 8./3.));
346 n = pow(2., ceil(log2(tau0*params->tSampling)));
347 signalFD1 = XLALCreateREAL4Vector(n);
348 signalFD2 = XLALCreateREAL4Vector(n);
349 signalTD1 = XLALCreateREAL4Vector(n);
350 signalTD2 = XLALCreateREAL4Vector(n);
351 if (!signalFD1 || !signalFD2 || !signalTD1 || !signalTD2) {
352 if (signalFD1) XLALDestroyREAL4Vector(signalFD1);
353 if (signalFD2) XLALDestroyREAL4Vector(signalFD2);
354 if (signalTD1) XLALDestroyREAL4Vector(signalTD1);
355 if (signalTD2) XLALDestroyREAL4Vector(signalTD2);
357 }
358
359 /* generate the phenomenological waveform in frequency domain */
360 switch (params->approximant) {
361 /* non-spinning binaries. Ref. http://arxiv.org/pdf/0710.2335 */
362 case IMRPhenomA:
363 XLALBBHPhenWaveFD (&phenParams, params, signalFD1);
364 params->startPhase += LAL_PI_2;
365 XLALBBHPhenWaveFD (&phenParams, params, signalFD2);
366 break;
367 /* aligned-spin binaries. Ref. http://arxiv.org/abs/0909.2867 */
368 case IMRPhenomB:
369 XLALBBHPhenWaveFD2 (&phenParams, params, signalFD1);
370 params->startPhase += LAL_PI_2;
371 XLALBBHPhenWaveFD2 (&phenParams, params, signalFD2);
372 break;
373 default:
374 XLALDestroyREAL4Vector(signalFD1);
375 XLALDestroyREAL4Vector(signalFD2);
376 XLALDestroyREAL4Vector(signalTD1);
377 XLALDestroyREAL4Vector(signalTD2);
378 XLALPrintError(LALINSPIRALH_MSGESWITCH);
380 }
381
382 /* apply the softening window function */
383 fRes = params->tSampling/n;
384
385 winFLo = (fLowerOrig + fLower)/2.;
386 // UNUSED!!: REAL8 winFHi = (fCut + phenParams.fCut)/2.;
387 sigLo = 4.;
388 sigHi = 4.;
389
390 signalFD1->data[0] = 0.;
391 signalFD2->data[0] = 0.;
392 for (k = 1; k <= n/2; k++) {
393 f = k*fRes;
394 softWin = (1+tanh((4.*(f-winFLo)/sigLo)))*(1-tanh(4.*(f-fCut)/sigHi))/4.;
395 signalFD1->data[k] *= softWin;
396 signalFD1->data[n-k] *= softWin;
397 signalFD2->data[k] *= softWin;
398 signalFD2->data[n-k] *= softWin;
399 }
400
401 /* Inverse Fourier transform */
402 revPlan = XLALCreateReverseREAL4FFTPlan(n, 0);
403 if (!revPlan) {
404 XLALDestroyREAL4Vector(signalFD1);
405 XLALDestroyREAL4Vector(signalFD2);
406 XLALDestroyREAL4Vector(signalTD1);
407 XLALDestroyREAL4Vector(signalTD2);
408 XLALPrintError(LALINSPIRALH_MSGEMEM);
410 }
411 XLALREAL4VectorFFT(signalTD1, signalFD1, revPlan);
412 XLALREAL4VectorFFT(signalTD2, signalFD2, revPlan);
413 XLALDestroyREAL4Vector(signalFD1);
414 XLALDestroyREAL4Vector(signalFD2);
416
417 /* FFT normalisation. The LAL implementation of the FFT omits the factor 1/n.*/
418 for (i = 0; i < n; i++) {
419 signalTD1->data[i] *= params->tSampling/n;
420 signalTD2->data[i] *= params->tSampling/n;
421 }
422
423 /* apply a linearly decresing window at the end
424 * of the waveform in order to avoid edge effects. */
425 windowLength = 10.*totalMass * LAL_MTSUN_SI*params->tSampling;
426 for (i=0; i< windowLength; i++){
427 signalTD1->data[n-i-1] *= i/windowLength;
428 signalTD2->data[n-i-1] *= i/windowLength;
429 }
430
431 /* compute the instantaneous frequency */
433 fVec = XLALCreateREAL4Vector(n);
435
436 XLALComputeInstantFreq(fVec, signalTD1, signalTD2, dt);
437 peakAmp = 0.;
438 peakAmpIdx = 0;
439
440 /* compute the amplitude and phase. find the peak amplitude */
441 for (i=0; i<n; i++){
442 a->data[i] = sqrt(pow(signalTD1->data[i],2.) + pow(signalTD2->data[i],2.));
443 phi->data[i] = -atan2(signalTD2->data[i], signalTD1->data[i]);
444
445 /* find the peak amplitude*/
446 if (a->data[i] > peakAmp) {
447 peakAmp = a->data[i];
448 peakAmpIdx = i;
449 }
450 }
451
452 /* If the instantaneous amplitude a(t) is too low, noise will corrupt the
453 * estimation of the instantaneous frequency f(t). Hence, if the ratio of
454 * a(t) with the peak amplitude is lower than a threshold, we set f(t) = 0
455 * for those time bins. Choosing of such a threshold is a tricky business!
456 * The current value is chosen using the following argument: a(t) is
457 * proportional to the square of the speed v(t) of the binary (the PN
458 * parameter). At the merger, v(t) ~ 1. At the start of the waveform
459 * v(t) = (pi M f_lower)^{1/3}. Thus, the expected ratio at the start of
460 * the waveform is a_0/a_peak = (pi M f_lower)^{-2/3}. Divide by 2 as a
461 * safety margin */
462 expectedAmplRatio = pow(LAL_PI*LAL_MTSUN_SI*totalMass*fLowerOrig, 2./3.)/2.;
463
464 /* Also find the index of the fVec corresponding to fLowerOrig */
465 iLower = 0;
466 for (i=0; i< fVec->length; i++) {
467 if (a->data[i] < expectedAmplRatio*peakAmp) {
468 fVec->data[i] = 0.0;
469 }
470 if ((iLower == 0) && (fVec->data[i] >= fLowerOrig)) iLower = i;
471 }
472
473 /* unwrap the phase */
475
476 /* reassign the original values of fLower, fCut and startPhase */
477 params->fLower = fLowerOrig;
478 params->fFinal = fCut;
479 params->startPhase = startPhaseOrig;
480
481 /* apply a phase-shift to the waveforms such that the phase at coalescence is
482 * equal to params->startPhase */
483 phaseShift = phi->data[peakAmpIdx] + params->startPhase;
484 phiC = phi->data[peakAmpIdx] - params->startPhase;
485 for (i=iLower; i < fVec->length; i++) {
486 sig1 = signalTD1->data[i]*cos(phaseShift) - signalTD2->data[i]*sin(phaseShift);
487 sig2 = signalTD1->data[i]*sin(phaseShift) + signalTD2->data[i]*cos(phaseShift);
488 signalTD1->data[i] = sig1;
489 signalTD2->data[i] = sig2;
490 phi->data[i] -= phiC;
491 }
492
493 /* copy the frequency components above fLower to the signal vector to be returned */
494 j = params->nStartPad;
495 tF0 = iLower*dt; /* time bin corresponding to f(t) = fLower */
496
497 sigLength = 0;
498 if (signalvec1) sigLength = signalvec1->length;
499 else if (signalvec2) sigLength = signalvec2->length;
500 else if (freqVec) sigLength = freqVec->length;
501 else if (phiVec) sigLength = phiVec->length;
502 else if (aVec) sigLength = aVec->length / 2;
503 else if (h) sigLength = h->length / 2;
504
505 /* inclination-weights on two polarizations */
506 z1 = -0.5*(1. + cosI*cosI);
507 z2 = -cosI;
508
509 for (i=iLower; i < fVec->length; i++) {
510 if (j<sigLength) {
511 k = 2*j;
512 l = k+1;
513
514 if (signalvec1) signalvec1->data[j] = signalTD1->data[i]; /* signal with phi_c = 0 */
515 if (signalvec2) signalvec2->data[j] = signalTD2->data[i]; /* signal with phic_c = pi/2 */
516 if (freqVec) freqVec->data[j] = fVec->data[i]; /* instant. freq */
517 if (phiVec) phiVec->data[j] = phi->data[i]; /* phase evolution */
518 if (aVec) {
519 aVec->data[k] = a->data[i]; /* inst. amplitude, assuming that h+ & hx have equal ...*/
520 aVec->data[l] = a->data[i]; /* ... amplitude. valid in the absence of precession */
521 }
522 if (h) {
523 h->data[k] = z1 * signalTD1->data[i]; /* h+ ~ -[ 1+cos^2(iota) ]/2 cos (phi) */
524 h->data[l] = z2 * signalTD2->data[i]; /* hx ~ -cos(iota) sin (phi) */
525 }
526 j++;
527 }
528 }
529 *countback = j;
530
531 /* free the memory */
535 XLALDestroyREAL4Vector(signalTD1);
536 XLALDestroyREAL4Vector(signalTD2);
537
538 /* store some parameters for record keeping */
539 params->vFinal = pow(LAL_PI*LAL_MTSUN_SI*totalMass*params->fFinal, 1./3.);
540 params->tC = peakAmpIdx*dt-tF0; /* time of coalescence. defined as the time
541 corresponding to the peak amplitude*/
542
543 return 0;
544}
545
547 CoherentGW *waveform, /**< allocated, but completely zeroed CoherentGW structure; this function allocates sub-structures that you must free */
548 InspiralTemplate *params, /**< UNDOCUMENTED */
549 PPNParamStruc *ppnParams) /**< UNDOCUMENTED */ {
550 REAL4Vector *a=NULL; /* amplitude data */
551 REAL4Vector *h=NULL; /* polarization data */
552 REAL4Vector *ff=NULL; /* frequency data */
553 REAL8Vector *phi=NULL; /* phase data */
554
555 UINT4 count, i;
556 REAL8 s, phiC; /* phase at coalescence */
557 CHAR message[256];
558 LIGOTimeGPS zero_time = {0, 0};
559 InspiralInit paramsInit;
560 int returnval = XLAL_SUCCESS;
561
562 /* check inputs */
563 if (!params || !waveform) {
564 XLALPrintError(LALINSPIRALH_MSGENULL);
566 }
567 if ((params->nStartPad < 0) || (params->nEndPad < 0)
568 || (params->fLower <= 0) || (params->tSampling <= 0)) {
569 XLALPrintError(LALINSPIRALH_MSGECHOICE);
571 }
572
573 /* Make sure waveform fields don't exist */
574 if (waveform->a || waveform->h || waveform->f
575 || waveform->phi || waveform->shift) {
576 XLALPrintError(LALINSPIRALH_MSGENULL);
578 }
579
580 if (params->ampOrder) {
581 params->ampOrder = (LALPNOrder) 0;
582 snprintf(message, 256, "WARNING: Amp Order has been reset to %d", params->ampOrder);
583 XLALPrintInfo("%s", message);
584 }
585
586 /* Compute some parameters */
587 if (XLALInspiralInit(params, &paramsInit))
589 if (paramsInit.nbins == 0) return XLAL_SUCCESS;
590
591 /* allocate temporary structures */
592 h = XLALCreateREAL4Vector(2 * paramsInit.nbins);
593 a = XLALCreateREAL4Vector(2 * paramsInit.nbins);
594 ff = XLALCreateREAL4Vector(paramsInit.nbins);
595 phi = XLALCreateREAL8Vector(paramsInit.nbins);
596 memset(h->data, 0, 2 * paramsInit.nbins * sizeof(REAL4));
597 memset(a->data, 0, 2 * paramsInit.nbins * sizeof(REAL4));
598 memset(ff->data, 0, paramsInit.nbins * sizeof(REAL4));
599 memset(phi->data, 0, paramsInit.nbins * sizeof(REAL8));
600 if (!h || !a || !ff || !phi) {
601 XLALError(__func__, __FILE__, __LINE__, XLAL_EFAULT);
602 returnval = XLAL_FAILURE;
603 goto done;
604 }
605
606 /* generate two orthogonal waveforms */
607 params->startPhase = ppnParams->phi;
608 if (XLALBBHPhenTimeDomEngine(NULL, NULL, h, a, ff, phi, &count, params)) {
609 XLALError(__func__, __FILE__, __LINE__, XLAL_EFUNC);
610 returnval = XLAL_FAILURE;
611 goto done;
612 }
613
614 /* Check an empty waveform hasn't been returned */
615 for (i = 0; i < phi->length; i++) {
616 if (phi->data[i] != 0.0) break;
617 }
618 if (i == phi->length) {
619 XLALError(__func__, __FILE__, __LINE__, XLAL_ERANGE);
620 returnval = XLAL_FAILURE;
621 goto done;
622 }
623
624 /* print some messages */
625 snprintf(message, 256, "fFinal = %f", params->fFinal);
626 XLALPrintInfo("%s", message);
627 s = 0.5 * (phi->data[count-1] - phi->data[0]);
628 snprintf(message, 256, "cycles = %f", s/3.14159);
629 XLALPrintInfo("%s", message);
630 if (s < LAL_TWOPI){
631 snprintf(message, 256, "The waveform has only %f cycles; we don't "
632 "keep waveform with less than 2 cycles.", (double) s/ (double)LAL_PI );
633 XLALPrintWarning("%s", message);
634 goto done;
635 }
636
637 /* compute phase */
638 phiC = phi->data[count-1] ;
639 for (i=0; i<count;i++) {
640 phi->data[i] = -phiC + phi->data[i] + ppnParams->phi;
641 }
642
643 /* Allocate the CoherentGW structures. */
646 if (!(waveform->h) || !(waveform->a)) {
647 XLALError(__func__, __FILE__, __LINE__, XLAL_ENOMEM);
648 returnval = XLAL_FAILURE;
649 goto done;
650 }
651 memset(waveform->h, 0, sizeof(REAL4TimeVectorSeries));
652 memset(waveform->a, 0, sizeof(REAL4TimeVectorSeries));
653 waveform->h->data = XLALCreateREAL4VectorSequence(count, 2);
654 waveform->a->data = XLALCreateREAL4VectorSequence(count, 2);
655 waveform->f = XLALCreateREAL4TimeSeries("Phenom inspiral frequency", &zero_time, 0, 1. / params->tSampling, &lalHertzUnit, count);
656 waveform->phi = XLALCreateREAL8TimeSeries("Phenom inspiral phase", &zero_time, 0, 1 / params->tSampling, &lalDimensionlessUnit, count);
657 if (!(waveform->h->data) || !(waveform->a->data) || !(waveform->f) || !(waveform->phi)) {
658 XLALError(__func__, __FILE__, __LINE__, XLAL_ENOMEM);
659 returnval = XLAL_FAILURE;
660 goto done;
661 }
662
663 /* copy the frequency, amplitude and phase data to the waveform structure */
664 /*
665 * Truncation (count < paramsInit.nbins in general) is necessary
666 * because tC is defined at the end of the vector in some codes and not
667 * others.
668 */
669 memcpy(waveform->h->data->data, h->data, 2*count*(sizeof(REAL4)));
670 memcpy(waveform->a->data->data, a->data, 2*count*(sizeof(REAL4)));
671 memcpy(waveform->f->data->data, ff->data, count*(sizeof(REAL4)));
672 memcpy(waveform->phi->data->data, phi->data, count*(sizeof(REAL8)));
673
674 /* also set other parameters in the waveform structure */
675 waveform->h->sampleUnits = waveform->a->sampleUnits = lalStrainUnit;
676 waveform->h->deltaT = waveform->a->deltaT = 1./params->tSampling;
677 waveform->position = ppnParams->position;
678 waveform->psi = ppnParams->psi;
679 snprintf(waveform->h->name, LALNameLength, "Phenom inspiral polarizations");
680 snprintf(waveform->a->name, LALNameLength, "Phenom inspiral amplitudes");
681 /* epoch is set elsewhere */
682
683 /* fill some output */
684 ppnParams->tc = (double)(count-1) / params->tSampling;
685 ppnParams->length = count;
686 ppnParams->dfdt = ((REAL4)(waveform->f->data->data[count-1]
687 - waveform->f->data->data[count-2]))* ppnParams->deltaT;
688 ppnParams->fStop = params->fFinal;
690 ppnParams->termDescription = GENERATEPPNINSPIRALH_MSGEFSTOP;
691 ppnParams->fStart = ppnParams->fStartIn;
692
693done:
694
695 if (ff) XLALDestroyREAL4Vector(ff);
697 if (phi) XLALDestroyREAL8Vector(phi);
698 if (h) XLALDestroyREAL4Vector(h);
699 if (returnval != XLAL_SUCCESS) {
700 if (waveform->h && waveform->h->data) XLALDestroyREAL4VectorSequence(waveform->h->data);
701 if (waveform->h) XLALFree(waveform->h);
702 if (waveform->a && waveform->a->data) XLALDestroyREAL4VectorSequence(waveform->a->data);
703 if (waveform->a) XLALFree(waveform->a);
704 if (waveform->f) XLALDestroyREAL4TimeSeries(waveform->f);
705 if (waveform->phi) XLALDestroyREAL8TimeSeries(waveform->phi);
706 }
707 return returnval;
708}
709
710
711/*
712 *
713 * Core numerical routines; the Phenom prescription is in FD
714 *
715 */
716
717/*********************************************************************/
718/* Compute phenomenological parameters for non-spinning binaries */
719/* Ref. Eq.(4.18) of http://arxiv.org/pdf/0710.2335 and */
720/* Table I of http://arxiv.org/pdf/0712.0343 */
721/*********************************************************************/
723 BBHPhenomParams *phenParams,
725 REAL8 totalMass, piM, eta, fMerg_a, fMerg_b, fMerg_c, fRing_a, fRing_b, etap2;
726 REAL8 fRing_c, sigma_a, sigma_b, sigma_c, fCut_a, fCut_b, fCut_c;
727 REAL8 psi0_a, psi0_b, psi0_c, psi2_a, psi2_b, psi2_c, psi3_a, psi3_b, psi3_c;
728 REAL8 psi4_a, psi4_b, psi4_c, psi6_a, psi6_b, psi6_c, psi7_a, psi7_b, psi7_c;
729
730 /* calculate the total mass and symmetric mass ratio */
731
732 if (params) {
733
734 totalMass = params->mass1+params->mass2;
735 eta = params->mass1*params->mass2/pow(totalMass,2.);
736 piM = totalMass*LAL_PI*LAL_MTSUN_SI;
737 }
738 else {
739 return;
740 }
741
742 fMerg_a = BBHPHENOMCOEFFSH_FMERG_A;
743 fMerg_b = BBHPHENOMCOEFFSH_FMERG_B;
744 fMerg_c = BBHPHENOMCOEFFSH_FMERG_C;
745
746 fRing_a = BBHPHENOMCOEFFSH_FRING_A;
747 fRing_b = BBHPHENOMCOEFFSH_FRING_B;
748 fRing_c = BBHPHENOMCOEFFSH_FRING_C;
749
750 sigma_a = BBHPHENOMCOEFFSH_SIGMA_A;
751 sigma_b = BBHPHENOMCOEFFSH_SIGMA_B;
752 sigma_c = BBHPHENOMCOEFFSH_SIGMA_C;
753
757
761
765
769
773
777
781
782 /* Evaluate the polynomials. See Eq. (4.18) of P. Ajith et al
783 * arXiv:0710.2335 [gr-qc] */
784 if (phenParams) {
785 etap2 = eta*eta;
786 phenParams->fCut = (fCut_a*etap2 + fCut_b*eta + fCut_c)/piM;
787 phenParams->fMerger = (fMerg_a*etap2 + fMerg_b*eta + fMerg_c)/piM;
788 phenParams->fRing = (fRing_a*etap2 + fRing_b*eta + fRing_c)/piM;
789 phenParams->sigma = (sigma_a*etap2 + sigma_b*eta + sigma_c)/piM;
790
791 phenParams->psi0 = (psi0_a*etap2 + psi0_b*eta + psi0_c)/(eta*pow(piM, 5./3.));
792 phenParams->psi1 = 0.;
793 phenParams->psi2 = (psi2_a*etap2 + psi2_b*eta + psi2_c)/(eta*pow(piM, 3./3.));
794 phenParams->psi3 = (psi3_a*etap2 + psi3_b*eta + psi3_c)/(eta*pow(piM, 2./3.));
795 phenParams->psi4 = (psi4_a*etap2 + psi4_b*eta + psi4_c)/(eta*pow(piM, 1./3.));
796 phenParams->psi5 = 0.;
797 phenParams->psi6 = (psi6_a*etap2 + psi6_b*eta + psi6_c)/(eta*pow(piM, -1./3.));
798 phenParams->psi7 = (psi7_a*etap2 + psi7_b*eta + psi7_c)/(eta*pow(piM, -2./3.));
799 }
800
801 return;
802
803}
804
805/*********************************************************************/
806/* Compute phenomenological parameters for aligned-spin binaries */
807/* Ref. Eq.(2) and Table I of http://arxiv.org/pdf/0909.2867 */
808/*********************************************************************/
810 BBHPhenomParams *phenParams,
812
813 REAL8 totalMass, piM, eta, chi, delta;
814 REAL8 etap2, chip2, etap3, etap2chi, etachip2, etachi;
815
816 if (!params || !phenParams) return;
817
818 /* calculate the total mass and symmetric mass ratio */
819 totalMass = params->mass1+params->mass2;
820 eta = params->mass1*params->mass2/pow(totalMass,2.);
821 piM = totalMass*LAL_PI*LAL_MTSUN_SI;
822 delta = sqrt(1.-4.*eta); /* asymmetry parameter */
823
824 /* spin parameter used for the search */
825 chi = 0.5*(params->spin1[2]*(1.+delta) + params->spin2[2]*(1.-delta));
826
827 /* spinning phenomenological waveforms */
828 etap2 = eta*eta;
829 chip2 = chi*chi;
830 etap3 = etap2*eta;
831 etap2chi = etap2*chi;
832 etachip2 = eta*chip2;
833 etachi = eta*chi;
834
835 phenParams->psi0 = 3./(128.*eta);
836
837 phenParams->psi2 = 3715./756. +
838 -9.2091e+02*eta + 4.9213e+02*etachi + 1.3503e+02*etachip2 +
839 6.7419e+03*etap2 + -1.0534e+03*etap2chi +
840 -1.3397e+04*etap3 ;
841
842 phenParams->psi3 = -16.*LAL_PI + 113.*chi/3. +
843 1.7022e+04*eta + -9.5659e+03*etachi + -2.1821e+03*etachip2 +
844 -1.2137e+05*etap2 + 2.0752e+04*etap2chi +
845 2.3859e+05*etap3 ;
846
847 phenParams->psi4 = 15293365./508032. - 405.*chip2/8. +
848 -1.2544e+05*eta + 7.5066e+04*etachi + 1.3382e+04*etachip2 +
849 8.7354e+05*etap2 + -1.6573e+05*etap2chi +
850 -1.6936e+06*etap3 ;
851
852 phenParams->psi6 = -8.8977e+05*eta + 6.3102e+05*etachi + 5.0676e+04*etachip2 +
853 5.9808e+06*etap2 + -1.4148e+06*etap2chi +
854 -1.1280e+07*etap3 ;
855
856 phenParams->psi7 = 8.6960e+05*eta + -6.7098e+05*etachi + -3.0082e+04*etachip2 +
857 -5.8379e+06*etap2 + 1.5145e+06*etap2chi +
858 1.0891e+07*etap3 ;
859
860 phenParams->psi8 = -3.6600e+05*eta + 3.0670e+05*etachi + 6.3176e+02*etachip2 +
861 2.4265e+06*etap2 + -7.2180e+05*etap2chi +
862 -4.5524e+06*etap3;
863
864 phenParams->fMerger = 1. - 4.4547*pow(1.-chi,0.217) + 3.521*pow(1.-chi,0.26) +
865 6.4365e-01*eta + 8.2696e-01*etachi + -2.7063e-01*etachip2 +
866 -5.8218e-02*etap2 + -3.9346e+00*etap2chi +
867 -7.0916e+00*etap3 ;
868
869 phenParams->fRing = (1. - 0.63*pow(1.-chi,0.3))/2. +
870 1.4690e-01*eta + -1.2281e-01*etachi + -2.6091e-02*etachip2 +
871 -2.4900e-02*etap2 + 1.7013e-01*etap2chi +
872 2.3252e+00*etap3 ;
873
874 phenParams->sigma = (1. - 0.63*pow(1.-chi,0.3))*pow(1.-chi,0.45)/4. +
875 -4.0979e-01*eta + -3.5226e-02*etachi + 1.0082e-01*etachip2 +
876 1.8286e+00*etap2 + -2.0169e-02*etap2chi +
877 -2.8698e+00*etap3 ;
878
879 phenParams->fCut = 3.2361e-01 + 4.8935e-02*chi + 1.3463e-02*chip2 +
880 -1.3313e-01*eta + -8.1719e-02*etachi + 1.4512e-01*etachip2 +
881 -2.7140e-01*etap2 + 1.2788e-01*etap2chi +
882 4.9220e+00*etap3 ;
883
884 phenParams->fCut /= piM;
885 phenParams->fMerger/= piM;
886 phenParams->fRing /= piM;
887 phenParams->sigma /= piM;
888
889 phenParams->psi1 = 0.;
890 phenParams->psi5 = 0.;
891}
892
893/*********************************************************************/
894/* Compute frequency-domain waveforms for non-spinning binaries */
895/* Ref. Eq.(4.13) and (4.16) of http://arxiv.org/pdf/0710.2335 */
896/*********************************************************************/
899 InspiralTemplate *insp_template,
900 REAL4Vector *signalvec) {
901
902 REAL8 df, shft, phi, amp0, ampEff=0, psiEff, fMerg, fNorm;
903 REAL8 f, fRing, sigma, totalMass, eta;
904 INT4 i, j, n, nby2;
905
906 /* freq resolution and the low-freq bin */
907 df = insp_template->tSampling/signalvec->length;
908 n = signalvec->length;
909
910 shft = LAL_TWOPI*(insp_template->nStartPad/insp_template->tSampling + insp_template->startTime);
911 phi = insp_template->startPhase;
912
913 /* phenomenological parameters*/
914 fMerg = params->fMerger;
915 fRing = params->fRing;
916 sigma = params->sigma;
917 totalMass = insp_template->mass1 + insp_template->mass2;
918 eta = insp_template->mass1 * insp_template->mass2 / pow(totalMass, 2.);
919
920 /* Now compute the amplitude. NOTE the params->distance is assumed to
921 * me in meters. This is, in principle, inconsistent with the LAL
922 * documentation (inspiral package). But this seems to be the convention
923 * employed in the injection codes */
924 amp0 = pow(LAL_MTSUN_SI*totalMass, 5./6.)*pow(fMerg,-7./6.)/pow(LAL_PI,2./3.);
925 amp0 *= pow(5.*eta/24., 1./2.)/(insp_template->distance/LAL_C_SI);
926
927 /* fill the zero and Nyquist frequency with zeros */
928 *(signalvec->data+0) = 0.;
929 *(signalvec->data+n/2) = 0.;
930
931 nby2 = n/2;
932
933 /* now generate the waveform at all frequency bins */
934 for (i=1; i<nby2; i++) {
935 /* this is the index of the imaginary part */
936 j = n-i;
937
938 /* fourier frequency corresponding to this bin */
939 f = i * df;
940 fNorm = f/fMerg;
941
942 /* compute the amplitude */
943 if ((f < insp_template->fLower) || (f > params->fCut))
944 ampEff = 0.;
945 else if (f <= fMerg)
946 ampEff = amp0*pow(fNorm, -7./6.);
947 else if ((f > fMerg) & (f <= fRing))
948 ampEff = amp0*pow(fNorm, -2./3.);
949 else if (f > fRing) {
950 ampEff = XLALLorentzianFn ( f, fRing, sigma);
951 ampEff *= amp0*LAL_PI_2*pow(fRing/fMerg,-2./3.)*sigma;
952 }
953
954 /* now compute the phase */
955 psiEff = shft*f + phi
956 + params->psi0*pow(f,-5./3.)
957 + params->psi1*pow(f,-4./3.)
958 + params->psi2*pow(f,-3./3.)
959 + params->psi3*pow(f,-2./3.)
960 + params->psi4*pow(f,-1./3.)
961 + params->psi5*pow(f,0.)
962 + params->psi6*pow(f,1./3.)
963 + params->psi7*pow(f,2./3.);
964
965 /* generate the waveform */
966 *(signalvec->data+i) = (REAL4) (ampEff * cos(psiEff)); /* real */
967 *(signalvec->data+j) = (REAL4) (ampEff * sin(psiEff)); /* imag */
968 }
969
970}
971
972
973/*********************************************************************/
974/* Compute frequency-domain waveforms for aligned-spin binaries */
975/* Ref. Eq.(1) of http://arxiv.org/pdf/0909.2867 */
976/*********************************************************************/
979 InspiralTemplate *insp_template,
980 REAL4Vector *signalvec) {
981
982 REAL8 df, shft, phi, amp0, ampEff=0, psiEff, fMerg, fNorm;
983 REAL8 f, fRing, sigma, totalMass, eta;
984 INT4 i, j, n, nby2;
985 REAL8 v, alpha2, alpha3, w1, vMerg, v2, v3, v4, v5, v6, v7, v8;
986 REAL8 epsilon_1, epsilon_2, w2, vRing, chi, mergPower, delta;
987
988 /* freq resolution and the low-freq bin */
989 df = insp_template->tSampling/signalvec->length;
990 n = signalvec->length;
991
992 shft = LAL_TWOPI*(insp_template->nStartPad/insp_template->tSampling + insp_template->startTime);
993
994 /* the negative sign is introduced such that the definition of the
995 * polarisations (phi = 0 for plus, phi = pi/2 for cross) is consistent
996 * with the IMRPhenomA waveforms */
997 phi = -insp_template->startPhase;
998
999 /* phenomenological parameters*/
1000 fMerg = params->fMerger;
1001 fRing = params->fRing;
1002 sigma = params->sigma;
1003
1004 /* physical parameters*/
1005 totalMass = insp_template->mass1 + insp_template->mass2;
1006 eta = insp_template->eta = insp_template->mass1 * insp_template->mass2 / pow(totalMass, 2.);
1007 delta = sqrt(1.-4.*eta);
1008 chi = 0.5*(insp_template->spin1[2]*(1.+delta) + insp_template->spin2[2]*(1.-delta));
1009
1010 /* Now compute the amplitude. NOTE the params->distance is assumed to
1011 * me in meters. This is, in principle, inconsistent with the LAL
1012 * documentation (inspiral package). But this seems to be the convention
1013 * employed in the injection codes */
1014 amp0 = pow(LAL_MTSUN_SI*totalMass, 5./6.)*pow(fMerg,-7./6.)/pow(LAL_PI,2./3.);
1015 amp0 *= pow(5.*eta/24., 1./2.)/(insp_template->distance/LAL_C_SI);
1016
1017 /* fill the zero and Nyquist frequency with zeros */
1018 *(signalvec->data+0) = 0.;
1019 *(signalvec->data+n/2) = 0.;
1020
1021 /***********************************************************************/
1022 /* these are the parameters required for the "new" phenomenological IMR
1023 * waveforms*/
1024 /***********************************************************************/
1025
1026 /* PN corrections to the frequency domain amplitude of the (2,2) mode */
1027 alpha2 = -323./224. + 451.*insp_template->eta/168.;
1028 alpha3 = (27./8. - 11.*insp_template->eta/6.)*chi;
1029
1030 /* leading order power law of the merger amplitude */
1031 mergPower = -2./3.;
1032
1033 /* spin-dependant corrections to the merger amplitude */
1034 epsilon_1 = 1.4547*chi - 1.8897;
1035 epsilon_2 = -1.8153*chi + 1.6557;
1036
1037 /* normalisation constant of the inspiral amplitude */
1038 vMerg = pow(LAL_PI*totalMass*LAL_MTSUN_SI*fMerg, 1./3.);
1039 vRing = pow(LAL_PI*totalMass*LAL_MTSUN_SI*fRing, 1./3.);
1040
1041 w1 = 1. + alpha2*pow(vMerg,2.) + alpha3*pow(vMerg,3.);
1042 w1 = w1/(1. + epsilon_1*vMerg + epsilon_2*vMerg*vMerg);
1043 w2 = w1*(LAL_PI*sigma/2.)*pow(fRing/fMerg, mergPower)*(1. + epsilon_1*vRing
1044 + epsilon_2*vRing*vRing);
1045
1046 /***********************************************************************/
1047 /* now generate the waveform at all frequency bins */
1048 /***********************************************************************/
1049 nby2 = n/2;
1050 for (i=1; i<nby2; i++) {
1051 /* this is the index of the imaginary part */
1052 j = n-i;
1053
1054 /* fourier frequency corresponding to this bin */
1055 f = i * df;
1056 fNorm = f/fMerg;
1057
1058 /* PN expansion parameter */
1059 v = pow(LAL_PI*totalMass*LAL_MTSUN_SI*f, 1./3.);
1060
1061 v2 = v*v; v3 = v2*v; v4 = v2*v2; v5 = v4*v; v6 = v3*v3; v7 = v6*v, v8 = v7*v;
1062
1063 /* compute the amplitude */
1064 if ((f < insp_template->fLower) || (f > params->fCut))
1065 ampEff = 0.;
1066 else if (f <= fMerg)
1067 ampEff = pow(fNorm, -7./6.)*(1. + alpha2*v2 + alpha3*v3);
1068 else if ((f > fMerg) & (f <= fRing))
1069 ampEff = w1*pow(fNorm, mergPower)*(1. + epsilon_1*v + epsilon_2*v2);
1070 else if (f > fRing)
1071 ampEff = w2*XLALLorentzianFn ( f, fRing, sigma);
1072
1073 /* now compute the phase */
1074 psiEff = shft*f + phi
1075 + 3./(128.*eta*v5)*(1 + params->psi2*v2
1076 + params->psi3*v3 + params->psi4*v4
1077 + params->psi5*v5 + params->psi6*v6
1078 + params->psi7*v7 + params->psi8*v8);
1079
1080 /* generate the waveform */
1081 *(signalvec->data+i) = (REAL4) (amp0 * ampEff * cos(psiEff)); /* real */
1082 *(signalvec->data+j) = (REAL4) (-amp0 * ampEff * sin(psiEff)); /* imag */
1083 }
1084}
1085
1086
1087/*
1088 * Private support functions
1089 */
1090
1092 REAL4Vector *Freq,
1093 REAL4Vector *hp,
1094 REAL4Vector *hc,
1095 REAL8 dt) {
1096 REAL8Vector *hpDot = NULL, *hcDot = NULL;
1097 UINT4 k, len;
1098 REAL8 fOfT;
1099
1100 len = hp->length;
1101
1102 hpDot= XLALCreateREAL8Vector(len);
1103 hcDot= XLALCreateREAL8Vector(len);
1104
1105 /* Construct the dot vectors (2nd order differencing) */
1106 hpDot->data[0] = 0.0;
1107 hpDot->data[len-1] = 0.0;
1108 hcDot->data[0] = 0.0;
1109 hcDot->data[len-1] = 0.0;
1110 for( k = 1; k < len-1; k++) {
1111 hpDot->data[k] = 1./(2.*dt)*(hp->data[k+1]-hp->data[k-1]);
1112 hcDot->data[k] = 1./(2.*dt)*(hc->data[k+1]-hc->data[k-1]);
1113 }
1114
1115 /* Compute frequency using the fact that */
1116 /*h(t) = A(t) e^(-i Phi) = Re(h) - i Im(h) = h_+ - i h_x */
1117 for( k = 0; k < len; k++) {
1118 fOfT = -hcDot->data[k] * hp->data[k] +hpDot->data[k] * hc->data[k];
1119 fOfT /= LAL_TWOPI;
1120 fOfT /= (pow(hp->data[k],2.) + pow(hc->data[k], 2.));
1121 Freq->data[k] = (REAL4) fOfT;
1122 }
1123
1124 /* free the memory allocated for the derivative vectors */
1127
1128 return;
1129
1130}
1131
1133 REAL8 freq,
1134 REAL8 fRing,
1135 REAL8 sigma) {
1136 return sigma / (2 * LAL_PI * ((freq - fRing)*(freq - fRing)
1137 + sigma*sigma / 4.0));
1138}
1139
1141 REAL4Vector *h,
1142 REAL4Vector *freq,
1143 REAL8 cutFreq,
1144 REAL8 UNUSED deltaT) {
1145 UINT4 k, k0, kMid, len;
1146 REAL4 currentFreq;
1147 len = freq->length;
1148
1149 /* Since the boundaries of this freq vector are likely to have */
1150 /* FFT crap, let's scan the freq values starting from the middle */
1151 kMid = len/2;
1152 currentFreq = freq->data[kMid];
1153 k = kMid;
1154
1155 /* freq is an increasing function of time */
1156 /* If we are above the cutFreq we move to the left; else to the right */
1157 if (currentFreq > cutFreq && k > 0) {
1158 while (currentFreq > cutFreq)
1159 currentFreq = freq->data[k--];
1160 k0 = k;
1161 }
1162 else {
1163 while (currentFreq < cutFreq && k < len)
1164 currentFreq = freq->data[k++];
1165 k0 = k;
1166 }
1167
1168 /* zero the waveform below the cutoff frequency */
1169 for (k = 0; k < k0; k++)
1170 h->data[k] = 0.0;
1171
1172 return h;
1173}
1174
1175
1176/*
1177 * Deprecated LAL wrapper functions
1178 */
1179
1182 REAL4Vector *signalvec,
1186 XLAL_PRINT_DEPRECATION_WARNING("XLALBBHPhenWaveFreqDom");
1187 switch (params->approximant) {
1188 case IMRPhenomA:
1190 break;
1191 case IMRPhenomB:
1193 break;
1194 default:
1195 ABORT(status, LALINSPIRALH_ESWITCH, LALINSPIRALH_MSGESWITCH);
1196 }
1198 RETURN(status);
1199}
1200
1203 REAL4Vector *signalvec1,
1204 REAL4Vector *signalvec2,
1208 XLAL_PRINT_DEPRECATION_WARNING("XLALBBHPhenWaveFreqDomTemplates");
1209 switch (params->approximant) {
1210 case IMRPhenomA:
1211 if (XLALBBHPhenWaveAFreqDomTemplates(signalvec1, signalvec2, params)) ABORTXLAL(status);
1212 break;
1213 case IMRPhenomB:
1214 if (XLALBBHPhenWaveAFreqDomTemplates(signalvec1, signalvec2, params)) ABORTXLAL(status);
1215 break;
1216 default:
1217 ABORT(status, LALINSPIRALH_ESWITCH, LALINSPIRALH_MSGESWITCH);
1218 }
1220 RETURN(status);
1221}
1222
1225 REAL4Vector *signalvec1,
1226 InspiralTemplate *insp_template) {
1229 XLAL_PRINT_DEPRECATION_WARNING("XLALBBHPhenWaveTimeDom");
1230 if (XLALBBHPhenWaveTimeDom(signalvec1, insp_template)) ABORTXLAL(status);
1233 RETURN (status);
1234}
1235
1238 REAL4Vector *signalvec1,
1239 REAL4Vector *signalvec2,
1240 InspiralTemplate *insp_template) {
1243 XLAL_PRINT_DEPRECATION_WARNING("XLALBBHPhenWaveTimeDomTemplates");
1244 if (XLALBBHPhenWaveTimeDomTemplates(signalvec1, signalvec2, insp_template))
1248 RETURN (status);
1249}
1250
1253 REAL4Vector *signalvec1,
1254 REAL4Vector *signalvec2,
1255 REAL4Vector *h,
1256 REAL4Vector *aVec,
1257 REAL4Vector *freqVec,
1258 REAL8Vector *phiVec,
1259 UINT4 *countback,
1263 XLAL_PRINT_DEPRECATION_WARNING("XLALBBHPhenTimeDomEngine");
1264 if (XLALBBHPhenTimeDomEngine(signalvec1, signalvec2, h, aVec,
1265 freqVec, phiVec, countback, params)) ABORTXLAL(status);
1267 RETURN (status);
1268}
1269
1272 CoherentGW *waveform,
1274 PPNParamStruc *ppnParams) {
1277 XLAL_PRINT_DEPRECATION_WARNING("XLALBBHPhenWaveTimeDomForInjection");
1278 if (XLALBBHPhenWaveTimeDomForInjection(waveform, params, ppnParams))
1281 RETURN(status);
1282}
int XLALInspiralInit(InspiralTemplate *params, InspiralInit *paramsInit)
static UNUSED REAL4Vector * XLALCutAtFreq(REAL4Vector *h, REAL4Vector *freq, REAL8 cutFreq, REAL8 deltaT)
int XLALBBHPhenWaveTimeDomTemplates(REAL4Vector *signalvec1, REAL4Vector *signalvec2, InspiralTemplate *insp_template)
static void XLALComputePhenomParams(BBHPhenomParams *phenParams, InspiralTemplate *params)
void LALBBHPhenWaveTimeDomTemplates(LALStatus *status, REAL4Vector *signalvec1, REAL4Vector *signalvec2, InspiralTemplate *insp_template)
static void XLALComputeInstantFreq(REAL4Vector *Freq, REAL4Vector *hp, REAL4Vector *hc, REAL8 dt)
int XLALBBHPhenWaveAFreqDomTemplates(REAL4Vector *signalvec1, REAL4Vector *signalvec2, InspiralTemplate *params)
int XLALBBHPhenWaveBFreqDom(REAL4Vector *signalvec, InspiralTemplate *params)
int XLALBBHPhenWaveBFreqDomTemplates(REAL4Vector *signalvec1, REAL4Vector *signalvec2, InspiralTemplate *params)
int XLALBBHPhenTimeDomEngine(REAL4Vector *signalvec1, REAL4Vector *signalvec2, REAL4Vector *h, REAL4Vector *aVec, REAL4Vector *freqVec, REAL8Vector *phiVec, UINT4 *countback, InspiralTemplate *params)
static void XLALBBHPhenWaveFD(BBHPhenomParams *params, InspiralTemplate *insp_template, REAL4Vector *signalvec)
void LALBBHPhenWaveFreqDomTemplates(LALStatus *status, REAL4Vector *signalvec1, REAL4Vector *signalvec2, InspiralTemplate *params)
void LALBBHPhenWaveTimeDom(LALStatus *status, REAL4Vector *signalvec1, InspiralTemplate *insp_template)
static void XLALComputePhenomParams2(BBHPhenomParams *phenParams, InspiralTemplate *params)
static void XLALBBHPhenWaveFD2(BBHPhenomParams *params, InspiralTemplate *insp_template, REAL4Vector *signalvec)
void LALBBHPhenTimeDomEngine(LALStatus *status, REAL4Vector *signalvec1, REAL4Vector *signalvec2, REAL4Vector *h, REAL4Vector *aVec, REAL4Vector *freqVec, REAL8Vector *phiVec, UINT4 *countback, InspiralTemplate *params)
int XLALBBHPhenWaveTimeDomForInjection(CoherentGW *waveform, InspiralTemplate *params, PPNParamStruc *ppnParams)
static REAL8 XLALLorentzianFn(REAL8 freq, REAL8 fRing, REAL8 sigma)
int XLALBBHPhenWaveAFreqDom(REAL4Vector *signalvec, InspiralTemplate *params)
void LALBBHPhenWaveFreqDom(LALStatus *status, REAL4Vector *signalvec, InspiralTemplate *params)
void LALBBHPhenWaveTimeDomForInjection(LALStatus *status, CoherentGW *waveform, InspiralTemplate *params, PPNParamStruc *ppnParams)
int XLALBBHPhenWaveTimeDom(REAL4Vector *signalvec1, InspiralTemplate *insp_template)
static double double delta
#define ABORT(statusptr, code, mesg)
#define CHECKSTATUSPTR(statusptr)
#define ATTATCHSTATUSPTR(statusptr)
#define DETATCHSTATUSPTR(statusptr)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
#define ABORTXLAL(sp)
int s
int l
double i
coeffs k0
const double w2
void XLALDestroyREAL4VectorSequence(REAL4VectorSequence *vecseq)
REAL4VectorSequence * XLALCreateREAL4VectorSequence(UINT4 length, UINT4 veclen)
#define BBHPHENOMCOEFFSH_PSI4_X
#define BBHPHENOMCOEFFSH_PSI3_Y
#define BBHPHENOMCOEFFSH_PSI2_Z
#define BBHPHENOMCOEFFSH_PSI6_Z
#define BBHPHENOMCOEFFSH_FRING_C
#define BBHPHENOMCOEFFSH_PSI2_Y
#define BBHPHENOMCOEFFSH_FRING_B
#define BBHPHENOMCOEFFSH_SIGMA_A
#define BBHPHENOMCOEFFSH_FRING_A
#define BBHPHENOMCOEFFSH_PSI0_X
#define BBHPHENOMCOEFFSH_PSI3_X
#define BBHPHENOMCOEFFSH_PSI3_Z
#define BBHPHENOMCOEFFSH_PSI6_X
#define BBHPHENOMCOEFFSH_PSI6_Y
#define BBHPHENOMCOEFFSH_FMERG_B
#define BBHPHENOMCOEFFSH_FCUT_B
#define BBHPHENOMCOEFFSH_PSI2_X
#define BBHPHENOMCOEFFSH_PSI7_X
#define BBHPHENOMCOEFFSH_PSI0_Z
#define BBHPHENOMCOEFFSH_PSI4_Y
#define BBHPHENOMCOEFFSH_SIGMA_B
#define BBHPHENOMCOEFFSH_SIGMA_C
#define BBHPHENOMCOEFFSH_FMERG_A
#define BBHPHENOMCOEFFSH_PSI4_Z
#define BBHPHENOMCOEFFSH_PSI7_Y
#define BBHPHENOMCOEFFSH_FMERG_C
#define BBHPHENOMCOEFFSH_FCUT_C
#define BBHPHENOMCOEFFSH_FCUT_A
#define BBHPHENOMCOEFFSH_PSI7_Z
#define BBHPHENOMCOEFFSH_PSI0_Y
#define GENERATEPPNINSPIRALH_EFSTOP
Reached requested termination frequency.
#define LAL_PI_2
#define LAL_C_SI
#define LAL_PI
#define LAL_TWOPI
#define LAL_MTSUN_SI
double REAL8
char CHAR
uint32_t UINT4
int32_t INT4
float REAL4
LALNameLength
#define LALINSPIRALH_ENULL
Arguments contained an unexpected null pointer.
Definition: LALInspiral.h:56
#define LALINSPIRALH_ESWITCH
Unknown case in switch.
Definition: LALInspiral.h:75
void * XLALMalloc(size_t n)
void XLALFree(void *p)
LALPNOrder
IMRPhenomA
IMRPhenomB
static const INT4 a
REAL4FFTPlan * XLALCreateReverseREAL4FFTPlan(UINT4 size, int measurelvl)
void XLALDestroyREAL4FFTPlan(REAL4FFTPlan *plan)
int XLALREAL4VectorFFT(REAL4Vector *_LAL_RESTRICT_ output, const REAL4Vector *_LAL_RESTRICT_ input, const REAL4FFTPlan *plan)
REAL4TimeSeries * XLALCreateREAL4TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
void XLALDestroyREAL4TimeSeries(REAL4TimeSeries *series)
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
const LALUnit lalStrainUnit
const LALUnit lalHertzUnit
const LALUnit lalDimensionlessUnit
REAL4Vector * XLALCreateREAL4Vector(UINT4 length)
void XLALDestroyREAL4Vector(REAL4Vector *vector)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
int XLALREAL8VectorUnwrapAngle(REAL8Vector *out, const REAL8Vector *in)
void XLALError(const char *func, const char *file, int line, int errnum)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_ERROR(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
int int XLALPrintWarning(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_PRINT_DEPRECATION_WARNING(replacement)
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFAULT
XLAL_ERANGE
XLAL_ENAME
XLAL_EFUNC
XLAL_EDOM
XLAL_FAILURE
int n
REAL4TimeSeries * shift
SkyPosition position
REAL4TimeVectorSeries * a
REAL8TimeSeries * phi
REAL4TimeVectorSeries * h
REAL4TimeSeries * f
The inspiral waveform parameter structure containing information about the waveform to be generated.
Definition: LALInspiral.h:205
REAL8 startTime
starting time of the waveform (in sec); if different from zero, the waveform will start with an insta...
Definition: LALInspiral.h:222
INT4 nStartPad
Number of leading elements in the signal generation to be set to zero (input) If template is requeste...
Definition: LALInspiral.h:307
REAL8 eta
symmetric mass ratio (input/output)
Definition: LALInspiral.h:291
REAL8 spin2[3]
Spin vector of the secondary.
Definition: LALInspiral.h:266
REAL8 mass1
Mass of the primary in solar mass (input/output)
Definition: LALInspiral.h:211
REAL8 distance
Distance to the binary in seconds.
Definition: LALInspiral.h:219
REAL8 startPhase
starting phase of the waveform in radians (input)
Definition: LALInspiral.h:221
REAL8 spin1[3]
Spin vector of the primary.
Definition: LALInspiral.h:265
REAL8 mass2
Mass of the secondary in solar mass (mass1 need not be larger than mass2 (input/output)
Definition: LALInspiral.h:212
REAL8 tSampling
Sampling rate in Hz (input)
Definition: LALInspiral.h:218
This structure stores the parameters for constructing a restricted post-Newtonian waveform.
const CHAR * termDescription
The termination code description (above)
INT4 termCode
The termination condition (above) that stopped computation of the waveform.
REAL4 dfdt
The maximum value of encountered over any timestep used in generating the waveform.
REAL4 phi
The phase at coalescence (or arbitrary reference phase for a post -Newtonian approximation),...
REAL4 psi
polarization angle (radians)
REAL4 fStart
The actual starting frequency of the waveform, in Hz (normally close but not identical to fStartIn)
SkyPosition position
location of source on sky
UINT4 length
The length of the generated waveform.
REAL4 fStop
The frequency at the termination of the waveform, in Hz.
REAL8 tc
The time from the start of the waveform to coalescence (in the point-mass approximation),...
REAL4 fStartIn
The requested starting frequency of the waveform, in Hz.
REAL8 deltaT
The requested sampling interval of the waveform, in s.
REAL4Sequence * data
CHAR name[LALNameLength]
REAL4VectorSequence * data
REAL4 * data
REAL8Sequence * data
REAL8 * data
double inclination
double df