Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInspiral 5.0.3.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALPSpinInspiralRingdownWave.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2010 R. Sturani, adapted from LALInspiralRingdownWave.c
3*
4* This program is free software; you can redistribute it and/or modify
5* it under the terms of the GNU General Public License as published by
6* the Free Software Foundation; either version 2 of the License, or
7* (at your option) any later version.
8*
9* This program is distributed in the hope that it will be useful,
10* but WITHOUT ANY WARRANTY; without even the implied warranty of
11* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12* GNU General Public License for more details.
13*
14* You should have received a copy of the GNU General Public License
15* along with with program; see the file COPYING. If not, write to the
16* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17* MA 02110-1301 USA
18*/
19
20
21/**
22 * \file
23 *
24 * \brief Module to compute the ring-down waveform as linear combination
25 * of quasi-normal-modes decaying waveforms, which can be attached to
26 * the phenomenological spin Taylor waveform.
27 *
28 * ### Prototypes ###
29 *
30 * <tt>XLALXLALPSpinInspiralRingdownWave()</tt>
31 * <ul>
32 * <li> <tt>rdwave,</tt> Output, the ring-down waveform
33 * </li><li> <tt>params,</tt> Input, the parameters where ring-down waveforms are computed
34 * </li><li> <tt>inspwave,</tt> Input, the inspiral waveform with given multiple
35 * </li><li> <tt>modefreqs,</tt> Input, the frequencies of the quasi-normal-modes
36 * </li><li> <tt>nmodes,</tt> Input, the number of quasi-normal-modes to be combined.
37 * </li></ul>
38 *
39 * <tt>XLALGenerateWaveDerivative()</tt>
40 * <ul>
41 * <li> <tt>dwave,</tt> Output, time derivative of the input waveform
42 * </li><li> <tt>wave,</tt> Input, waveform to be differentiated in time
43 * </li><li> <tt>params,</tt> Input, the parameters of the input waveform.
44 * </li></ul>
45 *
46 * <tt>XLALPSpinGenerateQNMFreq()</tt>
47 * <ul>
48 * <li> <tt>ptfwave,</tt> Output, the frequencies of the quasi-normal-modes
49 * </li><li> <tt>params,</tt> Input, the parameters of the binary system
50 * </li><li> <tt>l,</tt> Input, the l of the modes
51 * </li><li> <tt>m,</tt> Input, the m of the modes
52 * </li><li> <tt>nmodes,</tt> Input, the number of overtones considered.
53 * </li></ul>
54 *
55 * <tt>XLALPSpinFinalMassSpin()</tt>
56 * <ul>
57 * <li> <tt>finalMass,</tt> Output, the mass of the final Kerr black hole
58 * </li><li> <tt>finalSpin,</tt> Input, the spin of the final Kerr balck hole
59 * </li><li> <tt>params,</tt> Input, the parameters of the binary system.
60 * </li><li> <tt>energy,</tt> Input, the binding energy at the time final time.
61 * </li></ul>
62 *
63 * <tt>XLALPSpinInspiralAttachRingdownWave()</tt>
64 * <ul>
65 * <li> <tt>sigl,</tt> Output, the waveform filled with ring-down phase
66 * </li> <li> <tt>params,</tt> Input, inspiral parameters
67 * </li> <li> <tt>attpos,</tt> Input, position of the start of the ring-down
68 * </li> <li> <tt>nmodes,</tt> Input, number of ring-down modes
69 * </li> <li> <tt>l,</tt> Input, spherical harmonic l-number of the ring-down mode
70 * </li> <li> <tt>m,</tt> Input, spherical harmonic m-number of the ring-down mode
71 * </li> <li> <tt>finalMass,</tt> Input, estimated final mass of the black hole
72 * </li> <li> <tt>finalSpin,</tt> Input, estimated final spin of the black hole</li>
73 * </ul>
74 *
75 * ### Description ###
76 *
77 * This module generate ring-down waveforms.
78 *
79 * ### Algorithm ###
80 *
81 *
82 * ### Uses ###
83 *
84 * \code
85 * LALMalloc
86 * LALFree
87 * \endcode
88 *
89 * ### Notes ###
90 *
91 */
92
93#include <stdlib.h>
94#include <lal/LALStdlib.h>
95#include <lal/AVFactories.h>
96#include <lal/SeqFactories.h>
97#include <lal/LALInspiralBank.h>
98#include <lal/LALNoiseModels.h>
99#include <lal/LALConstants.h>
100#include <gsl/gsl_linalg.h>
101#include <gsl/gsl_interp.h>
102#include <gsl/gsl_spline.h>
103#include <lal/XLALGSL.h>
104
105#define XLAL_BEGINGSL \
106 { \
107 gsl_error_handler_t *saveGSLErrorHandler_; \
108 saveGSLErrorHandler_ = gsl_set_error_handler_off();
109
110#define XLAL_ENDGSL \
111 gsl_set_error_handler( saveGSLErrorHandler_ ); \
112 }
113
115 REAL8Vector *rdwave,
117 REAL8Vector *matchinspwave,
118 COMPLEX8Vector *modefreqs,
119 UINT4 nmodes
120 )
121
122{
123 /* XLAL error handling */
124 INT4 errcode = XLAL_SUCCESS;
125
126 /* Needed to check GSL return codes */
127 INT4 gslStatus;
128
129 UINT4 i, j;
130
131 /* Sampling rate from input */
132 REAL8 dt;
133 gsl_matrix *coef;
134 gsl_vector *hderivs;
135 gsl_vector *x;
136 gsl_permutation *p;
137 REAL8Vector *modeamps;
138
139 int s;
140 REAL8 tj;
141
142 dt = 1.0 / params -> tSampling;
143
144 if ( modefreqs->length != nmodes )
145 {
147 }
148
149 /* Solving the linear system for QNMs amplitude coefficients using gsl routine */
150 /* Initialize matrices and supporting variables */
151
153 coef = (gsl_matrix *) gsl_matrix_alloc(2 * nmodes, 2 * nmodes);
154 hderivs = (gsl_vector *) gsl_vector_alloc(2 * nmodes);
155 x = (gsl_vector *) gsl_vector_alloc(2 * nmodes);
156 p = (gsl_permutation *) gsl_permutation_alloc(2 * nmodes);
157
158 /* Check all matrices and variables were allocated */
159 if ( !coef || !hderivs || !x || !p )
160 {
161 if (coef) gsl_matrix_free(coef);
162 if (hderivs) gsl_vector_free(hderivs);
163 if (x) gsl_vector_free(x);
164 if (p) gsl_permutation_free(p);
166 }
167
168 /* Define the linear system Ax=y */
169 /* Matrix A (2*nmodes by 2*nmodes) has block symmetry. Define half of A here as "coef" */
170 /* Define y here as "hderivs" */
171
172 j=0;
173 while (j<nmodes) {
174 switch (j) {
175 case 0:
176 for (i = 0; i < nmodes; i++) {
177 gsl_matrix_set(coef, 2*j, i, 1.);
178 gsl_matrix_set(coef, 2*j, i+nmodes, 0.);
179 gsl_matrix_set(coef, 2*j+1, i, -cimagf(modefreqs->data[i]));
180 gsl_matrix_set(coef, 2*j+1, i+nmodes, crealf(modefreqs->data[i]));
181 }
182 break;
183 case 1:
184 for (i = 0; i < nmodes; i++) {
185 gsl_matrix_set(coef, 2*j, i, cimagf(modefreqs->data[i])*cimagf(modefreqs->data[i])-crealf(modefreqs->data[i])*crealf(modefreqs->data[i]));
186 gsl_matrix_set(coef, 2*j, i+nmodes, -2.*cimagf(modefreqs->data[i])*crealf(modefreqs->data[i]));
187 gsl_matrix_set(coef, 2*j+1, i, -cimagf(modefreqs->data[i])*cimagf(modefreqs->data[i])*cimagf(modefreqs->data[i])+3.*cimagf(modefreqs->data[i])*crealf(modefreqs->data[i])*crealf(modefreqs->data[i]));
188 gsl_matrix_set(coef, 2*j+1, i+nmodes, -crealf(modefreqs->data[i])*crealf(modefreqs->data[i])*crealf(modefreqs->data[i])+3.*crealf(modefreqs->data[i])*cimagf(modefreqs->data[i])*cimagf(modefreqs->data[i]));
189 }
190 break;
191 case 2:
192 for (i = 0; i < nmodes; i++) {
193 gsl_matrix_set(coef, 2*j, i, pow(cimagf(modefreqs->data[i]),4.)+pow(crealf(modefreqs->data[i]),4.)-6.*pow(crealf(modefreqs->data[i])*cimagf(modefreqs->data[i]),2.));
194 gsl_matrix_set(coef, 2*j, i+nmodes, -4.*pow(cimagf(modefreqs->data[i]),3.)*crealf(modefreqs->data[i])+4.*pow(crealf(modefreqs->data[i]),3.)*cimagf(modefreqs->data[i]));
195 gsl_matrix_set(coef, 2*j+1, i, -pow(cimagf(modefreqs->data[i]),5.)+10.*pow(cimagf(modefreqs->data[i]),3.)*pow(crealf(modefreqs->data[i]),2.)-5.*cimagf(modefreqs->data[i])*pow(crealf(modefreqs->data[i]),4.));
196 gsl_matrix_set(coef, 2*j+1, i+nmodes, 5.*pow(cimagf(modefreqs->data[i]),4.)*crealf(modefreqs->data[i])-10.*pow(cimagf(modefreqs->data[i]),2.)*pow(crealf(modefreqs->data[i]),3.)+pow(crealf(modefreqs->data[i]),5.));
197 }
198 break;
199 default:
200 XLALPrintError("*** LALPSpinInspiralRingDown ERROR ***: nmode must be <=2, %d selected\n",nmodes);
201 gsl_matrix_free(coef);
202 gsl_vector_free(hderivs);
203 gsl_vector_free(x);
204 gsl_permutation_free(p);
206 }
207 gsl_vector_set(hderivs, 2*j, matchinspwave->data[2*j]);
208 gsl_vector_set(hderivs, 2*j+1, matchinspwave->data[2*j+1]);
209 j++;
210 }
211
212 /* Call gsl LU decomposition to solve the linear system */
213 gslStatus = gsl_linalg_LU_decomp(coef, p, &s);
214 if ( gslStatus == GSL_SUCCESS )
215 {
216 gslStatus = gsl_linalg_LU_solve(coef, p, hderivs, x);
217 }
218
219 if ( gslStatus != GSL_SUCCESS )
220 {
221 gsl_matrix_free(coef);
222 gsl_vector_free(hderivs);
223 gsl_vector_free(x);
224 gsl_permutation_free(p);
226 }
227
228 /* Putting solution to an XLAL vector */
229 modeamps = XLALCreateREAL8Vector(2*nmodes);
230
231 if ( !modeamps )
232 {
233 gsl_matrix_free(coef);
234 gsl_vector_free(hderivs);
235 gsl_vector_free(x);
236 gsl_permutation_free(p);
238 }
239
240 for (i = 0; i < 2*nmodes; i++) {
241 modeamps->data[i] = gsl_vector_get(x, i);
242 }
243
244 /* Free all gsl linear algebra objects */
245 gsl_matrix_free(coef);
246 gsl_vector_free(hderivs);
247 gsl_vector_free(x);
248 gsl_permutation_free(p);
250
251 /* Build ring-down waveforms */
252 UINT4 Nrdwave=rdwave->length;
253 for (j = 0; j < Nrdwave; j++) {
254 tj = j * dt;
255 rdwave->data[j] = 0.;
256 for (i = 0; i < nmodes; i++) {
257 rdwave->data[j] += exp(- tj * cimagf(modefreqs->data[i]))
258 * ( modeamps->data[i] * cos(tj * crealf(modefreqs->data[i]))
259 + modeamps->data[i + nmodes] * sin(tj * crealf(modefreqs->data[i])) );
260 }
261 }
262
263 XLALDestroyREAL8Vector(modeamps);
264 return errcode;
265} /*End of XLALPSpinInspiralRingdownWave */
266
267
269 REAL8Vector *dwave,
270 REAL8Vector *wave,
271 REAL8 dt
272 )
273{
274 /* XLAL error handling */
275 INT4 errcode = XLAL_SUCCESS;
276
277 /* For checking GSL return codes */
278 INT4 gslStatus;
279
280 UINT4 j;
281 double *x, *y;
282 double dy;
283 gsl_interp_accel *acc;
284 gsl_spline *spline;
285
286 if (wave->length!=dwave->length)
288
289 /* Getting interpolation and derivatives of the waveform using gsl spline routine */
290 /* Initialize arrays and supporting variables for gsl */
291
292 x = (double *) LALMalloc(wave->length * sizeof(double));
293 y = (double *) LALMalloc(wave->length * sizeof(double));
294
295 if ( !x || !y )
296 {
297 if ( x ) LALFree (x);
298 if ( y ) LALFree (y);
300 }
301
302 for (j = 0; j < wave->length; ++j)
303 {
304 x[j] = j;
305 y[j] = wave->data[j];
306 }
307
308 XLAL_CALLGSL( acc = (gsl_interp_accel*) gsl_interp_accel_alloc() );
309 XLAL_CALLGSL( spline = (gsl_spline*) gsl_spline_alloc(gsl_interp_cspline, wave->length) );
310 if ( !acc || !spline )
311 {
312 if ( acc ) gsl_interp_accel_free(acc);
313 if ( spline ) gsl_spline_free(spline);
314 LALFree( x );
315 LALFree( y );
317 }
318
319 /* Gall gsl spline interpolation */
320 XLAL_CALLGSL( gslStatus = gsl_spline_init(spline, x, y, wave->length) );
321 if ( gslStatus != GSL_SUCCESS )
322 {
323 gsl_spline_free(spline);
324 gsl_interp_accel_free(acc);
325 LALFree( x );
326 LALFree( y );
328 }
329
330 /* Getting first and second order time derivatives from gsl interpolations */
331 for (j = 0; j < wave->length; ++j)
332 {
333 XLAL_CALLGSL(gslStatus = gsl_spline_eval_deriv_e( spline, j, acc, &dy ) );
334 if (gslStatus != GSL_SUCCESS )
335 {
336 gsl_spline_free(spline);
337 gsl_interp_accel_free(acc);
338 LALFree( x );
339 LALFree( y );
341 }
342 dwave->data[j] = (REAL8)(dy / dt);
343
344 }
345
346
347 /* Free gsl variables */
348 gsl_spline_free(spline);
349 gsl_interp_accel_free(acc);
350 LALFree(x);
351 LALFree(y);
352
353 return errcode;
354}
355
357 COMPLEX8Vector *modefreqs,
359 UINT4 l,
360 INT4 m,
361 UINT4 nmodes,
362 REAL8 finalMass,
363 REAL8 finalSpin
364 )
365
366{
367 /* XLAL error handling */
368 INT4 errcode = XLAL_SUCCESS;
369 UINT4 i;
370 REAL8 totalMass;
371 /* Fitting coefficients for QNM frequencies from PRD73, 064030, gr-qc/0512160, tables VIII and IX */
372 REAL4 BCW22re[3][3] = { {1.5251, -1.1568, 0.1292}, {1.3673, -1.0260, 0.1628}, { 1.3223, -1.0257, 0.1860} };
373 REAL4 BCW22im[3][3] = { {0.7000, 1.4187, -0.4990}, {0.1000, 0.5436, -0.4731}, {-0.1000, 0.4206, -0.4256} };
374
375 /*REAL4 BCW2m2re[3][3] = { {0.2938, 0.0782, 1.3546}, {0.2528, 0.0921, 1.3344}, { 0.1873, 0.1117, 1.3322} };
376 REAL4 BCW2m2im[3][3] = { {1.6700, 0.4192, 1.4700}, {0.4550, 0.1729, 1.3617}, { 0.1850, 0.1266, 1.3661} };*/
377
378 REAL4 BCW21re[3][3] = { {0.60000, -0.2339, 0.4175}, {0.5800, -0.2416, 0.4708}, { 0.5660, -0.2740, 0.4960} };
379 REAL4 BCW21im[3][3] = { {-0.30000, 2.3561, -0.2277}, {-0.3300, 0.9501, -0.2072}, { -0.1000, 0.4173, -0.2774} };
380
381 /*REAL4 BCW2m1re[3][3] = { {0.3441, 0.0293, 2.0010}, {0.3165, 0.0301, 2.3415}, {0.2696, 0.0315, 2.7755} };
382 REAL4 BCW2m1im[3][3] = { {2.0000, 0.1078, 5.0069}, {0.6100, 0.0276, 13.1683}, {0.2900, 0.0276, 6.4715} };*/
383
384 REAL4 BCW20re[3][3] = { {0.4437, -0.0739, 0.3350}, {0.4185, -0.0768, 0.4355}, { 0.3734, -0.0794, 0.6306} };
385 REAL4 BCW20im[3][3] = { {4.0000, -1.9550, 0.1420}, {1.2500, -0.6359, 0.1614}, {0.5600, -0.2589, -0.3034} };
386
387 REAL4 BCW33re[3][3] = { {1.8596, -1.3043, 0.1818}, {1.8566, -1.2818, 0.1934}, {1.8004, -1.2558, 0.2133} };
388 REAL4 BCW33im[3][3] = { {0.9000, 2.3430, -0.4810}, {0.2274, 0.8173, -0.4731}, {0.0400, 0.5445, -0.4539} };
389
390 /*REAL4 BCW3m3re[3][3] = { {0.4673, 0.1296, 1.3255}, {0.4413, 0.1387, 1.3178}, {0.3933, 0.1555, 1.3037} };
391 REAL4 BCW3m3im[3][3] = { {2.5500, 0.6576, 1.3378}, {0.7900, 0.2381, 1.3706}, {0.4070, 0.1637, 1.3819} };*/
392
393 REAL4 BCW32re[3][3] = { {1.1481, -0.5552, 0.3002}, {1.1226, -0.5471, 0.3264}, {1.0989, -0.5550, 0.3569} };
394 REAL4 BCW32im[3][3] = { {0.8313, 2.3773, -0.3655}, {0.2300, 0.8025, -0.3684}, {0.1000, 0.4804, -0.3784}};
395
396 /*REAL4 BCW3m2re[3][3] = { {0.5158, 0.8195, 1.408}, {0.4413, 0.1378, 1.3178}, {0.4567, 0.09300, 1.4469} };
397 REAL4 BCW3m2im[3][3] = { {2.9000, 0.3365, 2.3050}, {0.9000, 0.1295, 1.6142}, {0.4900, 0.0848, 1.9737} };*/
398
399 REAL4 BCW31re[3][3] = { {0.8345, -0.2405, 0.4095}, {0.8105, -0.2342, 0.4660}, {0.7684, -0.2252, 0.5805} };
400 REAL4 BCW31im[3][3] = { {23.8450, -20.724, 0.03837}, {8.8530, -7.8506, 0.03418}, {2.1800, -1.6273, 0.1163} };
401
402 /*REAL4 BCW3m1re[3][3] = { {0.5751, 0.02508, 3.1360}, {0.5584, 0.02514, 3.4154}, {0.5271, 0.02561, 3.8011} };
403 REAL4 BCW3m1im[3][3] = { {3.0464, 0.1162, -0.2812}, {1.2000, -0.1928, 0.1037}, {1.0000, -0.4424, 0.02467} };*/
404
405 REAL4 BCW30re[3][3] = { {0.6873, -0.09282, 0.3479}, {0.6687, -0.09155, 0.4021}, {0.6343, -0.08915, 0.5117} };
406 REAL4 BCW30im[3][3] = { {6.7841, -3.6112, 0.09480}, {2.0075, -0.9930, 0.1197}, {0.9000, -0.3409, 0.2679} };
407
408 REAL4 BCW44re[3][3] = { {2.3, -1.5056, 0.2244}, {2.3, -1.5173, 0.2271}, {2.3, -1.5397, 0.2321} };
409 REAL4 BCW44im[3][3] = { {1.1929, 3.1191, -0.4825}, {0.3, 1.1034, -0.4703}, {0.11, 0.6997, -0.4607} };
410
411 /*REAL4 BCW4m4re[3][3] = { {0.6256, 0.18, 1.3218}, {0.6061, 0.1869, 1.3168}, {0.5686, 0.2003, 1.3068} };
412 REAL4 BCW4m4im[3][3] = { {3.4, 0.8696, 1.4074}, {1.08, 0.3095, 1.3279}, {0.5980, 0.2015, 1.3765} };*/
413
414 REAL4 BCW43re[3][3] = { {1.6869, -0.8862, 0.2822}, {1.6722, -0.8843, 0.2923}, {1.6526, -0.8888, 0.3081} };
415 REAL4 BCW43im[3][3] = { {1.4812, 2.8096, -0.4271}, {0.4451, 0.9569, -0.425}, {0.22, 0.5904, -0.4236} };
416
417 /*REAL4 BCW4m3re[3][3] = { {0.6728, 0.1338, 1.3413}, {0.6562, 0.1377, 1.3456}, {0.6244, 0.1454, 1.3513} };
418 REAL4 BCW4m3im[3][3] = { {3.7, 0.5829, 1.6681}, {1.18, 0.2111, 1.4129}, {0.66, 0.1385, 1.3742} };*/
419
420 REAL4 BCW42re[3][3] = { {1.2702, -0.4685, 0.3835}, {1.2462, -0.4580, 0.4139}, {1.2025, -0.4401, 0.4769} };
421 REAL4 BCW42im[3][3] = { {-3.6, 7.7749, -0.1491}, {-1.5, 2.8601, -0.1392}, {-1.5, 2.2784, -0.1124}};
422
423 /*REAL4 BCW4m2re[3][3] = { {0.7294, 0.07842, 1.5646}, {0.7154, 0.07979, 1.5852}, {0.6885, 0.08259, 1.6136} };
424 REAL4 BCW4m2im[3][3] = { {4., 0.2777, 2.0647}, {1.32, 0.08694, 4.3255}, {0.75, 0.05803, 3.7971} };*/
425
426 REAL4 BCW41re[3][3] = { {1.0507, -0.2478, 0.4348}, {1.0337, -0.2439, 0.4695}, {1.0019, -0.2374, 0.5397} };
427 REAL4 BCW41im[3][3] = { {14., -9.8240, 0.09047}, {4.2, -2.8399, 0.1081}, {2.2, -1.4195, 0.1372} };
428
429 /*REAL4 BCW4m1re[3][3] = { {0.7908, 0.02024, 5.4628}, {0.7785, 0.02005, 5.8547}, {0.7549, 0.01985, 6.5272} };
430 REAL4 BCW4m1im[3][3] = { {4.6, -0.4038, 0.4629}, {1.6, -0.2323, 0.2306}, {1.6, -0.8136, 0.03163} };*/
431
432 REAL4 BCW40re[3][3] = { {0.9175, -0.1144, 0.3511}, {0.9028, -0.1127, 0.3843}, {0.8751, -0.1096, 0.4516} };
433 REAL4 BCW40im[3][3] = { {7.0, -2.7934, 0.1708}, {2.2, -0.8308, 0.2023}, {1.2, -0.4159, 0.2687} };
434
435
436 /* Get a local copy of the intrinstic parameters */
437 totalMass = params->totalMass;
438
439 /* QNM frequencies from the fitting given in PRD73, 064030 */
440
441 if ((l==2)&&(abs(m)==2)) {
442 for (i = 0; i < nmodes; i++)
443 {
444 REAL8 real_part = BCW22re[i][0] + BCW22re[i][1] * pow(1.- finalSpin, BCW22re[i][2]);
445 REAL8 imag_part = real_part / 2. / (BCW22im[i][0] + BCW22im[i][1] * pow(1.- finalSpin, BCW22im[i][2]));
446 modefreqs->data[i] = crectf(real_part, imag_part);
447 modefreqs->data[i] *= ((REAL4) 1./ finalMass / (totalMass * LAL_MTSUN_SI));
448 }
449 }
450 else {
451 if ((l==2)&&(m==0)) {
452 for (i = 0; i < nmodes; i++)
453 {
454 REAL8 real_part = BCW20re[i][0] + BCW20re[i][1] * pow(1.- finalSpin, BCW20re[i][2]);
455 REAL8 imag_part = real_part / 2. / (BCW20im[i][0] + BCW20im[i][1] * pow(1.- finalSpin, BCW20im[i][2]));
456 modefreqs->data[i] = crectf(real_part, imag_part);
457 modefreqs->data[i] /= ((REAL4) finalMass * totalMass * LAL_MTSUN_SI);
458 }
459 }
460 else {
461 if ((l==2)&&(abs(m)==1)) {
462 for (i = 0; i < nmodes; i++) {
463 REAL8 real_part = BCW21re[i][0] + BCW21re[i][1] * pow(1.- finalSpin, BCW21re[i][2]);
464 REAL8 imag_part = real_part / 2. / (BCW21im[i][0] + BCW21im[i][1] * pow(1.- finalSpin, BCW21im[i][2]));
465 modefreqs->data[i] = crectf(real_part, imag_part);
466 modefreqs->data[i] /= ((REAL4) finalMass * totalMass * LAL_MTSUN_SI);
467 }
468 }
469 else {
470 if ((l==3)&&(abs(m)==3)) {
471 for (i = 0; i < nmodes; i++) {
472 REAL8 real_part = BCW33re[i][0] + BCW33re[i][1] * pow(1.- finalSpin, BCW33re[i][2]);
473 REAL8 imag_part = real_part / 2. / (BCW33im[i][0] + BCW33im[i][1] * pow(1.- finalSpin, BCW33im[i][2]));
474 modefreqs->data[i] = crectf(real_part, imag_part);
475 modefreqs->data[i] /= ((REAL4) finalMass * totalMass * LAL_MTSUN_SI);
476 }
477 }
478 else
479 if ((l==3)&&(abs(m)==2)) {
480 for (i = 0; i < nmodes; i++) {
481 REAL8 real_part = BCW32re[i][0] + BCW32re[i][1] * pow(1.- finalSpin, BCW32re[i][2]);
482 REAL8 imag_part = real_part / 2. / (BCW32im[i][0] + BCW32im[i][1] * pow(1.- finalSpin, BCW32im[i][2]));
483 modefreqs->data[i] = crectf(real_part, imag_part);
484 modefreqs->data[i] /= ((REAL4) finalMass * totalMass * LAL_MTSUN_SI);
485 }
486 }
487 else {
488 if ((l==3)&&(abs(m)==1)) {
489 for (i = 0; i < nmodes; i++) {
490 REAL8 real_part = BCW31re[i][0] + BCW31re[i][1] * pow(1.- finalSpin, BCW31re[i][2]);
491 REAL8 imag_part = real_part / 2. / (BCW31im[i][0] + BCW31im[i][1] * pow(1.- finalSpin, BCW31im[i][2]));
492 modefreqs->data[i] = crectf(real_part, imag_part);
493 modefreqs->data[i] /= ((REAL4) finalMass * totalMass * LAL_MTSUN_SI);
494 }
495 }
496 else {
497 if ((l==3)&&(m==0)) {
498 for (i = 0; i < nmodes; i++) {
499 REAL8 real_part = BCW30re[i][0] + BCW30re[i][1] * pow(1.- finalSpin, BCW30re[i][2]);
500 REAL8 imag_part = real_part / 2. / (BCW30im[i][0] + BCW30im[i][1] * pow(1.- finalSpin, BCW30im[i][2]));
501 modefreqs->data[i] = crectf(real_part, imag_part);
502 modefreqs->data[i] /= ((REAL4) finalMass * totalMass * LAL_MTSUN_SI);
503 }
504 }
505 else {
506 if ((l==4)&&(abs(m)==4)) {
507 for (i = 0; i < nmodes; i++) {
508 REAL8 real_part = BCW44re[i][0] + BCW44re[i][1] * pow(1.- finalSpin, BCW44re[i][2]);
509 REAL8 imag_part = real_part / 2. / (BCW44im[i][0] + BCW44im[i][1] * pow(1.- finalSpin, BCW44im[i][2]));
510 modefreqs->data[i] = crectf(real_part, imag_part);
511 modefreqs->data[i] /= ((REAL4) finalMass * totalMass * LAL_MTSUN_SI);
512 }
513 }
514 else {
515 if ((l==4)&&(abs(m)==3)) {
516 for (i = 0; i < nmodes; i++) {
517 REAL8 real_part = BCW43re[i][0] + BCW43re[i][1] * pow(1.- finalSpin, BCW43re[i][2]);
518 REAL8 imag_part = real_part / 2. / (BCW43im[i][0] + BCW43im[i][1] * pow(1.- finalSpin, BCW43im[i][2]));
519 modefreqs->data[i] = crectf(real_part, imag_part);
520 modefreqs->data[i] /= ((REAL4) finalMass * totalMass * LAL_MTSUN_SI);
521 }
522 }
523 else {
524 if ((l==4)&&(abs(m)==2)) {
525 for (i = 0; i < nmodes; i++) {
526 REAL8 real_part = BCW42re[i][0] + BCW42re[i][1] * pow(1.- finalSpin, BCW42re[i][2]);
527 REAL8 imag_part = real_part / 2. / (BCW42im[i][0] + BCW42im[i][1] * pow(1.- finalSpin, BCW42im[i][2]));
528 modefreqs->data[i] = crectf(real_part, imag_part);
529 modefreqs->data[i] /= ((REAL4) finalMass * totalMass * LAL_MTSUN_SI);
530 }
531 }
532 else {
533 if ((l==4)&&(abs(m)==1)) {
534 for (i = 0; i < nmodes; i++) {
535 REAL8 real_part = BCW41re[i][0] + BCW41re[i][1] * pow(1.- finalSpin, BCW41re[i][2]);
536 REAL8 imag_part = real_part / 2. / (BCW41im[i][0] + BCW41im[i][1] * pow(1.- finalSpin, BCW41im[i][2]));
537 modefreqs->data[i] = crectf(real_part, imag_part);
538 modefreqs->data[i] /= ((REAL4) finalMass * totalMass * LAL_MTSUN_SI);
539 }
540 }
541 else {
542 if ((l==4)&&(m==0)) {
543 for (i = 0; i < nmodes; i++) {
544 REAL8 real_part = BCW40re[i][0] + BCW40re[i][1] * pow(1.- finalSpin, BCW40re[i][2]);
545 REAL8 imag_part = real_part / 2. / (BCW40im[i][0] + BCW40im[i][1] * pow(1.- finalSpin, BCW40im[i][2]));
546 modefreqs->data[i] = crectf(real_part, imag_part);
547 modefreqs->data[i] /= ((REAL4) finalMass * totalMass * LAL_MTSUN_SI);
548 }
549 }
550 else {
551 fprintf(stderr,"*** LALPSpinInspiralRingdownWave ERROR: Ringdown modes for l=%d m=%d not availbale\n",l,m);
553 }
554 }
555 }
556 }
557 }
558 }
559 }
560 }
561 }
562 }
563 }
564
565 return errcode;
566}
567
569 REAL8 *finalMass,
570 REAL8 *finalSpin,
572 REAL8 energy,
573 REAL8 *LNhvec
574 )
575{
576 /* XLAL error handling */
577 INT4 errcode = XLAL_SUCCESS;
578 REAL8 qq,ll,eta;
579
580 /* See eq.(6) in arXiv:0904.2577 */
581 REAL8 ma1,ma2,a12,a12l;
582 REAL8 cosa1=0.;
583 REAL8 cosa2=0.;
584 REAL8 cosa12=0.;
585
586 REAL8 t0=-2.9;
587 REAL8 t3=2.6;
588 REAL8 s4=-0.123;
589 REAL8 s5=0.45;
590 REAL8 t2=16.*(0.6865-t3/64.-sqrt(3.)/2.);
591
592 /* get a local copy of the intrinstic parameters */
593 qq=params->mass2/params->mass1;
594 eta = params->eta;
595 /* done */
596 ma1=sqrt( params->spin1[0]*params->spin1[0] + params->spin1[1]*params->spin1[1] + params->spin1[2]*params->spin1[2] );
597 ma2=sqrt( params->spin2[0]*params->spin2[0] + params->spin2[1]*params->spin2[1] + params->spin2[2]*params->spin2[2] );
598
599 if (ma1>0.) cosa1 = (params->spin1[0]*LNhvec[0]+params->spin1[1]*LNhvec[1]+params->spin1[2]*LNhvec[2])/ma1;
600 else cosa1=0.;
601 if (ma2>0.) cosa2 = (params->spin2[0]*LNhvec[0]+params->spin2[1]*LNhvec[1]+params->spin2[2]*LNhvec[2])/ma2;
602 else cosa2=0.;
603 if ((ma1>0.)&&(ma2>0.)) {
604 cosa12 = (params->spin1[0]*params->spin2[0] + params->spin1[1]*params->spin2[1] + params->spin1[2]*params->spin2[2])/ma1/ma2;
605 }
606 else cosa12=0.;
607
608 a12 = ma1*ma1 + ma2*ma2*qq*qq*qq*qq + 2.*ma1*ma2*qq*qq*cosa12 ;
609 a12l = ma1*cosa1 + ma2*cosa2*qq*qq ;
610 ll = 2.*sqrt(3.)+ t2*eta + t3*eta*eta + s4*a12/(1.+qq*qq)/(1.+qq*qq) + (s5*eta+t0+2.)/(1.+qq*qq)*a12l;
611
612 /* Estimate final mass by adding the negative binding energy to the rest mass*/
613 *finalMass = 1. + energy;
614 if (*finalMass < 0.) {
615 fprintf(stderr,"*** LALPSpinInspiralRingdownWave ERROR: Estimated final mass <0 : %12.6f\n ",*finalMass);
616 fprintf(stderr,"*** Final mass set to initial mass\n");
618 *finalMass = 1.;
619 }
620
621 *finalSpin = sqrt( a12 + 2.*ll*qq*a12l + ll*ll*qq*qq)/(1.+qq)/(1.+qq);
622 if ((*finalSpin > 1.)||(*finalSpin < 0.)) {
623 if ((*finalSpin>=1.)&&(*finalSpin<1.01)) {
624 fprintf(stderr,"*** LALPSpinInspiralRingdownWave WARNING: Estimated final Spin slightly >1 : %11.3e\n ",*finalSpin);
625 fprintf(stderr," (m1=%8.3f m2=%8.3f s1=(%8.3f,%8.3f,%8.3f) s2=(%8.3f,%8.3f,%8.3f) ) final spin set to 1 and code goes on\n",params->mass1,params->mass2,params->spin1[0],params->spin1[1],params->spin1[2],params->spin2[0],params->spin2[1],params->spin2[2]);
626 *finalSpin = .99999;
627 }
628 else {
629 fprintf(stderr,"*** LALPSpinInspiralRingdownWave ERROR: Unphysical estimation of final Spin : %11.3e\n ",*finalSpin);
630 fprintf(stderr," (m1=%8.3f m2=%8.3f s1=(%8.3f,%8.3f,%8.3f) s2=(%8.3f,%8.3f,%8.3f) )\n",params->mass1,params->mass2,params->spin1[0],params->spin1[1],params->spin1[2],params->spin2[0],params->spin2[1],params->spin2[2]);
631 fprintf(stderr,"*** Code aborts\n");
632 *finalSpin = 0.;
634 }
635 }
636
637 /*For reference these are the formula used in the EOBNR construction*/
638 //*finalMass = 1. - 0.057191 * eta - 0.498 * eta*eta;
639 //*finalSpin = 3.464102 * eta - 2.9 * eta*eta;
640
641 return errcode;
642}
643
645 REAL8Vector *sigl,
647 UINT4 *attpos,
648 UINT4 nmodes,
649 UINT4 l,
650 INT4 m,
651 REAL8 finalMass,
652 REAL8 finalSpin
653 )
654{
655 const UINT4 Npatch=40;
656 const UINT4 offsetAttch = 2;
657
658 COMPLEX8Vector *modefreqs;
659 UINT4 Nrdwave;
660
661 UINT4 i=0;
662 UINT4 j=0;
663 UINT4 k=0;
664 UINT4 atpos;
665 INT4 errcode;
666
667 REAL8Vector *rdwave;
668 REAL8Vector *inspwave,*dinspwave;
669 REAL8Vector *matchinspwave;
670 REAL8 dt;
671
672 dt = 1./params->tSampling;
673 atpos=(*attpos);
674
675 /* Create memory for the QNM frequencies */
676 modefreqs = XLALCreateCOMPLEX8Vector( nmodes );
677 if ( !modefreqs )
678 {
680 }
681 errcode = XLALPSpinGenerateQNMFreq( modefreqs, params, l, m, nmodes, finalMass, finalSpin);
682 if ( errcode != XLAL_SUCCESS )
683 {
684 XLALDestroyCOMPLEX8Vector( modefreqs );
686 }
687
688 /* Ringdown signal length: 10 times the decay time of the n=0 mode */
689 Nrdwave = (INT4) (10. / cimagf(modefreqs->data[0]) / dt);
690 /* Patch length, centered around the matching point "attpos" */
691
692 (*attpos)+=Nrdwave;
693
694 /* Check the value of attpos, to prevent memory access problems later */
695 if ( atpos < Npatch || atpos + Npatch >= sigl->length )
696 {
697 XLALPrintError( "Value of attpos inconsistent with given value of Npatch: atpos=%d Npatch=%d, sign->length=%d, m1=%11.5f m2=%11.5f s1z=%8.3f s2z=%8.3f fL=%11.3e\n",atpos,Npatch,sigl->length,params->mass1,params->mass2,params->spin1[2],params->spin2[2],params->fLower);
698 XLALDestroyCOMPLEX8Vector( modefreqs );
700 }
701
702 /* Create memory for the ring-down and full waveforms, derivatives of inspirals
703 and waveforms and its derivative values at the attach point */
704
705 rdwave = XLALCreateREAL8Vector( Nrdwave );
706 inspwave = XLALCreateREAL8Vector( Npatch );
707 dinspwave = XLALCreateREAL8Vector( Npatch );
708 matchinspwave = XLALCreateREAL8Vector( 2*nmodes );
709
710 /* Check memory was allocated */
711 if ( !rdwave || !inspwave || !dinspwave || !matchinspwave )
712 {
713 XLALDestroyCOMPLEX8Vector( modefreqs );
714 if (rdwave) XLALDestroyREAL8Vector( rdwave );
715 if (inspwave) XLALDestroyREAL8Vector( inspwave );
716 if (dinspwave) XLALDestroyREAL8Vector( dinspwave );
717 if (matchinspwave) XLALDestroyREAL8Vector( matchinspwave );
719 }
720
721 /* Generate derivatives of the last part of inspiral waves */
722 /* Take the last part of sigl1 */
723
724 for (i=0; i<2; i++) {
725 /* i=0(1) for real(imaginary) part */
726 for (j = 0; j < Npatch; j++) {
727 inspwave->data[j] = sigl->data[2*(atpos - Npatch + j)+i];
728 }
729
730 for (k=0;k<2*nmodes;k++) {
731 matchinspwave->data[k] = inspwave->data[Npatch-1-offsetAttch];
732 if ((k+1)<2*nmodes) {
733 errcode = XLALGenerateWaveDerivative( dinspwave, inspwave, dt);
734 if ( (errcode != XLAL_SUCCESS) ) {
735 XLALDestroyCOMPLEX8Vector( modefreqs );
736 XLALDestroyREAL8Vector( rdwave );
737 XLALDestroyREAL8Vector( inspwave );
738 XLALDestroyREAL8Vector( dinspwave );
739 XLALDestroyREAL8Vector( matchinspwave );
741 }
742 for (j=0; j<Npatch; j++) {
743 inspwave->data[j]=dinspwave->data[j];
744 }
745 }
746 }
747
748 errcode = XLALPSpinInspiralRingdownWave( rdwave, params, matchinspwave, modefreqs, nmodes );
749
750 if ( errcode != XLAL_SUCCESS ) {
751 XLALDestroyCOMPLEX8Vector( modefreqs );
752 XLALDestroyREAL8Vector( rdwave );
753 XLALDestroyREAL8Vector( inspwave );
754 XLALDestroyREAL8Vector( dinspwave );
755 XLALDestroyREAL8Vector( matchinspwave );
757 }
758 /* Generate full waveforms, by stitching inspiral and ring-down waveforms */
759
760 for (j = 0; j < Nrdwave; j++) {
761 sigl->data[2*j + 2*(atpos - 1 - offsetAttch) + i ] = rdwave->data[j];
762 }
763
764 }
765
766 /* Free memory */
767 XLALDestroyCOMPLEX8Vector( modefreqs );
768 XLALDestroyREAL8Vector( rdwave );
769 XLALDestroyREAL8Vector( inspwave );
770 XLALDestroyREAL8Vector( dinspwave );
771 XLALDestroyREAL8Vector( matchinspwave );
772
773 return errcode;
774}
#define LALMalloc(n)
#define LALFree(p)
INT4 XLALPSpinGenerateQNMFreq(COMPLEX8Vector *modefreqs, InspiralTemplate *params, UINT4 l, INT4 m, UINT4 nmodes, REAL8 finalMass, REAL8 finalSpin)
#define XLAL_ENDGSL
INT4 XLALPSpinFinalMassSpin(REAL8 *finalMass, REAL8 *finalSpin, InspiralTemplate *params, REAL8 energy, REAL8 *LNhvec)
INT4 XLALGenerateWaveDerivative(REAL8Vector *dwave, REAL8Vector *wave, REAL8 dt)
#define XLAL_BEGINGSL
INT4 XLALPSpinInspiralAttachRingdownWave(REAL8Vector *sigl, InspiralTemplate *params, UINT4 *attpos, UINT4 nmodes, UINT4 l, INT4 m, REAL8 finalMass, REAL8 finalSpin)
INT4 XLALPSpinInspiralRingdownWave(REAL8Vector *rdwave, InspiralTemplate *params, REAL8Vector *matchinspwave, COMPLEX8Vector *modefreqs, UINT4 nmodes)
#define fprintf
#define XLAL_CALLGSL(statement)
int s
int l
double i
const double qq
#define LAL_MTSUN_SI
double REAL8
#define crectf(re, im)
uint32_t UINT4
int32_t INT4
float REAL4
static const INT4 m
COMPLEX8Vector * XLALCreateCOMPLEX8Vector(UINT4 length)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
void XLALDestroyCOMPLEX8Vector(COMPLEX8Vector *vector)
#define XLAL_ERROR(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_EBADLEN
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_ERANGE
XLAL_EFUNC
XLAL_EDOM
XLAL_EFAILED
list y
p
x
double t0
COMPLEX8 * data
The inspiral waveform parameter structure containing information about the waveform to be generated.
Definition: LALInspiral.h:205
REAL8 * data