Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-5e288d3
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimFindAttachTime.c
Go to the documentation of this file.
1#include <stdlib.h>
2#include <gsl/gsl_vector.h>
3#include <gsl/gsl_deriv.h>
4
5
6
7/*#include <lal/LALDatatypes.h>
8#include <lal/LALSimInspiral.h>
9#include <lal/TimeSeries.h>
10#include <lal/LALStdlib.h>
11#include <lal/LALConstants.h>
12#include <lal/TimeSeries.h>
13#include <lal/Units.h>
14
15#include <lal/LALSimIMR.h>
16#include <lal/Date.h>
17
18#include <lal/SeqFactories.h>*/
19
20#include "LALSimIMREOBNRv2.h"
21#include "LALSimIMRSpinEOB.h"
22#include "LALSimFindAttachTime.h"
23#include "LALSimIMRSpinPrecEOBHcapExactDerivative.c"
24#include "LALSimIMRSpinEOBHamiltonianPrec.c"
25#include "LALSimIMREOBHybridRingdownPrec.c"
26#include "LALSpinPrecHcapRvecDerivative_v3opt.c"
27
28
29double XLALSimLocateOmegaTime(
30 REAL8Array *dynamicsHi,
31 unsigned int numdynvars,
32 unsigned int retLenHi,
33 SpinEOBParams seobParams,
34 SpinEOBHCoeffs seobCoeffs,
35 REAL8 m1,
36 REAL8 m2,
37 REAL8Vector *radiusVec,
38 int *found,
39 REAL8* tMaxOmega,
40 INT4 use_optimized
41 )
42{
43 /* check that retLenHi is at least 2 */
44 XLAL_CHECK_ABORT( retLenHi > 2 && "retLenHi must be greater than 2" );
45
46 *tMaxOmega = 0; //Zach E: Fixes Heisenbug with ICC 16 and 17 compilers (5181); removing this line will result in segfaults with both compilers.
47 /*
48 * Locate merger point (max omega),
49 * WaveStep 1.1: locate merger point */
50 int debugPK = 0;
51 int debugRD = 0;
52 FILE *out = NULL;
53 gsl_spline *spline = NULL;
54 gsl_interp_accel *acc = NULL;
55
56 if (debugPK) {debugRD = 0;}
57
58 unsigned int peakIdx, i, j;
59 REAL8Vector *values = NULL;
60 REAL8Vector *dvalues = NULL;
61 REAL8Vector *omegaHi = NULL;
62
63
64 if ( !(values = XLALCreateREAL8Vector( numdynvars )) )
65 {
66 XLAL_ERROR( XLAL_ENOMEM );
67 }
68 if ( !(dvalues = XLALCreateREAL8Vector( numdynvars )) )
69 {
70 XLAL_ERROR( XLAL_ENOMEM );
71 }
72 if ( !(omegaHi = XLALCreateREAL8Vector( retLenHi )) )
73 {
74 XLAL_ERROR( XLAL_ENOMEM );
75 }
76 REAL8 rdotvec[3] = {0,0,0};
77 REAL8 rvec[3] = {0,0,0};
78 REAL8 rcrossrdot[3] = {0,0,0};
79 REAL8Vector timeHi;
80
81 timeHi.length = retLenHi;
82 timeHi.data = dynamicsHi->data;
83
84 double dt = timeHi.data[1] - timeHi.data[0];
85 double XLAL_INIT_DECL(ddradiusVec, [timeHi.length - 1]);
86 unsigned int k;
87 for (k = 1; k < timeHi.length-1; k++) {
88 ddradiusVec[k] = (radiusVec->data[k+1] - 2.*radiusVec->data[k] + radiusVec->data[k-1])/dt/dt;
89 // XLAL_PRINT_INFO("%3.10f %3.10f\n", timeHi->data[k], ddradiusVec[k]);
90 }
91 // for (k = timeHi->length-3; k>=1; k--) {
92 // XLAL_PRINT_INFO("%3.10f %3.10f\n", timeHi->data[k], ddradiusVec[k]);
93 // if (ddradiusVec[k] < 0) {
94 // break;
95 // }
96 // }
97 ddradiusVec[0] = ddradiusVec[1];
98 for (k = 0; k < timeHi.length-2; k++) {
99 // XLAL_PRINT_INFO("%3.10f %3.10f\n", timeHi->data[k], ddradiusVec[k]);
100 if (dt*k > dt*( timeHi.length-2)-20 && ddradiusVec[k] > 0) {
101 break;
102 }
103 }
104 double minoff = dt*( timeHi.length-2 - k) > 0.2 ? dt*( timeHi.length-2 - k) : 0.2;
105 if (debugPK) {
106 XLAL_PRINT_INFO("Change of sign in ddr %3.10f M before the end\n", minoff );
107 }
108 // First we search for the maximum (extremum) of amplitude
109 double maxoff = 20.0;
110 unsigned int Nps = timeHi.length;
111 // this definesthe search interval for maximum (we might use min0ff= 0.051 instead)
112 double tMin = timeHi.data[Nps-1] - maxoff;
113 // FIXME
114 // minoff = 0.0;
115 double tMax = timeHi.data[Nps-1] - minoff;
116 *tMaxOmega = tMax;
117 tMin = tMax - 20.;
118 if ( debugPK ) {
119 XLAL_PRINT_INFO("tMin, tMax = %3.10f %3.10f\n", tMin, tMax);
120 }
121
122 double omega = 0.0;
123
124 double magR;
125 double time1, time2, omegaDeriv, omegaDerivMid, tPeakOmega;
126
127 if(debugPK) {
128 out = fopen( "omegaHi.dat", "w" );
129 XLAL_PRINT_INFO("length of values = %d, retLenHi = %d\n", values->length, retLenHi);
130 fflush(NULL);
131 }
132 if(debugRD) {
133 out = fopen( "omegaHi.dat", "w" );
134 }
135
136 for ( i = 0; i < retLenHi; i++ )
137 {
138 for ( j = 0; j < values->length; j++ )
139 { values->data[j] = *(dynamicsHi->data+(j+1)*retLenHi+i); }
140
141 /* Calculate dr/dt */
142 memset( dvalues->data, 0, numdynvars*sizeof(dvalues->data[0]));
143 if(use_optimized){
144 if( XLALSpinPrecHcapRvecDerivative_exact( 0, values->data, dvalues->data,
145 &seobParams) != XLAL_SUCCESS )
146 {
147 XLAL_PRINT_INFO(
148 " Calculation of dr/dt failed while computing omegaHi time series\n");
149 XLAL_ERROR( XLAL_EFUNC );
150 }
151 } else {
152 if( XLALSpinPrecHcapRvecDerivative( 0, values->data, dvalues->data,
153 &seobParams) != XLAL_SUCCESS )
154 {
155 XLAL_PRINT_INFO(
156 " Calculation of dr/dt failed while computing omegaHi time series\n");
157 XLAL_ERROR( XLAL_EFUNC );
158 }
159 }
160
161 /* Calculare r x dr/dt */
162 for (j=0; j<3; j++){
163 rvec[j] = values->data[j];
164 rdotvec[j] = dvalues->data[j];
165 }
166
167 //memcpy(rdotvec, dvalues->data, 3*sizeof(REAL8));
168 //rvec[0] = posVecxHi.data[i]; rvec[1] = posVecyHi.data[i];
169 //rvec[2] = posVeczHi.data[i];
170 cross_product( rvec, rdotvec, rcrossrdot );
171
172 /* Calculate omega = |r x dr/dt| / r*r */
173 magR = sqrt(inner_product(rvec, rvec));
174 omegaHi->data[i] = sqrt(inner_product(rcrossrdot, rcrossrdot)) / (magR*magR);
175 if(debugPK || debugRD){
176 fprintf( out, "%.16e\t%.16e\n", timeHi.data[i], omegaHi->data[i]);
177 }
178 }
179
180
181 // Searching for crude omega_max (extremum)
182 peakIdx = 0;
183 *found = 0;
184 for ( i = 1, peakIdx = 0; i < retLenHi-1; i++ ){
185 omega = omegaHi->data[i];
186
187 if (omega >= omegaHi->data[i-1] && omega > omegaHi->data[i+1] && tMax>=timeHi.data[i] && timeHi.data[i]>=tMin){
188 peakIdx = i;
189 *found = 1;
190 if (debugPK){
191 XLAL_PRINT_INFO("PK: Crude peak of Omega is at idx = %d. t = %f, OmegaPeak = %.16e\n",
192 peakIdx, timeHi.data[peakIdx], omega);
193 fflush(NULL);
194
195 }
196 }
197 }
198
199 if(debugPK) {
200 fclose(out);
201 if (peakIdx ==0){
202 XLAL_PRINT_INFO("Stas: peak of orbital frequency was not found. peakIdx = %d, retLenHi = %d, i at exit = %d\n", peakIdx, retLenHi, i);
203 fflush(NULL);
204 }
205 }
206 if(debugRD) {
207 fclose(out);
208 }
209
210 // refining the omega_max search (if it is found)
211 tPeakOmega = 0.0;
212 if(peakIdx != 0){
213 spline = gsl_spline_alloc( gsl_interp_cspline, retLenHi );
214 acc = gsl_interp_accel_alloc();
215 time1 = timeHi.data[peakIdx-2];
216 gsl_spline_init( spline, timeHi.data, omegaHi->data, retLenHi );
217 omegaDeriv = gsl_spline_eval_deriv( spline, time1, acc );
218
219 if ( omegaDeriv > 0. ) { time2 = timeHi.data[peakIdx+2]; }
220 else{
221 time2 = time1;
222 peakIdx = peakIdx-2;
223 time1 = timeHi.data[peakIdx-2];
224 omegaDeriv = gsl_spline_eval_deriv( spline, time1, acc );
225 }
226
227 do
228 {
229 tPeakOmega = ( time1 + time2 ) / 2.;
230 omegaDerivMid = gsl_spline_eval_deriv( spline, tPeakOmega, acc );
231
232 if ( omegaDerivMid * omegaDeriv < 0.0 ) { time2 = tPeakOmega; }
233 else
234 {
235 omegaDeriv = omegaDerivMid;
236 time1 = tPeakOmega;
237 }
238 if (debugPK){
239 XLAL_PRINT_INFO("Stas: searching for orbital max: %f, %f, %f, %f \n", time1, time2, omegaDeriv, omegaDerivMid);
240 }
241 } while ( time2 - time1 > 1.0e-5 );
242 if(debugPK) {
243 XLAL_PRINT_INFO( "Estimation of the orbital peak is now at time %.16e \n", tPeakOmega);
244 fflush(NULL);
245 }
246 }
247
248 if(*found == 0 || debugRD || debugPK){
249 if(debugPK){
250 XLAL_PRINT_INFO("Stas: We couldn't find the maximum of orbital frequency, search for maximum of A(r)/r^2 \n");
251 }
252 REAL8 rad, rad2, m1PlusetaKK, bulk, logTerms, deltaU, u, u2, u3, u4, u5;
253 REAL8 listAOverr2[retLenHi];
254 REAL8 Aoverr2;
255 REAL8Vector *sigmaStar = NULL;
256 REAL8Vector *sigmaKerr = NULL;
257 if ( !(sigmaStar = XLALCreateREAL8Vector( 3 )) )
258 {
259 XLALDestroyREAL8Vector( sigmaStar );
260 XLAL_ERROR( XLAL_ENOMEM );
261 }
262 if ( !(sigmaKerr = XLALCreateREAL8Vector( 3 )) )
263 {
264 XLALDestroyREAL8Vector( sigmaStar );
265 XLAL_ERROR( XLAL_ENOMEM );
266 }
267 REAL8Vector s1Vec, s2Vec;
268 s1Vec.length = s2Vec.length = 3;
269 REAL8 s1Data[3], s2Data[3];
270 REAL8 mTotal = m1 + m2;
271 REAL8 a;
272 REAL8 eta = m1*m2/(mTotal*mTotal);
273
274 if(debugPK || debugRD){
275 out = fopen( "OutAofR.dat", "w" );
276 }
277 for ( i = 0; i < retLenHi; i++ )
278 {
279 for ( j = 0; j < values->length; j++ )
280 {
281 values->data[j] = *(dynamicsHi->data+(j+1)*retLenHi+i);
282 }
283 for( j = 0; j < 3; j++ )
284 {
285 //s1DataNorm[k] = values->data[k+6];
286 //s2DataNorm[k] = values->data[k+9];
287 s1Data[j] = values->data[j+6] * mTotal * mTotal;
288 s2Data[j] = values->data[j+9] * mTotal * mTotal;
289 }
290 s1Vec.data = s1Data;
291 s2Vec.data = s2Data;
292 XLALSimIMRSpinEOBCalculateSigmaStar( sigmaStar, m1, m2, &s1Vec, &s2Vec );
293 XLALSimIMRSpinEOBCalculateSigmaKerr( sigmaKerr, m1, m2, &s1Vec, &s2Vec );
294
295 seobParams.a = a = sqrt(inner_product(sigmaKerr->data, sigmaKerr->data));
296 m1PlusetaKK = -1. + eta * seobCoeffs.KK;
297 rad2 = values->data[0]*values->data[0] + values->data[1]*values->data[1] + values->data[2]*values->data[2];
298 rad = sqrt(rad2);
299 u = 1./rad;
300 u2 = u*u;
301 u3 = u2*u;
302 u4 = u2*u2;
303 u5 = u4*u;
304 bulk = 1./(m1PlusetaKK*m1PlusetaKK) + (2.*u)/m1PlusetaKK + a*a*u2;
305 logTerms = 1. + eta*seobCoeffs.k0 + eta*log(1. + seobCoeffs.k1*u + seobCoeffs.k2*u2 + seobCoeffs.k3*u3 + seobCoeffs.k4*u4 + seobCoeffs.k5*u5 + seobCoeffs.k5l*u5*log(u));
306 deltaU = bulk*logTerms;
307 listAOverr2[i] = deltaU / rad2;
308 if(debugPK || debugRD){
309 fprintf(out, "%3.10f %3.10f\n", timeHi.data[i], listAOverr2[i]);
310 }
311
312 }
313 if(debugPK || debugRD ) fclose(out);
314 if (*found == 0){
315 // searching formaximum of A(r)/r^2
316 peakIdx = 0;
317 //*found = 0;
318 for ( i = 1, peakIdx = 0; i < retLenHi-1; i++ ){
319 Aoverr2 = listAOverr2[i];
320 if (Aoverr2 >= listAOverr2[i-1] && Aoverr2 > listAOverr2[i+1]){
321 if (timeHi.data[i] > tMin){
322 peakIdx = i;
323 tPeakOmega = timeHi.data[i];
324 *found = 1;
325 if (debugPK){
326 XLAL_PRINT_INFO("PK: Peak of A(r)/r^2 is at idx = %d. t = %f, Peak ampl. = %.16e\n",
327 peakIdx, timeHi.data[peakIdx], Aoverr2);
328 fflush(NULL);
329 }
330 break;
331 }
332 }
333 }
334 }
335
336 if(debugPK) {
337 if (peakIdx ==0){
338 XLAL_PRINT_INFO("Stas: peak of A(r)/r^2 was not found. \
339 peakIdx = %d, retLenHi = %d, i at exit = %d\n", peakIdx, retLenHi, i);
340 fflush(NULL);
341 }
342 }
343 XLALDestroyREAL8Vector(sigmaStar);
344 XLALDestroyREAL8Vector(sigmaKerr);
345 }
346 if (spline != NULL)
347 gsl_spline_free(spline);
348 if (acc != NULL)
349 gsl_interp_accel_free(acc);
350 XLALDestroyREAL8Vector( values );
351 XLALDestroyREAL8Vector( dvalues );
352 XLALDestroyREAL8Vector( omegaHi );
353 if (*found == 0){
354 return(timeHi.data[retLenHi-1]);
355 }
356 else{
357 return(tPeakOmega);
358 }
359}
360
361double XLALSimLocateMaxAmplTime(
362 REAL8Vector *timeHi,
363 COMPLEX16Vector *hP22,
364 int *found)
365{
366 int debug = 0;
367 FILE *out = NULL;
368 unsigned int i, peakIdx;
369 unsigned int NpsSmall = 0;
370 double tAmpMax = 0.;
371 double AmpMax = 0.;
372 double Ampl = 0.;
373 double AmplN = 0.0;
374 double AmplO = 0.0;
375 *found = 0;
376
377 if(debug) {
378 out = fopen( "AmpPHi.dat", "w" );
379 }
380 int iMin = 0;
381 int iMax = timeHi->length-1;
382 NpsSmall = iMax - iMin + 1;
383 AmplO = sqrt(creal(hP22->data[iMin + 0])*creal(hP22->data[iMin + 0]) + cimag(hP22->data[iMin + 0])*cimag(hP22->data[iMin + 0]));
384 Ampl = AmplO;
385 tAmpMax = timeHi->data[iMin];
386
387 peakIdx = iMin;
388 for (i=0; i<NpsSmall-1; i++){
389 //tSeries[i] = timeHi->data[iMin + i];
390 AmplN = sqrt(creal(hP22->data[iMin + i+1])*creal(hP22->data[iMin + i+1]) + cimag(hP22->data[iMin + i+1])*cimag(hP22->data[iMin + i+1]));
391 //Ampl = sqrt(hreP22->data[i]*hreP22->data[i] + himP22->data[i]*himP22->data[i]);
392 if(debug){
393 fprintf(out, "%3.10f %3.10f\n", timeHi->data[iMin+i], Ampl);
394 }
395 if (Ampl >= AmplO && Ampl >AmplN){
396 if (*found !=1){
397 *found = 1;
398 tAmpMax = timeHi->data[iMin + i];
399 AmpMax = Ampl;
400 peakIdx = iMin + i;
401 }else{
402 if (debug){
403 XLAL_PRINT_INFO("Stas dismissing time %3.10f outside limits %3.10f, %3.10f \n",
404 timeHi->data[iMin + i], timeHi->data[iMin], timeHi->data[iMax]);
405 }
406 }
407 }
408 AmplO = Ampl;
409 Ampl = AmplN;
410 }
411 if (debug)
412 {
413 fclose(out);
414 if (*found ==0){
415 XLAL_PRINT_INFO("Peak of 2,2 mode in P-frame was not found. peakIdx = %d, retLenHi = %d, i at exit = %d\n", peakIdx, timeHi->length, i);
416 fflush(NULL);
417 }else{
418 XLAL_PRINT_INFO("We have found maximum of amplitude %3.10f at t = %3.10f \n", AmpMax, tAmpMax);
419 }
420 }
421
422 return (tAmpMax);
423
424
425
426}
427
428double XLALSimLocateAmplTime(
429 REAL8Vector *timeHi,
430 COMPLEX16Vector *hP22,
431 REAL8Vector *radiusVec,
432 int *found,
433 REAL8* tMaxAmp)
434{
435 int debugPK = 0;
436 int debugRD = 0;
437 FILE *out = NULL;
438 gsl_spline *spline = NULL;
439 gsl_interp_accel *acc = NULL;
440 if (debugPK) {debugRD = 0;}
441
442 double dt = timeHi->data[1] - timeHi->data[0];
443 double XLAL_INIT_DECL(ddradiusVec, [timeHi->length - 1]);
444 unsigned int k;
445
446 /* check that timeHi->length is at least 2 */
447 XLAL_CHECK_ABORT( timeHi->length > 2 && "timeHi->length must be greater than 2" );
448
449 for (k = 1; k < timeHi->length-1; k++) {
450 ddradiusVec[k] = (radiusVec->data[k+1] - 2.*radiusVec->data[k] + radiusVec->data[k-1])/dt/dt;
451// XLAL_PRINT_INFO("%3.10f %3.10f\n", timeHi->data[k], ddradiusVec[k]);
452 }
453// for (k = timeHi->length-3; k>=1; k--) {
454// XLAL_PRINT_INFO("%3.10f %3.10f\n", timeHi->data[k], ddradiusVec[k]);
455// if (ddradiusVec[k] < 0) {
456// break;
457// }
458// }
459 ddradiusVec[0]=ddradiusVec[1];
460 for (k = 0; k < timeHi->length-2; k++) {
461// XLAL_PRINT_INFO("%3.10f %3.10f\n", timeHi->data[k], ddradiusVec[k]);
462 if (dt*k > dt*( timeHi->length-2)-20 && ddradiusVec[k] > 0) {
463 break;
464 }
465 }
466 double minoff = dt*( timeHi->length-2 - k) > 0.2 ? dt*( timeHi->length-2 - k) : 0.2;
467 if (debugPK) {
468 XLAL_PRINT_INFO("Change of sign in ddr %3.10f M before the end\n", minoff );
469 }
470 // First we search for the maximum (extremum) of amplitude
471 unsigned int i, peakIdx;
472 double maxoff = 20.0;
473 unsigned int Nps = timeHi->length;
474 // this definesthe search interval for maximum (we might use min0ff= 0.051 instead)
475
476 //FIXME
477 //minoff = 0.0;
478 double tMin = timeHi->data[Nps-1] - maxoff;
479 double tMax = timeHi->data[Nps-1] - minoff;
480 *tMaxAmp = tMax;
481 tMin = tMax - maxoff;
482 if ( debugPK ) {
483 XLAL_PRINT_INFO("tMin, tMax = %3.10f %3.10f \n", tMin, tMax);
484 }
485 unsigned int iMin = ceil((tMin-timeHi->data[0])/dt);
486 unsigned int iMax = floor((tMax-timeHi->data[0])/dt);
487 unsigned int NpsSmall = iMax - iMin + 1;
488
489 double AmplN, AmplO;
490 double tAmpMax, tAmp;
491 double AmpMax = 0;
492 tAmpMax = 0.;
493 REAL8 tSeries[NpsSmall], Ampl[NpsSmall];
494
495 if(debugPK || debugRD) {
496 out = fopen( "AmpPHi.dat", "w" );
497 }
498 AmplO = sqrt(creal(hP22->data[iMin + 0])*creal(hP22->data[iMin + 0]) + cimag(hP22->data[iMin + 0])*cimag(hP22->data[iMin + 0]));
499 Ampl[0] = AmplO;
500 peakIdx = 0;
501 for (i=0; i<NpsSmall-1; i++){
502 tSeries[i] = timeHi->data[iMin + i];
503 AmplN = sqrt(creal(hP22->data[iMin + i+1])*creal(hP22->data[iMin + i+1]) + cimag(hP22->data[iMin + i+1])*cimag(hP22->data[iMin + i+1]));
504 //Ampl = sqrt(hreP22->data[i]*hreP22->data[i] + himP22->data[i]*himP22->data[i]);
505 if(debugPK || debugRD){
506 fprintf(out, "%3.10f %3.10f\n", tSeries[i], Ampl[i]);
507 }
508 if (Ampl[i] >= AmplO && Ampl[i] >AmplN){
509 if (*found !=1){
510 tAmp = timeHi->data[iMin + i];
511 if (tAmp >=tMin && tAmp <= tMax ){
512 *found = 1;
513 tAmpMax = tAmp;
514 AmpMax = Ampl[i];
515 peakIdx = iMin + i;
516 }else{
517 if (debugPK){
518 XLAL_PRINT_INFO("Stas dismissing time %3.10f outside limits %3.10f, %3.10f \n",
519 tAmp, tMin, tMax);
520 }
521 }
522 }
523 }
524 AmplO = Ampl[i];
525 Ampl[i+1] = AmplN;
526 }
527
528 if (debugPK)
529 {
530 fclose(out);
531 if (*found ==0){
532 XLAL_PRINT_INFO("Stas: peak of 2,2 mode in P-frame was not found. peakIdx = %d, retLenHi = %d, i at exit = %d\n", peakIdx, Nps, i);
533 fflush(NULL);
534 }else{
535 XLAL_PRINT_INFO("Stas: we have found maximum of amplitude %3.10f at t = %3.10f \n", AmpMax, tAmpMax);
536 }
537 }
538 if (debugRD)
539 {
540 fclose(out);
541 }
542
543 if (*found ==0 || debugRD || debugPK){
544 // we haven't found the maximum of amplitude -> search for minimum of derivative (extremum)
545// spline = gsl_spline_alloc( gsl_interp_cspline, NpsSmall );
546// acc = gsl_interp_accel_alloc();
547// gsl_spline_init( spline, tSeries, Ampl, NpsSmall );
548
549
550 REAL8 AmpDot[NpsSmall];
551 REAL8 AmpDDot[NpsSmall];
552
553 for (i=1; i<NpsSmall-2; i++){
554// AmpDot[i] = gsl_spline_eval_deriv(spline, tSeries[i] , acc);
555 AmpDot[i] = (Ampl[i+1] - Ampl[i-1])/2./dt;
556 AmpDDot[i] = (Ampl[i+1] - 2.0*Ampl[i] + Ampl[i-1])/(dt*dt);
557 }
558 AmpDot[0] = AmpDot[1];
559 AmpDot[NpsSmall -2] = AmpDot[NpsSmall-3];
560 AmpDot[NpsSmall -1] = AmpDot[NpsSmall-2];
561 AmpDDot[0] = AmpDDot[1];
562 AmpDDot[NpsSmall -2] = AmpDDot[NpsSmall-3];
563 AmpDDot[NpsSmall -1] = AmpDDot[NpsSmall-2];
564 //XLAL_PRINT_INFO("Stas, check AmDot %f, %f, %f \n", AmpDot[NpsSmall -3], AmpDot[NpsSmall -2], AmpDot[NpsSmall -1]);
565 //XLAL_PRINT_INFO("Stas, check AmDDot %f, %f, %f \n", AmpDDot[NpsSmall -3], AmpDDot[NpsSmall -2], AmpDDot[NpsSmall -1]);
566
567
568
569 REAL8 AmpDotSmooth[NpsSmall];
570 // Compute moving average over 7 points
571 unsigned int win = 3;
572 // unsigned int win = 5;
573 //int j;
574 double norm = 1.0/(2.0*win+1.0);
575
576 AmpDotSmooth[win] = 0;
577 for (i=0; i<(2*win +1); i++){
578 AmpDotSmooth[win] += AmpDot[i];
579 }
580 AmpDotSmooth[win] *= norm;
581 for (i=0; i<win; i++){
582 AmpDotSmooth[i] = AmpDotSmooth[win];
583 }
584
585 for (i=win+1; i<NpsSmall-1 -win; i++){
586 AmpDotSmooth[i] = AmpDotSmooth[i-1] + norm*(AmpDot[i+win] - AmpDot[i-win-1]);
587 }
588 for (i=0; i<win; i++){
589 AmpDotSmooth[NpsSmall-win-1+i] = AmpDotSmooth[NpsSmall-win-2];
590 }
591
592 // second deriv (in case)
593 REAL8 AmpDDotSmooth[NpsSmall];
594 unsigned int win2 = 100;
595 // unsigned int win = 5;
596 //int j;
597 norm = 1.0/(2.0*win2+1.0);
598
599 AmpDDotSmooth[win2] = 0;
600 for (i=0; i<(2*win2 +1); i++){
601 AmpDDotSmooth[win2] += AmpDDot[i];
602 }
603 AmpDDotSmooth[win2] *= norm;
604 for (i=0; i<win2; i++){
605 AmpDDotSmooth[i] = AmpDDotSmooth[win2];
606 }
607
608 for (i=win2+1; i<NpsSmall-1 -win2; i++){
609 AmpDDotSmooth[i] = AmpDDotSmooth[i-1] + norm*(AmpDDot[i+win2] - AmpDDot[i-win2-1]);
610 }
611 for (i=0; i<win2; i++){
612 AmpDDotSmooth[NpsSmall-win2-1+i] = AmpDDotSmooth[NpsSmall-win2-2];
613 }
614
615
616 if(debugPK || debugRD) {
617 out = fopen( "DotAmpPHi.dat", "w" );
618 for (i=0; i<NpsSmall - 1; i++){
619 fprintf(out, "%3.10f %3.10f %3.10f %3.10f %3.10f\n", tSeries[i], AmpDot[i], AmpDotSmooth[i], AmpDDot[i], AmpDDotSmooth[i]);
620 }
621 fclose(out);
622 }
623 if (*found ==0){
624 if (debugPK || debugRD){
625 XLAL_PRINT_INFO("Max of Amplitude is not found, looking for min of dot{Ampl} %d \n", iMin);
626 }
627 for (i=1; i<NpsSmall-1-win; i++){
628 if (AmpDotSmooth[i] < AmpDotSmooth[i-1] && AmpDotSmooth[i] < AmpDotSmooth[i+1]){
629 tAmp = tSeries[i];
630 //tAmp = tSeries[i-iMin];
631 //XLAL_PRINT_INFO("Stas check i = %d, tSeries = %3.10f, tAmp = %3.10f \n", i, tSeries[i], tAmp);
632 if (tAmp >=tMin && tAmp <= tMax && *found==0){
633 *found = 1;
634 tAmpMax = tAmp;
635 AmpMax = AmpDotSmooth[i];
636 //AmpMax = AmpDotSmooth[i-iMin];
637 peakIdx = i;
638 if (debugPK || debugRD){
639 XLAL_PRINT_INFO("we have found min of Adot at t= %f\n", tAmpMax);
640 }
641 //break;
642 }else{
643 if (debugPK){
644 XLAL_PRINT_INFO("Stas, AmplDot - dismissing time %3.10f outside limits %3.10f, %3.10f \n",
645 tAmp, tMin, tMax);
646 }
647 }
648
649 }
650 }
651 }
652 if (*found ==0){
653 if (debugPK || debugRD)
654 XLAL_PRINT_INFO("Min of Adot is not found, looking for min of ddot{Ampl} \n");
655 for (i=win2*2; i<NpsSmall-1-win2; i++){
656 if (AmpDDotSmooth[i] < AmpDDotSmooth[i-1] && AmpDDotSmooth[i] < AmpDDotSmooth[i+1]){
657 tAmp = tSeries[i];
658 //XLAL_PRINT_INFO("Stas check i = %d, tSeries = %3.10f, tAmp = %3.10f, %3.10f, %3.10f, %3.10f \n", i, tSeries[i], tAmp, AmpDDotSmooth[i-1], AmpDDotSmooth[i], AmpDDotSmooth[i+1]);
659 //if (tAmp >=tMin && tAmp <= tMax && *found==0)
660 if (tAmp >=tMin && tAmp <= tMax){
661 *found = 1;
662 tAmpMax = tAmp;
663 AmpMax = AmpDDotSmooth[i];
664 peakIdx = i;
665 //break;
666 }else{
667 if (debugPK){
668 XLAL_PRINT_INFO("Stas, AmplDDot - dismissing time %3.10f outside limits %3.10f, %3.10f \n",
669 tAmp, tMin, tMax);
670 }
671 }
672 }
673 }
674 }
675
676// if(*found==0){
677// REAL8 hRe[NpsSmall], hIm[NpsSmall];
678// REAL8 dhRe[NpsSmall], dhIm[NpsSmall];
679// for (i=0; i<NpsSmall - 1; i++) {
680// hRe[i] = creal(hP22->data[iMin + i]);
681// hIm[i] = cimag(hP22->data[iMin + i]);
682// }
683// for (i=1; i<NpsSmall - 2; i++) {
684// dhRe[i] = (hRe[i+1] - hRe[i-1])/2./dt;
685// dhIm[i] = (hIm[i+1] - hIm[i-1])/2./dt;
686// }
687// dhRe[0]=dhRe[1];
688// dhIm[0]=dhIm[1];
689// dhRe[NpsSmall-1]=dhRe[NpsSmall-2];
690// dhIm[NpsSmall-1]=dhIm[NpsSmall-2];
691//
692// REAL8 OmegaWave[NpsSmall], dOmegaWave[NpsSmall];
693// double hNorm2;
694// for (i=0; i<NpsSmall - 1; i++) {
695// hNorm2 = hRe[i]*hRe[i] + hIm[i]*hIm[i];
696// OmegaWave[i] = (-hRe[i]*dhIm[i] + hIm[i]*dhRe[i])/hNorm2;
697// }
698// for (i=1; i<NpsSmall - 2; i++) {
699// hNorm2 = hRe[i]*hRe[i] + hIm[i]*hIm[i];
700// dOmegaWave[i] = (OmegaWave[i+1] - OmegaWave[i-1])/2./dt;
701// }
702// dOmegaWave[0]=dOmegaWave[1];
703// dOmegaWave[NpsSmall-1]=dOmegaWave[NpsSmall-2];
704//
705// if(debugPK || debugRD) {
706// out = fopen( "OmegaW.dat", "w" );
707// for (i=0; i<NpsSmall - 1; i++){
708// fprintf(out, "%3.10f %3.10f %3.10f\n", tSeries[i], OmegaWave[i], dOmegaWave[i]);
709// }
710// }
711// fclose(out);
712//
713// }
714
715 /*
716 for (i=1; i<Nps-3; i++){
717 AmplDerivN2 = gsl_spline_eval_deriv(spline, timeHi->data[i+2], acc);
718 if(debugPK || debugRD){
719 fprintf(out, "%3.10f %3.10f\n", timeHi->data[i], AmplDeriv);
720 }
721
722 if (*found == 0){
723 if ((AmplDeriv <= AmplDerivO1 && AmplDeriv < AmplDerivN1) && (AmplDerivO2 > AmplDerivO1 && AmplDerivN2>=AmplDerivN1) ){
724 XLAL_PRINT_INFO("check %.16e, %.16e, %.16e, %.16e, %.16e \n", AmplDerivO2, AmplDerivO1, AmplDeriv, AmplDerivN1, AmplDerivN2);
725 tAmp = timeHi->data[i];
726 if (tAmp >=tMin && tAmp <= tMax && *found==0){
727 *found = 1;
728 tAmpMax = tAmp;
729 AmpMax = AmplDeriv;
730 peakIdx = i;
731 //break;
732 }else{
733 if (debugPK){
734 XLAL_PRINT_INFO("Stas dismissing time %3.10f outside limits %3.10f, %3.10f \n",
735 tAmp, tMin, tMax);
736 }
737 }
738 }
739 }
740 AmplDerivO2 = AmplDerivO1;
741 AmplDerivO1 = AmplDeriv;
742 AmplDeriv = AmplDerivN1;
743 AmplDerivN1 = AmplDerivN2;
744 }
745
746 if (debugPK)
747 {
748 fclose(out);
749 if (*found ==0){
750 XLAL_PRINT_INFO("Stas: peak of 2,2 mode in P-frame was not found. peakIdx = %d, retLenHi = %d, i at exit = %d\n", peakIdx, Nps, i);
751 fflush(NULL);
752 }else{
753 XLAL_PRINT_INFO("Stas: we have found maximum of amplitude %3.10f at t = %3.10f \n", AmpMax, tAmpMax);
754 }
755 }
756 if (debugRD)
757 {
758 fclose(out);
759 }
760 */
761
762 }
763
764
765
766 if (spline != NULL)
767 gsl_spline_free(spline);
768 if (acc != NULL)
769 gsl_interp_accel_free(acc);
770 if (*found == 0){
771 return(timeHi->data[Nps-1]);
772 }
773 else{
774 return(tAmpMax);
775 }
776
777}
778
779
780INT4 XLALSimCheckRDattachment(
781 REAL8Vector * signal1, /**<< Real of inspiral waveform to which we attach ringdown */
782 REAL8Vector * signal2, /**<< Imag of inspiral waveform to which we attach ringdown */
783 REAL8* ratio, /**<< output ratio */
784 const REAL8 tAtt, /**<< time of RD attachment */
785 const INT4 l, /**<< Current mode l */
786 const INT4 m, /**<< Current mode m */
787 const REAL8 dt, /**<< Sample time step (in seconds) */
788 const REAL8 mass1, /**<< First component mass (in Solar masses) */
789 const REAL8 mass2, /**<< Second component mass (in Solar masses) */
790 const REAL8 spin1x, /**<<The spin of the first object; only needed for spin waveforms */
791 const REAL8 spin1y, /**<<The spin of the first object; only needed for spin waveforms */
792 const REAL8 spin1z, /**<<The spin of the first object; only needed for spin waveforms */
793 const REAL8 spin2x, /**<<The spin of the second object; only needed for spin waveforms */
794 const REAL8 spin2y, /**<<The spin of the second object; only needed for spin waveforms */
795 const REAL8 spin2z, /**<<The spin of the second object; only needed for spin waveforms */
796 REAL8Vector * timeVec, /**<< Vector containing the time values */
797 REAL8Vector * matchrange, /**<< Time values chosen as points for performing comb matching */
798 Approximant approximant, /**<<The waveform approximant being used */
799 const REAL8 JLN, /**<< cosine of the angle between J and LN at the light ring */
800 REAL8 * timediff /**<< Time diff b/w peaks */
801 )
802{
803 int debugPK = 0;
804 unsigned int i;
805 unsigned int i_att = 0;
806 REAL8 Amp[signal1->length];
807 // sanity check
808 int ind_att = (int) matchrange->data[1]*(((mass1 + mass2) * LAL_MTSUN_SI / dt)) + 1;
809 if (debugPK){
810 XLAL_PRINT_INFO("attach_ind = %d, t =%f, %f \n", ind_att, matchrange->data[1], timeVec->data[ind_att]); fflush(NULL);
811 }
812 if (signal1->data[ind_att] == 0.0 && signal2->data[ind_att] == 0.0){
813 XLAL_PRINT_INFO("Opyat' signal = 0 \n");fflush(NULL);
814 //FILE *out1 = fopen( "Andrea1.dat","w");
815 //for (i = 0; i < timeVec->length; i++) {
816 // fprintf(out1, "%.16e %.16e %.16e\n", timeVec->data[i], signal1->data[i], signal2->data[i]);
817 //XLAL_PRINT_INFO("%.16e %.16e\n", timeVec->data[j], y[j]);
818 //}
819 //fclose(out1);
820 //exit(0);
821 XLAL_ERROR(XLAL_EFAILED);
822
823
824 }
825
826
827 if ( XLALSimIMREOBHybridAttachRingdownPrec( signal1, signal2, l, m,
828 dt, mass1, mass2, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
829 timeVec, matchrange, approximant, JLN )
830 == XLAL_FAILURE )
831 {
832 XLAL_ERROR( XLAL_EFUNC );
833 }
834 Amp[0] = sqrt(signal1->data[0]*signal1->data[0] + signal2->data[0]*signal2->data[0]);
835 for (i=1; i<signal1->length; i++){
836 Amp[i] = sqrt(signal1->data[i]*signal1->data[i] + signal2->data[i]*signal2->data[i]);
837 if ( i < timeVec->length){
838 if (timeVec->data[i-1] <= tAtt && timeVec->data[i] > tAtt){
839 i_att = i;
840 //if (debugPK){
841 // XLAL_PRINT_INFO(" the attachment index = %d, time = %f, %f \n", i_att, timeVec->data[i_att], tAtt);
842 //}
843 }
844 }
845 }
846 REAL8 timemaxR = timeVec->data[0], timemaxL = 0.;
847
848 REAL8 maxL = Amp[0];
849 UINT4 maxLidx = i_att - 1;
850 for (i=0; i<i_att; i++){
851 if (Amp[i] >= maxL){
852 maxL = Amp[i];
853 timemaxL = timeVec->data[i];
854 maxLidx = i;
855 if (debugPK){
856 XLAL_PRINT_INFO("i, timemaxL, maxL = %d %f %f\n",i,timemaxL,maxL);fflush(NULL);
857 }
858 }
859 }
860
861 timemaxR = timeVec->data[i_att];
862 REAL8 timemaxRtmp = timemaxR;
863 REAL8 maxR = Amp[i_att];
864 for (i=i_att; i<signal1->length; i++){
865 timemaxRtmp = timemaxRtmp + dt;
866 if (Amp[i] >= maxR){
867 maxR = Amp[i];
868 timemaxR = timemaxRtmp;
869 if (debugPK){
870 XLAL_PRINT_INFO("i, timemaxR, maxR = %d %f %f\n",i, timemaxR,maxR);fflush(NULL);
871 }
872 }
873 }
874 if ( maxLidx == i_att - 1 ) {
875 *timediff = timemaxR - timemaxL;
876 }
877 else {
878 *timediff = 0.;
879 }
880 if (debugPK){
881 XLAL_PRINT_INFO(" the ratio of amplitudes = %f , ampls = %f, %f \n", maxR/maxL, maxR, maxL);fflush(NULL);
882 XLAL_PRINT_INFO(" timemaxR, timemaxL = %f %f \n", timemaxR,timemaxL);fflush(NULL);
883 }
884
885 *ratio = maxR/maxL;
886 if (maxR/maxL != maxR/maxL){
887 //this is nan
888 *ratio = 1000.0;
889 }
890
891
892 return XLAL_SUCCESS;
893
894}
895
896
897int XLALSimAdjustRDattachmentTime(
898 REAL8Vector * signal1, /**<< Output Real of inspiral waveform to which we attach ringdown */
899 REAL8Vector * signal2, /**<< Output Imag of inspiral waveform to which we attach ringdown */
900 COMPLEX16TimeSeries* h22, /**<< input time series (inspiral) */
901 COMPLEX16TimeSeries* h2m2, /**<< input time series (inspiral) */
902 REAL8* ratio22, /**<< output ratio for 2,2 mode */
903 REAL8* ratio2m2, /**<< output ratio for 2,-2 mode*/
904 REAL8* tAttach, /**<< output/input time of RD attachment */
905 const REAL8 thr, /**<< threshold on the ratio */
906 const REAL8 dt, /**<< Sample time step (in seconds) */
907 const REAL8 m1, /**<< First component mass (in Solar masses) */
908 const REAL8 m2, /**<< Second component mass (in Solar masses) */
909 const REAL8 spin1x, /**<<The spin of the first object; only needed for spin waveforms */
910 const REAL8 spin1y, /**<<The spin of the first object; only needed for spin waveforms */
911 const REAL8 spin1z, /**<<The spin of the first object; only needed for spin waveforms */
912 const REAL8 spin2x, /**<<The spin of the second object; only needed for spin waveforms */
913 const REAL8 spin2y, /**<<The spin of the second object; only needed for spin waveforms */
914 const REAL8 spin2z, /**<<The spin of the second object; only needed for spin waveforms */
915 REAL8Vector * timeVec, /**<< Vector containing the time values */
916 REAL8Vector * matchrange, /**<< Time values chosen as points for performing comb matching */
917 Approximant approximant, /**<<The waveform approximant being used */
918 const REAL8 JLN, /**<< cosine of the angle between J and LN at the light ring */
919 const REAL8 combSize, /**<< combsize for RD attachment */
920 const REAL8 tMaxOmega, /**<< Time up to which we can trust dynamics */
921 const REAL8 tMaxAmp /**<< Time up to which we can trust waveform */
922 )
923{
924 int debugPK = 0;
925 REAL8 timediff = 0.;
926 REAL8 timediffStore = 0.;
927 unsigned int retLenHi = h22->data->length;
928 unsigned int i;
929 int pass = 0;
930 REAL8 tAtt;
931 tAtt = *tAttach; // replace with the loop
932 REAL8 maxDeltaT = 10.0;
933 REAL8 thrStore22L = 0., thrStore2m2L = 0., thrStore22R = 0., thrStore2m2R = 0., tLBest = *tAttach, tRBest = *tAttach;
934
935 REAL8 mTScaled = (retLenHi-1)*dt/matchrange->data[2];
936 REAL8 tMax = timeVec->data[retLenHi - 2] - 0.5 ;
937 double hNorm2, dsignal1, dsignal2; double omegaVec[retLenHi - 1];
938
939// XLAL_PRINT_INFO("tMaxOmega, tMaxAmp %f %f %f\n", tMaxOmega, tMaxAmp, tMax);
940
941 if ( tMaxAmp < tMax) {
942 tMax = tMaxAmp;
943 }
944 if ( tMaxOmega < tMax) {
945 tMax = tMaxOmega;
946 }
947 tMax = tMax - 3.;
948 if(tMax > tAtt + 5.0){
949 tMax = tAtt + 5.0;
950 }
951 REAL8 eta = m1*m2/(m1+m2)/(m1+m2);
952 REAL8 chiS = 0.5*(spin1z + spin2z);
953 REAL8 chiA = 0.5*(spin1z - spin2z);
954 REAL8 chi = chiS + chiA*sqrt(fabs(1. - 4.*eta))/(1. - 2.*eta);
955 if ( chi >= 0.8 ) {
956 if (tMaxAmp < tMaxOmega) {
957 tMax = tMaxAmp;
958 }
959 else {
960 tMax = tMaxOmega;
961 }
962 }
963 if (debugPK){
964 XLAL_PRINT_INFO("tmax = %f, tAtt = %f, tmaxAmp = %f, tmaxOm = %f\n", tMax, tAtt, tMaxAmp, tMaxOmega);
965 }
966// XLAL_PRINT_INFO("tAtt, tMax = %f %f\n", tAtt, tMax);
967 while(pass == 0 && (tAtt >= *tAttach - maxDeltaT)){
968 tAtt = tAtt - 0.5;
969 memset( signal1->data, 0, signal1->length * sizeof( signal1->data[0] ));
970 memset( signal2->data, 0, signal2->length * sizeof( signal2->data[0] ));
971 for ( i = 0; i < retLenHi; i++ )
972 {
973 signal1->data[i] = creal(h22->data->data[i]);
974 signal2->data[i] = cimag(h22->data->data[i]);
975 }
976
977 matchrange->data[0] = combSize < tAtt ? tAtt - combSize : 0;
978 matchrange->data[1] = tAtt;
979 matchrange->data[0] -= fmod( matchrange->data[0], dt/mTScaled );
980 matchrange->data[1] -= fmod( matchrange->data[1], dt/mTScaled );
981 if (debugPK) XLAL_PRINT_INFO("left 2,2 mode tAtt = %f ", tAtt);
982 if( XLALSimCheckRDattachment(signal1, signal2, ratio22, tAtt, 2, 2,
983 dt, m1, m2, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
984 timeVec, matchrange, approximant, JLN, &timediff ) == XLAL_FAILURE )
985 {
986 XLAL_ERROR( XLAL_EFUNC );
987 }
988
989 memset( signal1->data, 0, signal1->length * sizeof( signal1->data[0] ));
990 memset( signal2->data, 0, signal2->length * sizeof( signal2->data[0] ));
991 for ( i = 0; i < retLenHi; i++ )
992 {
993 signal1->data[i] = creal(h2m2->data->data[i]);
994 signal2->data[i] = cimag(h2m2->data->data[i]);
995 }
996
997 if (debugPK) XLAL_PRINT_INFO("left 2,-2 mode tAtt = %f ", tAtt);
998 if( XLALSimCheckRDattachment(signal1, signal2, ratio2m2, tAtt, 2, -2,
999 dt, m1, m2, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
1000 timeVec, matchrange, approximant, JLN, &timediff ) == XLAL_FAILURE )
1001 {
1002 XLAL_ERROR( XLAL_EFUNC );
1003 }
1004 if (debugPK) XLAL_PRINT_INFO("Quality %f\n", (*ratio22 - thr)*(*ratio22 - thr) + (*ratio2m2 - thr)*(*ratio2m2 - thr));
1005 if (debugPK) XLAL_PRINT_INFO("timediff %f\n", timediff);
1006
1007 if ( thrStore22L != 0. && thrStore2m2L != 0. ) {
1008 if ( (*ratio22 - thr)*(*ratio22 - thr) + (*ratio2m2 - thr)*(*ratio2m2 - thr) + timediff < (thrStore22L - thr)*(thrStore22L - thr) + (thrStore2m2L - thr)*(thrStore2m2L - thr) + timediffStore ) {
1009 thrStore22L = *ratio22;
1010 thrStore2m2L = *ratio2m2;
1011 timediffStore = timediff;
1012 tLBest = tAtt;
1013 if(debugPK)XLAL_PRINT_INFO("tLBest is now %f %f %f %f\n", tLBest, *ratio22 ,*ratio2m2, (*ratio22 - thr)*(*ratio22 - thr) + (*ratio2m2 - thr)*(*ratio2m2 - thr));
1014 }
1015 }
1016 else {
1017 thrStore22L = *ratio22;
1018 thrStore2m2L = *ratio2m2;
1019 timediffStore = timediff;
1020 tLBest = tAtt;
1021 }
1022
1023// if (*ratio22 <= thr && *ratio2m2 <= thr){
1024// pass = 1;
1025// }
1026
1027
1028 }
1029 memset( signal1->data, 0, signal1->length * sizeof( signal1->data[0] ));
1030 memset( signal2->data, 0, signal2->length * sizeof( signal2->data[0] ));
1031 if(debugPK){
1032 if (pass == 1){
1033 XLAL_PRINT_INFO("Going left, we have found better attachment point: new tAtt = %f, old = %f, ratios = %f, %f \n", tAtt, *tAttach, *ratio22, *ratio2m2);
1034 }else{
1035 XLAL_PRINT_INFO("Going left did nto find the best attachment point\n");
1036 }
1037
1038 }
1039 //REAL8 left_r22 = *ratio22;
1040 //REAL8 left_r2m2 = *ratio2m2;
1041 //REAL8 left_tAtt = tLBest;
1042 int pass_left = pass;
1043 int iBad = retLenHi - 1;
1044
1045 pass = 0;
1046 tAtt = *tAttach;
1047 timediffStore = 0.;
1048 while(pass == 0 && (tAtt < tMax-0.5)){
1049// XLAL_PRINT_INFO("tAtt tMax %f %f\n", tAtt, tMax);
1050 tAtt = tAtt + 0.5;
1051 memset( signal1->data, 0, signal1->length * sizeof( signal1->data[0] ));
1052 memset( signal2->data, 0, signal2->length * sizeof( signal2->data[0] ));
1053 matchrange->data[0] = combSize < tAtt ? tAtt - combSize : 0;
1054 matchrange->data[1] = tAtt;
1055 matchrange->data[0] -= fmod( matchrange->data[0], dt/mTScaled );
1056 matchrange->data[1] -= fmod( matchrange->data[1], dt/mTScaled );
1057 for ( i = 0; i < retLenHi; i++ )
1058 {
1059 signal1->data[i] = creal(h22->data->data[i]);
1060 signal2->data[i] = cimag(h22->data->data[i]);
1061 }
1062 for ( i = 0; i < retLenHi-1; i++ )
1063 {
1064 dsignal1 = (signal1->data[i+1] - signal1->data[i]);
1065 dsignal2 = (signal2->data[i+1] - signal2->data[i]);
1066 hNorm2 = signal1->data[i]*signal1->data[i] + signal2->data[i]*signal2->data[i];
1067 omegaVec[i] = (-dsignal1*signal2->data[i] + dsignal2*signal1->data[i])/hNorm2;
1068 }
1069 for ( i = 0; i < retLenHi-2; i++ )
1070 {
1071 if (omegaVec[i]*omegaVec[i+1] < 0) {
1072 iBad = i;
1073 break;
1074 }
1075 }
1076 //for ( i = iBad; i < retLenHi; i++ )
1077 //{
1078 // signal1->data[i] = 0.;
1079 // signal2->data[i] = 0.;
1080 //}
1081
1082
1083 if (debugPK) XLAL_PRINT_INFO("right 2,2 mode tAtt = %f ", tAtt);
1084 if (matchrange->data[1] < iBad){
1085 if( XLALSimCheckRDattachment(signal1, signal2, ratio22, tAtt, 2, 2,
1086 dt, m1, m2, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
1087 timeVec, matchrange, approximant, JLN, &timediff ) == XLAL_FAILURE )
1088 {
1089 XLAL_ERROR( XLAL_EFUNC );
1090 }
1091
1092
1093 memset( signal1->data, 0, signal1->length * sizeof( signal1->data[0] ));
1094 memset( signal2->data, 0, signal2->length * sizeof( signal2->data[0] ));
1095 for ( i = 0; i < retLenHi; i++ )
1096 {
1097 signal1->data[i] = creal(h2m2->data->data[i]);
1098 signal2->data[i] = cimag(h2m2->data->data[i]);
1099 }
1100
1101 if (debugPK) XLAL_PRINT_INFO("right 2,-2 mode tAtt = %f ", tAtt);
1102 if( XLALSimCheckRDattachment(signal1, signal2, ratio2m2, tAtt, 2, -2,
1103 dt, m1, m2, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
1104 timeVec, matchrange, approximant, JLN, &timediff ) == XLAL_FAILURE )
1105 {
1106 XLAL_ERROR( XLAL_EFUNC );
1107 }
1108 if (debugPK) XLAL_PRINT_INFO("Quality %f\n", (*ratio22 - thr)*(*ratio22 - thr) + (*ratio2m2 - thr)*(*ratio2m2 - thr));
1109 if (debugPK) XLAL_PRINT_INFO("timediff %f\n", timediff);
1110
1111
1112 if ( thrStore22R != 0. && thrStore2m2R != 0. ) {
1113 if ( tRBest < timeVec->data[iBad] && (*ratio22 - thr)*(*ratio22 - thr) + (*ratio2m2 - thr)*(*ratio2m2 - thr) + timediff < (thrStore22R - thr)*(thrStore22R - thr) + (thrStore2m2R - thr)*(thrStore2m2R - thr) + timediffStore ) {
1114 thrStore22R = *ratio22;
1115 thrStore2m2R = *ratio2m2;
1116 timediffStore = timediff;
1117 tRBest = tAtt;
1118 if(debugPK)XLAL_PRINT_INFO("tRBest is now %f %f %f %f\n", tRBest, *ratio22 , *ratio2m2, (*ratio22 - thr)*(*ratio22 - thr) + (*ratio2m2 - thr)*(*ratio2m2 - thr));
1119 }
1120 }
1121 else {
1122 thrStore22R = *ratio22;
1123 thrStore2m2R = *ratio2m2;
1124 timediffStore = timediff;
1125 tRBest = tAtt;
1126 }
1127 }
1128
1129
1130// if (*ratio22 <= thr && *ratio2m2 <= thr){
1131// pass = 1;
1132// }
1133
1134 }
1135 if(debugPK){
1136 if (pass == 1){
1137 XLAL_PRINT_INFO("Going right, we have found better attachment point: new tAtt = %f, old = %f, ratios = %f, %f \n", tAtt, *tAttach, *ratio22, *ratio2m2);
1138 }else{
1139 XLAL_PRINT_INFO("Going right did nto find the best attachment point\n");
1140 }
1141
1142 }
1143 int pass_right = pass;
1144
1145 if (1==1 || (pass_right == 0 && pass_left == 0) ) {
1146 if ( debugPK ) {
1147 XLAL_PRINT_INFO("Cannot go below required threshold on RD/insp amplitude\n");
1148 }
1149 if ( (thrStore22L - thr)*(thrStore22L - thr) + (thrStore2m2L - thr)*(thrStore2m2L - thr) < (thrStore22R - thr)*(thrStore22R - thr) + (thrStore2m2R - thr)*(thrStore2m2R - thr)) {
1150 *tAttach = tLBest;
1151 if ( debugPK ) {
1152 XLAL_PRINT_INFO("tLBest %f\n", tLBest);
1153 }
1154 }
1155 else {
1156 *tAttach = tRBest;
1157 if ( debugPK ) {
1158 XLAL_PRINT_INFO("tRBest %f\n", tRBest);
1159 }
1160
1161 }
1162 return(2);
1163 }
1164
1165 /*if( pass_right == 1 && pass_left == 0){
1166 *tAttach = tAtt;
1167 return(1);
1168 }
1169 if( pass_left == 1 && pass_right == 0){
1170 *tAttach =left_tAtt;
1171 *ratio22 = left_r22;
1172 *ratio2m2 = left_r2m2;
1173 return(1);
1174 }
1175 if (pass_left == 1 && pass_right == 1){
1176
1177 // hard choice
1178 if (left_r22 <= 1.0 || left_r2m2 <= 1.0){
1179 if (*ratio22 <= 1.0 || *ratio2m2 <= 1.0){
1180 *tAttach =left_tAtt;
1181 *ratio22 = left_r22;
1182 *ratio2m2 = left_r2m2;
1183 return(1); // choose left
1184 }else{
1185 *tAttach = tAtt;
1186 return(1); // choose right
1187 }
1188 }else{
1189 if (*ratio22 <= 1.0 || *ratio2m2 <= 1.0){
1190 *tAttach =left_tAtt;
1191 *ratio22 = left_r22;
1192 *ratio2m2 = left_r2m2;
1193 return(1); // choose left
1194 }else{ // booth looks ok so far
1195 if (left_r22*left_r2m2 < (*ratio22)*(*ratio2m2)){
1196 *tAttach =left_tAtt;
1197 *ratio22 = left_r22;
1198 *ratio2m2 = left_r2m2;
1199 return(1); // choose left
1200 }else{
1201 *tAttach = tAtt;
1202 return(1); // choose right
1203 }
1204
1205 }
1206 }
1207 }*/
1208 return(0);
1209
1210}
XLALSimLocateAmplTime
double XLALSimLocateAmplTime(REAL8Vector *timeHi, COMPLEX16Vector *hP22, REAL8Vector *radiusVec, int *found, REAL8 *tMaxAmp)
Definition: LALSimFindAttachTime.c:428
XLALSimAdjustRDattachmentTime
int XLALSimAdjustRDattachmentTime(REAL8Vector *signal1, REAL8Vector *signal2, COMPLEX16TimeSeries *h22, COMPLEX16TimeSeries *h2m2, REAL8 *ratio22, REAL8 *ratio2m2, REAL8 *tAttach, const REAL8 thr, const REAL8 dt, const REAL8 m1, const REAL8 m2, 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, const REAL8 combSize, const REAL8 tMaxOmega, const REAL8 tMaxAmp)
Definition: LALSimFindAttachTime.c:897
XLALSimLocateMaxAmplTime
double XLALSimLocateMaxAmplTime(REAL8Vector *timeHi, COMPLEX16Vector *hP22, int *found)
Definition: LALSimFindAttachTime.c:361
XLALSimLocateOmegaTime
double XLALSimLocateOmegaTime(REAL8Array *dynamicsHi, unsigned int numdynvars, unsigned int retLenHi, SpinEOBParams seobParams, SpinEOBHCoeffs seobCoeffs, REAL8 m1, REAL8 m2, REAL8Vector *radiusVec, int *found, REAL8 *tMaxOmega, INT4 use_optimized)
Definition: LALSimFindAttachTime.c:29
XLALSimCheckRDattachment
INT4 XLALSimCheckRDattachment(REAL8Vector *signal1, REAL8Vector *signal2, REAL8 *ratio, const REAL8 tAtt, 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, REAL8 *timediff)
Definition: LALSimFindAttachTime.c:780
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
XLALSimIMRSpinEOBCalculateSigmaKerr
static int XLALSimIMRSpinEOBCalculateSigmaKerr(REAL8Vector *sigmaKerr, REAL8 mass1, REAL8 mass2, REAL8Vector *s1, REAL8Vector *s2)
Function to calculate normalized spin of the deformed-Kerr background in SEOBNRv1.
Definition: LALSimIMRSpinEOBAuxFuncs.c:67
XLALSimIMRSpinEOBCalculateSigmaStar
static int XLALSimIMRSpinEOBCalculateSigmaStar(REAL8Vector *sigmaStar, REAL8 mass1, REAL8 mass2, REAL8Vector *s1, REAL8Vector *s2)
Function to calculate normalized spin of the test particle in SEOBNRv1.
Definition: LALSimIMRSpinEOBAuxFuncs.c:90
cross_product
static REAL8 * cross_product(const REAL8 values1[], const REAL8 values2[], REAL8 result[])
Definition: LALSimIMRSpinEOBHamiltonianPrec.c:694
XLALSpinPrecHcapRvecDerivative
static int XLALSpinPrecHcapRvecDerivative(double UNUSED t, const REAL8 values[], REAL8 dvalues[], void *funcParams)
Function to calculate numerical derivatives of the spin EOB Hamiltonian, to get , as decribed in Eqs.
Definition: LALSimIMRSpinEOBHamiltonianPrec.c:1026
inner_product
static REAL8 inner_product(const REAL8 values1[], const REAL8 values2[])
Functions to compute the inner product and cross products between vectors.
Definition: LALSimIMRSpinEOBHamiltonianPrec.c:684
LALSpinPrecHcapRvecDerivative_v3opt.c
XLALSpinPrecHcapRvecDerivative_exact
static int XLALSpinPrecHcapRvecDerivative_exact(double UNUSED t, const REAL8 values[], REAL8 dvalues[], void *funcParams)
Function to calculate exact derivatives of the spin EOB Hamiltonian, to get , as decribed in Eqs.
Definition: LALSpinPrecHcapRvecDerivative_v3opt.c:18
fprintf
#define fprintf
l
int l
Definition: bh_qnmode.c:135
dt
double dt
Definition: bh_ringdown.c:113
i
double i
Definition: bh_ringdown.c:118
deltaU
const double deltaU
Definition: exact_derivatives-HrealHeader.c:63
u
const double u
Definition: exact_derivatives-HrealHeader.c:17
u3
const double u3
Definition: exact_derivatives-HrealHeader.c:19
bulk
const double bulk
Definition: exact_derivatives-HrealHeader.c:57
u5
const double u5
Definition: exact_derivatives-HrealHeader.c:21
data
sigmaKerr data[0]
Definition: exact_derivatives-HrealHeader.c:4
m1PlusetaKK
const double m1PlusetaKK
Definition: exact_derivatives-HrealHeader.c:36
u2
const double u2
Definition: exact_derivatives-HrealHeader.c:18
u4
const double u4
Definition: exact_derivatives-HrealHeader.c:20
logTerms
const double logTerms
Definition: exact_derivatives-HrealHeader.c:62
LAL_MTSUN_SI
#define LAL_MTSUN_SI
REAL8
double REAL8
XLAL_INIT_DECL
#define XLAL_INIT_DECL(var,...)
UINT4
uint32_t UINT4
INT4
int32_t INT4
Approximant
Approximant
Enum that specifies the PN approximant to be used in computing the waveform.
Definition: LALSimInspiral.h:229
m
static const INT4 m
a
static const INT4 a
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLAL_PRINT_INFO
#define XLAL_PRINT_INFO(...)
XLAL_ERROR
#define XLAL_ERROR(...)
XLAL_CHECK_ABORT
#define XLAL_CHECK_ABORT(assertion)
XLAL_ENOMEM
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EFUNC
XLAL_EFAILED
XLAL_EFAILED
XLAL_FAILURE
XLAL_FAILURE
COMPLEX16TimeSeries::data
COMPLEX16Sequence * data
COMPLEX16Vector::data
COMPLEX16 * data
REAL8Array::data
REAL8 * data
REAL8Vector::data
REAL8 * data
SpinEOBHCoeffs
Parameters for the spinning EOB model, used in calculating the Hamiltonian.
Definition: LALSimIMRSpinEOB.h:57
SpinEOBParams
Parameters for the spinning EOB model.
Definition: LALSimIMRSpinEOB.h:116