LALSimulation  5.4.0.1-fe68b98
LALSimIMREOBHybridRingdownPrec.c
Go to the documentation of this file.
1 /**
2  * \author Yi Pan, Craig Robinson, Prayush Kumar, Stas Babak, Andrea Taracchini
3 
4  *
5  * \brief Module to compute the ring-down waveform as linear combination
6  * of quasi-normal-modes decaying waveforms, which can be attached to
7  * the inspiral part of the compat binary coalescing waveform.
8  * The method is describe in Sec. II C of Pan et al. PRD 84, 124052 (2011),
9  * specifically Eqs. 30 - 32.
10  * Eqs. 30 and 31 are written in explicity linear equation systems in
11  * DCC document T1100433.
12  * This method is currently used for EOBNRv2 and SEOBNRv1 models. The only difference
13  * between the two models in ring-down waveform is the pseudo-QNM introduced
14  * in the latter (see Taracchini et al. PRD 86, 024011 (2012) for more details).
15  */
16 #ifndef _LALSIMIMREOBHYBRIDRINGDOWNPREC_C
17 #define _LALSIMIMREOBHYBRIDRINGDOWNPREC_C
18 
19 #include <math.h>
20 #include <complex.h>
21 #include <stdlib.h>
22 #include <lal/LALStdlib.h>
23 #include <lal/AVFactories.h>
24 #include <lal/SeqFactories.h>
25 #include <gsl/gsl_linalg.h>
26 #include <gsl/gsl_interp.h>
27 #include <gsl/gsl_spline.h>
28 #include <gsl/gsl_blas.h>
29 
30 //#include "LALSimIMREOBNRv2.h"
31 //#include "LALSimBlackHoleRingdown.h"
32 #include "LALSimBlackHoleRingdownPrec.h"
33 #include "LALSimIMREOBNQCCorrection.c"
34 //#include "LALSimIMREOBHybridRingdown.c"
35 
36 #define LMax 6
37 #define MMax 6
38 
39 #ifdef __GNUC__
40 #define UNUSED __attribute__ ((unused))
41 #else
42 #define UNUSED
43 #endif
44 
45 //////////////////////////////////////////////////////////////////////////////////////////////////////////
46 /* Defining matrices to store the coefficients that enter amplitude and phase of the RD signal of SEOBNRv4 and SEOBNRv4HM */
47 
48 
49 static REAL8 XLALCalculateRDAmplitudeCoefficient1( UINT4 modeL, UINT4 modeM, REAL8 eta, REAL8 chi){
50  REAL8 A1coeff00[LMax][MMax]={{0}};
51  REAL8 A1coeff01[LMax][MMax]={{0}};
52  REAL8 A1coeff02[LMax][MMax]={{0}};
53  REAL8 A1coeff10[LMax][MMax]={{0}};
54  REAL8 A1coeff11[LMax][MMax]={{0}};
55  REAL8 A1coeff20[LMax][MMax]={{0}};
56  REAL8 A1coeff21[LMax][MMax]={{0}};
57  REAL8 A1coeff31[LMax][MMax]={{0}};
58 
59  //////////////////////////////////////////////////////////////////////////////////////////////////////////
60  /* Fit coefficients that enter amplitude of the RD signal for SEOBNRv4 and SEOBNRv4HM*/
61  /* Eq. (54) of https://dcc.ligo.org/T1600383 for SEOBNRv4 and 22 mode SEOBNRv4HM*/
62  /* Eqs. (C1, C3, C5, C7) of https://arxiv.org/pdf/1803.10701.pdf for SEOBNRv4HM*/
63 
64  A1coeff00[2][2] = 0.0830664;
65  A1coeff01[2][2] = -0.0196758;
66  A1coeff02[2][2] = -0.0136459;
67  A1coeff10[2][2] = 0.0612892;
68  A1coeff11[2][2] = 0.00146142;
69  A1coeff20[2][2] = -0.0893454;
70 
71  A1coeff00[2][1] = 0.07780330893915006;
72  A1coeff01[2][1] = -0.05070638166864379;
73  A1coeff10[2][1] = 0.24091006920322164;
74  A1coeff11[2][1] = 0.38582622780596576;
75  A1coeff20[2][1] = -0.7456327888190485;
76  A1coeff21[2][1] = -0.9695534075470388;
77 
78  A1coeff00[3][3] = 0.07638733045623343;
79  A1coeff01[3][3] = -0.030993441267953236;
80  A1coeff10[3][3] = 0.2543447497371546;
81  A1coeff11[3][3] = 0.2516879591102584;
82  A1coeff20[3][3] = -1.0892686061231245;
83  A1coeff21[3][3] = -0.7980907313033606;
84 
85  A1coeff00[4][4] = -0.06392710223439678;
86  A1coeff01[4][4] = -0.03646167590514318;
87  A1coeff10[4][4] = 0.345195277237925;
88  A1coeff11[4][4] = 1.2777441574118558;
89  A1coeff20[4][4] = -1.764352185878576;
90  A1coeff21[4][4] = -14.825262897834696;
91  A1coeff31[4][4] = 40.67135475479875;
92 
93  A1coeff00[5][5] = -0.06704614393611373;
94  A1coeff01[5][5] = 0.021905949257025503;
95  A1coeff10[5][5] = -0.24754936787743445;
96  A1coeff11[5][5] = -0.0943771022698497;
97  A1coeff20[5][5] = 0.7588042862093705;
98  A1coeff21[5][5] = 0.4357768883690394;
99 
100 
101 
102  return A1coeff00[modeL][modeM] + A1coeff01[modeL][modeM] * chi + A1coeff02[modeL][modeM] * chi * chi + A1coeff10[modeL][modeM] * eta +
103  A1coeff11[modeL][modeM] * eta * chi + A1coeff20[modeL][modeM] * eta * eta + A1coeff21[modeL][modeM]*eta*eta*chi + A1coeff31[modeL][modeM]*eta*eta*eta*chi;
104 
105 }
106 
107 
108 static REAL8 XLALCalculateRDAmplitudeCoefficient2( UINT4 modeL, UINT4 modeM, REAL8 eta, REAL8 chi){
109 
110  REAL8 A2coeff00[LMax][MMax]={{0}};
111  REAL8 A2coeff01[LMax][MMax]={{0}};
112  REAL8 A2coeff10[LMax][MMax]={{0}};
113  REAL8 A2coeff11[LMax][MMax]={{0}};
114  REAL8 A2coeff20[LMax][MMax]={{0}};
115  REAL8 A2coeff21[LMax][MMax]={{0}};
116  REAL8 A2coeff31[LMax][MMax]={{0}};
117 
118 
119  //////////////////////////////////////////////////////////////////////////////////////////////////////////
120  /* Fit coefficients that enter amplitude of the RD signal of SEOBNRv4 and SEOBNRv4HM*/
121  /* Eq. (55) of https://dcc.ligo.org/T1600383 for SEOBNRv4 and 22 mode SEOBNRv4HM*/
122  /* Eqs. (C1, C3, C5, C7) of https://arxiv.org/pdf/1803.10701.pdf for SEOBNRv4HM*/
123  A2coeff00[2][2] = -0.623953;
124  A2coeff01[2][2] = -0.371365;
125  A2coeff10[2][2] = 1.39777;
126  A2coeff11[2][2] = 2.40203;
127  A2coeff20[2][2] = -1.82173;
128  A2coeff21[2][2] = -5.25339;
129 
130  A2coeff00[2][1] = -1.2451852641667298;
131  A2coeff01[2][1] = -1.195786238319961;
132  A2coeff10[2][1] = 6.134202504312409;
133  A2coeff11[2][1] = 15.66696619631313;
134  A2coeff20[2][1] = -14.67251127485556;
135  A2coeff21[2][1] = -44.41982201432511;
136 
137  A2coeff00[3][3] = -0.8325292359346013;
138  A2coeff01[3][3] = -0.598880303198448;
139  A2coeff10[3][3] = 2.767989795032018;
140  A2coeff11[3][3] = 5.904371617277156;
141  A2coeff20[3][3] = -7.028151926115957;
142  A2coeff21[3][3] = -18.232606706124482;
143 
144  A2coeff00[4][4] = 0.7813275473485185;
145  A2coeff01[4][4] = 0.8094706044462984;
146  A2coeff10[4][4] = -5.18689829943586;
147  A2coeff11[4][4] = -5.38343327318501;
148  A2coeff20[4][4] = 14.026415859369477;
149  A2coeff21[4][4] = 0.1051625997942299;
150  A2coeff31[4][4] = 46.978434956814006;
151 
152  A2coeff00[5][5] = 1.6763424265367357;
153  A2coeff01[5][5] = 0.4925695499534606;
154  A2coeff10[5][5] = -5.604559311983177;
155  A2coeff11[5][5] = -6.209095657439377;
156  A2coeff20[5][5] = 16.751319143123386;
157  A2coeff21[5][5] = 16.778452555342554;
158 
159  return A2coeff00[modeL][modeM] + A2coeff01[modeL][modeM] * chi + A2coeff10[modeL][modeM] * eta +
160  A2coeff11[modeL][modeM] * eta * chi + A2coeff20[modeL][modeM] * eta * eta + A2coeff21[modeL][modeM] * eta * eta * chi + A2coeff31[modeL][modeM]*eta*eta*eta*chi;
161 }
162 
163 static REAL8 XLALCalculateRDPhaseCoefficient1( UINT4 modeL, UINT4 modeM, REAL8 eta, REAL8 chi){
164 
165  REAL8 P1coeff00[LMax][MMax]={{0}};
166  REAL8 P1coeff01[LMax][MMax]={{0}};
167  REAL8 P1coeff02[LMax][MMax]={{0}};
168  REAL8 P1coeff10[LMax][MMax]={{0}};
169  REAL8 P1coeff11[LMax][MMax]={{0}};
170  REAL8 P1coeff20[LMax][MMax]={{0}};
171  REAL8 P1coeff21[LMax][MMax]={{0}};
172  REAL8 P1coeff31[LMax][MMax]={{0}};
173 
174  //////////////////////////////////////////////////////////////////////////////////////////////////////////
175  /* Fit coefficients that enter the phase of the RD signal of SEOBNRv4 and SEOBNRv4HM*/
176  /* Eq. (62) of https://dcc.ligo.org/T1600383 for SEOBNRv4 and 22 mode SEOBNRv4HM*/
177  /* Eqs. (C9, C11, C13, C15) of https://arxiv.org/pdf/1803.10701.pdf for SEOBNRv4HM*/
178 
179  P1coeff00[2][2] = 0.147584;
180  P1coeff01[2][2] = 0.00779176;
181  P1coeff02[2][2] = -0.0244358;
182  P1coeff10[2][2] = 0.263456;
183  P1coeff11[2][2] = -0.120853;
184  P1coeff20[2][2] = -0.808987;
185 
186  P1coeff00[2][1] = 0.15601401627613815;
187  P1coeff01[2][1] = 0.10219957917717622;
188  P1coeff10[2][1] = 0.023346852452720928;
189  P1coeff11[2][1] = -0.9435308286367039;
190  P1coeff20[2][1] = 0.15326558178697175;
191  P1coeff21[2][1] = 1.7979082057513565;
192 
193  P1coeff00[3][3] = 0.11085299117493969;
194  P1coeff01[3][3] = 0.018959099261813613;
195  P1coeff10[3][3] = 0.9999800463662053;
196  P1coeff11[3][3] = -0.729149797691242;
197  P1coeff20[3][3] = -3.3983315694441125;
198  P1coeff21[3][3] = 2.5192011762934037;
199 
200  P1coeff00[4][4] = 0.11498976343440313;
201  P1coeff01[4][4] = 0.008389519706605305;
202  P1coeff10[4][4] = 1.6126522800609633;
203  P1coeff11[4][4] = -0.8069979888526699;
204  P1coeff20[4][4] = -6.255895564079467;
205  P1coeff21[4][4] = 7.595651881827078;
206  P1coeff31[4][4] = -19.32367406125053;
207 
208  P1coeff00[5][5] = 0.16465380962882128;
209  P1coeff01[5][5] = -0.026574817803812007;
210  P1coeff10[5][5] = -0.19184523794508765;
211  P1coeff11[5][5] = -0.05519618962738479;
212  P1coeff20[5][5] = 0.33328424135336965;
213  P1coeff21[5][5] = 0.3194274548351241;
214 
215 
216  return P1coeff00[modeL][modeM] + P1coeff01[modeL][modeM] * chi + P1coeff02[modeL][modeM] * chi * chi + P1coeff10[modeL][modeM] * eta +
217  P1coeff11[modeL][modeM] * eta * chi + P1coeff20[modeL][modeM] * eta * eta + P1coeff21[modeL][modeM]*eta*eta*chi + P1coeff31[modeL][modeM]*eta*eta*eta*chi;
218 
219 }
220 
221 static REAL8 XLALCalculateRDPhaseCoefficient2( UINT4 modeL, UINT4 modeM, REAL8 eta, REAL8 chi){
222 
223  REAL8 P2coeff00[LMax][MMax]={{0}};
224  REAL8 P2coeff01[LMax][MMax]={{0}};
225  REAL8 P2coeff02[LMax][MMax]={{0}};
226  REAL8 P2coeff12[LMax][MMax]={{0}};
227  REAL8 P2coeff10[LMax][MMax]={{0}};
228  REAL8 P2coeff11[LMax][MMax]={{0}};
229  REAL8 P2coeff20[LMax][MMax]={{0}};
230  REAL8 P2coeff21[LMax][MMax]={{0}};
231  REAL8 P2coeff31[LMax][MMax]={{0}};
232 
233  //////////////////////////////////////////////////////////////////////////////////////////////////////////
234  /* Fit coefficients that enter the phase of the RD signal of SEOBNRv4 and SEOBNRv4HM*/
235  /* Eq. (63) of https://dcc.ligo.org/T1600383 for SEOBNRv4 and 22 mode SEOBNRv4HM*/
236  /* Eqs. (C10, C12, C14, C16) of https://arxiv.org/pdf/1803.10701.pdf for SEOBNRv4HM*/
237  P2coeff00[2][2] = 2.46654;
238  P2coeff01[2][2] = 3.13067;
239  P2coeff02[2][2] = 0.581626;
240  P2coeff10[2][2] = -6.99396;
241  P2coeff11[2][2] = -9.61861;
242  P2coeff20[2][2] = 17.5646;
243 
244  P2coeff00[2][1] = 2.7886287922318105;
245  P2coeff01[2][1] = 4.29290053494256;
246  P2coeff10[2][1] = -0.8145406685320334;
247  P2coeff11[2][1] = -15.93796979597706;
248  P2coeff20[2][1] = 5.549338798935832;
249  P2coeff21[2][1] = 12.649775582333442;
250  P2coeff02[2][1] = 2.5582321247274726;
251  P2coeff12[2][1] = -10.232928498772893;
252 
253  P2coeff00[3][3] = 2.7825237371542735;
254  P2coeff01[3][3] = 2.8796835808075003;
255  P2coeff10[3][3] = -7.844741660437831;
256  P2coeff11[3][3] = -34.7670039322078;
257  P2coeff20[3][3] = 27.181024362399302;
258  P2coeff21[3][3] = 127.13948436435182;
259 
260  P2coeff00[4][4] = 3.111817347262856;
261  P2coeff01[4][4] = 5.399341180960216;
262  P2coeff10[4][4] = 15.885333959709488;
263  P2coeff11[4][4] = -87.92421137153823;
264  P2coeff20[4][4] = -79.64931908155609;
265  P2coeff21[4][4] = 657.7156442271963;
266  P2coeff31[4][4] = -1555.2968529739226;
267  P2coeff02[4][4] = 2.3832321567874686;
268  P2coeff12[4][4] = -9.532928476043567;
269 
270  P2coeff00[5][5] = 11.102447263357977;
271  P2coeff01[5][5] = 6.015112119742853;
272  P2coeff10[5][5] = -58.605776859097084;
273  P2coeff11[5][5] = -81.68032025902797;
274  P2coeff20[5][5] = 176.60643662729498;
275  P2coeff21[5][5] = 266.47345742836745;
276 
277 
278  return P2coeff00[modeL][modeM] + P2coeff01[modeL][modeM] * chi + P2coeff02[modeL][modeM] * chi * chi + P2coeff12[modeL][modeM] * chi * chi * eta + P2coeff10[modeL][modeM] * eta +
279  P2coeff11[modeL][modeM] * eta * chi + P2coeff20[modeL][modeM] * eta * eta + P2coeff21[modeL][modeM]*eta*eta*chi + P2coeff31[modeL][modeM]*eta*eta*eta*chi;
280 }
281 
282 
283 static INT4 XLALSimFindIndexMaxAmpli( UINT4 * indAmax, REAL8Vector * timeVec, REAL8Vector * ampWave, REAL8 * valAmax, REAL8 tofAmax );
284 
285 //static INT4 XLALComputePhenomFitRD( REAL8* ampcf1, REAL8* ampcf2, REAL8* phasecf1, REAL8* phasecf2, REAL8 chi, REAL8 eta, INT4 l, INT4 m);
286 
287 /**
288  * Generates the ringdown wave associated with the given real
289  * and imaginary parts of the inspiral waveform. The parameters of
290  * the ringdown, such as amplitude and phase offsets, are determined
291  * by solving the linear equations defined in the DCC document T1100433.
292  * In the linear equations Ax=y,
293  * A is a 16-by-16 matrix depending on QNM (complex) frequencies,
294  * x is a 16-d vector of the 8 unknown complex QNM amplitudes,
295  * y is a 16-d vector depending on inspiral-plunge waveforms and their derivatives near merger.
296  */
297 static INT4 XLALSimIMREOBHybridRingdownWave(
298  REAL8Vector *rdwave1, /**<< OUTPUT, Real part of ringdown waveform */
299  REAL8Vector *rdwave2, /**<< OUTPUT, Imag part of ringdown waveform */
300  const REAL8 dt, /**<< Sampling interval */
301  const REAL8 mass1, /**<< First component mass (in Solar masses) */
302  const REAL8 mass2, /**<< Second component mass (in Solar masses) */
303  REAL8VectorSequence *inspwave1, /**<< Values and derivs of real part inspiral waveform */
304  REAL8VectorSequence *inspwave2, /**<< Values and derivs of imag part inspiral waveform */
305  COMPLEX16Vector *modefreqs, /**<< Complex freqs of ringdown (scaled by total mass) */
306  REAL8Vector *matchrange /**<< Times which determine the comb of ringdown attachment */
307 )
308 {
309 
310  /* XLAL error handling */
311  INT4 errcode = XLAL_SUCCESS;
312 
313  /* For checking GSL return codes */
314  INT4 gslStatus;
315  INT4 debugSB = 0;
316 
317  UINT4 i, j, k, nmodes = 8;
318 
319  /* Sampling rate from input */
320  REAL8 t1, t2, t3, t4, t5, rt;
321  gsl_matrix *coef;
322  gsl_vector *hderivs;
323  gsl_vector *x;
324  gsl_permutation *p;
325  REAL8Vector *modeamps;
326  int s;
327  REAL8 tj;
328  REAL8 m;
329 
330  /* mass in geometric units */
331  m = (mass1 + mass2) * LAL_MTSUN_SI;
332  t5 = (matchrange->data[0] - matchrange->data[1]) * m;
333  rt = -t5 / 5.;
334 
335  t4 = t5 + rt;
336  t3 = t4 + rt;
337  t2 = t3 + rt;
338  t1 = t2 + rt;
339 
340  if ( inspwave1->length != 3 || inspwave2->length != 3 ||
341  modefreqs->length != nmodes )
342  {
343  XLAL_ERROR( XLAL_EBADLEN );
344  }
345 
346  /* Solving the linear system for QNMs amplitude coefficients using gsl routine */
347  /* Initiate matrices and supporting variables */
348  XLAL_CALLGSL( coef = (gsl_matrix *) gsl_matrix_alloc(2 * nmodes, 2 * nmodes) );
349  XLAL_CALLGSL( hderivs = (gsl_vector *) gsl_vector_alloc(2 * nmodes) );
350  XLAL_CALLGSL( x = (gsl_vector *) gsl_vector_alloc(2 * nmodes) );
351  XLAL_CALLGSL( p = (gsl_permutation *) gsl_permutation_alloc(2 * nmodes) );
352 
353  /* Check all matrices and variables were allocated */
354  if ( !coef || !hderivs || !x || !p )
355  {
356  if (coef) gsl_matrix_free(coef);
357  if (hderivs) gsl_vector_free(hderivs);
358  if (x) gsl_vector_free(x);
359  if (p) gsl_permutation_free(p);
360 
361  XLAL_ERROR( XLAL_ENOMEM );
362  }
363 
364  /* Define the linear system Ax=y */
365  /* Matrix A (2*n by 2*n) has block symmetry. Define half of A here as "coef" */
366  /* The half of A defined here corresponds to matrices M1 and -M2 in the DCC document T1100433 */
367  /* Define y here as "hderivs" */
368  for (i = 0; i < nmodes; ++i)
369  {
370  gsl_matrix_set(coef, 0, i, 1);
371  gsl_matrix_set(coef, 1, i, - cimag(modefreqs->data[i]));
372  gsl_matrix_set(coef, 2, i, exp(-cimag(modefreqs->data[i])*t1) * cos(creal(modefreqs->data[i])*t1));
373  gsl_matrix_set(coef, 3, i, exp(-cimag(modefreqs->data[i])*t2) * cos(creal(modefreqs->data[i])*t2));
374  gsl_matrix_set(coef, 4, i, exp(-cimag(modefreqs->data[i])*t3) * cos(creal(modefreqs->data[i])*t3));
375  gsl_matrix_set(coef, 5, i, exp(-cimag(modefreqs->data[i])*t4) * cos(creal(modefreqs->data[i])*t4));
376  gsl_matrix_set(coef, 6, i, exp(-cimag(modefreqs->data[i])*t5) * cos(creal(modefreqs->data[i])*t5));
377  gsl_matrix_set(coef, 7, i, exp(-cimag(modefreqs->data[i])*t5) *
378  (-cimag(modefreqs->data[i]) * cos(creal(modefreqs->data[i])*t5)
379  -creal(modefreqs->data[i]) * sin(creal(modefreqs->data[i])*t5)));
380  gsl_matrix_set(coef, 8, i, 0);
381  gsl_matrix_set(coef, 9, i, - creal(modefreqs->data[i]));
382  gsl_matrix_set(coef, 10, i, -exp(-cimag(modefreqs->data[i])*t1) * sin(creal(modefreqs->data[i])*t1));
383  gsl_matrix_set(coef, 11, i, -exp(-cimag(modefreqs->data[i])*t2) * sin(creal(modefreqs->data[i])*t2));
384  gsl_matrix_set(coef, 12, i, -exp(-cimag(modefreqs->data[i])*t3) * sin(creal(modefreqs->data[i])*t3));
385  gsl_matrix_set(coef, 13, i, -exp(-cimag(modefreqs->data[i])*t4) * sin(creal(modefreqs->data[i])*t4));
386  gsl_matrix_set(coef, 14, i, -exp(-cimag(modefreqs->data[i])*t5) * sin(creal(modefreqs->data[i])*t5));
387  gsl_matrix_set(coef, 15, i, exp(-cimag(modefreqs->data[i])*t5) *
388  ( cimag(modefreqs->data[i]) * sin(creal(modefreqs->data[i])*t5)
389  -creal(modefreqs->data[i]) * cos(creal(modefreqs->data[i])*t5)));
390  }
391  for (i = 0; i < 2; ++i)
392  {
393  gsl_vector_set(hderivs, i, inspwave1->data[(i + 1) * inspwave1->vectorLength - 1]);
394  gsl_vector_set(hderivs, i + nmodes, inspwave2->data[(i + 1) * inspwave2->vectorLength - 1]);
395  gsl_vector_set(hderivs, i + 6, inspwave1->data[i * inspwave1->vectorLength]);
396  gsl_vector_set(hderivs, i + 6 + nmodes, inspwave2->data[i * inspwave2->vectorLength]);
397  }
398  gsl_vector_set(hderivs, 2, inspwave1->data[4]);
399  gsl_vector_set(hderivs, 2 + nmodes, inspwave2->data[4]);
400  gsl_vector_set(hderivs, 3, inspwave1->data[3]);
401  gsl_vector_set(hderivs, 3 + nmodes, inspwave2->data[3]);
402  gsl_vector_set(hderivs, 4, inspwave1->data[2]);
403  gsl_vector_set(hderivs, 4 + nmodes, inspwave2->data[2]);
404  gsl_vector_set(hderivs, 5, inspwave1->data[1]);
405  gsl_vector_set(hderivs, 5 + nmodes, inspwave2->data[1]);
406 
407  /* Complete the definition for the rest half of A */
408  for (i = 0; i < nmodes; ++i)
409  {
410  for (k = 0; k < nmodes; ++k)
411  {
412  gsl_matrix_set(coef, i, k + nmodes, - gsl_matrix_get(coef, i + nmodes, k));
413  gsl_matrix_set(coef, i + nmodes, k + nmodes, gsl_matrix_get(coef, i, k));
414  }
415  }
416 
417 #if 0
418  /* print ringdown-matching linear system: coefficient matrix and RHS vector */
419  XLAL_PRINT_INFO("\nRingdown matching matrix:\n");
420  for (i = 0; i < 16; ++i)
421  {
422  for (j = 0; j < 16; ++j)
423  {
424  XLAL_PRINT_INFO("%.12e ",gsl_matrix_get(coef,i,j));
425  }
426  XLAL_PRINT_INFO("\n");
427  }
428  XLAL_PRINT_INFO("RHS: ");
429  for (i = 0; i < 16; ++i)
430  {
431  XLAL_PRINT_INFO("%.12e ",gsl_vector_get(hderivs,i));
432  }
433  XLAL_PRINT_INFO("\n");
434 #endif
435 
436  /* Call gsl LU decomposition to solve the linear system */
437  XLAL_CALLGSL( gslStatus = gsl_linalg_LU_decomp(coef, p, &s) );
438  if ( gslStatus == GSL_SUCCESS )
439  {
440  XLAL_CALLGSL( gslStatus = gsl_linalg_LU_solve(coef, p, hderivs, x) );
441  }
442  if ( gslStatus != GSL_SUCCESS )
443  {
444  gsl_matrix_free(coef);
445  gsl_vector_free(hderivs);
446  gsl_vector_free(x);
447  gsl_permutation_free(p);
448  XLAL_ERROR( XLAL_EFUNC );
449  }
450 
451  /* Putting solution to an XLAL vector */
452  modeamps = XLALCreateREAL8Vector(2 * nmodes);
453 
454  if ( !modeamps )
455  {
456  gsl_matrix_free(coef);
457  gsl_vector_free(hderivs);
458  gsl_vector_free(x);
459  gsl_permutation_free(p);
460  XLAL_ERROR( XLAL_ENOMEM );
461  }
462 
463  for (i = 0; i < nmodes; ++i)
464  {
465  modeamps->data[i] = gsl_vector_get(x, i);
466  modeamps->data[i + nmodes] = gsl_vector_get(x, i + nmodes);
467  }
468 
469  /* Free all gsl linear algebra objects */
470  gsl_matrix_free(coef);
471  gsl_vector_free(hderivs);
472  gsl_vector_free(x);
473  gsl_permutation_free(p);
474 
475  /* Build ring-down waveforms */
476 
477  REAL8 timeOffset = fmod( matchrange->data[1], dt/m) * dt;
478 
479  if (debugSB){
480  /// Print the solution:
481  for (i=0; i<nmodes; i++){
482  printf("RD info: QNM: (re) %.16e, (im) %.16e, Amp %16e \n", creal(modefreqs->data[i]),
483  cimag(modefreqs->data[i]), modeamps->data[i]);
484  }
485 
486  }
487 
488  for (j = 0; j < rdwave1->length; ++j)
489  {
490  tj = j * dt - timeOffset;
491  rdwave1->data[j] = 0;
492  rdwave2->data[j] = 0;
493  for (i = 0; i < nmodes; ++i)
494  {
495  rdwave1->data[j] += exp(- tj * cimag(modefreqs->data[i]))
496  * ( modeamps->data[i] * cos(tj * creal(modefreqs->data[i]))
497  + modeamps->data[i + nmodes] * sin(tj * creal(modefreqs->data[i])) );
498  rdwave2->data[j] += exp(- tj * cimag(modefreqs->data[i]))
499  * (- modeamps->data[i] * sin(tj * creal(modefreqs->data[i]))
500  + modeamps->data[i + nmodes] * cos(tj * creal(modefreqs->data[i])) );
501  }
502  }
503 
504  XLALDestroyREAL8Vector(modeamps);
505  return errcode;
506 }
507 
508 /**
509  * Function which calculates the value of the waveform, plus its
510  * first and second derivatives, for the points which will be required
511  * in the hybrid comb attachment of the ringdown.
512  */
513 static INT4 XLALGenerateHybridWaveDerivatives (
514  REAL8Vector *rwave, /**<< OUTPUT, values of the waveform at comb points */
515  REAL8Vector *dwave, /**<< OUTPUT, 1st deriv of the waveform at comb points */
516  REAL8Vector *ddwave, /**<< OUTPUT, 2nd deriv of the waveform at comb points */
517  REAL8Vector *timeVec, /**<< Vector containing the time */
518  REAL8Vector *wave, /**<< Last part of inspiral waveform */
519  REAL8Vector *matchrange, /**<< Times which determine the size of the comb */
520  REAL8 dt, /**<< Sample time step */
521  REAL8 mass1, /**<< First component mass (in Solar masses) */
522  REAL8 mass2 /**<< Second component mass (in Solar masses) */
523 )
524 {
525 
526  /* XLAL error handling */
527  INT4 errcode = XLAL_SUCCESS;
528 
529  /* For checking GSL return codes */
530  INT4 gslStatus;
531 
532  UINT4 j;
533  UINT4 vecLength;
534  REAL8 m;
535  double *y;
536  double ry, dy, dy2;
537  double rt;
538  double *tlist;
539  gsl_interp_accel *acc;
540  gsl_spline *spline;
541 
542  /* Total mass in geometric units */
543  m = (mass1 + mass2) * LAL_MTSUN_SI;
544 
545  tlist = (double *) LALMalloc(6 * sizeof(double));
546  rt = (matchrange->data[1] - matchrange->data[0]) / 5.;
547  tlist[0] = matchrange->data[0];
548  tlist[1] = tlist[0] + rt;
549  tlist[2] = tlist[1] + rt;
550  tlist[3] = tlist[2] + rt;
551  tlist[4] = tlist[3] + rt;
552  tlist[5] = matchrange->data[1];
553 
554  /* Set the length of the interpolation vectors */
555  vecLength = ( m * matchrange->data[2] / dt ) + 1;
556 
557  /* Getting interpolation and derivatives of the waveform using gsl spline routine */
558  /* Initiate arrays and supporting variables for gsl */
559  y = (double *) LALMalloc(vecLength * sizeof(double));
560 
561  if ( !y )
562  {
563  XLAL_ERROR( XLAL_ENOMEM );
564  }
565 
566  for (j = 0; j < vecLength; ++j)
567  {
568  y[j] = wave->data[j];
569  }
570 
571 
572  XLAL_CALLGSL( acc = (gsl_interp_accel*) gsl_interp_accel_alloc() );
573  XLAL_CALLGSL( spline = (gsl_spline*) gsl_spline_alloc(gsl_interp_cspline, vecLength) );
574  if ( !acc || !spline )
575  {
576  if ( acc ) gsl_interp_accel_free(acc);
577  if ( spline ) gsl_spline_free(spline);
578  LALFree( y );
579  XLAL_ERROR( XLAL_ENOMEM );
580  }
581 
582  /* Gall gsl spline interpolation */
583  gslStatus = gsl_spline_init(spline, timeVec->data, y, vecLength);
584  if ( gslStatus != GSL_SUCCESS )
585  {
586  gsl_spline_free(spline);
587  gsl_interp_accel_free(acc);
588  LALFree( y );
589  XLAL_ERROR( XLAL_EFUNC );
590  }
591 
592  /* Getting first and second order time derivatives from gsl interpolations */
593  for (j = 0; j < 6; ++j)
594  {
595  gslStatus = gsl_spline_eval_e( spline, tlist[j], acc, &ry );
596  if ( gslStatus == GSL_SUCCESS )
597  {
598  gslStatus = gsl_spline_eval_deriv_e(spline, tlist[j], acc, &dy );
599  gslStatus = gsl_spline_eval_deriv2_e(spline, tlist[j], acc, &dy2 );
600  }
601  if (gslStatus != GSL_SUCCESS )
602  {
603  gsl_spline_free(spline);
604  gsl_interp_accel_free(acc);
605  LALFree( y );
606  XLAL_ERROR( XLAL_EFUNC );
607  }
608  rwave->data[j] = (REAL8)(ry);
609  dwave->data[j] = (REAL8)(dy/m);
610  ddwave->data[j] = (REAL8)(dy2/m/m);
611 
612  }
613 
614  /* Free gsl variables */
615  gsl_spline_free(spline);
616  gsl_interp_accel_free(acc);
617  LALFree( tlist );
618  LALFree(y);
619 
620  return errcode;
621 }
622 
623 /**
624  * Translates the major approximant version into an unsigned integer
625  * to minimize if and or statements throughout the function
626  * XLALSimIMREOBHybridAttachRingdownPrec. Translation into an
627  * unsigned integer prevents issues with passing the approximant
628  * parameter to other functions.
629  */
630 static INT4 get_major_approximant( Approximant approximant /**<<The waveform approximant being used */ )
631 {
632  INT4 major_approximant = 0;
633 
634  if(approximant == SEOBNRv1) major_approximant = 1;
635  if(approximant == SEOBNRv2) major_approximant = 2;
636  if(approximant == SEOBNRv3 || approximant == SEOBNRv3_opt) major_approximant = 3;
637 
638  if(major_approximant == 0) XLAL_PRINT_ERROR("In LALSimIMREOBHybridRingdownPrec.c major_approximant not assigned.\n");
639 
640  return major_approximant;
641 }
642 
643 /**
644  * The main workhorse function for performing the ringdown attachment for EOB
645  * models EOBNRv2 and SEOBNRv1. This is the function which gets called by the
646  * code generating the full IMR waveform once generation of the inspiral part
647  * has been completed.
648  * The ringdown is attached using the hybrid comb matching detailed in
649  * The method is describe in Sec. II C of Pan et al. PRD 84, 124052 (2011),
650  * specifically Eqs. 30 - 32.. Further details of the
651  * implementation of the found in the DCC document T1100433.
652  * In SEOBNRv1, the last physical overtone is replace by a pseudoQNM. See
653  * Taracchini et al. PRD 86, 024011 (2012) for details.
654  * STEP 1) Get mass and spin of the final black hole and the complex ringdown frequencies
655  * STEP 2) Based on least-damped-mode decay time, allocate memory for rigndown waveform
656  * STEP 3) Get values and derivatives of inspiral waveforms at matching comb points
657  * STEP 4) Solve QNM coefficients and generate ringdown waveforms
658  * STEP 5) Stitch inspiral and ringdown waveoforms
659  * SEOBNRv2 prescriptions can be found in https://dcc.ligo.org/T1400476
660  * SEONBRv3 prescriptions can be found in XXX
661  */
662 UNUSED static INT4
663 XLALSimIMREOBHybridAttachRingdownPrec(
664  REAL8Vector * signal1, /**<< OUTPUT, Real of inspiral waveform to which we attach ringdown */
665  REAL8Vector * signal2, /**<< OUTPUT, Imag of inspiral waveform to which we attach ringdown */
666  const INT4 l, /**<< Current mode l */
667  const INT4 m, /**<< Current mode m */
668  const REAL8 dt, /**<< Sample time step (in seconds) */
669  const REAL8 mass1, /**<< First component mass (in Solar masses) */
670  const REAL8 mass2, /**<< Second component mass (in Solar masses) */
671  const REAL8 spin1x, /**<<The spin of the first object; only needed for spin waveforms */
672  const REAL8 spin1y, /**<<The spin of the first object; only needed for spin waveforms */
673  const REAL8 spin1z, /**<<The spin of the first object; only needed for spin waveforms */
674  const REAL8 spin2x, /**<<The spin of the second object; only needed for spin waveforms */
675  const REAL8 spin2y, /**<<The spin of the second object; only needed for spin waveforms */
676  const REAL8 spin2z, /**<<The spin of the second object; only needed for spin waveforms */
677  REAL8Vector * timeVec, /**<< Vector containing the time values */
678  REAL8Vector * matchrange, /**<< Time values chosen as points for performing comb matching */
679  Approximant approximant, /**<<The waveform approximant being used */
680  const REAL8 JLN /**<< cosine of the angle between J and LN at the light ring */
681 )
682 {
683  INT4 debugout = 0;
684 
685  COMPLEX16Vector *modefreqs;
686  //COMPLEX16 freq7sav;
687  UINT4 Nrdwave;
688  UINT4 j;
689 
690  UINT4 nmodes;
691  REAL8Vector *rdwave1;
692  REAL8Vector *rdwave2;
693  REAL8Vector *rinspwave;
694  REAL8Vector *dinspwave;
695  REAL8Vector *ddinspwave;
696  REAL8VectorSequence *inspwaves1;
697  REAL8VectorSequence *inspwaves2;
698  REAL8 eta , a, chi;
699  REAL8 NRPeakOmega22 = 0; /* To generate pQNM
700  * frequency */
701  REAL8 sh , kk, kt1, kt2; /* To generate pQNM frequency */
702  REAL8 mTot; /* In geometric units */
703  REAL8 spin1 [3] = {spin1x, spin1y, spin1z};
704  REAL8 spin2 [3] = {spin2x, spin2y, spin2z};
705  REAL8 finalMass, finalSpin;
706  REAL8 chi1 , chi2, theta1, theta2;
707  INT4 mHere = 0;
708 
709  INT4 major_approximant = get_major_approximant(approximant);
710 
711  mTot = (mass1 + mass2) * LAL_MTSUN_SI;
712  eta = mass1 * mass2 / ((mass1 + mass2) * (mass1 + mass2));
713 
714  /*
715  * STEP 1) Get mass and spin of the final black hole and the complex
716  * ringdown frequencies
717  */
718 
719  /* Create memory for the QNM frequencies */
720  nmodes = 8;
721  modefreqs = XLALCreateCOMPLEX16Vector(nmodes);
722  if (!modefreqs) {
723  XLAL_ERROR(XLAL_ENOMEM);
724  }
725  if (XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs, mass1, mass2, spin1, spin2, l, m, nmodes, approximant) == XLAL_FAILURE) {
726  XLALDestroyCOMPLEX16Vector(modefreqs);
727  XLAL_ERROR(XLAL_EFUNC);
728  }
729  if(major_approximant == 3) {
730  if (JLN > 0){
731  mHere = (int)fabs((REAL8) m);
732  }else{
733  mHere = -(int)fabs((REAL8) m);
734  }
735 
736  if (XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs, mass1, mass2, spin1, spin2, l, (int)fabs((REAL8) m), nmodes, approximant) == XLAL_FAILURE) {
737  XLALDestroyCOMPLEX16Vector(modefreqs);
738  XLAL_ERROR(XLAL_EFUNC);
739  }
740  // FIXME
741  if (m < 0 && JLN >= 0) {
742  for (j = 0; j < nmodes; j++) {
743  modefreqs->data[j] = conjl(-1.0 * modefreqs->data[j]);
744  }
745  }
746  if (m > 0 && JLN < 0) {
747  for (j = 0; j < nmodes; j++) {
748  modefreqs->data[j] = conjl(-1.0 * modefreqs->data[j]);
749  }
750  }
751  if (m == 0) {
752  for (j = 5; j < nmodes; j++) {
753  modefreqs->data[j] = conjl(-1.0 * modefreqs->data[j - 5]);
754  }
755  }
756  }
757 
758  /*
759  * Call XLALSimIMREOBFinalMassSpinPrec() to get mass and spin of the
760  * final black hole
761  */
762  if (XLALSimIMREOBFinalMassSpinPrec(&finalMass, &finalSpin, mass1, mass2, spin1, spin2, approximant) == XLAL_FAILURE) {
763  XLAL_ERROR(XLAL_EFUNC);
764  }
765  if (debugout) {
766  XLAL_PRINT_INFO("Modes:\n");
767  for (j = 0; j < nmodes; j++) {
768  XLAL_PRINT_INFO("%d, %e %e \n", j, creal(modefreqs->data[j]), cimag(modefreqs->data[j]) );
769  }
770  XLAL_PRINT_INFO("spin1 [%e, %e, %e], spin2 [%e, %e, %e] \n", spin1x, spin1y, spin1z, spin2x, spin2y, spin2z);
771  XLAL_PRINT_INFO("(Mf, af) = (%3.10f, %3.10f)\n", finalMass, finalSpin);
772  }
773  if (major_approximant == 1) {
774  /* Replace the last QNM with pQNM */
775  /* We assume aligned/antialigned spins here */
776  a = (spin1[2] + spin2[2]) / 2. * (1.0 - 2.0 * eta) + (spin1[2] - spin2[2]) / 2. * (mass1 - mass2) / (mass1 + mass2);
777  NRPeakOmega22 = XLALSimIMREOBGetNRSpinPeakOmega(l, m, eta, a) / mTot;
778  /*
779  * XLAL_PRINT_INFO("a and NRomega in QNM freq: %.16e %.16e %.16e %.16e
780  * %.16e\n",spin1[2],spin2[2],
781  * mTot/LAL_MTSUN_SI,a,NRPeakOmega22*mTot);
782  */
783  modefreqs->data[7] = (NRPeakOmega22 / finalMass + creal(modefreqs->data[0])) / 2.;
784  modefreqs->data[7] += I * 10. / 3. * cimag(modefreqs->data[0]);
785  }
786  if (major_approximant == 2)
787  //See pages 6 to 12 of the dcc document T1400476 - v3 for expressions
788  // in this block.
789  {
790  /* Replace the last two QNMs with pQNMs */
791  /* We assume aligned/antialigned spins here */
792  /*
793  * Definitions of a, chi and NRPeakOmega22, where the
794  * last one is an approximation of \phi'[tmatch] in
795  * T1400476-v3.
796  */
797  a = (spin1[2] + spin2[2]) / 2. * (1.0 - 2.0 * eta) + (spin1[2] - spin2[2]) / 2. * (mass1 - mass2) / (mass1 + mass2);
798  NRPeakOmega22 = XLALSimIMREOBGetNRSpinPeakOmegav2(l, m, eta, a) / mTot;
799 
800  /* Define chi */
801  chi = (spin1[2] + spin2[2]) / 2. + (spin1[2] - spin2[2]) / 2. * ((mass1 - mass2) / (mass1 + mass2)) / (1. - 2. * eta);
802 
803  /*
804  * For extreme chi (>= 0.8), there are scale factors
805  * in both complex pseudo-QNM frequencies. kk, kt1,
806  * kt2 describe those factors.
807  */
808  //Below definitions of kk, kt1 and kt2 are likely obsolete
809  kk = kt1 = kt2 = 1.;
810  if (chi >= 0.8) {
811  kk = 0.7 + 0.3 * exp(100. * (eta - 0.25));
812  kt1 = 0.5 * sqrt(1. + 800.0 * eta * eta / 3.0) - 0.125;
813  kt2 = 0.5 * pow(1. + 0.5 * eta * sqrt(eta) / 0.0225, 2. / 3.) - 0.2;
814  }
815  //Above definitions of kk, kt1 and kt2 are likely obsolete
816  /*
817  * XLAL_PRINT_INFO("a, chi and NRomega in QNM freq: %.16e
818  * %.16e %.16e %.16e %.16e %.16e\n",
819  * spin1[2],spin2[2],mTot/LAL_MTSUN_SI,a,chi,NRPeakOme
820  * ga22*mTot);
821  */
822  sh = 0.;
823  //freq7sav = modefreqs->data[7];
824 
825  /*
826  * Cases 1, 2 and 3 in T1400476-v3. Note that the
827  * difference between the chi1=chi2=0 case and the
828  * chi<0.7 cases is only in Dtcomb, which is not
829  * specified or used in this file.
830  */
831  modefreqs->data[7] = (2. / 3. * NRPeakOmega22 / finalMass) + (1. / 3. * creal(modefreqs->data[0]));
832  modefreqs->data[7] += I * 3.5 / 0.9 * cimag(modefreqs->data[0]);
833  modefreqs->data[6] = (3. / 4. * NRPeakOmega22 / finalMass) + (1. / 4. * creal(modefreqs->data[0]));
834  modefreqs->data[6] += I * 3.5 * cimag(modefreqs->data[0]);
835 
836  /*
837  * For extreme chi (>= 0.8), the ringdown attachment
838  * should end slightly before the value given passed
839  * to this function. sh gives exactly how long
840  * before.
841  */
842  if (chi >= 0.7 && chi < 0.8) {
843  sh = -9. * (eta - 0.25);
844  }
845  if ((eta > 30. / 31. / 31. && eta <= 10. / 121. && chi >= 0.8) || (eta <= 30. / 31. / 31. && chi >= 0.8 && chi < 0.9)) {
846  //This is case 4 in T1400476 - v3
847  sh = -9. * (eta - 0.25) * (1. + 2. * exp(-(chi - 0.85) * (chi - 0.85) / 0.05 / 0.05)) * (1. + 1. / (1. + exp((eta - 0.01) / 0.001)));
848  kk = 0.7 + 0.3 * exp(100. * (eta - 0.25));
849  kt1 = 0.5 * sqrt(1. + 800.0 * eta * eta / 3.0) - 0.125;
850  kt2 = 0.5 * pow(1. + 0.5 * eta * sqrt(eta) / 0.0225, 2. / 3.) - 0.2;
851  modefreqs->data[4] = 0.4 * (1. + kk) * creal(modefreqs->data[6])
852  + I * cimag(modefreqs->data[6]) / (2.5 * kt2 * exp(-(eta - 0.005) / 0.03));
853  modefreqs->data[5] = 0.4 * (1. + kk) * creal(modefreqs->data[7])
854  + I * cimag(modefreqs->data[7]) / (1.5 * kt1 * exp(-(eta - 0.005) / 0.03));
855  modefreqs->data[6] = kk * creal(modefreqs->data[6]) + I * cimag(modefreqs->data[6]) / kt2;
856  modefreqs->data[7] = kk * creal(modefreqs->data[7]) + I * cimag(modefreqs->data[7]) / kt1;
857  }
858  if (eta < 30. / 31. / 31. && chi >= 0.9) {
859  //This is case 5 in T1400476 - v3
860  sh = 0.55 - 9. * (eta - 0.25) * (1. + 2. * exp(-(chi - 0.85) * (chi - 0.85) / 0.05 / 0.05)) * (1. + 1. / (1. + exp((eta - 0.01) / 0.001)));
861  kk = 0.7 + 0.3 * exp(100. * (eta - 0.25));
862  kt1 = 0.5 * sqrt(1. + 800.0 * eta * eta / 3.0) - 0.125;
863  kt2 = 0.5 * pow(1. + 0.5 * eta * sqrt(eta) / 0.0225, 2. / 3.) - 0.2;
864  modefreqs->data[4] = 1.1 * 0.4 * (1. + kk) * creal(modefreqs->data[6])
865  + I * cimag(modefreqs->data[6]) / (1.05 * 2.5 * kt2 * exp(-(eta - 0.005) / 0.03));
866  modefreqs->data[5] = 0.4 * (1. + kk) * creal(modefreqs->data[7])
867  + I * cimag(modefreqs->data[7]) / (1.05 * 1.5 * kt1 * exp(-(eta - 0.005) / 0.03));
868  modefreqs->data[6] = kk * creal(modefreqs->data[6]) + I * cimag(modefreqs->data[6]) / kt2;
869  modefreqs->data[7] = kk * creal(modefreqs->data[7]) + I * cimag(modefreqs->data[7]) / kt1;
870  }
871  if (eta > 10. / 121. && chi >= 0.8) {
872  //This is case 6 in T1400476 - v3
873  sh = 1. - 9. * (eta - 0.25) * (1. + 2. * exp(-(chi - 0.85) * (chi - 0.85) / 0.05 / 0.05)) * (1. + 1. / (1. + exp((eta - 0.01) / 0.001)));
874  kk = 0.7 + 0.3 * exp(100. * (eta - 0.25));
875  kt1 = 0.45 * sqrt(1. + 200.0 * eta * eta / 3.0) - 0.125;
876  kt2 = 0.5 * pow(1. + 0.5 * eta * sqrt(eta) / 0.0225, 2. / 3.) - 0.2;
877  modefreqs->data[6] = kk * creal(modefreqs->data[6]) + I * cimag(modefreqs->data[6]) / 0.95 / kt2;
878  modefreqs->data[7] = kk * creal(modefreqs->data[7]) + I * cimag(modefreqs->data[7]) / kt1;
879  }
880  //The last line of T1400476 - v3
881  matchrange->data[0] -= sh;
882  matchrange->data[1] -= sh;
883  /*
884  * modefreqs->data[7] = 0.38068371/mTot +
885  * I/1.4677128/mTot; modefreqs->data[6] =
886  * 0.37007703/mTot + I/1.3359367/mTot;
887  * modefreqs->data[5] = 0.36980703/mTot +
888  * I/1.7791212/mTot; modefreqs->data[4] =
889  * 0.3595034/mTot + I/2.6989764/mTot; XLAL_PRINT_INFO("sh =
890  * %f\n",sh); XLAL_PRINT_INFO("PeakOmega = %f, mTot =
891  * %f\n",NRPeakOmega22,mTot); XLAL_PRINT_INFO("w0 = %f, t0 =
892  * %f\n",creal(modefreqs->data[0])*mTot,
893  * 1./cimag(modefreqs->data[0])/mTot); XLAL_PRINT_INFO("w1 =
894  * %f, t1 = %f\n",creal(modefreqs->data[6])*mTot,
895  * 1./cimag(modefreqs->data[6])/mTot); XLAL_PRINT_INFO("w2 =
896  * %f, t2 = %f\n",creal(modefreqs->data[7])*mTot,
897  * 1./cimag(modefreqs->data[7])/mTot); XLAL_PRINT_INFO("w3 =
898  * %f, t3 = %f\n",creal(modefreqs->data[4])*mTot,
899  * 1./cimag(modefreqs->data[4])/mTot); XLAL_PRINT_INFO("w4 =
900  * %f, t4 = %f\n",creal(modefreqs->data[5])*mTot,
901  * 1./cimag(modefreqs->data[5])/mTot);
902  */
903  }
904  if (major_approximant == 3) {
905  REAL8 kappa_thr = 0.175;
906  //REAL8 eJL_thr = 7.5e-3;
907  REAL8 eJL_thr = 5.0e-2;
908  //REAL8 eJL_thr = 1.0e-2;
909  chi1 = sqrt(spin1[0] * spin1[0] + spin1[1] * spin1[1] + spin1[2] * spin1[2]);
910  chi2 = sqrt(spin2[0] * spin2[0] + spin2[1] * spin2[1] + spin2[2] * spin2[2]);
911  if (chi1 < 1.0e-15) {
912  theta1 = 0.;
913  } else {
914  theta1 = acos(spin1[2] / chi1);
915  }
916  if (chi2 < 1.0e-15) {
917  theta2 = 0.;
918  } else {
919  theta2 = acos(spin2[2] / chi2);
920  }
921  chi1 = chi1 * cos(theta1);
922  chi2 = chi2 * cos(theta2);
923 
924  /*Compute wf freq at matching point*/
925  double *y;
926  double ry, dy;
927  gsl_interp_accel *acc;
928  gsl_spline *spline;
929  //UINT4 vecLength = ((mass1 + mass2) * LAL_MTSUN_SI * matchrange->data[2] / dt) + 1;
930  //y = (double *)LALMalloc(vecLength * sizeof(double));
931  y = (double *)LALMalloc(timeVec->length * sizeof(double));
932  double hRe, dhRe, hIm, dhIm;
933 
934  if (!y) {
935  XLAL_ERROR(XLAL_ENOMEM);
936  }
937  //FILE *out1 = fopen( "Andrea1.dat","w");
938  for (j = 0; j < timeVec->length; ++j) {
939  y[j] = signal1->data[j];
940  //fprintf(out1, "%.16e %.16e\n", timeVec->data[j], y[j]);
941  //XLAL_PRINT_INFO("%.16e %.16e\n", timeVec->data[j], y[j]);
942  }
943  //fclose(out1);
944  //exit(0);
945 
946  XLAL_CALLGSL(acc = (gsl_interp_accel *) gsl_interp_accel_alloc());
947  //XLAL_CALLGSL(spline = (gsl_spline *) gsl_spline_alloc(gsl_interp_cspline, vecLength));
948  XLAL_CALLGSL(spline = (gsl_spline *) gsl_spline_alloc(gsl_interp_cspline, timeVec->length));
949  if (!acc || !spline) {
950  if (acc)
951  gsl_interp_accel_free(acc);
952  if (spline)
953  gsl_spline_free(spline);
954  LALFree(y);
955  XLAL_ERROR(XLAL_ENOMEM);
956  }
957 
958  //XLAL_PRINT_INFO("here1.... last point %f, match_t %f \n", timeVec->data[timeVec->length-1], matchrange->data[1]);
959  //INT4 gslStatus = gsl_spline_init(spline, timeVec->data, y, vecLength);
960  INT4 gslStatus = gsl_spline_init(spline, timeVec->data, y, timeVec->length);
961  if (gslStatus != GSL_SUCCESS) {
962  gsl_spline_free(spline);
963  gsl_interp_accel_free(acc);
964  LALFree(y);
965  XLAL_ERROR(XLAL_EFUNC);
966  }
967  gslStatus = gsl_spline_eval_e(spline, matchrange->data[1], acc, &ry);
968  //int ind_stas = (int) matchrange->data[1]*(((mass1 + mass2) * LAL_MTSUN_SI / dt));
969  //XLAL_PRINT_INFO("here2.... %f, %f, %f, %f \n", ry, timeVec->data[ind_stas], matchrange->data[1], y[ind_stas]);
970  if (gslStatus == GSL_SUCCESS) {
971  gslStatus = gsl_spline_eval_deriv_e(spline, matchrange->data[1], acc, &dy);
972  //XLAL_PRINT_INFO("here2.... ");
973  }
974  if (gslStatus != GSL_SUCCESS) {
975  gsl_spline_free(spline);
976  gsl_interp_accel_free(acc);
977  LALFree(y);
978  XLAL_ERROR(XLAL_EFUNC);
979  }
980  hRe = (REAL8) (ry);
981  dhRe = (REAL8) (dy);
982 
983  //XLAL_PRINT_INFO("here hRe = %f, dhRe = %f \n", hRe, dhRe);
984 
985  //FILE *out2 = fopen( "Andrea2.dat","w");
986  //for (j = 0; j < vecLength; ++j) {
987  for (j = 0; j < timeVec->length; ++j) {
988  y[j] = signal2->data[j];
989  //fprintf(out2, "%.16e %.16e\n", timeVec->data[j], y[j]);
990  }
991  //fclose(out2);
992  //XLAL_CALLGSL(acc = (gsl_interp_accel *) gsl_interp_accel_alloc());
993  //XLAL_CALLGSL(spline = (gsl_spline *) gsl_spline_alloc(gsl_interp_cspline, vecLength));
994  //XLAL_CALLGSL(spline = (gsl_spline *) gsl_spline_alloc(gsl_interp_cspline, signal2->length));
995  if (!acc || !spline) {
996  if (acc)
997  gsl_interp_accel_free(acc);
998  if (spline)
999  gsl_spline_free(spline);
1000  LALFree(y);
1001  XLAL_ERROR(XLAL_ENOMEM);
1002  }
1003  //gslStatus = gsl_spline_init(spline, timeVec->data, y, vecLength);
1004  gslStatus = gsl_spline_init(spline, timeVec->data, y, timeVec->length);
1005  if (gslStatus != GSL_SUCCESS) {
1006  gsl_spline_free(spline);
1007  gsl_interp_accel_free(acc);
1008  LALFree(y);
1009  XLAL_ERROR(XLAL_EFUNC);
1010  }
1011  gslStatus = gsl_spline_eval_e(spline, matchrange->data[1], acc, &ry);
1012  if (gslStatus == GSL_SUCCESS) {
1013  gslStatus = gsl_spline_eval_deriv_e(spline, matchrange->data[1], acc, &dy);
1014  }
1015  if (gslStatus != GSL_SUCCESS) {
1016  gsl_spline_free(spline);
1017  gsl_interp_accel_free(acc);
1018  LALFree(y);
1019  XLAL_ERROR(XLAL_EFUNC);
1020  }
1021  hIm = (REAL8) (ry);
1022  dhIm = (REAL8) (dy);
1023 
1024  //XLAL_PRINT_INFO("Stas, check hRe = %f, dhRe = %f, hIm = %f, dhIm = %f \n", hRe, dhRe, hIm, dhIm);
1025  double hNorm2 = hRe*hRe + hIm*hIm;
1026  double omegaWavePeak = creal(modefreqs->data[0]);
1027  if (hNorm2 != 0.0){
1028  omegaWavePeak = ((-dhRe*hIm + dhIm*hRe)/hNorm2) / mTot;
1029  }
1030  else{
1031  if(debugout) {XLAL_PRINT_INFO("hNorm=0.0 for (l,m)=(%d,%d)\n",l,m);}
1032 // XLAL_ERROR(XLAL_EFAILED);
1033  }
1034 
1035 
1036  a = (chi1 + chi2) / 2. * (1.0 - 2.0 * eta) + (chi1 - chi2) / 2. * (mass1 - mass2) / (mass1 + mass2);
1037  NRPeakOmega22 = fabs(omegaWavePeak);
1038 // NRPeakOmega22 = XLALSimIMREOBGetNRSpinPeakOmegav2(l, m, eta, a) / mTot;
1039  gsl_spline_free(spline);
1040  gsl_interp_accel_free(acc);
1041  LALFree(y);
1042  // FIXME
1043  //NRPeakOmega22 = 0.3;
1044  //NRPeakOmega22 = XLALSimIMREOBGetNRSpinPeakOmegav2(l, m, eta, a) / mTot;
1045 // NRPeakOmega22 = omegaWavePeak/mTot;
1046 // XLAL_PRINT_INFO("(hRe, hIm, dhRe, dhIm)=(%.16e, %.16e, %.16e, %.16e)\n", hRe, hIm, dhRe, dhIm);
1047 // XLAL_PRINT_INFO("hNorm %.16e, tmatch %.16e\n",sqrt(hNorm2), matchrange->data[1]);
1048 // XLAL_PRINT_INFO("NRPeakOmega22 %.16e, omegaWavePeak %.16e\n",NRPeakOmega22,omegaWavePeak);
1049  chi = (chi1 + chi2) / 2. + (chi1 - chi2) / 2. * ((mass1 - mass2) / (mass1 + mass2)) / (1. - 2. * eta);
1050 
1051  sh = 0.;
1052  if (chi >= 0.7 && chi < 0.8) {
1053  sh = -9. * (eta - 0.25);
1054  }
1055  if ((eta > 30. / 31. / 31. && eta <= 10. / 121. && chi >= 0.8) || (eta <= 30. / 31. / 31. && chi >= 0.8 && chi < 0.9)) {
1056  //This is case 4 in T1400476 - v3
1057  sh = -9. * (eta - 0.25) * (1. + 2. * exp(-(chi - 0.85) * (chi - 0.85) / 0.05 / 0.05)) * (1. + 1. / (1. + exp((eta - 0.01) / 0.001)));
1058  }
1059  if (eta < 30. / 31. / 31. && chi >= 0.9) {
1060  //This is case 5 in T1400476 - v3
1061  sh = 0.55 - 9. * (eta - 0.25) * (1. + 2. * exp(-(chi - 0.85) * (chi - 0.85) / 0.05 / 0.05)) * (1. + 1. / (1. + exp((eta - 0.01) / 0.001)));
1062  }
1063  if (eta > 10. / 121. && chi >= 0.8) {
1064  //This is case 6 in T1400476 - v3
1065  sh = 1. - 9. * (eta - 0.25) * (1. + 2. * exp(-(chi - 0.85) * (chi - 0.85) / 0.05 / 0.05)) * (1. + 1. / (1. + exp((eta - 0.01) / 0.001)));
1066  }
1067  if ((int)fabs((REAL8) m) == 2 && l == 2) {
1068  /* if (m == 2 && l ==2){ */
1069 
1070  /*
1071  * For extreme chi (>= 0.8), there are scale factors
1072  * in both complex pseudo-QNM frequencies. kk, kt1,
1073  * kt2 describe those factors.
1074  */
1075  /*
1076  * XLAL_PRINT_INFO("Stas: a, chi and NRomega in QNM freq:
1077  * %.16e %.16e %.16e %.16e %.16e %.16e\n",
1078  * spin1[2],spin2[2],mTot/LAL_MTSUN_SI,a,chi,NRPeakOme
1079  * ga22*mTot);
1080  */
1081  kk = kt1 = kt2 = 1.;
1082  if (chi >= 0.8) {
1083  kk = 0.7 + 0.3 * exp(100. * (eta - 0.25));
1084  kt1 = -0.125 + sqrt(1. + 200. * pow(eta, 2) / 3.) / 2.;
1085  kt2 = -0.2 + pow(1. + 200. * pow(eta, 3. / 2.) / 9., 2. / 3.) / 2.;
1086  }
1087  //Computing pQNMs
1088  // FIXME
1089  //if (m < 0) {
1090  //if (m > 0) {
1091  //printf("Stas: J.L = %f, J.L*m = %f Re(mode) = %f \n", JLN, JLN*m, creal(modefreqs->data[0]));
1092  if (JLN*m < 0.0){
1093  modefreqs->data[6] = (-3. / 4. * NRPeakOmega22 / finalMass) + (1. / 4. * creal(modefreqs->data[0]));
1094  modefreqs->data[7] = (-2. / 3. * NRPeakOmega22 / finalMass) + (1. / 3. * creal(modefreqs->data[0]));
1095  } else {
1096  modefreqs->data[6] = (3. / 4. * NRPeakOmega22 / finalMass) + (1. / 4. * creal(modefreqs->data[0]));
1097  modefreqs->data[7] = (2. / 3. * NRPeakOmega22 / finalMass) + (1. / 3. * creal(modefreqs->data[0]));
1098  //XLAL_PRINT_INFO("Stas , here m = %d modes: 0 = %4.10f, 6 = %4.10f, 7 = %4.10f \n", m, creal(modefreqs->data[0])*mTot, creal(modefreqs->data[6])*mTot, creal(modefreqs->data[7])*mTot);
1099  }
1100 
1101  modefreqs->data[7] += I * 3.5 / 0.9 * cimag(modefreqs->data[0]);
1102  modefreqs->data[6] += I * 3.5 * cimag(modefreqs->data[0]);
1103 
1104  if ((eta > 30. / 31. / 31. && eta <= 10. / 121. && chi >= 0.8) || (eta <= 30. / 31. / 31. && chi >= 0.8 && chi < 0.9)) {
1105  //This is case 4 in T1400476 - v3
1106  if (debugout)
1107  XLAL_PRINT_INFO("Stas, this is case4 for pQNM eta = %f, chi = %f \n", eta, chi);
1108  sh = -9. * (eta - 0.25) * (1. + 2. * exp(-(chi - 0.85) * (chi - 0.85) / 0.05 / 0.05)) * (1. + 1. / (1. + exp((eta - 0.01) / 0.001)));
1109  kk = 0.7 + 0.3 * exp(100. * (eta - 0.25));
1110  kt1 = 0.5 * sqrt(1. + 800.0 * eta * eta / 3.0) - 0.125;
1111  kt2 = 0.5 * pow(1. + 0.5 * eta * sqrt(eta) / 0.0225, 2. / 3.) - 0.2;
1112 
1113  modefreqs->data[4] = 0.4 * (1. + kk) * creal(modefreqs->data[6])
1114  + I * cimag(modefreqs->data[6]) / (2.5 * kt2 * exp(-(eta - 0.005) / 0.03));
1115  modefreqs->data[5] = 0.4 * (1. + kk) * creal(modefreqs->data[7])
1116  + I * cimag(modefreqs->data[7]) / (1.5 * kt1 * exp(-(eta - 0.005) / 0.03));
1117  modefreqs->data[6] = kk * creal(modefreqs->data[6]) + I * cimag(modefreqs->data[6]) / kt2;
1118  modefreqs->data[7] = kk * creal(modefreqs->data[7]) + I * cimag(modefreqs->data[7]) / kt1;
1119  }
1120  if (eta < 30. / 31. / 31. && chi >= 0.9) {
1121  //This is case 5 in T1400476 - v3
1122  if (debugout)
1123  XLAL_PRINT_INFO("Stas, this is case5 for pQNM eta = %f, chi = %f \n", eta, chi);
1124  sh = 0.55 - 9. * (eta - 0.25) * (1. + 2. * exp(-(chi - 0.85) * (chi - 0.85) / 0.05 / 0.05)) * (1. + 1. / (1. + exp((eta - 0.01) / 0.001)));
1125  kk = 0.7 + 0.3 * exp(100. * (eta - 0.25));
1126  kt1 = 0.5 * sqrt(1. + 800.0 * eta * eta / 3.0) - 0.125;
1127  kt2 = 0.5 * pow(1. + 0.5 * eta * sqrt(eta) / 0.0225, 2. / 3.) - 0.2;
1128  modefreqs->data[4] = 1.1 * 0.4 * (1. + kk) * creal(modefreqs->data[6])
1129  + I * cimag(modefreqs->data[6]) / (1.05 * 2.5 * kt2 * exp(-(eta - 0.005) / 0.03));
1130  modefreqs->data[5] = 0.4 * (1. + kk) * creal(modefreqs->data[7])
1131  + I * cimag(modefreqs->data[7]) / (1.05 * 1.5 * kt1 * exp(-(eta - 0.005) / 0.03));
1132  modefreqs->data[6] = kk * creal(modefreqs->data[6]) + I * cimag(modefreqs->data[6]) / kt2;
1133  modefreqs->data[7] = kk * creal(modefreqs->data[7]) + I * cimag(modefreqs->data[7]) / kt1;
1134  }
1135  if (eta > 10. / 121. && chi >= 0.8) {
1136  //This is case 6 in T1400476 - v3
1137  if (debugout)
1138  XLAL_PRINT_INFO("Stas, this is case6 for pQNM eta = %f, chi = %f \n", eta, chi);
1139  sh = 1. - 9. * (eta - 0.25) * (1. + 2. * exp(-(chi - 0.85) * (chi - 0.85) / 0.05 / 0.05)) * (1. + 1. / (1. + exp((eta - 0.01) / 0.001)));
1140  kk = 0.7 + 0.3 * exp(100. * (eta - 0.25));
1141  kt1 = 0.45 * sqrt(1. + 200.0 * eta * eta / 3.0) - 0.125;
1142  kt2 = 0.5 * pow(1. + 0.5 * eta * sqrt(eta) / 0.0225, 2. / 3.) - 0.2;
1143  modefreqs->data[6] = kk * creal(modefreqs->data[6]) + I * cimag(modefreqs->data[6]) / 0.95 / kt2;
1144  modefreqs->data[7] = kk * creal(modefreqs->data[7]) + I * cimag(modefreqs->data[7]) / kt1;
1145  }
1146  if (debugout)
1147  XLAL_PRINT_INFO("Stas, default (if nothing specified above) for pQNM eta = %f, chi = %f chi1 = %f \n", eta, chi, chi1);
1148 
1149  //FIXME It is not compatible with v2!!!!
1150 // if (chi1 >= 0.96){// && eta < 10.0/11./11.){
1151 // // modefreqs->data[7] += I * 3.5 / 0.9 * cimag(modefreqs->data[0]) - I * cimag(modefreqs->data[7]);
1152 // if (debugout){
1153 // XLAL_PRINT_INFO("Stas, the decay time will be modified for pQNM chi1 = %f, from %f to %f \n", chi1, 1.0/cimag(modefreqs->data[7])/mTot, 1.0/(5.0 * cimag(modefreqs->data[0]))/mTot );
1154 // }
1155 // modefreqs->data[7] += I * 5.0 * cimag(modefreqs->data[0]) - I * cimag(modefreqs->data[7]);
1156 // }
1157  //The last line of T1400476 - v3
1158 // matchrange->data[0] -= sh;
1159 // matchrange->data[1] -= sh;
1160 
1161  /*
1162  * This is a v1 of pQNM in RD attachment if (m<0){
1163  * modefreqs->data[7] = (-NRPeakOmega22/finalMass +
1164  * creal(modefreqs->data[0])) / 2.; } else {
1165  * modefreqs->data[7] = (NRPeakOmega22/finalMass +
1166  * creal(modefreqs->data[0])) / 2.; }
1167  * modefreqs->data[7] += I * 10./3. *
1168  * cimag(modefreqs->data[0]);
1169  */
1170 // if(m==-2 || m==-2) {
1171 // modefreqs->data[6] = 0.495264/mTot + I/mTot/3.716525;
1172 // modefreqs->data[7] = 0.514602/mTot + I/mTot/3.344873;
1173 // }
1174 // else if (m==-2){
1175 // modefreqs->data[6] = -0.495264/mTot + I/mTot/3.716525;
1176 // modefreqs->data[7] = -0.514602/mTot + I/mTot/3.344873;
1177 // }
1178 
1179 // XLAL_PRINT_INFO("finalSpin = %f\n",finalSpin);
1180 // XLAL_PRINT_INFO("finalMass = %f\n",finalMass);
1181 // XLAL_PRINT_INFO("PeakOmega = %f\n",NRPeakOmega22*mTot);
1182 // XLAL_PRINT_INFO("w0 = %f, t0 = %f\n",creal(modefreqs->data[0])*mTot, 1./cimag(modefreqs->data[0])/mTot);
1183 // XLAL_PRINT_INFO("w1 = %f, t1 = %f\n",creal(modefreqs->data[1])*mTot, 1./cimag(modefreqs->data[1])/mTot);
1184 // XLAL_PRINT_INFO("w2 = %f, t2 = %f\n",creal(modefreqs->data[2])*mTot, 1./cimag(modefreqs->data[2])/mTot);
1185 // XLAL_PRINT_INFO("w3 = %f, t3 = %f\n",creal(modefreqs->data[3])*mTot, 1./cimag(modefreqs->data[3])/mTot);
1186 // XLAL_PRINT_INFO("w4 = %f, t4 = %f\n",creal(modefreqs->data[4])*mTot, 1./cimag(modefreqs->data[4])/mTot);
1187 // XLAL_PRINT_INFO("w5 = %f, t5 = %f\n",creal(modefreqs->data[5])*mTot, 1./cimag(modefreqs->data[5])/mTot);
1188 // XLAL_PRINT_INFO("w6 = %f, t6 = %f\n",creal(modefreqs->data[6])*mTot, 1./cimag(modefreqs->data[6])/mTot);
1189 // XLAL_PRINT_INFO("w7 = %f, t7 = %f\n",creal(modefreqs->data[7])*mTot, 1./cimag(modefreqs->data[7])/mTot);
1190  // FIXME
1191  // {{{
1192 
1193  if (debugout) {
1194  XLAL_PRINT_INFO("Stas: JLN*eta = %f, \n", JLN*eta);
1195  XLAL_PRINT_INFO("NRPeakOmega22 = %3.10f, %3.10f, %3.10f, %3.10f\n", NRPeakOmega22 * mTot /finalMass, NRPeakOmega22 / finalMass, finalMass, mTot);
1196  for (j = 0; j < nmodes; j++) {
1197  XLAL_PRINT_INFO("QNM frequencies: %d %d %d %3.10f %3.10f\n", l, m, j, creal(modefreqs->data[j]) * mTot, 1. / cimag(modefreqs->data[j]) / mTot);
1198  }
1199  }
1200 
1201  /*if (JLN*m < 0.0){
1202  modefreqs->data[5] = (-3. / 4. * NRPeakOmega22 / finalMass) + (1. / 4. * creal(modefreqs->data[0]));
1203  modefreqs->data[4] = (-2. / 3. * NRPeakOmega22 / finalMass) + (1. / 3. * creal(modefreqs->data[0]));
1204  } else {
1205  modefreqs->data[5] = (3. / 4. * NRPeakOmega22 / finalMass) + (1. / 4. * creal(modefreqs->data[0]));
1206  modefreqs->data[4] = (2. / 3. * NRPeakOmega22 / finalMass) + (1. / 3. * creal(modefreqs->data[0]));
1207  //XLAL_PRINT_INFO("Stas , here m = %d modes: 0 = %4.10f, 6 = %4.10f, 7 = %4.10f \n", m, creal(modefreqs->data[0])*mTot, creal(modefreqs->data[6])*mTot, creal(modefreqs->data[7])*mTot);
1208  }
1209  modefreqs->data[5] += I * 3.5 / 0.9 * cimag(modefreqs->data[0]);
1210  modefreqs->data[4] += I * 3.5 * cimag(modefreqs->data[0]);*/
1211 
1212  COMPLEX16Vector *modefreqs_xtr;
1213  modefreqs_xtr = XLALCreateCOMPLEX16Vector(nmodes);
1214  if (!modefreqs_xtr) {
1215  XLAL_ERROR(XLAL_ENOMEM);
1216  }
1217  //modefreqs->data[4] = ;
1218  //modefreqs->data[5] = 0.0;
1219  //modefreqs->data[6] = 0.0;
1220  //modefreqs->data[7] = 0.0;
1221  //XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, 1, nmodes, approximant);
1222  //if (m == -2){
1223  // XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, -1, nmodes, approximant);
1224  //}
1225  //modefreqs->data[4] = modefreqs_xtr->data[0];
1226  //modefreqs->data[5] = modefreqs_xtr->data[1];
1227 
1228  /*
1229  XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, -2, nmodes, approximant);
1230  if ((JLN > 0.0 && JLN < 0.98) || (eta*JLN < eJL_thr && eta*JLN>0.0)){
1231  modefreqs->data[5] = modefreqs_xtr->data[0];
1232  if (JLN < kappa_thr || eta*JLN < eJL_thr){
1233  modefreqs->data[6] = modefreqs_xtr->data[1];
1234  }
1235  if (m == -2){
1236  modefreqs->data[5] = conjl(-1.0 * modefreqs_xtr->data[0]);
1237  if (JLN < kappa_thr || eta*JLN < eJL_thr){
1238  modefreqs->data[6] = conjl(-1.0 * modefreqs_xtr->data[1]);
1239  }
1240  }
1241  }
1242  */
1243  if ((JLN > 0.0 && JLN < 0.98) || (eta*JLN < eJL_thr && eta*JLN>0.0)){
1244 
1245  XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, -2, nmodes, approximant);
1246  modefreqs->data[5] = modefreqs_xtr->data[0];
1247  if (m == -2){
1248  modefreqs->data[5] = conjl(-1.0 * modefreqs_xtr->data[0]);
1249  }
1250  }
1251 
1252  if ((JLN < 0.0 && JLN > -0.98) || (eta*JLN > -eJL_thr && eta*JLN<0.0) ){
1253  spin1[0] *= -1;
1254  spin1[1] *= -1;
1255  spin1[2] *= -1;
1256  spin2[0] *= -1;
1257  spin2[1] *= -1;
1258  spin2[2] *= -1;
1259  XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, -mHere, nmodes, approximant);
1260  spin1[0] *= -1;
1261  spin1[1] *= -1;
1262  spin1[2] *= -1;
1263  spin2[0] *= -1;
1264  spin2[1] *= -1;
1265  spin2[2] *= -1;
1266  modefreqs->data[5] = modefreqs_xtr->data[0];
1267 // if (JLN > -1.0*kappa_thr || fabs(eta*JLN) < eJL_thr){
1268 // modefreqs->data[6] = modefreqs_xtr->data[1];
1269 // }
1270  if (m == -2){
1271  modefreqs->data[5] = conjl(-1.0 * modefreqs->data[5]);
1272 // if (JLN > -1.0*kappa_thr || fabs(eta*JLN) < eJL_thr){
1273 // modefreqs->data[6] = conjl(-1.0 * modefreqs_xtr->data[1]);
1274 // }
1275  }
1276  }
1277  XLALDestroyCOMPLEX16Vector(modefreqs_xtr);
1278  if (debugout) {
1279  XLAL_PRINT_INFO("l,m = %d %d\n",l,m);
1280  XLAL_PRINT_INFO("finalSpin = %f\n",finalSpin);
1281  XLAL_PRINT_INFO("finalMass = %f\n",finalMass);
1282  XLAL_PRINT_INFO("PeakOmega = %f\n",NRPeakOmega22*mTot);
1283  XLAL_PRINT_INFO("w0 = %f, t0 = %f\n",creal(modefreqs->data[0])*mTot, 1./cimag(modefreqs->data[0])/mTot);
1284  XLAL_PRINT_INFO("w1 = %f, t1 = %f\n",creal(modefreqs->data[1])*mTot, 1./cimag(modefreqs->data[1])/mTot);
1285  XLAL_PRINT_INFO("w2 = %f, t2 = %f\n",creal(modefreqs->data[2])*mTot, 1./cimag(modefreqs->data[2])/mTot);
1286  XLAL_PRINT_INFO("w3 = %f, t3 = %f\n",creal(modefreqs->data[3])*mTot, 1./cimag(modefreqs->data[3])/mTot);
1287  XLAL_PRINT_INFO("w4 = %f, t4 = %f\n",creal(modefreqs->data[4])*mTot, 1./cimag(modefreqs->data[4])/mTot);
1288  XLAL_PRINT_INFO("w5 = %f, t5 = %f\n",creal(modefreqs->data[5])*mTot, 1./cimag(modefreqs->data[5])/mTot);
1289  XLAL_PRINT_INFO("w6 = %f, t6 = %f\n",creal(modefreqs->data[6])*mTot, 1./cimag(modefreqs->data[6])/mTot);
1290  XLAL_PRINT_INFO("w7 = %f, t7 = %f\n",creal(modefreqs->data[7])*mTot, 1./cimag(modefreqs->data[7])/mTot);
1291  XLAL_PRINT_INFO("matchrange %.16e,%.16e,%.16e\n", matchrange->data[0] * mTot / dt, matchrange->data[1] * mTot / dt, matchrange->data[2] * mTot / dt - 2);
1292  }
1293  /*XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, 1, nmodes, approximant);
1294 
1295 
1296 
1297  modefreqs->data[6] = modefreqs_xtr->data[0];
1298  modefreqs->data[7] = modefreqs_xtr->data[1];
1299  //modefreqs->data[6] = conjl(-1.0 * modefreqs->data[0]);
1300  //modefreqs->data[7] = conjl(-1.0 * modefreqs->data[1]);
1301  //modefreqs->data[6] = conjl(-1.0 * modefreqs_xtr->data[0]);
1302  //modefreqs->data[7] = conjl(-1.0 * modefreqs_xtr->data[1]);
1303  if (m <0) {
1304  //modefreqs->data[6] = modefreqs_xtr->data[0];
1305  //modefreqs->data[7] = modefreqs_xtr->data[1];
1306  modefreqs->data[6] = conjl(-1.0 * modefreqs_xtr->data[0]);
1307  modefreqs->data[7] = conjl(-1.0 * modefreqs_xtr->data[1]);
1308  }
1309  XLALDestroyCOMPLEX16Vector(modefreqs_xtr);*/
1310 
1311 
1312 
1313  // }}}
1314 
1315  } //end if m= 2, l = 2
1316 
1317  if ((int)fabs((REAL8) m) == 1 && l == 2) {
1318  if (debugout)
1319  XLAL_PRINT_INFO("Stas, default (if nothing specified above) for pQNM eta = %f, chi = %f chi1 = %f \n", eta, chi, chi1);
1320  if (debugout) {
1321  XLAL_PRINT_INFO("Stas: JLN*eta = %f, \n", JLN*eta);
1322  XLAL_PRINT_INFO("NRPeakOmega22 = %3.10f, %3.10f, %3.10f, %3.10f\n", NRPeakOmega22 * mTot /finalMass, NRPeakOmega22 / finalMass, finalMass, mTot);
1323  for (j = 0; j < nmodes; j++) {
1324  XLAL_PRINT_INFO("QNM frequencies: %d %d %d %3.10f %3.10f\n", l, m, j, creal(modefreqs->data[j]) * mTot, 1. / cimag(modefreqs->data[j]) / mTot);
1325  }
1326  }
1327 
1328  COMPLEX16Vector *modefreqs_xtr;
1329  modefreqs_xtr = XLALCreateCOMPLEX16Vector(nmodes);
1330  if (!modefreqs_xtr) {
1331  XLAL_ERROR(XLAL_ENOMEM);
1332  }
1333 
1334  if ( 1==1 || (JLN < 0.0 && eta <= 0.1)){
1335  spin1[0] *= -1;
1336  spin1[1] *= -1;
1337  spin1[2] *= -1;
1338  spin2[0] *= -1;
1339  spin2[1] *= -1;
1340  spin2[2] *= -1;
1341  XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, -mHere, nmodes, approximant);
1342  spin1[0] *= -1;
1343  spin1[1] *= -1;
1344  spin1[2] *= -1;
1345  spin2[0] *= -1;
1346  spin2[1] *= -1;
1347  spin2[2] *= -1;
1348  modefreqs->data[7] = modefreqs_xtr->data[0];
1349  modefreqs->data[6] = NRPeakOmega22 + I/mTot/((1./cimag(modefreqs->data[0])/mTot)/2.);
1350 
1351  if (m == -1){
1352  modefreqs->data[6] = conjl(-1.0 * modefreqs->data[6]);
1353  modefreqs->data[7] = conjl(-1.0 * modefreqs->data[7]);
1354  }
1355  }
1356 // modefreqs->data[7] = 0.5*(NRPeakOmega22 + modefreqs->data[0]) + I/mTot/((1./cimag(modefreqs->data[0])/mTot)/3.);
1357  XLALDestroyCOMPLEX16Vector(modefreqs_xtr);
1358  if (debugout) {
1359  XLAL_PRINT_INFO("l,m = %d %d\n",l,m);
1360  XLAL_PRINT_INFO("finalSpin = %f\n",finalSpin);
1361  XLAL_PRINT_INFO("finalMass = %f\n",finalMass);
1362  XLAL_PRINT_INFO("PeakOmega = %f\n",NRPeakOmega22*mTot);
1363  XLAL_PRINT_INFO("w0 = %f, t0 = %f\n",creal(modefreqs->data[0])*mTot, 1./cimag(modefreqs->data[0])/mTot);
1364  XLAL_PRINT_INFO("w1 = %f, t1 = %f\n",creal(modefreqs->data[1])*mTot, 1./cimag(modefreqs->data[1])/mTot);
1365  XLAL_PRINT_INFO("w2 = %f, t2 = %f\n",creal(modefreqs->data[2])*mTot, 1./cimag(modefreqs->data[2])/mTot);
1366  XLAL_PRINT_INFO("w3 = %f, t3 = %f\n",creal(modefreqs->data[3])*mTot, 1./cimag(modefreqs->data[3])/mTot);
1367  XLAL_PRINT_INFO("w4 = %f, t4 = %f\n",creal(modefreqs->data[4])*mTot, 1./cimag(modefreqs->data[4])/mTot);
1368  XLAL_PRINT_INFO("w5 = %f, t5 = %f\n",creal(modefreqs->data[5])*mTot, 1./cimag(modefreqs->data[5])/mTot);
1369  XLAL_PRINT_INFO("w6 = %f, t6 = %f\n",creal(modefreqs->data[6])*mTot, 1./cimag(modefreqs->data[6])/mTot);
1370  XLAL_PRINT_INFO("w7 = %f, t7 = %f\n",creal(modefreqs->data[7])*mTot, 1./cimag(modefreqs->data[7])/mTot);
1371  XLAL_PRINT_INFO("matchrange %.16e,%.16e,%.16e\n", matchrange->data[0] * mTot / dt, matchrange->data[1] * mTot / dt, matchrange->data[2] * mTot / dt - 2);
1372  }
1373  }
1374  // FIXME
1375  if (1==0 &&(m==1 || m==-1)){
1376  /*NRPeakOmega22 *= 0.5;
1377  if (JLN*m < 0.0){
1378  modefreqs->data[5] = (-3. / 4. * NRPeakOmega22 / finalMass) + (1. / 4. * creal(modefreqs->data[0]));
1379  modefreqs->data[4] = (-2. / 3. * NRPeakOmega22 / finalMass) + (1. / 3. * creal(modefreqs->data[0]));
1380  } else {
1381  modefreqs->data[5] = (3. / 4. * NRPeakOmega22 / finalMass) + (1. / 4. * creal(modefreqs->data[0]));
1382  modefreqs->data[4] = (2. / 3. * NRPeakOmega22 / finalMass) + (1. / 3. * creal(modefreqs->data[0]));
1383  //XLAL_PRINT_INFO("Stas , here m = %d modes: 0 = %4.10f, 6 = %4.10f, 7 = %4.10f \n", m, creal(modefreqs->data[0])*mTot, creal(modefreqs->data[6])*mTot, creal(modefreqs->data[7])*mTot);
1384  }
1385  modefreqs->data[5] += I * 3.5 / 0.9 * cimag(modefreqs->data[0]);
1386  modefreqs->data[4] += I * 3.5 * cimag(modefreqs->data[0]);*/
1387 
1388  COMPLEX16Vector *modefreqs_xtr;
1389  modefreqs_xtr = XLALCreateCOMPLEX16Vector(nmodes);
1390  if (!modefreqs_xtr) {
1391  XLAL_ERROR(XLAL_ENOMEM);
1392  }
1393  //XLAL_PRINT_INFO("Stas, generating QNM freqs... ");
1394  if (JLN > 0.0){
1395  XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, -1, nmodes, approximant);
1396  //XLAL_PRINT_INFO("done\n");
1397  modefreqs->data[5] = modefreqs_xtr->data[0];
1398  //modefreqs->data[6] = modefreqs_xtr->data[1];
1399  //modefreqs->data[7] = modefreqs_xtr->data[2];
1400  if (JLN < kappa_thr || eta*JLN < eJL_thr){
1401  modefreqs->data[6] = modefreqs_xtr->data[1];
1402  modefreqs->data[7] = modefreqs_xtr->data[2];
1403  }
1404  if (m == -1){
1405  //XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, -1, nmodes, approximant);
1406  modefreqs->data[5] = conjl(-1.0 * modefreqs_xtr->data[0]);
1407  //modefreqs->data[6] = conjl(-1.0 * modefreqs_xtr->data[1]);
1408  //modefreqs->data[7] = conjl(-1.0 * modefreqs_xtr->data[2]);
1409  if (JLN < kappa_thr || eta*JLN < eJL_thr){
1410  modefreqs->data[6] = conjl(-1.0 * modefreqs_xtr->data[1]);
1411  modefreqs->data[7] = conjl(-1.0 * modefreqs_xtr->data[2]);
1412  }
1413  }
1414  }
1415  //XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, -2, nmodes, approximant);
1416  //if (JLN > 0.0){
1417  // modefreqs->data[6] = modefreqs_xtr->data[0];
1418  // modefreqs->data[5] = modefreqs_xtr->data[1];
1419  // if (m == -1){
1420  // XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, 2, nmodes, approximant);
1421  // modefreqs->data[6] = conjl(-1.0 * modefreqs_xtr->data[0]);
1422  // modefreqs->data[5] = conjl(-1.0 * modefreqs_xtr->data[1]);
1423  // }
1424  //}
1425  if (JLN < 0.0){
1426  //modefreqs->data[4] = conjl(-1.0 * modefreqs->data[4]);
1427  //modefreqs->data[5] = conjl(-1.0 * modefreqs->data[5]);
1428  XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, 1, nmodes, approximant);
1429  //XLAL_PRINT_INFO("done\n");
1430  modefreqs->data[5] = modefreqs_xtr->data[0];
1431  //modefreqs->data[6] = modefreqs_xtr->data[1];
1432  //modefreqs->data[7] = modefreqs_xtr->data[2];
1433  if (JLN > -1.0*kappa_thr || fabs(eta*JLN) < eJL_thr){
1434  modefreqs->data[6] = modefreqs_xtr->data[1];
1435  modefreqs->data[7] = modefreqs_xtr->data[2];
1436  }
1437  if (m == -1){
1438  //XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, -1, nmodes, approximant);
1439  modefreqs->data[5] = conjl(-1.0 * modefreqs_xtr->data[0]);
1440  //modefreqs->data[6] = conjl(-1.0 * modefreqs_xtr->data[1]);
1441  //modefreqs->data[7] = conjl(-1.0 * modefreqs_xtr->data[2]);
1442  if (JLN > -1.0*kappa_thr || fabs(eta*JLN) < eJL_thr){
1443  modefreqs->data[6] = conjl(-1.0 * modefreqs_xtr->data[1]);
1444  modefreqs->data[7] = conjl(-1.0 * modefreqs_xtr->data[2]);
1445  }
1446  }
1447  //XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, 2, nmodes, approximant);
1448  //modefreqs->data[6] = modefreqs_xtr->data[0];
1449  //modefreqs->data[5] = modefreqs_xtr->data[1];
1450  //if (m == -1){
1451  //XLALSimIMREOBGenerateQNMFreqV2Prec(modefreqs_xtr, mass1, mass2, spin1, spin2, l, 2, nmodes, approximant);
1452  // modefreqs->data[6] = conjl(-1.0 * modefreqs_xtr->data[0]);
1453  //modefreqs->data[5] = conjl(-1.0 * modefreqs_xtr->data[1]);
1454  //}
1455  }
1456 
1457 
1458  XLALDestroyCOMPLEX16Vector(modefreqs_xtr);
1459  }
1460  } //v3
1461 
1462  // FIXME !!!!
1463  //for (j =0; j<nmodes; j++){
1464  //if (m > 0) {
1465  // for (j = 0; j < nmodes; j++) {
1466  // modefreqs->data[j] = conjl(-1.0 * modefreqs->data[j]);
1467  // }
1468  //}
1469 
1470 
1471  // Move ringdown comb boundaries to sampling points to avoid numerical artifacts.
1472 // matchrange->data[0] -= fmod(matchrange->data[0], dt / mTot);
1473 // matchrange->data[1] -= fmod(matchrange->data[1], dt / mTot);
1474  if (debugout) {
1475  XLAL_PRINT_INFO("NRPeakOmega22 = %3.10f, %3.10f, %3.10f, %3.10f\n", NRPeakOmega22 * mTot /finalMass, NRPeakOmega22 / finalMass, finalMass, mTot);
1476  for (j = 0; j < nmodes; j++) {
1477  XLAL_PRINT_INFO("QNM frequencies: %d %d %d %3.10f %3.10f\n", l, m, j, creal(modefreqs->data[j]) * mTot, 1. / cimag(modefreqs->data[j]) / mTot);
1478  }
1479  }
1480  /*
1481  * Ringdown signal length: 10 times the decay time of the n=0
1482  * mode
1483  */
1484  Nrdwave = (INT4) (EOB_RD_EFOLDS / cimag(modefreqs->data[0]) / dt);
1485  //XLAL_PRINT_INFO("STas, EOB_RD_EFOLDS = %f, freq = %f \n", EOB_RD_EFOLDS, cimag(modefreqs->data[0])*dt);
1486 
1487  /*
1488  * Check the value of attpos, to prevent memory access
1489  * problems later
1490  */
1491  if (matchrange->data[0] * mTot / dt < 5 || matchrange->data[1] * mTot / dt > matchrange->data[2] * mTot / dt - 2) {
1492  XLALPrintError("More inspiral points needed for ringdown matching.\n");
1493  XLAL_PRINT_INFO("%.16e,%.16e,%.16e\n", matchrange->data[0] * mTot / dt, matchrange->data[1] * mTot / dt, matchrange->data[2] * mTot / dt - 2);
1494  XLALDestroyCOMPLEX16Vector(modefreqs);
1495  XLAL_ERROR(XLAL_EFAILED);
1496  }
1497  /*
1498  * STEP 2) Based on least-damped-mode decay time, allocate
1499  * memory for rigndown waveform
1500  */
1501 
1502  /*
1503  * Create memory for the ring-down and full waveforms, and
1504  * derivatives of inspirals
1505  */
1506 
1507  rdwave1 = XLALCreateREAL8Vector(Nrdwave);
1508  rdwave2 = XLALCreateREAL8Vector(Nrdwave);
1509  rinspwave = XLALCreateREAL8Vector(6);
1510  dinspwave = XLALCreateREAL8Vector(6);
1511  ddinspwave = XLALCreateREAL8Vector(6);
1512  inspwaves1 = XLALCreateREAL8VectorSequence(3, 6);
1513  inspwaves2 = XLALCreateREAL8VectorSequence(3, 6);
1514 
1515  /* Check memory was allocated */
1516  if (!rdwave1 || !rdwave2 || !rinspwave || !dinspwave
1517  || !ddinspwave || !inspwaves1 || !inspwaves2) {
1518  XLALDestroyCOMPLEX16Vector(modefreqs);
1519  if (rdwave1)
1520  XLALDestroyREAL8Vector(rdwave1);
1521  if (rdwave2)
1522  XLALDestroyREAL8Vector(rdwave2);
1523  if (rinspwave)
1524  XLALDestroyREAL8Vector(rinspwave);
1525  if (dinspwave)
1526  XLALDestroyREAL8Vector(dinspwave);
1527  if (ddinspwave)
1528  XLALDestroyREAL8Vector(ddinspwave);
1529  if (inspwaves1)
1530  XLALDestroyREAL8VectorSequence(inspwaves1);
1531  if (inspwaves2)
1532  XLALDestroyREAL8VectorSequence(inspwaves2);
1533  XLAL_ERROR(XLAL_ENOMEM);
1534  }
1535  memset(rdwave1->data, 0, rdwave1->length * sizeof(REAL8));
1536  memset(rdwave2->data, 0, rdwave2->length * sizeof(REAL8));
1537 
1538  /*
1539  * STEP 3) Get values and derivatives of inspiral waveforms
1540  * at matching comb points
1541  */
1542 
1543  /* Generate derivatives of the last part of inspiral waves */
1544  /* Get derivatives of signal1 */
1545  if (XLALGenerateHybridWaveDerivatives(rinspwave, dinspwave, ddinspwave, timeVec, signal1,
1546  matchrange, dt, mass1, mass2) == XLAL_FAILURE) {
1547  XLALDestroyCOMPLEX16Vector(modefreqs);
1548  XLALDestroyREAL8Vector(rdwave1);
1549  XLALDestroyREAL8Vector(rdwave2);
1550  XLALDestroyREAL8Vector(rinspwave);
1551  XLALDestroyREAL8Vector(dinspwave);
1552  XLALDestroyREAL8Vector(ddinspwave);
1553  XLALDestroyREAL8VectorSequence(inspwaves1);
1554  XLALDestroyREAL8VectorSequence(inspwaves2);
1555  XLAL_ERROR(XLAL_EFUNC);
1556  }
1557  for (j = 0; j < 6; j++) {
1558  inspwaves1->data[j] = rinspwave->data[j];
1559  inspwaves1->data[j + 6] = dinspwave->data[j];
1560  inspwaves1->data[j + 12] = ddinspwave->data[j];
1561  }
1562 
1563  /* Get derivatives of signal2 */
1564  if (XLALGenerateHybridWaveDerivatives(rinspwave, dinspwave, ddinspwave, timeVec, signal2,
1565  matchrange, dt, mass1, mass2) == XLAL_FAILURE) {
1566  XLALDestroyCOMPLEX16Vector(modefreqs);
1567  XLALDestroyREAL8Vector(rdwave1);
1568  XLALDestroyREAL8Vector(rdwave2);
1569  XLALDestroyREAL8Vector(rinspwave);
1570  XLALDestroyREAL8Vector(dinspwave);
1571  XLALDestroyREAL8Vector(ddinspwave);
1572  XLALDestroyREAL8VectorSequence(inspwaves1);
1573  XLALDestroyREAL8VectorSequence(inspwaves2);
1574  XLAL_ERROR(XLAL_EFUNC);
1575  }
1576  for (j = 0; j < 6; j++) {
1577  inspwaves2->data[j] = rinspwave->data[j];
1578  inspwaves2->data[j + 6] = dinspwave->data[j];
1579  inspwaves2->data[j + 12] = ddinspwave->data[j];
1580  }
1581 
1582 
1583  /*
1584  * STEP 4) Solve QNM coefficients and generate ringdown
1585  * waveforms
1586  */
1587 
1588  /* Generate ring-down waveforms */
1589  if (XLALSimIMREOBHybridRingdownWave(rdwave1, rdwave2, dt, mass1, mass2, inspwaves1, inspwaves2,
1590  modefreqs, matchrange) == XLAL_FAILURE) {
1591  XLALDestroyCOMPLEX16Vector(modefreqs);
1592  XLALDestroyREAL8Vector(rdwave1);
1593  XLALDestroyREAL8Vector(rdwave2);
1594  XLALDestroyREAL8Vector(rinspwave);
1595  XLALDestroyREAL8Vector(dinspwave);
1596  XLALDestroyREAL8Vector(ddinspwave);
1597  XLALDestroyREAL8VectorSequence(inspwaves1);
1598  XLALDestroyREAL8VectorSequence(inspwaves2);
1599  XLAL_ERROR(XLAL_EFUNC);
1600  }
1601 
1602  //XLAL_PRINT_INFO("Stas we're down with RD solving, attaching it now....");
1603  /*
1604  * STEP 5) Stitch inspiral and ringdown waveoforms
1605  */
1606 
1607  /*
1608  * Generate full waveforms, by stitching inspiral and
1609  * ring-down waveforms
1610  */
1611  //UINT4 attachIdx = matchrange->data[1] * mTot / dt;
1612  UINT4 attachIdx = round(matchrange->data[1] * mTot / dt);
1613  for (j = 1; j < Nrdwave; ++j) {
1614  signal1->data[j + attachIdx] = rdwave1->data[j];
1615  signal2->data[j + attachIdx] = rdwave2->data[j];
1616  }
1617  //XLAL_PRINT_INFO(" Stas, done 1 length = %d, Nrdwave = %d, attachIdx = %d \n", signal1->length, Nrdwave, attachIdx);
1618 
1619  memset(signal1->data + Nrdwave + attachIdx, 0, (signal1->length - Nrdwave - attachIdx) * sizeof(REAL8));
1620  memset(signal2->data + Nrdwave + attachIdx, 0, (signal2->length - Nrdwave - attachIdx) * sizeof(REAL8));
1621 
1622  //XLAL_PRINT_INFO(" Stas, done 2 \n");
1623 
1624  /* Free memory */
1625  XLALDestroyCOMPLEX16Vector(modefreqs);
1626  XLALDestroyREAL8Vector(rdwave1);
1627  XLALDestroyREAL8Vector(rdwave2);
1628  XLALDestroyREAL8Vector(rinspwave);
1629  XLALDestroyREAL8Vector(dinspwave);
1630  XLALDestroyREAL8Vector(ddinspwave);
1631  XLALDestroyREAL8VectorSequence(inspwaves1);
1632  XLALDestroyREAL8VectorSequence(inspwaves2);
1633 
1634  //XLAL_PRINT_INFO(" Stas, done 3 \n");
1635 // matchrange->data[0] += fmod(matchrange->data[0], dt / mTot);
1636 // matchrange->data[1] += fmod(matchrange->data[1], dt / mTot);
1637 
1638  return XLAL_SUCCESS;
1639  }
1640 
1641 
1642 /* Search for index at which the maximum of the amplitude occurs */
1643 static INT4 XLALSimFindIndexMaxAmpli( UINT4 * indAmax, REAL8Vector * timeVec, REAL8Vector * ampWave, REAL8 * valAmax, REAL8 tofAmax ) {
1644  if ( indAmax == NULL || timeVec == NULL || ampWave == NULL || valAmax == NULL ) {
1645  return XLAL_FAILURE;
1646  }
1647  INT4 debugSB = 0;
1648  *indAmax = 0;
1649  INT4 found = 0;
1650  for (UINT4 i = 0; i < timeVec->length - 1; i++) {
1651  if (timeVec->data[i] == tofAmax) {
1652  found = 1;
1653  *indAmax = i;
1654  *valAmax = ampWave->data[i];
1655  }
1656  }
1657  if (found == 0) {
1658  return XLAL_FAILURE;
1659  }
1660  else {
1661  if (debugSB) {
1662  printf(" found maximum times: %f, %f \n", timeVec->data[*indAmax], tofAmax );
1663  }
1664  return XLAL_SUCCESS;
1665  }
1666 }
1667 /**
1668  * The main function for performing the ringdown attachment for SEOBNRv4 (and beyond)
1669  * This is the function which gets called by the code generating the full IMR waveform once
1670  * generation of the inspiral part has been completed.
1671  * The ringdown is attached by factoring the less damped harmonics and apply tanh fit to the rest.
1672  * STEP 1) Get mass and spin of the final black hole and the complex ringdown frequencies
1673  * STEP 2) Apply the fit function from the attachment time
1674  * STEP 3) Construct the full RD by applying (factor) of less damped 220 mode
1675  * STEP 4) Constructing the RD stitched to the inspiral-merger
1676  * See Section 5.2 of https://dcc.ligo.org/T1600383
1677  */
1678 static UNUSED INT4 XLALSimIMREOBAttachFitRingdown(
1679  REAL8Vector * signal1, /**<< OUTPUT, Real of inspiral waveform to which we attach ringdown */
1680  REAL8Vector * signal2, /**<< OUTPUT, Imag of inspiral waveform to which we attach ringdown */
1681  const INT4 l, /**<< Current mode l */
1682  const INT4 m, /**<< Current mode m */
1683  const REAL8 dt, /**<< Sample time step (in seconds) */
1684  const REAL8 mass1, /**<< First component mass (in Solar masses) */
1685  const REAL8 mass2, /**<< Second component mass (in Solar masses) */
1686  const REAL8 spin1x, /**<<The spin of the first object; */
1687  const REAL8 spin1y, /**<<The spin of the first object; */
1688  const REAL8 spin1z, /**<<The spin of the first object; */
1689  const REAL8 spin2x, /**<<The spin of the second object; */
1690  const REAL8 spin2y, /**<<The spin of the second object; */
1691  const REAL8 spin2z, /**<<The spin of the second object; */
1692  const REAL8 finalM, /**<<The spin of the final BH; */
1693  const REAL8 finalS, /**<< The magnitude of the final spin */
1694  REAL8Vector * timeVec, /**<< Vector containing the time values */
1695  REAL8Vector * matchrange,
1696  /**<< Time values chosen as points for performing comb matching */
1697  Approximant approximant,/**<<The waveform approximant being used */
1698  UINT4 * indAmpMax /**<<This is used only for SEOBNRv4HM, it is needed to remember the attaching point of the 22 mode */
1699  ) {
1700  if ( mass1 <= 0. || mass2 <= 0.
1701  || fabs(spin1x) > 1. || fabs(spin1y) > 1. || fabs(spin1z) > 1.
1702  || fabs(spin2x) > 1. || fabs(spin2y) > 1. || fabs(spin2z) > 1.
1703  || spin1x*spin1x + spin1y*spin1y + spin1z*spin1z > 1.
1704  || spin2x*spin2x + spin2y*spin2y + spin2z*spin2z > 1.
1705  ) {
1706  return XLAL_FAILURE;
1707  }
1708  INT4 debugSB = 0;
1709  UINT4 i;
1710  UNUSED INT4 phasecount;
1711  REAL8Vector *ampWave;
1712  REAL8Vector *phWave;
1713  REAL8Vector *rdtime;
1714  COMPLEX16Vector *modefreqs;
1715  COMPLEX16Vector *modefreqs22; //RC: we use this variable to compute the damping time of the 22 mode which will be used to set the lenght of the ringdown for all the modes
1716 
1717  ampWave = XLALCreateREAL8Vector(signal1->length);
1718  phWave = XLALCreateREAL8Vector(signal1->length);
1719 
1720  REAL8 mtot = mass1 + mass2;
1721  REAL8 eta = mass1 * mass2 / (mtot * mtot);
1722 
1723  /* Here I assume that the spins were properly projected (is precessing) and only spin1z, spin2z
1724  * are relevant, if not we need to add extra function to define what we call chi1, chi2 */
1725  REAL8 chi1 = spin1z;
1726  REAL8 chi2 = spin2z;
1727  //REAL8 spin1[3] = { spin1x, spin1y, spin1z };
1728  //REAL8 spin2[3] = { spin2x, spin2y, spin2z };
1729 
1730  Approximant appr;
1731  appr = approximant;
1732 
1733  gsl_spline *spline0 = NULL;
1734  gsl_interp_accel *acc0 = NULL;
1735 
1736 
1737 
1738  if (debugSB) {
1739  printf("RDfit: we use spin1 %f, spin2 %f, and it should be dimensionless [-1,1] \n", chi1, chi2);
1740  printf("We use approximant = %d \n", appr);
1741  }
1742  REAL8 chi = 0.5 * (chi1 + chi2) + 0.5 * (chi1 - chi2) * (mass1 - mass2)/(mass1 + mass2) / (1.0 - 2.0 * eta);
1743 
1744 
1745 
1746  /*********************************************************************************************/
1747  /* QNMs of the remnant */
1748  /*********************************************************************************************/
1749  /* Getting QNMs */
1750  modefreqs = XLALCreateCOMPLEX16Vector(1);
1751  if (XLALSimIMREOBGenerateQNMFreqV2FromFinalPrec(modefreqs, mass1, mass2, finalM, finalS , l, m, 1) == XLAL_FAILURE) {
1752  XLALDestroyCOMPLEX16Vector(modefreqs);
1753  XLAL_ERROR(XLAL_EFUNC);
1754  }
1755 
1756  //RC: we use this variable to compute the damping time of the 22 mode which will be used to set the lenght of the ringdown for all the modes
1757  modefreqs22 = XLALCreateCOMPLEX16Vector(1);
1758  if (XLALSimIMREOBGenerateQNMFreqV2FromFinalPrec(modefreqs22, mass1, mass2, finalM, finalS, 2, 2, 1) == XLAL_FAILURE) {
1759  XLALDestroyCOMPLEX16Vector(modefreqs);
1760  XLAL_ERROR(XLAL_EFUNC);
1761  }
1762  /* Least-damped QNM */
1763  COMPLEX16 sigmalm0;
1764  sigmalm0 = (-cimag(modefreqs->data[0]) - I * creal(modefreqs->data[0])) * (mtot * LAL_MTSUN_SI);
1765  // sigmalm0 = -0.08119703120243188 - I*0.5040714995216017;
1766 
1767  if (debugSB) {
1768  printf("Mode %d %d \n",l, m);
1769  printf("matchpoints are: %.16f, %.16f, %.16f \n", matchrange->data[0], matchrange->data[1], matchrange->data[2]);
1770  printf("the 0-overtone is: %.16f + i %.16f \n", creal(sigmalm0), cimag(sigmalm0));
1771  }
1772 
1773 
1774 
1775  /*********************************************************************************************/
1776  /* Prepare inspiral signal */
1777  /* This is needed to guarantee continuity with RD signal */
1778  /*********************************************************************************************/
1779  /* Compute amplitude and the unwrapped phase of the inspiral */
1780  phasecount = 0;
1781  for (i = 0; i < signal1->length; i++) {
1782  ampWave->data[i] = 0.0;
1783  phWave->data[i] = 0.0;
1784  }
1785  REAL8 prev, now, corph;
1786  prev = atan2(signal2->data[0], signal1->data[0]);
1787  phWave->data[0] = prev;
1788  for (i = 0; i < timeVec->length; i++) {
1789  ampWave->data[i] = sqrt(signal1->data[i] * signal1->data[i] + signal2->data[i] * signal2->data[i]);
1790  now = atan2(signal2->data[i], signal1->data[i]);
1791  if (i > 0) {
1792  /* Unwrapping */
1793  corph = now - prev;
1794  corph = corph > LAL_PI ? corph - LAL_TWOPI : (corph < -LAL_PI ? corph + LAL_TWOPI : corph);
1795 
1796  phWave->data[i] = phWave->data[i - 1] + corph;
1797  prev = now;
1798  }
1799  //phWave->data[i] = now;
1800  }
1801  FILE *fout = NULL;
1802  if (debugSB) {
1803  char filename2[sizeof "CheckStasAmplPhaseXX.dat"];
1804  sprintf(filename2,"CheckStasAmplPhase%01d%01d.dat",l,m);
1805  fout = fopen(filename2, "w");
1806  printf("Check the length: time %d, signal %d \n", timeVec->length, signal1->length);
1807  for (i = 0; i < timeVec->length; i++) {
1808  fprintf(fout, "%f %.16e %.16e %.16e %.16e \n", timeVec->data[i], ampWave->data[i], phWave->data[i], signal1->data[i], signal2->data[i]);
1809  }
1810 
1811  fclose(fout);
1812  }
1813  /* Search for index at which the maximum of the amplitude occurs */
1814  REAL8 valAmax = ampWave->data[0];
1815  REAL8 valdAmax = 0;
1816  REAL8 tofAmax = matchrange->data[1];
1817  UINT4 indAmax;
1818  if (l == 5 && m==5) {
1819  tofAmax = matchrange->data[3];
1820  }
1821  //RC we should find this only for the mode 22, for the other one the attaching point is the same, with the exception of the 55
1822  //For the 55 mode we need to find tpeak22 -10M, which in that case is fed as matchrange->data[1]
1823  if((l == 2 && m == 2) || (l == 5 && m == 5)){
1824  if ( XLALSimFindIndexMaxAmpli( &indAmax, timeVec, ampWave, &valAmax, tofAmax ) == XLAL_FAILURE )
1825  {
1826  XLALPrintError("Time of maximum amplitude is not found .\n");
1827  XLALDestroyCOMPLEX16Vector(modefreqs);
1828  XLAL_ERROR(XLAL_EFUNC);
1829  }
1830  *indAmpMax = indAmax;
1831  //RC for the 55 mode we also need to find the derivative at the attachment point, since we are not making the attachment
1832  //at the peak of the mode where the derivative is zero
1833  if(l == 5 && m == 5){
1834  spline0 = gsl_spline_alloc (gsl_interp_cspline, timeVec->length);
1835  acc0 = gsl_interp_accel_alloc ();
1836  gsl_spline_init (spline0, timeVec->data, ampWave->data, timeVec->length);
1837  gsl_interp_accel_reset (acc0);
1838  valdAmax = gsl_spline_eval_deriv (spline0, tofAmax, acc0);
1839  }
1840  }
1841  else{
1842  indAmax = *indAmpMax;
1843  valAmax = ampWave->data[indAmax]; //For the modes different from 22, this IS NOT the value of the maximum of the amplitude of the mode, but is the value of the amplitude at the peak of the 22 mode
1844  //RC: for the SEOBNRv4HM model we need to compute also the derivative at the attaching point. This is not zero because we are attaching the HMs not at the paek of the lm mode, but at the peak of the 22 mode
1845  spline0 = gsl_spline_alloc (gsl_interp_cspline, timeVec->length);
1846  acc0 = gsl_interp_accel_alloc ();
1847  gsl_spline_init (spline0, timeVec->data, ampWave->data, timeVec->length);
1848  gsl_interp_accel_reset (acc0);
1849  valdAmax = gsl_spline_eval_deriv (spline0, tofAmax, acc0);
1850  // valAmax = 0.019835031225903733*AMP0;
1851  // valdAmax = 0.00020163171965464758*AMP0;
1852  }
1853 
1854 
1855  if (debugSB) {
1856  printf("Check: The maximum of amplitude is %.16e found at t=%f, index = %d (out of %d) \n", valAmax, tofAmax, indAmax, timeVec->length);
1857  printf("Amp@Attachment = %.16f \n", valAmax);
1858  printf("dAmp@Attachment = %.16f \n", valdAmax);
1859  FILE *indMax = NULL;
1860  char filename2[sizeof "indexMaxXXHi.dat"];
1861  sprintf(filename2,"indexMax%01d%01dHi.dat",l,m);
1862  indMax = fopen (filename2, "w");
1863  fprintf(indMax, "%d \n",indAmax);
1864  fclose(indMax);
1865  }
1866 
1867 
1868 
1869  /*********************************************************************************************/
1870  /* Constant coefficients entering RD fitting formulas */
1871  /*********************************************************************************************/
1872  /* Computing fit coefficients that enter amplitude and phase of the RD signal */
1873  /* Numeircal constants are defined at the top of this file */
1874  /* Eq. (28) of https://dcc.ligo.org/T1600383 */
1875  REAL8 ampcf1;
1876  // ampcf1 = A1coeff00[l][m] + A1coeff01[l][m] * chi + A1coeff02[l][m] * chi * chi + A1coeff10[l][m] * eta + A1coeff11[l][m] * eta * chi + A1coeff20[l][m] * eta * eta;
1877  ampcf1 = XLALCalculateRDAmplitudeCoefficient1(l, m, eta, chi);
1878  //ampcf1 = 0.08953902308302486;
1879  if(debugSB){
1880  printf("ampcf1 = %.16f \n", ampcf1);
1881  }
1882 
1883  /* Eq. (29) of https://dcc.ligo.org/T1600383 */
1884  REAL8 ampcf2;
1885  // ampcf2 = A2coeff00[l][m] + A2coeff01[l][m] * chi + A2coeff10[l][m] * eta + A2coeff11[l][m] * eta * chi + A2coeff20[l][m] * eta * eta + A2coeff21[l][m] * eta * eta * chi;
1886  ampcf2 = XLALCalculateRDAmplitudeCoefficient2(l, m, eta, chi);
1887  //ampcf2 = -0.506368869538912;
1888  if(debugSB){
1889  printf("ampcf2 = %.16f \n", ampcf2);
1890  }
1891 // printf("creal(sigma220), 2.*ampcf1*tanh(ampcf2) = %.16e %.16e\n",1./creal(sigma220), 2.*ampcf1*tanh(ampcf2));
1892  /* Eqs. (31)-(32) of https://dcc.ligo.org/T1600383 */
1893  if( l == 2 && m == 2){
1894  if (creal(sigmalm0) > 2. * ampcf1 * tanh(ampcf2)) {
1895  ampcf1 = creal(sigmalm0) / (2. * tanh(ampcf2));
1896  }
1897  }
1898  /* Eq. (36) of https://dcc.ligo.org/T1600383 */
1899  REAL8 phasecf1;
1900  // phasecf1 = P1coeff00[l][m] + P1coeff01[l][m] * chi + P1coeff02[l][m] * chi * chi + P1coeff10[l][m] * eta + P1coeff11[l][m] * eta * chi + P1coeff20[l][m] * eta * eta;
1901  phasecf1 = XLALCalculateRDPhaseCoefficient1(l, m, eta, chi);
1902  if(debugSB){
1903  printf("phasecf1 = %.16f \n", phasecf1);
1904  }
1905 
1906  /* Eq. (37) of https://dcc.ligo.org/T1600383 */
1907  REAL8 phasecf2;
1908  // phasecf2 = P2coeff00[l][m] + P2coeff01[l][m] * chi + P2coeff02[l][m] * chi * chi + P2coeff10[l][m] * eta + P2coeff11[l][m] * eta * chi + P2coeff20[l][m] * eta * eta;
1909  phasecf2 = XLALCalculateRDPhaseCoefficient2(l, m, eta, chi);
1910  if(debugSB){
1911  printf("phasecf2 = %.16f \n", phasecf2);
1912  }
1913 
1914  /*********************************************************************************************/
1915  /* RD fitting formulas */
1916  /*********************************************************************************************/
1917  /* Ringdown signal length: 10 times the decay time of the n=0 mode */
1918  UINT4 Nrdwave = (INT4) (EOB_RD_EFOLDS / cimag(modefreqs22->data[0]) / dt);
1919  //printf("Stas Nrdwave %d, dt = %f", Nrdwave, dt);
1920  REAL8 dtM = dt / (mtot * LAL_MTSUN_SI); // go to geometric units
1921  rdtime = XLALCreateREAL8Vector(Nrdwave);
1922  for (i = 0; i < Nrdwave; i++) {
1923  rdtime->data[i] = i * dtM; // this time for RD and it starts with 0
1924  }
1925  REAL8 tcons = 0.; //attachment point relative to peak; possibility for non zero values not fully implemented
1926 
1927  /* Rescalings to guarantee continuity with inspiral */
1928  REAL8 Arescaledtcons = valAmax * exp(-creal(sigmalm0) * tcons) / eta;
1929  REAL8 dtArescaledtcons = (valdAmax - creal(sigmalm0) * valAmax) * exp(-creal(sigmalm0) * tcons) / eta; // valdAmax = 0 - assumes extermum (max) of peak amplitude //RC: valdAmax = 0 for 22 mode
1930  REAL8 ampcc1 = dtArescaledtcons * pow(cosh(ampcf2), 2) / ampcf1;
1931  REAL8 ampcc2 = (Arescaledtcons * ampcf1 - dtArescaledtcons * cosh(ampcf2) * sinh(ampcf2)) / ampcf1;
1932 
1933 
1934  REAL8Vector *ampRD;
1935  REAL8Vector *phRD;
1936  ampRD = XLALCreateREAL8Vector(Nrdwave);
1937  phRD = XLALCreateREAL8Vector(Nrdwave);
1938 
1939  /* Construct RD amplitude */
1940  /* Eqs. (22)-(23) of https://dcc.ligo.org/T1600383 */
1941  for (i = 0; i < Nrdwave; i++) {
1942  ampRD->data[i] = eta * exp(creal(sigmalm0) * rdtime->data[i]) * (ampcc1 * tanh(ampcf1 * rdtime->data[i] + ampcf2) + ampcc2);
1943  }
1944 
1945  if (debugSB) {
1946  char filename[sizeof "StasAmpRD_fullXX.dat"];
1947  sprintf(filename,"StasAmpRD_full%01d%01d.dat",l,m);
1948  fout = fopen (filename, "w");
1949  // for (i = 0; i < indAmax; i++) {
1950  // fprintf(fout, "%.16e %.16e \n", timeVec->data[i], ampWave->data[i]/AMP0);
1951  // }
1952  for (i = 1; i < Nrdwave; i++) {
1953  fprintf(fout, "%.16e %.16e \n", rdtime->data[i] + tofAmax*0, ampRD->data[i]);
1954  }
1955  fclose(fout);
1956  }
1957 
1958  /* Construct RD phase */
1959  REAL8 omegarescaledtcons = (phWave->data[indAmax] - phWave->data[indAmax - 1]) / dtM - cimag(sigmalm0);
1960  // omegarescaledtcons = -0.2354501785436996 - cimag(sigmalm0); 21 mode
1961  // omegarescaledtcons = -0.6076660197512114 - cimag(sigmalm0); 33 mode
1962  // omegarescaledtcons = -0.8063013198270685 - cimag(sigmalm0); 44 mode
1963  // omegarescaledtcons = -0.7813933465833801 - cimag(sigmalm0); 55 mode
1964 
1965  if(debugSB){
1966  printf("omega0 = %.16f\n", phWave->data[0]);
1967  printf("omega0fit = %.16f, omega0numder = %.16f \n", XLALSimIMREOBGetNRSpinPeakOmegaV4 (l, m, eta, chi), (phWave->data[indAmax] - phWave->data[indAmax - 1]) / dtM);
1968  }
1969  REAL8 phasecc1 = omegarescaledtcons * (phasecf2 + 1.0) / phasecf2 / phasecf1;
1970  REAL8 phi0 = phWave->data[indAmax];
1971  // phi0 = -165.37581501842033; 21 mode
1972  // phi0 = -486.98433862793974; 33 mode
1973  // phi0 = -651.2335864724689; 44 mode
1974  // phi0 = -805.2637568966843; 55 mode
1975  if(debugSB){
1976  printf("phi_0 = %.16f \n", phi0);
1977  }
1978  REAL8 logargnum, logargden;
1979  /* Eq. (33)-(34) of https://dcc.ligo.org/T1600383 */
1980  for (i = 0; i < Nrdwave; i++) {
1981  logargnum = 1. + phasecf2 * exp(-phasecf1 * rdtime->data[i]);
1982  logargden = 1. + phasecf2;
1983  phRD->data[i] = phi0 - phasecc1 * log(logargnum / logargden) + cimag(sigmalm0) * rdtime->data[i];
1984  }
1985  if (debugSB) {
1986  char filename[sizeof "StasPhRD_fullXX.dat"];
1987  sprintf(filename,"StasPhRD_full%01d%01d.dat",l,m);
1988  fout = fopen (filename, "w");
1989  // for (i = 0; i < indAmax; i++) {
1990  // fprintf(fout, "%.16e %.16e \n", timeVec->data[i], phWave->data[i]);
1991  // }
1992  for (i = 1; i < Nrdwave; i++) {
1993  fprintf(fout, "%.16e %.16e \n", rdtime->data[i] + tofAmax*0, phRD->data[i]);
1994  }
1995  fclose(fout);
1996 
1997  // Compute the frequency of the ful signal
1998  UINT4 totSz = indAmax + Nrdwave;
1999  REAL8Vector *PhFull;
2000  PhFull = XLALCreateREAL8Vector(totSz);
2001  REAL8Vector *tFull;
2002  tFull = XLALCreateREAL8Vector(totSz);
2003  REAL8Vector *frFull;
2004  frFull = XLALCreateREAL8Vector(totSz);
2005 
2006  for (i = 0; i < indAmax; i++) {
2007  tFull->data[i] = timeVec->data[i];
2008  PhFull->data[i] = phWave->data[i];
2009  }
2010  for (i = 0; i < Nrdwave; i++) {
2011  tFull->data[i + indAmax] = rdtime->data[i] + tofAmax;
2012  PhFull->data[i + indAmax] = phRD->data[i];
2013  }
2014  fout = fopen("StasPhRD_full2.dat", "w");
2015  for (i = 0; i < totSz; i++) {
2016  fprintf(fout, "%.16e %.16e \n", tFull->data[i], PhFull->data[i]);
2017  }
2018  fclose(fout);
2019 
2020  gsl_spline *spline = NULL;
2021  gsl_interp_accel *acc = NULL;
2022  spline = gsl_spline_alloc(gsl_interp_cspline, totSz);
2023  acc = gsl_interp_accel_alloc();
2024  gsl_spline_init(spline, tFull->data, PhFull->data, totSz);
2025  for (i = 0; i < totSz; i++) {
2026  frFull->data[i] = gsl_spline_eval_deriv(spline, tFull->data[i], acc);
2027  }
2028  fout = fopen("StasFrRD_full.dat", "w");
2029  for (i = 0; i < totSz; i++) {
2030  fprintf(fout, "%.16e %.16e \n", tFull->data[i], frFull->data[i]);
2031  }
2032  fclose(fout);
2033 
2034  gsl_spline_free(spline);
2035  gsl_interp_accel_free(acc);
2036 
2037  XLALDestroyREAL8Vector(PhFull);
2038  XLALDestroyREAL8Vector(tFull);
2039  XLALDestroyREAL8Vector(frFull);
2040 
2041  }
2042 
2043  /* Construct RD signal */
2044  for (i = 0; i < Nrdwave; i++) {
2045  signal1->data[i + indAmax] = ampRD->data[i] * cos(phRD->data[i]);
2046  signal2->data[i + indAmax] = ampRD->data[i] * sin(phRD->data[i]);
2047  }
2048 
2049  gsl_spline_free (spline0);
2050  gsl_interp_accel_free (acc0);
2051 
2052  XLALDestroyREAL8Vector(ampWave);
2053  XLALDestroyREAL8Vector(phWave);
2054  XLALDestroyCOMPLEX16Vector(modefreqs);
2055  XLALDestroyCOMPLEX16Vector(modefreqs22);
2056  XLALDestroyREAL8Vector(rdtime);
2057  XLALDestroyREAL8Vector(ampRD);
2058  XLALDestroyREAL8Vector(phRD);
2059 
2060  return XLAL_SUCCESS;
2061 
2062 }
2063 
2064 #endif /* _LALSIMIMREOBHYBRIDRINGDOWNPREC_C */
cosh
double cosh(double)
LALMalloc
#define LALMalloc(n)
LALFree
#define LALFree(p)
XLALSimIMREOBFinalMassSpinPrec
INT4 XLALSimIMREOBFinalMassSpinPrec(REAL8 *finalMass, REAL8 *finalSpin, const REAL8 mass1, const REAL8 mass2, const REAL8 spin1[3], const REAL8 spin2[3], Approximant approximant)
Computes the final mass and spin of the black hole resulting from merger.
Definition: LALSimBlackHoleRingdownPrec.c:39
XLALSimIMREOBGenerateQNMFreqV2FromFinalPrec
INT4 XLALSimIMREOBGenerateQNMFreqV2FromFinalPrec(COMPLEX16Vector *modefreqs, const REAL8 mass1, const REAL8 mass2, const REAL8 finalMass, REAL8 finalSpin, UINT4 l, INT4 m, UINT4 nmodes)
This function generates the quasinormal mode frequencies for a black hole ringdown.
Definition: LALSimBlackHoleRingdownPrec.c:1926
XLALSimIMREOBGenerateQNMFreqV2Prec
INT4 XLALSimIMREOBGenerateQNMFreqV2Prec(COMPLEX16Vector *modefreqs, const REAL8 mass1, const REAL8 mass2, const REAL8 spin1[3], const REAL8 spin2[3], UINT4 l, INT4 m, UINT4 nmodes, Approximant approximant)
This function generates the quasinormal mode frequencies for a black hole ringdown.
Definition: LALSimBlackHoleRingdownPrec.c:329
XLALCalculateRDPhaseCoefficient1
static REAL8 XLALCalculateRDPhaseCoefficient1(UINT4 modeL, UINT4 modeM, REAL8 eta, REAL8 chi)
Definition: LALSimIMREOBHybridRingdownPrec.c:163
XLALSimIMREOBHybridRingdownWave
static INT4 XLALSimIMREOBHybridRingdownWave(REAL8Vector *rdwave1, REAL8Vector *rdwave2, const REAL8 dt, const REAL8 mass1, const REAL8 mass2, REAL8VectorSequence *inspwave1, REAL8VectorSequence *inspwave2, COMPLEX16Vector *modefreqs, REAL8Vector *matchrange)
Generates the ringdown wave associated with the given real and imaginary parts of the inspiral wavefo...
Definition: LALSimIMREOBHybridRingdownPrec.c:297
XLALCalculateRDPhaseCoefficient2
static REAL8 XLALCalculateRDPhaseCoefficient2(UINT4 modeL, UINT4 modeM, REAL8 eta, REAL8 chi)
Definition: LALSimIMREOBHybridRingdownPrec.c:221
XLALGenerateHybridWaveDerivatives
static INT4 XLALGenerateHybridWaveDerivatives(REAL8Vector *rwave, REAL8Vector *dwave, REAL8Vector *ddwave, REAL8Vector *timeVec, REAL8Vector *wave, REAL8Vector *matchrange, REAL8 dt, REAL8 mass1, REAL8 mass2)
Function which calculates the value of the waveform, plus its first and second derivatives,...
Definition: LALSimIMREOBHybridRingdownPrec.c:513
XLALSimFindIndexMaxAmpli
static INT4 XLALSimFindIndexMaxAmpli(UINT4 *indAmax, REAL8Vector *timeVec, REAL8Vector *ampWave, REAL8 *valAmax, REAL8 tofAmax)
Definition: LALSimIMREOBHybridRingdownPrec.c:1643
XLALCalculateRDAmplitudeCoefficient2
static REAL8 XLALCalculateRDAmplitudeCoefficient2(UINT4 modeL, UINT4 modeM, REAL8 eta, REAL8 chi)
Definition: LALSimIMREOBHybridRingdownPrec.c:108
XLALCalculateRDAmplitudeCoefficient1
static REAL8 XLALCalculateRDAmplitudeCoefficient1(UINT4 modeL, UINT4 modeM, REAL8 eta, REAL8 chi)
Definition: LALSimIMREOBHybridRingdownPrec.c:49
get_major_approximant
static INT4 get_major_approximant(Approximant approximant)
Translates the major approximant version into an unsigned integer to minimize if and or statements th...
Definition: LALSimIMREOBHybridRingdownPrec.c:630
XLALSimIMREOBHybridAttachRingdownPrec
static UNUSED INT4 XLALSimIMREOBHybridAttachRingdownPrec(REAL8Vector *signal1, REAL8Vector *signal2, const INT4 l, const INT4 m, const REAL8 dt, const REAL8 mass1, const REAL8 mass2, const REAL8 spin1x, const REAL8 spin1y, const REAL8 spin1z, const REAL8 spin2x, const REAL8 spin2y, const REAL8 spin2z, REAL8Vector *timeVec, REAL8Vector *matchrange, Approximant approximant, const REAL8 JLN)
The main workhorse function for performing the ringdown attachment for EOB models EOBNRv2 and SEOBNRv...
Definition: LALSimIMREOBHybridRingdownPrec.c:663
XLALSimIMREOBAttachFitRingdown
static UNUSED INT4 XLALSimIMREOBAttachFitRingdown(REAL8Vector *signal1, REAL8Vector *signal2, const INT4 l, const INT4 m, const REAL8 dt, const REAL8 mass1, const REAL8 mass2, const REAL8 spin1x, const REAL8 spin1y, const REAL8 spin1z, const REAL8 spin2x, const REAL8 spin2y, const REAL8 spin2z, const REAL8 finalM, const REAL8 finalS, REAL8Vector *timeVec, REAL8Vector *matchrange, Approximant approximant, UINT4 *indAmpMax)
The main function for performing the ringdown attachment for SEOBNRv4 (and beyond) This is the functi...
Definition: LALSimIMREOBHybridRingdownPrec.c:1678
XLALSimIMREOBGetNRSpinPeakOmega
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakOmega(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 a)
Peak frequency predicted by fitting NR results (currently only 2,2 available).
Definition: LALSimIMREOBNQCCorrection.c:1616
XLALSimIMREOBGetNRSpinPeakOmegav2
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakOmegav2(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 a)
Peak frequency predicted by fitting NR results (currently only 2,2 available).
Definition: LALSimIMREOBNQCCorrection.c:1756
XLALSimIMREOBGetNRSpinPeakOmegaV4
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakOmegaV4(INT4 modeL, INT4 modeM, REAL8 UNUSED eta, REAL8 a)
Peak frequency predicted by fitting NR results The coefficients for the (2,2) are in Eq.
Definition: LALSimIMREOBNQCCorrection.c:1336
fprintf
#define fprintf
XLAL_CALLGSL
#define XLAL_CALLGSL(statement)
s
int s
Definition: bh_qnmode.c:137
l
int l
Definition: bh_qnmode.c:135
dt
double dt
Definition: bh_ringdown.c:113
i
double i
Definition: bh_ringdown.c:118
data
sigmaKerr data[0]
Definition: exact_derivatives-HrealHeader.c:4
XLALCreateREAL8VectorSequence
REAL8VectorSequence * XLALCreateREAL8VectorSequence(UINT4 length, UINT4 veclen)
XLALDestroyREAL8VectorSequence
void XLALDestroyREAL8VectorSequence(REAL8VectorSequence *vecseq)
LAL_PI
#define LAL_PI
LAL_TWOPI
#define LAL_TWOPI
LAL_MTSUN_SI
#define LAL_MTSUN_SI
COMPLEX16
double complex COMPLEX16
REAL8
double REAL8
UINT4
uint32_t UINT4
INT4
int32_t INT4
EOB_RD_EFOLDS
#define EOB_RD_EFOLDS
The number of e-folds of ringdown which should be attached for EOBNR models.
Definition: LALSimIMR.h:71
Approximant
Approximant
Enum that specifies the PN approximant to be used in computing the waveform.
Definition: LALSimInspiral.h:229
SEOBNRv2
@ SEOBNRv2
Spin-aligned EOBNR model v2.
Definition: LALSimInspiral.h:326
SEOBNRv1
@ SEOBNRv1
Spin-aligned EOBNR model.
Definition: LALSimInspiral.h:324
SEOBNRv3
@ SEOBNRv3
Spin precessing EOBNR model v3.
Definition: LALSimInspiral.h:330
SEOBNRv3_opt
@ SEOBNRv3_opt
Optimized Spin precessing EOBNR model v3.
Definition: LALSimInspiral.h:334
m
static const INT4 m
a
static const INT4 a
XLALCreateCOMPLEX16Vector
COMPLEX16Vector * XLALCreateCOMPLEX16Vector(UINT4 length)
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLALDestroyCOMPLEX16Vector
void XLALDestroyCOMPLEX16Vector(COMPLEX16Vector *vector)
XLAL_PRINT_INFO
#define XLAL_PRINT_INFO(...)
XLAL_ERROR
#define XLAL_ERROR(...)
XLALPrintError
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_PRINT_ERROR
#define XLAL_PRINT_ERROR(...)
XLAL_EBADLEN
XLAL_EBADLEN
XLAL_ENOMEM
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EFUNC
XLAL_EFAILED
XLAL_EFAILED
XLAL_FAILURE
XLAL_FAILURE
x
list x
p
list p
y
list y
now
def now
filename
filename
test_gwsignal.approximant
string approximant
Definition: test_gwsignal.py:64
COMPLEX16Vector::data
COMPLEX16 * data
REAL8Vector::data
REAL8 * data