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>
105#define XLAL_BEGINGSL \
107 gsl_error_handler_t *saveGSLErrorHandler_; \
108 saveGSLErrorHandler_ = gsl_set_error_handler_off();
111 gsl_set_error_handler( saveGSLErrorHandler_ ); \
144 if ( modefreqs->
length != nmodes )
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);
159 if ( !coef || !hderivs || !
x || !
p )
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);
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]));
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]));
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.));
200 XLALPrintError(
"*** LALPSpinInspiralRingDown ERROR ***: nmode must be <=2, %d selected\n",nmodes);
201 gsl_matrix_free(coef);
202 gsl_vector_free(hderivs);
204 gsl_permutation_free(
p);
207 gsl_vector_set(hderivs, 2*j, matchinspwave->
data[2*j]);
208 gsl_vector_set(hderivs, 2*j+1, matchinspwave->
data[2*j+1]);
213 gslStatus = gsl_linalg_LU_decomp(coef,
p, &
s);
214 if ( gslStatus == GSL_SUCCESS )
216 gslStatus = gsl_linalg_LU_solve(coef,
p, hderivs,
x);
219 if ( gslStatus != GSL_SUCCESS )
221 gsl_matrix_free(coef);
222 gsl_vector_free(hderivs);
224 gsl_permutation_free(
p);
233 gsl_matrix_free(coef);
234 gsl_vector_free(hderivs);
236 gsl_permutation_free(
p);
240 for (
i = 0;
i < 2*nmodes;
i++) {
241 modeamps->
data[
i] = gsl_vector_get(
x,
i);
245 gsl_matrix_free(coef);
246 gsl_vector_free(hderivs);
248 gsl_permutation_free(
p);
253 for (j = 0; j < Nrdwave; j++) {
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])) );
283 gsl_interp_accel *acc;
302 for (j = 0; j < wave->
length; ++j)
305 y[j] = wave->
data[j];
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 )
312 if ( acc ) gsl_interp_accel_free(acc);
313 if ( spline ) gsl_spline_free(spline);
321 if ( gslStatus != GSL_SUCCESS )
323 gsl_spline_free(spline);
324 gsl_interp_accel_free(acc);
331 for (j = 0; j < wave->
length; ++j)
333 XLAL_CALLGSL(gslStatus = gsl_spline_eval_deriv_e( spline, j, acc, &dy ) );
334 if (gslStatus != GSL_SUCCESS )
336 gsl_spline_free(spline);
337 gsl_interp_accel_free(acc);
348 gsl_spline_free(spline);
349 gsl_interp_accel_free(acc);
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} };
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} };
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} };
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} };
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}};
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} };
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} };
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} };
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} };
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}};
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} };
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} };
437 totalMass =
params->totalMass;
441 if ((
l==2)&&(abs(
m)==2)) {
442 for (
i = 0;
i < nmodes;
i++)
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]));
451 if ((
l==2)&&(
m==0)) {
452 for (
i = 0;
i < nmodes;
i++)
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]));
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]));
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]));
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]));
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]));
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]));
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]));
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]));
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]));
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]));
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]));
551 fprintf(stderr,
"*** LALPSpinInspiralRingdownWave ERROR: Ringdown modes for l=%d m=%d not availbale\n",
l,
m);
581 REAL8 ma1,ma2,a12,a12l;
590 REAL8 t2=16.*(0.6865-t3/64.-sqrt(3.)/2.);
599 if (ma1>0.) cosa1 = (
params->spin1[0]*LNhvec[0]+
params->spin1[1]*LNhvec[1]+
params->spin1[2]*LNhvec[2])/ma1;
601 if (ma2>0.) cosa2 = (
params->spin2[0]*LNhvec[0]+
params->spin2[1]*LNhvec[1]+
params->spin2[2]*LNhvec[2])/ma2;
603 if ((ma1>0.)&&(ma2>0.)) {
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;
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");
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]);
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");
655 const UINT4 Npatch=40;
656 const UINT4 offsetAttch = 2;
689 Nrdwave = (
INT4) (10. / cimagf(modefreqs->
data[0]) /
dt);
695 if ( atpos < Npatch || atpos + Npatch >= sigl->
length )
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);
711 if ( !rdwave || !inspwave || !dinspwave || !matchinspwave )
724 for (
i=0;
i<2;
i++) {
726 for (j = 0; j < Npatch; j++) {
727 inspwave->
data[j] = sigl->
data[2*(atpos - Npatch + j)+
i];
730 for (k=0;k<2*nmodes;k++) {
731 matchinspwave->
data[k] = inspwave->
data[Npatch-1-offsetAttch];
732 if ((k+1)<2*nmodes) {
742 for (j=0; j<Npatch; j++) {
743 inspwave->
data[j]=dinspwave->
data[j];
760 for (j = 0; j < Nrdwave; j++) {
761 sigl->
data[2*j + 2*(atpos - 1 - offsetAttch) +
i ] = rdwave->
data[j];
INT4 XLALPSpinGenerateQNMFreq(COMPLEX8Vector *modefreqs, InspiralTemplate *params, UINT4 l, INT4 m, UINT4 nmodes, REAL8 finalMass, REAL8 finalSpin)
INT4 XLALPSpinFinalMassSpin(REAL8 *finalMass, REAL8 *finalSpin, InspiralTemplate *params, REAL8 energy, REAL8 *LNhvec)
INT4 XLALGenerateWaveDerivative(REAL8Vector *dwave, REAL8Vector *wave, REAL8 dt)
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 XLAL_CALLGSL(statement)
COMPLEX8Vector * XLALCreateCOMPLEX8Vector(UINT4 length)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
void XLALDestroyCOMPLEX8Vector(COMPLEX8Vector *vector)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
The inspiral waveform parameter structure containing information about the waveform to be generated.