Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimIMRPhenomXHM_internals.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2019 Marta Colleoni, Cecilio Garcia Quiros
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// LALSimIMRPhenomXHM_internals.c
21//
22// Created by Marta on 06/02/2019.
23//
24
25#include <gsl/gsl_linalg.h>
26
28
31
34
35
37
41
42
43/*Equations referenced in this file come from arXiv:2001.10914 [gr-qc], see also dcc:LIGO-P2000011 */
44
45/* Intialize Quasi Normal Modes' ringdown and damping frequencies */
47
48 // pointers to fring fits
53
54 //pointers to fdamp fits
59}
60
61/* Initialize Mixing Coefficients (complex). We assume mixing with just one mode. Only 32 mode has mixing. */
63{
68
69 // Adjust conventions so that they match the ones used for the hybrids
70 wf->mixingCoeffs[2]=-1.* wf->mixingCoeffs[2];
71 wf->mixingCoeffs[3]=-1.* wf->mixingCoeffs[3];
72
73 #if DEBUG == 1
74 printf("\nMixing coefficients:\n");
75 printf("re(a222)=%.16e,im(a222)=%.16e \n",evaluate_QNMfit_re_l2m2lp2(wf22->afinal),evaluate_QNMfit_im_l2m2lp2(wf22->afinal));
76 printf("re(a223)=%.16e,im(a223)=%.16e \n",evaluate_QNMfit_re_l2m2lp3(wf22->afinal),evaluate_QNMfit_im_l2m2lp3(wf22->afinal));
77 printf("re(a322)=%.16e,im(a322)=%.16e \n",evaluate_QNMfit_re_l3m2lp2(wf22->afinal),evaluate_QNMfit_im_l3m2lp2(wf22->afinal));
78 printf("re(a323)=%.16e,im(a323)=%.16e \n",evaluate_QNMfit_re_l3m2lp3(wf22->afinal),evaluate_QNMfit_im_l3m2lp3(wf22->afinal));
79 #endif
80}
81
82/* Store useful parameters specific of higher modes (thus not in IMRPhenomXWaveformStruct) in IMRPhenomXHMWaveformStruct */
83/* E.g.: MECO, ringdown and damping frequencies, the version of the fits, etc. */
85 int ell,
86 int emm,
89 QNMFits *qnms,
90 LALDict *LALParams
91)
92{
93
94 // read in which mode is being generated
95 wf->ell=ell;
96 wf->emm=emm;
97 wf->modeTag=ell*10+emm; //21, 33, 32, 44
98 wf->ampNorm = wf22->ampNorm;
99 wf->fMECOlm = wf22->fMECO*emm*0.5;
100 wf->Ampzero = 0; // Ampzero = 1 (true) for odd modes and equal black holes
101 wf->Amp0 = wf22->amp0;// * wf22->ampNorm;
102 wf->useFAmpPN = 0; // Only true for the 21, this mode has a different inspiral ansatz
103 wf->AmpEMR = 0; // Only one intermediate region
104 wf->InspiralAmpVeto=0;
106 wf->RingdownAmpVeto=0;
107
108 switch(wf->modeTag){
109 case 21:{
110 wf->modeInt=0;
111 wf->MixingOn=0;
112 if(wf22->q == 1. && wf22->chi1L == wf22->chi2L){ // Odd mode for equal mass, equal spin is zero
113 wf->Ampzero = 1;
114 }
115 break;
116 }
117 case 33:{
118 wf->modeInt=1;
119 wf->MixingOn=0;
120 if(wf22->q == 1. && wf22->chi1L == wf22->chi2L){ // Odd mode for equal mass, equal spin is zero
121 wf->Ampzero = 1;
122 }
123 break;
124 }
125 case 32:{ //Mode with Mixing
126 wf->modeInt=2;
127 wf->MixingOn=1;
128 break;
129 }
130 case 44:{
131 wf->modeInt=3;
132 wf->MixingOn=0;
133 break;
134 }
135 default:
136 {XLALPrintError("Error in IMRPhenomXHM_SetHMWaveformVariables: mode (%i,%i) selected is not currently available. Modes available are ((2,|2|),(2,|1|),(3,|2|),(3,|3|),(4,|4|)).\n", wf->ell, wf->emm);}
137 }
138
139 /* Here we select the version of the fits and of the reconstruction that will be used in the code.
140 Through XLALSimInspiralWaveformParamsLookupPhenomXHM(Inspiral/Intermediate/Ringdown)(Amp/Phase)Version we call the version of the fits used by the phase in each region.
141 Currently there is only one version available and is tagged by the release date in the format mmyyyy (122019)*/
143
147
151
155
156 /* Reconstruction version for the amplitude */
160
161 /* HM tuning is currently only setup for he l=m=3 moment. NOTE that default values for the parameters below are zero (IMRPhenomXPNRUseTunedCoprec) */
162 // NOTE that default values are all ZERO. We manually impose this here FOR ALL MODES
163 wf->MU1 = 0;
164 wf->MU2 = 0;
165 wf->MU3 = 0;
166 wf->MU4 = 0;
167 wf->NU0 = 0;
168 wf->NU4 = 0;
169 wf->NU5 = 0;
170 wf->NU6 = 0;
171 wf->ZETA1 = 0;
172 wf->ZETA2 = 0;
173 if(wf->modeTag==33){
174
175 /* Set parameters for coprecessing frame deviations. NOTE that we add a factor of delta so that all 33 deviations turn off at equal mass ratio -- i.e. there is NO calibration for equal mass ratio cases */
176 wf->PNR_DEV_PARAMETER = wf22->delta * (wf22->PNR_DEV_PARAMETER);
177
178 // IF PNR's coprecessing model is wanted, use its fits for the deviation parameters ELSE use the input values with defaults of zeros
180 {
181
182 /* ------------------------------------------------------ >>
183 Get them from the stored model fits that define PhenomXCP
184 within PhenomXPNR
185 << ------------------------------------------------------ */
186
187 /* MU1 modifies pAmp->lambda */
188 wf->MU1 = XLALSimIMRPhenomXCP_MU1_l3m3( wf22->theta_LS, wf22->eta, wf22->a1 );
189
190 // NOTE that the function for MU2 is not defined in the model
191 /* MU2 would modify pAmp->gamma2 */
192
193 /* MU2 */
194 wf->MU2 = XLALSimIMRPhenomXCP_MU2_l3m3( wf22->theta_LS, wf22->eta, wf22->a1 );
195
196 /* MU3 modifies pAmp->gamma3 */
197 wf->MU3 = XLALSimIMRPhenomXCP_MU3_l3m3( wf22->theta_LS, wf22->eta, wf22->a1 );
198
199 /* MU4 modifies V2 or V3 for the intermediate amplitude
200 for the DEFAULT value of IMRPhenomXIntermediateAmpVersion
201 use in IMRPhenomXPHM */
202 wf->MU4 = XLALSimIMRPhenomXCP_MU4_l3m3( wf22->theta_LS, wf22->eta, wf22->a1 );
203
204 // NOTE that we choose to disable time-shift tuning
205 // /* NU0 modifies pPhase->c0 */
206 // wf->NU0 = XLALSimIMRPhenomXCP_NU0_l3m3( wf22->theta_LS, wf22->eta, wf22->a1 );
207
208 /* NU4 modifies pPhase->cL */
209 wf->NU4 = XLALSimIMRPhenomXCP_NU4_l3m3( wf22->theta_LS, wf22->eta, wf22->a1 );
210
211 /* NU5 modifies wf->fRING [EXTRAP-PASS-TRUE] */
212 wf->NU5 = XLALSimIMRPhenomXCP_NU5_l3m3( wf22->theta_LS, wf22->eta, wf22->a1 );
213
214 /* NU6 modifies wf->fDAMP [EXTRAP-PASS-TRUE] */
215 wf->NU6 = XLALSimIMRPhenomXCP_NU6_l3m3( wf22->theta_LS, wf22->eta, wf22->a1 );
216
217 /* ZETA1 modifies pPhase->b4 */
218 wf->ZETA1 = XLALSimIMRPhenomXCP_ZETA1_l3m3( wf22->theta_LS, wf22->eta, wf22->a1 );
219
220 /* ZETA2 modifies pPhase->b1 */
221 wf->ZETA2 = XLALSimIMRPhenomXCP_ZETA2_l3m3( wf22->theta_LS, wf22->eta, wf22->a1 );
222
223 /* DEBUGGING: Turn off select deviations */
224 // wf->MU1 = 0;
225 // wf->MU2 = 0;
226 // wf->MU3 = 0;
227 // wf->MU4 = 0;
228
229 } else {
230
231 // NOTE that all default values are all ZERO
242
243 }
244
245 }
246
247
251
255
256 wf->nCollocPtsRDPhase = 0;
257
258 if (wf->IMRPhenomXHMReleaseVersion == 122019){
259 if(wf22->eta < 0.013886133703630232 && wf22->chi1L<=0.9){ // For q>70 and chi1<0.9 use two intermediate regions.
260 wf->AmpEMR = 1; // These cases have a more pronounced drop off in the amplitude at the end of the inspiral and need two intermediate regions to model it.
261 }
262 switch(wf->modeTag){
263 case 21:{
264 if(wf22->q < 8.){
265 wf->InspiralAmpVeto=1;
267 }
268 wf->RingdownAmpVeto=1;
269 if(wf22->eta >= 0.0237954){ //for EMR (q>40.)
270 wf->useFAmpPN = 1;
271 }
272 break;
273 }
274 case 33:{
275 if(wf22->q == 1. && wf22->chi1L == wf22->chi2L){ // Odd mode for equal mass, equal spin is zero
276 wf->Ampzero = 1;
277 }
278 break;
279 }
280 case 32:{ //Mode with Mixing
281 wf->RingdownAmpVeto=1;
282 wf->nCollocPtsRDPhase = 4;
283 break;
284 }
285 case 44:{
286 break;
287 }
288 default:
289 {XLALPrintError("Error in IMRPhenomXHM_SetHMWaveformVariables: mode (%i,%i) selected is not currently available. Modes available are ((2,|2|),(2,|1|),(3,|2|),(3,|3|),(4,|4|)).\n", wf->ell, wf->emm);}
290 }
294 }
295 else if (wf->IMRPhenomXHMReleaseVersion == 122022)
296 {
297 /* Here we change the versions for the collocation point frequencies, parameter space fits and ansatzaes for both amplitude and phase.
298 It is better to do it here instead of changing the default values in SimInspiralWaveformParams.c.
299 For a new release, copy this else if block and use it as a template. */
300
301 /* Amplitude */
302
303 // Versions for the cutting and collocation points frequencies
307 // Versions for the parameter space fits of coefficients or collocation points
311 // Versions for the ansatzaes
315
316 /* Notes for the new release and beyond:
317 - IMRPhenomXHMInspiralAmpVersion:
318 assuming that we have 3 collocation points, it indicates which collocation points are going to be used in the reconstruction.
319 e.g.: 123 means we use the three, 13 means we only use the 1st and the 3rd; 2 means we only use the 2nd, etc.
320 - IMRPhenomXHMIntermediateAmpVersion:
321 the number of digits indicate the number of collocation points in the intermediate part.
322 in this case we have 6 digits which relates to the 4 fitted collocation points + 2 boundaries
323 if the digit is 0 then we discard that collocation point
324 if the digit is 1 we use the value of the amplitude at that frequency
325 if the digit is 2 we use the value and the derivative of the derivative amplitude at that frequency
326 the digit 2 can only really go in the boundaries since the inner points are fixed parameter space fits of the amplitude value
327 at the boundaries we can use the inspiral or ringdown ansatz to compute the derivative
328 e.g.: 211112 means we will use the 4 parameter space fits + 2 value + 2 derivative at the boundaries = 8 degrees of freedom
329 e.g.: 110012 means we skip the left derivative and the two inner most collocation points = 5 degrees of freedom
330 - IMRPhenomXHMRingdownAmpVersion:
331 this is just an index for the different ansatz versions
332 currently it only takes 0 (old release) and 1 (122022 release)
333 */
334
335
336 // Mode specific tweaks for amplitude
337 if (wf->modeTag == 21){
338 /* The 21 behaves better if we skip the left the derivative.
339 The drawback is that for cases with drop the inspiral-intermediate transition will not be derivable.
340 We also drop two collocation points in the middle to avoid wavy behaviour in the extrapolation region */
342 }
343 else if (wf->modeTag == 32){
344 // Use fit of coefficients instead of collocation points
346 }
347
348 // As mentioned above, the number of digits in IMRPhenomXHMInsp(Intermediate)AmpVersion gives the number of collocation points for the inspiral(intermediate) part
349 wf->nCollocPtsInspAmp = snprintf(NULL, 0, "%i", wf->IMRPhenomXHMInspiralAmpVersion);
350 wf->nCollocPtsInterAmp = snprintf(NULL, 0, "%i", wf->IMRPhenomXHMIntermediateAmpVersion);
351
352
353 /* Phase */
354 /* HM Phase is unchanged respect to the old release */
355 if (wf->modeTag == 32){wf->nCollocPtsRDPhase = 4;}
356
357 /* Uncommenting the code below will employ a recalibration of the spheroidal phase for the 32.
358 This includes the collocation points for the derivative and the extra time and phase shift respect the 22 mode.
359 This can introduce introduce significant mismatch in the multimode waveform in the region without NR, so it is left for the next update. */
360 // if (wf->modeTag == 32){
361 // wf->IMRPhenomXHMIntermediatePhaseFreqsVersion = 122022;
362 // wf->IMRPhenomXHMRingdownPhaseFreqsVersion = 122022;
363 // wf->IMRPhenomXHMRingdownPhaseFitsVersion = 122022;
364 // wf->IMRPhenomXHMRingdownPhaseVersion = 122022;
365 // wf->nCollocPtsRDPhase = 5;
366 // }
367
368
369
370 }
371 else{
372 XLAL_ERROR_VOID(XLAL_EDOM, "Error in IMRPhenomXHM_SetHMWaveformVariables: IMRPhenomXHMReleaseVersion=%i is not valid.\n", wf->IMRPhenomXHMReleaseVersion);
373 }
374
375 /* Common for all releases */
376
377 if(wf->modeTag==32){
379 }
380 else{
382 }
383
384
385 /* Limit between comparable and extreme mass ratios for the phase */
386 wf->etaEMR=0.05;
387
388 /* Spin parameterisations: add here spin pars used for the higher modes that are not included in the 22 struct */
389 wf->chi_s = (wf22->chi1L + wf22->chi2L)*0.5;
390 wf->chi_a = (wf22->chi1L - wf22->chi2L)*0.5;
391
392 /* Ringdown and damping frequencies*/
393 wf->fRING = (qnms->fring_lm[wf->modeInt](wf22->afinal))/wf22->Mfinal;
394 wf->fDAMP = (qnms->fdamp_lm[wf->modeInt](wf22->afinal))/wf22->Mfinal;
395
396 /* Ringdown and damping frequencies*/
397 /* (IMRPhenomXPNRUseTunedCoprec) The EZH effective ringdown perscription uses the precessing final spin while the calibrated XPNR uses the non-precessing final spin to calculate its base waveform to which precession effect are added. */
398 #if DEBUG == 1
399 printf("\n>>> wf22->IMRPhenomXPNRUseTunedCoprec=%i ***\n",wf22->IMRPhenomXPNRUseTunedCoprec);
400 printf("*** PNR Co-precessing model in use for (l,m)=(%i,%i) ***\n\n",wf->ell,wf->emm);
401 printf("wf22->afinal : %e\n",wf22->afinal);
402 printf("wf22->afinal_prec : %e\n",wf22->afinal_prec);
403 #endif
404
405 /* We wish to use the EZH formula for the non-tuned HMs. This formula requires the precessing final spin,
406 * while tuning was performed relative to the non-precessing model and its non-precessing final spin
407 * (as would be conferred through wf22->afinal -- see IMRPhenomX_precession.c around about line 566).
408 *
409 * Just a note: when IMRPhenomXPNRUseTunedCoprec is false, afinal = afinal_prec and the code just above
410 * assigns the appropriate ringdown frequency and damping time. See LALSimIMRPhenomX_precession.c for
411 * the definition of afinal. */
413 {
414 wf->fRING = (qnms->fring_lm[wf->modeInt](wf22->afinal_prec))/wf22->Mfinal;
415 wf->fDAMP = (qnms->fdamp_lm[wf->modeInt](wf22->afinal_prec))/wf22->Mfinal;
416 #if DEBUG == 1
417 printf("\n** ell, emm ** : %i,%i\n",wf->ell,wf->emm);
418 printf("fring : %e\n",wf->fRING);
419 #endif
420 wf->fRING = wf->fRING - emm * wf22->fRINGEffShiftDividedByEmm;
421 #if DEBUG == 1
422 printf("fring shift : %e\n",- emm * wf22->fRINGEffShiftDividedByEmm);
423 printf("fring (coprec) : %e\n",wf->fRING);
424 printf(">> Note that fring has now been set to fringEff <<\n");
425 #endif
426 }
427
428 if ( wf22->IMRPhenomXPNRUseTunedCoprec33 ) {
429 if (wf->modeTag==33) {
430
431 #if DEBUG == 1
432 printf("fring : %e\n",wf->fRING);
433 #endif
434
435 // Apply PNR CoPrec deviations. NOTE that the are OFF when wf22->IMRPhenomXPNRUseTunedCoprec33 is false (see code above)
436 wf->fRING = wf->fRING - (wf->PNR_DEV_PARAMETER * wf->NU5);
437 wf->fDAMP = wf->fDAMP + (wf->PNR_DEV_PARAMETER * wf->NU6);
438
439 #if DEBUG == 1
440 printf("fring shift : %e\n",- (wf->PNR_DEV_PARAMETER * wf->NU5));
441 printf("fring (coprec) : %e\n",wf->fRING);
442 printf("fring prec : %e\n",(qnms->fring_lm[wf->modeInt](wf22->afinal_prec))/wf22->Mfinal);
443 printf("fring eff shift : %e\n",- emm * wf22->fRINGEffShiftDividedByEmm);
444 #endif
445
446 if( wf22->IMRPhenomXPNRUseTunedCoprec ){
447 // Note that we transition to the EZH effective ringdown frequency outside of the coprecessing calibration region, ie where PNR_DEV_PARAMETER and pnr_window become zero
448 wf->fRING = wf->fRING - (1.0-wf22->pnr_window) * emm * wf22->fRINGEffShiftDividedByEmm;
449 }
450
451 }
452 }
453
454 /* If (l,m)=(3,2), load the mixing coeffs to transform the spheroidal-harmonic ringdown ansatz back to spherical-harmonic */
455 if(wf->modeTag==32)
456 {
458 }
459
460 /* Linear part of the 22 phase. timeshift * ff + phaseshift. */
461 wf->timeshift =0;//XLALSimIMRPhenomXLinb(wf22->eta, wf22->chiPNHat, wf22->dchi, wf22->delta);
462 wf->phaseshift=0;//XLALSimIMRPhenomXLina(wf22->eta, wf22->chiPNHat, wf22->dchi, wf22->delta);
463 // current time-alignment of the hybrids
464 REAL8 psi4tostrain=XLALSimIMRPhenomXPsi4ToStrain(wf22->eta, wf22->STotR, wf22->dchi);
465 wf->DeltaT= -2.*LAL_PI*(500+psi4tostrain);
466
467 wf->fPhaseRDflat = 0.0;
468 wf->fAmpRDfalloff = 0.0;
469
470}
471
472/* Wrapper function to return ringdown frequency */
474 UINT4 ell,
475 UINT4 emm,
477 )
478{
479 /* emm is guaranteed to be positive */
480 UINT4 modeTag = ell * 10 + emm;
481
482 REAL8 fRING = 0.0;
483
484 /* if the tuned coprecessing tuning is activated, use the precessing final spin
485 * (the final spin is still precessing in the other case when called through XPHM,
486 * but the value is assigned to wf22->afinal instead) */
487 REAL8 afinal = (wf22->IMRPhenomXPNRUseTunedCoprec) ? wf22->afinal_prec : wf22->afinal;
488
489 switch(modeTag)
490 {
491 case 21:{
492 fRING = (evaluate_QNMfit_fring21(afinal))/wf22->Mfinal;
493 break;
494 }
495 case 22:{
496 fRING = wf22->fRING;
497 break;
498 }
499 case 33:{
500 fRING = (evaluate_QNMfit_fring33(afinal))/wf22->Mfinal;
501 break;
502 }
503 case 32:{
504 fRING = (evaluate_QNMfit_fring32(afinal))/wf22->Mfinal;
505 break;
506 }
507 case 44:{
508 fRING = (evaluate_QNMfit_fring44(afinal))/wf22->Mfinal;
509 break;
510 }
511 default:
512 {XLAL_ERROR_REAL8(XLAL_EDOM, "Error in IMRPhenomXHM_GenerateRingdownFrequency: mode (%i,%i) selected is not currently available. Modes available are ((2,|2|),(2,|1|),(3,|2|),(3,|3|),(4,|4|)).\n", ell, emm);}
513 }
514
515 /* if the coprecessing tuning is activated, return Effective RD frequency */
516 if( wf22->IMRPhenomXPNRUseTunedCoprec && ((ell!=2)||(emm!=2)) )
517 {
518 fRING -= emm * wf22->fRINGEffShiftDividedByEmm;
519 }
520
521 return fRING;
522}
523
524/* Store function names containing phase coefficient/collocation point fits in pPhase->[Inspiral|Intermediate|Ringdown]PhaseFits */
526
527 /* Here we fill some vectors of functions. Each function is the fit for one collocation point/coefficient.
528 There are three vectors, one for each region: inspiral, intermediate and ringdown.
529 The explicit fits can be found in the files of the corresponding region: IMRPhenomXHM_inspiral.c, _intermediate.c, _ringdown.c
530 */
531
532 //21
537
538 //21
545
546 //33
553
554 //32
561
562 //44
569
570 //32 Spheroidal
576
577}
578
579/* Store function names containing amplitude coeff./collocation point fits in pAmp->[Inspiral|Intermediate|Ringdown]AmpFits */
581
582 /* Here we fill some vectors of functions. Each function is the fit for one collocation point/coefficient.
583 There are three vectors, one for each region: inspiral, intermediate and ringdown.
584 The explicit fits can be found in the files of the corresponding region: IMRPhenomXHM_inspiral.c, _intermediate.c, _ringdown.c
585 */
586
587 //******Inspiral Fits for collocation points******/
588
589 //21 //Frequency of the collocation point
590 pAmp->InspiralAmpFits[0] = IMRPhenomXHM_Insp_Amp_21_iv1; //fcutInsp
591 pAmp->InspiralAmpFits[1] = IMRPhenomXHM_Insp_Amp_21_iv2; //fcutInsp*0.75
592 pAmp->InspiralAmpFits[2] = IMRPhenomXHM_Insp_Amp_21_iv3; //fcutInsp*0.5
593
594 //33
595 pAmp->InspiralAmpFits[3] = IMRPhenomXHM_Insp_Amp_33_iv1; //fcutInsp
596 pAmp->InspiralAmpFits[4] = IMRPhenomXHM_Insp_Amp_33_iv2; //fcutInsp*0.75
597 pAmp->InspiralAmpFits[5] = IMRPhenomXHM_Insp_Amp_33_iv3; //fcutInsp*0.5
598
599 //32
600 pAmp->InspiralAmpFits[6] = IMRPhenomXHM_Insp_Amp_32_iv1; //fcutInsp
601 pAmp->InspiralAmpFits[7] = IMRPhenomXHM_Insp_Amp_32_iv2; //fcutInsp*0.75
602 pAmp->InspiralAmpFits[8] = IMRPhenomXHM_Insp_Amp_32_iv3; //fcutInsp*0.5
603
604 //44
605 pAmp->InspiralAmpFits[9] = IMRPhenomXHM_Insp_Amp_44_iv1; //fcutInsp
606 pAmp->InspiralAmpFits[10] = IMRPhenomXHM_Insp_Amp_44_iv2; //fcutInsp*0.75
607 pAmp->InspiralAmpFits[11] = IMRPhenomXHM_Insp_Amp_44_iv3; //fcutInsp*0.5
608
609
610 /*****Intermediate Fits for EMR collocation points, 2 Intermediate regions*****/
611
612 //21 //Frequency of the collocation point
613 pAmp->IntermediateAmpFits[0] = IMRPhenomXHM_Inter_Amp_21_int1; //fcutInsp + (fcutRD-fcutInsp)/3
614 pAmp->IntermediateAmpFits[1] = IMRPhenomXHM_Inter_Amp_21_int2; //fcutInsp + 2(fcutRD-fcutInsp)/3
615
616 //33
617 pAmp->IntermediateAmpFits[2] = IMRPhenomXHM_Inter_Amp_33_int1; //fcutInsp + (fcutRD-fcutInsp)/3
618 pAmp->IntermediateAmpFits[3] = IMRPhenomXHM_Inter_Amp_33_int2; //fcutInsp + 2(fcutRD-fcutInsp)/3
619
620 //32
621 pAmp->IntermediateAmpFits[4] = IMRPhenomXHM_Inter_Amp_32_int1; //fcutInsp + (fcutRD-fcutInsp)/3
622 pAmp->IntermediateAmpFits[5] = IMRPhenomXHM_Inter_Amp_32_int2; //fcutInsp + 2(fcutRD-fcutInsp)/3
623
624 //44
625 pAmp->IntermediateAmpFits[6] = IMRPhenomXHM_Inter_Amp_44_int1; //fcutInsp + (fcutRD-fcutInsp)/3
626 pAmp->IntermediateAmpFits[7] = IMRPhenomXHM_Inter_Amp_44_int2; //fcutInsp + 2(fcutRD-fcutInsp)/3
627
628 //21 //fInt1 = fcutInsp + (fcutRD-fcutInsp)/3
629 pAmp->IntermediateAmpFits[8] = IMRPhenomXHM_Inter_Amp_21_int0; //fcutInsp + (fInt1 - fcutInsp)/3
630 pAmp->IntermediateAmpFits[9] = IMRPhenomXHM_Inter_Amp_21_dint0;//fcutInsp + (fInt1 - fcutInsp)/3
631
632 //33
633 pAmp->IntermediateAmpFits[10] = IMRPhenomXHM_Inter_Amp_33_int0; //fcutInsp + (fInt1 - fcutInsp)/3
634 pAmp->IntermediateAmpFits[11] = IMRPhenomXHM_Inter_Amp_33_dint0;//fcutInsp + (fInt1 - fcutInsp)/3
635
636 //32
637 pAmp->IntermediateAmpFits[12] = IMRPhenomXHM_Inter_Amp_32_int0; //fcutInsp + (fInt1 - fcutInsp)/3
638 pAmp->IntermediateAmpFits[13] = IMRPhenomXHM_Inter_Amp_32_dint0;//fcutInsp + (fInt1 - fcutInsp)/3
639
640 //44
641 pAmp->IntermediateAmpFits[14] = IMRPhenomXHM_Inter_Amp_44_int0; //fcutInsp + (fInt1 - fcutInsp)/3
642 pAmp->IntermediateAmpFits[15] = IMRPhenomXHM_Inter_Amp_44_dint0;//fcutInsp + (fInt1 - fcutInsp)/3
643
644 //21
647
648 //33
651
652 //32
655
656 //44
659
660 /****Ringdown Fits for coefficients*****/
661
662 //21
666
667 //33
670 pAmp->RingdownAmpFits[5] = IMRPhenomXHM_RD_Amp_33_sigma; //currently constant
671
672 //32
675 pAmp->RingdownAmpFits[8] = IMRPhenomXHM_RD_Amp_32_sigma; //currently constant
676
677 //44
680 pAmp->RingdownAmpFits[11] = IMRPhenomXHM_RD_Amp_44_sigma; //currently constant
681
682
683 /****Ringdown Fits for Collocation Points*****/
684
685 //21
689
690 //33
694
695 //32
699
700 //44
704
705 //32
708
709}
710
711
712/***************************************************/
713/* */
714/* Amplitude Cutting Frequencies */
715/* */
716/***************************************************/
717
718/* Inspiral cutting frequency for the Amplitude */
720
721 //Return the end frequency of the inspiral region and the beginning of the intermediate for the amplitude of one mode.
722
724 double fcut = 0.; //Cutting frequency for comparable mass ratios
725 double fMECO = pWFHM->fMECOlm;
726 double emm = 1.*(pWFHM->emm);
727 double eta = pWF22->eta;
728 double chi1 = pWF22->chi1L;
729 //fcutEMR is the cutting frequency for extreme mass ratios that is given by a fit to the frequncy of a particular geometrical structure of the amplitude
730 double fcutEMR = 1.25*emm*((0.011671068725758493 - 0.0000858396080377194*chi1 + 0.000316707064291237*pow(chi1,2))*(0.8447212540381764 + 6.2873167352395125*eta))/(1.2857082764038923 - 0.9977728883419751*chi1);
731
732 switch(version){
733 case 122018: // default version
734 {
735
736 double fring = pWFHM->fRING;
737 double chieff = pWF22->chiEff;
738 double fISCO = (pWF22->fISCO)*emm*0.5;
739
740 switch(pWFHM->modeTag)
741 {
742 case 21:{
743 if(eta < 0.023795359904818562){ //for EMR (q>40.)
744 fcut = fcutEMR;
745 }
746 else{ //for comparable q
747 fcut = fMECO + (0.75-0.235*chieff - 5./6.*chieff*chieff)*fabs(fISCO-fMECO);
748 }
749 break;
750 }
751 case 33:{
752 if(eta < 0.04535147392290249){ //for EMR (q>20.)
753 fcut = fcutEMR;
754 }
755 else{ //for comparable q
756 fcut = fMECO + (0.75-0.235*chieff-5./6.*chieff)*fabs(fISCO-fMECO);
757 }
758 break;
759 }
760 case 32:{
761 if(eta < 0.04535147392290249){ //for extreme mass ratios (q>20)
762 fcut = fcutEMR;
763 }
764 else{ //for comparable mass ratios
765 fcut = fMECO + (0.75-0.235*fabs(chieff))*fabs(fISCO-fMECO);
766 fcut = fcut*fring/pWF22->fRING;
767 }
768 break;
769 }
770 case 44:{
771 if(eta < 0.04535147392290249){ //for EMR (q>20)
772 fcut = fcutEMR;
773 }
774 else{ //for comparable q
775 fcut = fMECO + (0.75-0.235*chieff)*fabs(fISCO-fMECO);
776 }
777 break;
778 }
779 }
780 break;
781 }
782 case 122022:
783 {
784 if (pWF22->q < 20.){
785 fcut = fMECO;
786 }
787 else{
788 REAL8 transition_eta = 0.0192234; // q=50
789 REAL8 sharpness = 0.004;
790 REAL8 funcs = 0.5 + 0.5 * tanh((eta - transition_eta)/sharpness);
791 fcut = funcs * fMECO + (1 - funcs) * fcutEMR;
792 }
793 break;
794 }
795 default: {XLALPrintError("Error in IMRPhenomXHM_Amplitude_fcutInsp: version %i is not valid.", version);}
796 }
797 return fcut;
798}
799
800/* Ringdown cutting frequency for the amplitude */
802
803 //Returns the end of the intermediate region and the beginning of the ringdown for the amplitude of one mode
804
805 double fring = pWFHM->fRING, fdamp=pWFHM->fDAMP;
807 double eta = pWF22->eta;
808 double chi1 = pWF22->chi1L;
809 double fcut = 0.; //This is the cutting frequency
810
811 switch(version){
812 case 122018: // default version
813 {
814 switch(pWFHM->modeTag)
815 {
816 case 21:{
817 fcut = 0.75*fring;
818 break;
819 }
820 case 33:{
821 fcut = 0.95*fring;
822 break;
823 }
824 case 32:{
825 double fRD22 = pWF22->fRING;
826 double c = 0.5, r=5.;
827 if(eta < 0.0453515){
828 //for extreme mass ratios (q>20)
829 fcut = (fring*exp(c*r) + fRD22*exp(r*chi1))/(exp(c*r) + exp(r*chi1)) - fdamp; //This is a smooth step function between fRING (preferred by negative spins) and fRD22 (preferred by positive)
830 }else{
831 //for comparable mass ratios
832 fcut = fRD22;
833 }
834 if(0.02126654064272212<eta && eta<0.12244897959183673 && chi1>0.95) // for 6 < q < 45
835 {
836 fcut = fring - 2.*fdamp;
837 }
838 break;
839 }
840 case 44:{
841 fcut = 0.9*fring;
842 break;
843 }
844 }
845 break;
846 }
847 case 122022:
848 {
849 if(pWFHM->MixingOn == 1)
850 fcut = pWF22->fRING - 0.5*pWF22->fDAMP; //v8
851 else
852 fcut = fring - fdamp; //v2
853 break;
854 }
855 default: {XLALPrintError("Error in IMRPhenomXHM_Amplitude_fcutRD: version %i is not valid.", version);}
856 }
857 return fcut;
858}
859
860/***************************************************/
861/* */
862/* Phase Cutting & Collocation Points Frequencies */
863/* */
864/***************************************************/
865
866/* GetfcutInsp */
868 double MECOf=pWF22->fMECO;
869 double eta_factor=1.0+0.001*(0.25/pWF22->eta-1.);
870
871 return (eta_factor*pWFHM->emm*0.5*MECOf);
872}
873
874
875/********************** INTERMEDIATE PHASE COLLOCATION POINTS ******************/
877
878 //The frequencies are stored in the struct pPhase, which is initialized with new values as the code starts to work on each mode.
879
880 double fring = pWFHM->fRING, fdamp=pWFHM->fDAMP;
882
883 switch(version){
884 case 122019: // default version
885 case 122022:
886 {
887 double fcut= GetfcutInsp(pWF22,pWFHM);
888 pPhase->CollocationPointsFreqsPhaseInter[0]=fcut;
889
890 if(pWFHM->modeTag==32){
891
892
893 double fRD22 = pWF22->fRING, fdamp22 = pWF22->fDAMP;
894 double fEnd= fRD22- 0.5*fdamp22;
895 pPhase->CollocationPointsFreqsPhaseInter[1]=(sqrt(3)*(fcut - fEnd) + 2*(fcut + fEnd))/4.;
896 pPhase->CollocationPointsFreqsPhaseInter[2]=(3*fcut + fEnd)/4.;
897 pPhase->CollocationPointsFreqsPhaseInter[3]=(fcut + fEnd)/2.;
898 // we use first and second derivative at fEnd, so this frequency is duplicated here
899 pPhase->CollocationPointsFreqsPhaseInter[4]= fEnd;
900 pPhase->CollocationPointsFreqsPhaseInter[5]= fEnd;
901 pPhase->fPhaseMatchIM = fEnd;
902 // correct cutting frequency for EMR with negative spins
903 if(pWF22->eta<0.01&&pWF22->chi1L<0 && version==122019){
904 pPhase->fPhaseMatchIM=pPhase->fPhaseMatchIM*(1.2-0.25*pWF22->chi1L);
905 }
906
907
908 }
909
910 else{
911
912 pPhase->CollocationPointsFreqsPhaseInter[1]=(sqrt(3)*(fcut - fring) + 2*(fcut + fring))/4.;
913 pPhase->CollocationPointsFreqsPhaseInter[2]=(3*fcut + fring)/4.;
914 pPhase->CollocationPointsFreqsPhaseInter[3]=(fcut + fring)/2.;
915 pPhase->CollocationPointsFreqsPhaseInter[4]=(fcut + 3*fring)/4.;
916 pPhase->CollocationPointsFreqsPhaseInter[5]=(fcut + 7*fring)/8.;
917 pPhase->fPhaseMatchIM = fring-fdamp;
918
919 }
920
921 break;
922 }
923 default: {XLALPrintError("Error in IMRPhenomXHM_Intermediate_CollocPtsFreqs: version is not valid. Version recommended is 122019.");}
924 }
925
926 pPhase->fPhaseMatchIN = pWFHM->fMECOlm;
927}
928
929/********************** RINGDOWN PHASE COLLOCATION POINTS ******************/
930
931// this function initializes the frequencies of the collocation points for the spheroidal ringdown reconstruction
933
934 //The frequencies are stored in the struct pPhase, which is initialized with new values as the code starts to work on each mode.
935
936 double fringlm = pWFHM->fRING, fdamplm = pWFHM->fDAMP;
937 double fring22 = pWF22->fRING;
938
940 case 122019:{
941 pPhase->CollocationPointsFreqsPhaseRD[0] = fring22;
942 pPhase->CollocationPointsFreqsPhaseRD[2] = fringlm - 0.5*fdamplm;
943 pPhase->CollocationPointsFreqsPhaseRD[1] = fringlm - 1.5*fdamplm;
944 pPhase->CollocationPointsFreqsPhaseRD[3] = fringlm + 0.5*fdamplm;
945 break;
946 }
947 case 122022:{
948 double fdamp22 = pWF22->fDAMP;
949 pPhase->CollocationPointsFreqsPhaseRD[0] = fring22 - fdamp22;
950 pPhase->CollocationPointsFreqsPhaseRD[1] = fring22;
951 pPhase->CollocationPointsFreqsPhaseRD[2] = (fring22 + fringlm) * 0.5;
952 pPhase->CollocationPointsFreqsPhaseRD[3] = fringlm;
953 pPhase->CollocationPointsFreqsPhaseRD[4] = fringlm + fdamplm;
954 break;
955 }
956 default:{XLAL_ERROR_VOID(XLAL_EINVAL,"Error in IMRPhenomXHM_Ringdown_CollocPtsFreqs: version %i is not valid.", pWFHM->IMRPhenomXHMRingdownPhaseFreqsVersion);}
957 }
958
959}
960
961
962/*******************************************/
963/* */
964/* COMPUTE AMPLITUDE COEFFICIENTS */
965/* */
966/*******************************************/
967
968/*** Post-Newtonian Inspiral Ansatz Coefficients ***/
970
971 /* Fill pAmp with the coefficients of the power series in frequency of the Forier Domain Post-Newtonian Inspiral Ansatz.
972 The 21 mode by default does not used the power series because it breaks down before the end of the inspiral, but it corresponding power series is available here. */
973
974 /* The ansatz in Fourier Domain is built as follows: we multiply the Time-domain Post-Newtonian series up to 3PN by the phasing factor given by the Stationary-Phase-Approximation,
975 and the we reexpand in powers of f up to 3PN.
976 The only difference with the 21 is that this does not perform the last reexpansion, so it is not a real power series but it is a quantity better behaved. */
977
978 /* The coefficients below correspond to those in eqs E10-E14 in arXiv:2001.10914 */
979
980 //int inspversion = pWFHM->IMRPhenomXHMInspiralAmpFitsVersion;
981 double chiA = pWFHM->chi_a;
982 double chiS = pWFHM->chi_s;
983 double eta = pWF22->eta, delta = pWF22->delta;
984 double PI = powers_of_lalpiHM.itself;
985
986 const double prefactors[] = {sqrt(2)/3., 0.75*sqrt(5/7.), sqrt(5/7.)/3., 4*sqrt(2)/9*sqrt(5/7.)}; //Global factors of each PN hlm
987 pAmp->PNglobalfactor = pow(2./(pWFHM->emm),-7/6.)*prefactors[pWFHM->modeInt]; //This is to compensate that we rescale data with the leading order of the 22
988
989 /*switch(inspversion){ FIXMEE
990 case 122018: // default version
991 {*/
992 switch(pWFHM->modeTag)
993 {
994 case 21:{
995 if(pWFHM->useFAmpPN == 1){
996 Get21PNAmplitudeCoefficients(pAmp, pWF22);
997 pAmp->pnInitial = 0.;
998 pAmp->pnOneThird = 0.;
999 pAmp->pnTwoThirds = 0.;
1000 pAmp->pnThreeThirds = 0.;
1001 pAmp->pnFourThirds = 0.;
1002 pAmp->pnFiveThirds = 0.;
1003 pAmp->pnSixThirds = 0.;
1004 }
1005 else{
1006 IMRPhenomX_UsefulPowers powers_of_2d1;
1007 IMRPhenomX_Initialize_Powers(&powers_of_2d1, 2.);
1008 pAmp->pnInitial = 0.;
1010 pAmp->pnTwoThirds = (-3*(chiA + chiS*delta))/2.*powers_of_lalpiHM.two_thirds*powers_of_2d1.two_thirds;
1011 pAmp->pnThreeThirds = (335*delta + 1404*delta*eta)/672.*powers_of_lalpiHM.itself*powers_of_2d1.itself;
1012 pAmp->pnFourThirds = (3427*chiA - (I*672)*delta + 3427*chiS*delta - 8404*chiA*eta - 3860*chiS*delta*eta - 1344*delta*PI - (I*672)*delta*log(16))/1344.*powers_of_lalpiHM.four_thirds*powers_of_2d1.four_thirds;
1013 pAmp->pnFiveThirds = (-155965824*chiA*chiS - 964357*delta + 432843264*chiA*chiS*eta - 23670792*delta*eta + 24385536*chiA*PI + 24385536*chiS*delta*PI - 77982912*delta*chiA*chiA + 81285120*delta*eta*chiA*chiA - 77982912*delta*chiS*chiS + 39626496*delta*eta*chiS*chiS + 21535920*delta*eta*eta)/8.128512e6*powers_of_lalpiHM.five_thirds*powers_of_2d1.five_thirds;
1014 pAmp->pnSixThirds = (143063173*chiA - (I*1350720)*delta + 143063173*chiS*delta - 546199608*chiA*eta - (I*72043776)*delta*eta - 169191096*chiS*delta*eta - 9898560*delta*PI + 20176128*delta*eta*PI - (I*5402880)*delta*log(2) - (I*17224704)*delta*eta*log(2) + 61725888*chiS*delta*chiA*chiA - 81285120*chiS*delta*eta*chiA*chiA + 20575296*pow(chiA,3) - 81285120*eta*pow(chiA,3) + 61725888*chiA*chiS*chiS - 165618432*chiA*eta*chiS*chiS + 20575296*delta*pow(chiS,3) - 1016064*delta*eta*chiS*chiS*chiS + 128873808*chiA*eta*eta - 3859632*chiS*delta*eta*eta)/5.419008e6*powers_of_lalpiHM.two*powers_of_2d1.two;
1015 }
1016 break;
1017 }
1018 case 33:{
1019 IMRPhenomX_UsefulPowers powers_of_2d3;
1020 IMRPhenomX_Initialize_Powers(&powers_of_2d3, 2./3.);
1021 pAmp->pnInitial = 0.;
1023 pAmp->pnTwoThirds = 0.;
1024 pAmp->pnThreeThirds = (-1945*delta + 2268*delta*eta)/672.*powers_of_lalpiHM.itself*powers_of_2d3.itself;
1025 pAmp->pnFourThirds = (325*chiA - (I*504)*delta + 325*chiS*delta - 1120*chiA*eta - 80*chiS*delta*eta + 120*delta*PI + (I*720)*delta*log(1.5))/120.*powers_of_lalpiHM.four_thirds*powers_of_2d3.four_thirds;
1026 pAmp->pnFiveThirds = (-2263282560*chiA*chiS - 1077664867*delta + 9053130240*chiA*chiS*eta - 5926068792*delta*eta - 1131641280*delta*chiA*chiA + 4470681600*delta*eta*chiA*chiA - 1131641280*delta*chiS*chiS + 55883520*delta*eta*chiS*chiS + 2966264784*delta*eta*eta)/4.4706816e8*powers_of_lalpiHM.five_thirds*powers_of_2d3.five_thirds;
1027 pAmp->pnSixThirds = (22007835*chiA + (I*26467560)*delta + 22007835*chiS*delta - 80190540*chiA*eta - (I*98774368)*delta*eta - 31722300*chiS*delta*eta - 9193500*delta*PI + 17826480*delta*eta*PI - (I*37810800)*delta*log(1.5) + (I*37558080)*delta*eta*log(1.5) - 12428640*chiA*eta*eta - 6078240*chiS*delta*eta*eta)/2.17728e6*powers_of_lalpiHM.two*powers_of_2d3.two;
1028 break;
1029 }
1030 case 32:{
1031 pAmp->pnInitial = 0.;
1032 pAmp->pnOneThird = 0.;
1033 pAmp->pnTwoThirds = (-1 + 3*eta)*powers_of_lalpiHM.two_thirds;
1034 pAmp->pnThreeThirds = -4*chiS*eta*powers_of_lalpiHM.itself;
1035 pAmp->pnFourThirds = (10471 - 61625*eta + 82460*eta*eta)/10080.*powers_of_lalpiHM.four_thirds;
1036 pAmp->pnFiveThirds = ((I*2520) - 3955*chiS - 3955*chiA*delta - (I*11088)*eta + 10810*chiS*eta + 11865*chiA*delta*eta - 12600*chiS*eta*eta)/840.*powers_of_lalpiHM.five_thirds;
1037 pAmp->pnSixThirds = (824173699 + 2263282560*chiA*chiS*delta - 26069649*eta - 15209631360*chiA*chiS*delta*eta + 3576545280*chiS*eta*PI + 1131641280*chiA*chiA - 7865605440*eta*chiA*chiA + 1131641280*chiS*chiS - 11870591040*eta*chiS*chiS - 13202119896*eta*eta + 13412044800*chiA*chiA*eta*eta + 5830513920*chiS*chiS*eta*eta + 5907445488*pow(eta,3))/4.4706816e8*powers_of_lalpiHM.two;
1038 break;
1039 }
1040 case 44:{
1041 IMRPhenomX_UsefulPowers powers_of_2d4;
1042 IMRPhenomX_Initialize_Powers(&powers_of_2d4, 0.5);
1043 pAmp->pnInitial = 0.;
1044 pAmp->pnOneThird = 0.;
1045 pAmp->pnTwoThirds = (1 - 3*eta)*powers_of_lalpiHM.two_thirds*powers_of_2d4.two_thirds;
1046 pAmp->pnThreeThirds = 0.;
1047 pAmp->pnFourThirds = (-158383 + 641105*eta - 446460*eta*eta)/36960.*powers_of_lalpiHM.four_thirds*powers_of_2d4.four_thirds;
1048 pAmp->pnFiveThirds = ((I*-1008) + 565*chiS + 565*chiA*delta + (I*3579)*eta - 2075*chiS*eta - 1695*chiA*delta*eta + 240*PI - 720*eta*PI + (I*960)*log(2) - (I*2880)*eta*log(2) + 1140*chiS*eta*eta)/120.*powers_of_lalpiHM.five_thirds*powers_of_2d4.five_thirds;
1049 pAmp->pnSixThirds = (7888301437 - 147113366400*chiA*chiS*delta - 745140957231*eta + 441340099200*chiA*chiS*delta*eta - 73556683200*chiA*chiA + 511264353600*eta*chiA*chiA - 73556683200*chiS*chiS + 224302478400*eta*chiS*chiS + 2271682065240*eta*eta - 871782912000*chiA*chiA*eta*eta - 10897286400*chiS*chiS*eta*eta - 805075876080*pow(eta,3))/2.90594304e10*powers_of_lalpiHM.two*powers_of_2d4.two;
1050 break;
1051 }
1052 }/*
1053 break;
1054 }
1055 default: {XLALPrintError("Error in IMRPhenomXHM_GetPNAmplitudeCoefficients: version is not valid. Version recommended is 122018. ");}
1056}*/
1057}
1058
1059/*** Post-Newtonian Inspiral Ansatz Coefficients for the 21 mode ***/
1061
1062 /* The 21 ansatz in Fourier Domain is built multiplying the Time-domain Post-Newtonian series up to 3PN by the phasing factor given by the Stationary-Phase-Approximatio */
1063
1064 double m1 = pWF22->m1;
1065 double m2 = pWF22->m2;
1066 double m12 = m1*m1;
1067 double m22 = m2*m2;
1068 double m13 = m12*m1;
1069 double m23 = m22*m2;
1070 double m14 = m13*m1;
1071 double m24 = m23*m2;
1072 double m15 = m14*m1;
1073 double m25 = m24*m2;
1074 double m16 = m15*m1;
1075
1076 double chi1 = pWF22->chi1L;
1077 double chi2 = pWF22->chi2L;
1078 double chiS = (chi1 + chi2)*0.5;
1079 double chiA = (chi1 - chi2)*0.5;
1080 double delta = pWF22->delta;
1081 double eta = pWF22->eta;
1082 double chi12 = chi1*chi1;
1083 double chi22 = chi2*chi2;
1084 double Sc = m12*chi1 + m22*chi2;
1085 double Sigmac = m2*chi2 - m1*chi1;
1086
1087 IMRPhenomX_UsefulPowers powers_of_2d1;
1088 IMRPhenomX_Initialize_Powers(&powers_of_2d1, 2.);
1089 double logof2 = powers_of_2d1.log;
1090 double log4 = 1.3862943611198906;
1091 double EulerGamma = 0.5772156649015329;
1092
1093 /* Complex coefficients of the Time-Domain Post-Newtonian expansion. */
1094
1095 double factor = 8*pWF22->eta*powers_of_lalpiHM.two_thirds*powers_of_2d1.two_thirds*powers_of_lalpiHM.sqrt/sqrt(5.);
1096 pAmp->PNTDfactor = factor;
1097 pAmp->x05 = I*delta/3*powers_of_2d1.one_third*powers_of_lalpiHM.one_third;
1098 pAmp->x1 = -I*0.5*(chiA + chiS*delta)*powers_of_2d1.two_thirds*powers_of_lalpiHM.two_thirds;
1099 pAmp->x15 = I*delta*(-17./28 + 5*eta/7.)/3.*powers_of_2d1.itself*powers_of_lalpiHM.itself;
1100 pAmp->x2 = ( I*(-43./21*delta*Sc + (-79.+139*eta)/42.*Sigmac) + I/3*delta*(powers_of_lalpiHM.itself + I*(-0.5-2*logof2)))*powers_of_2d1.four_thirds*powers_of_lalpiHM.four_thirds;
1101 pAmp->x25 = I*delta*((-43-509*eta)/126. + 79*eta*eta/168.)/3.*powers_of_2d1.five_thirds*powers_of_lalpiHM.five_thirds;
1102 pAmp->x3 = I*delta*( (-17. + 6.*eta)/28.*powers_of_lalpiHM.itself + I*(17/56. + eta*(-353/28. - 3.*logof2/7.) + 17.*logof2/14.))/3.*powers_of_2d1.two*powers_of_lalpiHM.two;
1103
1104
1105 /* Coefficients of the phasing factor expansion */
1106 /* These coefficients correspond to those in equation E4 of arXiv:2001.10914. */
1107 pAmp->xdot5 = -(m1*m2*(-838252800*m1*m2 - 419126400*m12 - 419126400*m22)/3.274425e7)*powers_of_2d1.five_thirds*powers_of_2d1.five_thirds*powers_of_lalpiHM.five_thirds*powers_of_lalpiHM.five_thirds;
1108 pAmp->xdot6 = -(m1*m2*(1152597600*m2*m13 + 926818200*m22 + 2494800*m1*m2*(743 + 462*m22) + 1247400*m12*(743 + 1848*m22))/3.274425e7)*powers_of_2d1.four*powers_of_lalpiHM.four;
1109 pAmp->xdot65 = -(m1*m2*(-34927200*m1*m2*(-(m2*(75*chi1 + 376*chi2*m2)) + 96*powers_of_lalpiHM.itself) - 34927200*(-(m2*(75*chi2 + 188*(chi1 + chi2)*m2)) + 48*powers_of_lalpiHM.itself)*m12 - 2619540000*chi1*m13 + 13132627200*chi1*m2*m13 +
1110 6566313600*chi1*m14 - 34927200*(chi2*(75 - 188*m2)*m2 + 48*powers_of_lalpiHM.itself)*m22)/3.274425e7)*powers_of_2d1.eight_thirds*powers_of_2d1.five_thirds*powers_of_lalpiHM.eight_thirds*powers_of_lalpiHM.five_thirds;
1111 pAmp->xdot7 = -(m1*m2*(207900*m2*(-13661 - 19908*chi1*chi2 + 10206*chi12 + 10206*chi22)*m13 - 23100*(34103 + 91854*chi22)*m22 - 1373803200*m14*m22 +
1112 23100*m1*m2*(-2*(34103 + 45927*chi12 + 45927*chi22) + 9*(-13661 - 19908*chi1*chi2 + 10206*chi12 + 10206*chi22)*m22) - 2747606400*m13*m23 -
1113 23100*m12*(34103 + 91854*chi12 - 18*(-13661 - 19908*chi1*chi2 + 10206*chi12 + 10206*chi22)*m22 + 59472*m24))/3.274425e7)*powers_of_2d1.seven_thirds*powers_of_2d1.seven_thirds*powers_of_lalpiHM.seven_thirds*powers_of_lalpiHM.seven_thirds;
1114 pAmp->xdot75 = -(m1*m2*(-4036586400*chi1*m13 + 5821200*m2*(5861*chi1 + 1701*powers_of_lalpiHM.itself)*m13 + 17059026600*chi1*m14 + 14721814800*chi1*m2*m14 - 34962127200*chi1*m2*m15 +
1115 207900*(2*chi2*m2*(-9708 + 41027*m2) + 12477*powers_of_lalpiHM.itself)*m22 - 14721814800*chi2*m13*m22 - 69924254400*chi1*m14*m22 - 34962127200*(chi1 + chi2)*m13*m23 +
1116 207900*m12*(3*powers_of_lalpiHM.itself*(4159 + 31752*m22) + 2*m2*(9708*chi2 + 41027*(chi1 + chi2)*m2 - 35406*chi1*m22 - 168168*chi2*m23)) +
1117 415800*m1*m2*(9708*chi1*m2 + 82054*chi2*m22 + 3*powers_of_lalpiHM.itself*(4159 + 7938*m22) + 35406*chi2*m23 - 84084*chi2*m24))/3.274425e7)*powers_of_2d1.eight_thirds*powers_of_2d1.seven_thirds*powers_of_lalpiHM.eight_thirds*powers_of_lalpiHM.seven_thirds;
1118 pAmp->xdot8 = -(m1*m2*(- 10548014400*chi1*powers_of_lalpiHM.itself*m13 - 63392868000*chi1*chi2*m2*m14 + 34927200*chi1*(-375*chi1 + 752*powers_of_lalpiHM.itself)*m14 + 63392868000*chi12*m15-
1119 153213984000*m2*chi12*m15 - 76606992000*chi12*m16 - 63392868000*chi1*(chi1 - chi2)*m13*m22 -
1120 51975*(4869 + 2711352*chi1*chi2 + 1702428*chi12 + 228508*chi22)*m14*m22 - 103950*(4869 + 2711352*chi1*chi2 + 228508*chi12 + 228508*chi22)*m13*m23 +
1121 906328500*m15*m23 + 1812657000*m14*m24 + 906328500*m13*m25 +
1122 1925*m2*m13*(56198689 + 13635864*chi1*chi2 + 27288576*chi1*powers_of_lalpiHM.itself + 30746952*chi12 + 3617892*chi22 - 2045736*powers_of_lalpiHM.two) -
1123 3*m22*(16447322263 - 2277918720*EulerGamma - 23284800*chi2*m2*(-151 + 376*m2)*powers_of_lalpiHM.itself - 2277918720*log4 + 2321480700*chi22 + 4365900000*chi22*m22 -
1124 21130956000*chi22*m23 + 25535664000*chi22*m24 + 745113600*powers_of_lalpiHM.two) +
1125 m12*(6833756160*EulerGamma + 10548014400*chi2*m2*powers_of_lalpiHM.itself + 63392868000*(chi1 - chi2)*chi2*m23 - 51975*(4869 + 2711352*chi1*chi2 + 228508*chi12 + 1702428*chi22)*m24 -
1126 3850*m22*(-56198689 + 13580136*chi1*chi2 - 6822144*(chi1 + chi2)*powers_of_lalpiHM.itself - 6976422*chi12 - 6976422*chi22 + 2045736*powers_of_lalpiHM.two) -
1127 3*(16447322263 - 2277918720*log4 + 2321480700*chi12 + 745113600*powers_of_lalpiHM.two)) +
1128 m1*m2*(13667512320*EulerGamma + 10548014400*chi1*m2*powers_of_lalpiHM.itself - 63392868000*chi1*chi2*m23 - 153213984000*chi22*m24 -
1129 1925*m22*(-56198689 - 13635864*chi1*chi2 - 27288576*chi2*powers_of_lalpiHM.itself - 3617892*chi12 - 30746952*chi22 + 2045736*powers_of_lalpiHM.two) -
1130 6*(16447322263 - 2277918720*log4 + 1160740350*chi12 + 1160740350*chi22 + 745113600*powers_of_lalpiHM.two)))/3.274425e7)*powers_of_2d1.eight_thirds*powers_of_2d1.eight_thirds*powers_of_lalpiHM.eight_thirds*powers_of_lalpiHM.eight_thirds;
1131 pAmp->xdot8Log = -(m1*m2*3416878080)/3.274425e7*powers_of_2d1.eight_thirds*powers_of_2d1.eight_thirds*powers_of_lalpiHM.eight_thirds*powers_of_lalpiHM.eight_thirds;
1132 pAmp->xdot85 = -(m1*m2*(-14891068500*chi1*m13 + 1925*m2*(97151928*chi1 + 6613488*chi2 - 12912300*powers_of_lalpiHM.itself)*m13 + 87143248500*chi1*m14 + 33313480200*chi1*m2*m14 - 198816225300*chi1*m2*m15 +
1133 57750*(2*chi2*m2*(-128927 + 754487*m2) + 7947*powers_of_lalpiHM.itself)*m22 - 33313480200*chi2*m13*m22 - 138600*(3399633*chi1 + 530712*chi2 + 182990*powers_of_lalpiHM.itself)*m14*m22 -
1134 35665037100*chi1*m15*m22 + 84184254000*chi1*m16*m22 -
1135 23100*m1*m2*(15*powers_of_lalpiHM.itself*(-2649 + 71735*m22) + m2*(-(chi1*(644635 + 551124*m2)) + chi2*m2*(-8095994 - 1442142*m2 + 8606763*m22))) -
1136 69300*(4991769*(chi1 + chi2) + 731960*powers_of_lalpiHM.itself)*m13*m23 + 35665037100*chi2*m14*m23 + 170726094000*chi1*m15*m23 + 2357586000*chi2*m15*m23 +
1137 35665037100*chi1*m13*m24 + 88899426000*(chi1 + chi2)*m14*m24 + 9702000*(243*chi1 + 17597*chi2)*m13*m25 -
1138 11550*m12*(15*powers_of_lalpiHM.itself*(-2649 + 286940*m22 + 146392*m24) + 2*m2*
1139 (-644635*chi2 - 4874683*(chi1 + chi2)*m2 + 1442142*chi1*m22 + 54*(58968*chi1 + 377737*chi2)*m23 + 1543941*chi2*m24 - 3644340*chi2*m25)))/3.274425e7)*powers_of_2d1.eight_thirds*powers_of_2d1.eight_thirds*powers_of_2d1.one_third*powers_of_lalpiHM.eight_thirds*powers_of_lalpiHM.eight_thirds*powers_of_lalpiHM.one_third;
1140
1141 pAmp->log2pi_two_thirds = log(powers_of_lalpiHM.two_thirds * powers_of_2d1.two_thirds);
1142
1143 }
1144
1145/*** PHENOM AMPLITUDE COEFFICIENTS ***/
1147
1148 /* Take all the phenom coefficients accross the three regions (inspiral, intermeidate and ringdown) and all the needed parameters to reconstruct the amplitude (including mode-mixing). */
1149 pAmp->ampNorm = pWF22->ampNorm;
1150 pAmp->PNdominant = pWF22->ampNorm * pow(2./pWFHM->emm, -7/6.); // = Pi * Sqrt(2 eta/3) (2Pi /m)^(-7/6). Miss the f^(-7/6). The pi power included in ampNorm
1151 switch(pWFHM->modeTag){
1152 case 21:{
1153 if(pWF22->q == 1){
1154 pAmp->PNdominantlmpower = 2;
1155 pAmp->PNdominantlm = sqrt(2.) / 3. * 1.5 * pWF22->dchi * 0.5 * pow(2 * LAL_PI / pWFHM->emm, 2/3.);
1156 }
1157 else{
1158 pAmp->PNdominantlmpower = 1;
1159 pAmp->PNdominantlm = sqrt(2.) / 3. * pWF22->delta * pow(2 * LAL_PI / pWFHM->emm, 1/3.);
1160 }
1161 break;
1162 }
1163 case 33:{
1164 if(pWF22->q == 1){
1165 pAmp->PNdominantlmpower = 4;
1166 pAmp->PNdominantlm = 0.75 * sqrt(5./7) * pWF22->dchi * 0.5 * (65/24. - 28/3. * pWF22->eta) * pow(2 * LAL_PI / pWFHM->emm, 4/3.);
1167 }
1168 else{
1169 pAmp->PNdominantlmpower = 1;
1170 pAmp->PNdominantlm = 0.75 * sqrt(5./7) * pWF22->delta * pow(2 * LAL_PI / pWFHM->emm, 1/3.);
1171 }
1172 break;
1173 }
1174 case 32:{
1175 pAmp->PNdominantlmpower = 2;
1176 pAmp->PNdominantlm = 0.75 * sqrt(5./7) * pow(2 * LAL_PI / pWFHM->emm, 1/3.);
1177 break;
1178 }
1179 case 44:{
1180 pAmp->PNdominantlmpower = 2;
1181 pAmp->PNdominantlm = 4/9. * sqrt(10/7.) * (1 - 3 * pWF22->eta) * pow(2 * LAL_PI / pWFHM->emm, 1/3.);
1182 break;
1183 }
1184 default:
1185 {XLALPrintError("Error in IMRPhenomXHM_GetAmplitudeCoefficients: mode selected is not currently available. Modes available are ((2,|2|),(2,|1|),(3,|2|),(3,|3|),(4,|4|)).\n");}
1186 }
1187 pAmp->PNdominantlm = fabs(pAmp->PNdominantlm);
1188 pWFHM->fAmpRDfalloff = 0.;
1189 pAmp->nCoefficientsInter = 0;
1190
1191 /*** Proceed region by region ***/
1192 if(pWFHM->IMRPhenomXHMReleaseVersion != 122019){
1193 pAmp->InspRescaleFactor = 0;
1194 pAmp->RDRescaleFactor = 0;
1195 pAmp->InterRescaleFactor = 0;
1196
1197
1198 /* Transform IMRPhenomXHMIntermediateAmpVersion number to int array defining what to do for each collocation point */
1199 /* 0: don't use coll point, 1: use point, 2: use point and derivative (this only applies for boundaries) */
1200 // e.g. pAmp->VersionCollocPtsInter = {1, 1, 1, 1, 0, 2} -> use the two boundaries, add derivative to the right one, skip third collocation point
1202 for(UINT2 i = 0; i < pWFHM->nCollocPtsInterAmp; i++){
1203 pAmp->VersionCollocPtsInter[pWFHM->nCollocPtsInterAmp - i - 1] = num % 10;
1204 num/=10;
1205 }
1206
1207 pAmp->nCoefficientsInter = 0;
1208 for (UINT2 i = 0; i < pWFHM->nCollocPtsInterAmp; i++) pAmp->nCoefficientsInter += pAmp->VersionCollocPtsInter[i];
1209 /* The number of coefficients in the intermediate ansatz cannot be larger than the number of available collocation points in IMRPhenomXHMIntermediateAmpVersion
1210 If e.g. IMRPhenomXHMIntermediateAmpVersion has 6 slots, the maximum number of coefficients would be 6 + 2 because we count for the derivatives at the boundaries. */
1211 if (pAmp->nCoefficientsInter > pWFHM->nCollocPtsInterAmp + 2)
1212 XLAL_ERROR_VOID(XLAL_EFUNC, "IMRPhenomXHM_GetAmplitudeCoefficients failed. Inconsistent number of collocation points (%i) and free parameters (%i).", pWFHM->nCollocPtsInterAmp + 2, pAmp->nCoefficientsInter);
1213
1214 pAmp->nCoefficientsRDAux = 0;
1215 if (pWFHM->MixingOn){
1216 pAmp->nCollocPtsRDAux = 2;
1217 pAmp->nCoefficientsRDAux = 4;
1218 pAmp->fRDAux = pWFHM->fRING - pWFHM->fDAMP; //v2
1219 }
1220
1221 // Take the cutting frequencies at the inspiral and ringdown
1222 pAmp->fAmpMatchIN = IMRPhenomXHM_Amplitude_fcutInsp(pWFHM, pWF22);
1223 pAmp->fAmpMatchIM = IMRPhenomXHM_Amplitude_fcutRD(pWFHM, pWF22);
1224 pWFHM->fAmpRDfalloff = pWFHM->fRING + 2 * pWFHM->fDAMP;
1225
1226 /* Take Frequency Domain Post-Newtonian Amplitude Coefficients */
1227 IMRPhenomXHM_GetPNAmplitudeCoefficients(pAmp, pWFHM, pWF22);
1228
1230
1231 IMRPhenomXHM_RD_Amp_Coefficients(pWF22, pWFHM, pAmp);
1232
1233 IMRPhenomXHM_Intermediate_Amp_Coefficients(pAmp, pWFHM, pWF22, pPhase, pAmp22, pPhase22);
1234
1235 // printf("\nInsp Coll points\n");
1236 // for(UINT2 i = 0; i < 3; i++){
1237 // printf("%.16f %.16e\n", pAmp->CollocationPointsFreqsAmplitudeInsp[i], pAmp->CollocationPointsValuesAmplitudeInsp[i]);
1238 // }
1239 // printf("\nInter Coll points\n");
1240 // for(UINT2 i = 0; i < pAmp->nCoefficientsInter; i++){
1241 // printf("%.16f %.16e\n", pAmp->CollocationPointsFreqsAmplitudeInter[i], pAmp->CollocationPointsValuesAmplitudeInter[i]);
1242 // }
1243 // printf("\nRD Coll points\n");
1244 // for(UINT2 i = 0; i < 3; i++){
1245 // printf("%.16f %.16e\n", pAmp->CollocationPointsFreqsAmplitudeRD[i], pAmp->CollocationPointsValuesAmplitudeRD[i]);
1246 // }
1247 //
1248 // printf("\nInsp Coefficients\n");
1249 // for(UINT2 i = 0; i < 3; i++){
1250 // printf("%.16e\n", pAmp->InspiralCoefficient[i]);
1251 // }
1252 // printf("\nInter Coefficients %i\n", pAmp->nCoefficientsInter);
1253 // for(UINT2 i = 0; i < pAmp->nCoefficientsInter; i++){
1254 // printf("%.16e\n", pAmp->InterCoefficient[i]);
1255 // }
1256 // printf("\nRD Coefficients\n");
1257 // for(UINT2 i = 0; i < 5; i++){
1258 // printf("%.16e\n", pAmp->RDCoefficient[i]);
1259 // }
1260 // printf("\nRDAux Coll points\n");
1261 // for(UINT2 i = 0; i < pAmp->nCoefficientsRDAux; i++){
1262 // printf("%.16f %.16e\n", pAmp->CollocationPointsFreqsAmplitudeRDAux[i], pAmp->CollocationPointsValuesAmplitudeRDAux[i]);
1263 // }
1264 // printf("\nRDAux Coefficients\n");
1265 // for(UINT2 i = 0; i < pAmp->nCoefficientsRDAux; i++){
1266 // printf("%.16e\n", pAmp->RDAuxCoefficient[i]);
1267 // }
1268
1269
1270 /* Set Rescale Factors to build Strain (=0) */
1271 pAmp->InspRescaleFactor = 0;
1272 pAmp->RDRescaleFactor = 0;
1273 pAmp->InterRescaleFactor = 0;
1274 }
1275 else{
1276 pAmp->InspRescaleFactor = 1;
1277 pAmp->InterRescaleFactor = -1;
1278 pAmp->RDRescaleFactor = 1;
1279
1280 // Options for the extrapolation of the model outside the calibration region
1281 if(((pWFHM->modeTag==44) || (pWFHM->modeTag==33)) && pWF22->q>7. && pWF22->chi1L > 0.95){
1282 pAmp->useInspAnsatzRingdown = 1;
1283 }
1284 else{
1285 pAmp->useInspAnsatzRingdown = 0;
1286 }
1287 pAmp->WavyInsp = 0;
1288 pAmp->WavyInt = 0;
1289 if(pWFHM->modeTag==21){
1290 pAmp->WavyInsp = 1;
1291 pAmp->WavyInt = 1;
1292 }
1293 if(pWFHM->modeTag==32){
1294 pAmp->WavyInsp = 1;
1295 }
1296
1297 // Take the cutting frequencies at the inspiral and ringdown
1298 pAmp->fAmpMatchIN = IMRPhenomXHM_Amplitude_fcutInsp(pWFHM, pWF22);
1299 pAmp->fAmpMatchIM = IMRPhenomXHM_Amplitude_fcutRD(pWFHM, pWF22);
1300
1301 #if DEBUG == 1
1302 printf("\n\n*** IMRPhenomXHM_GetAmplitudeCoefficients ***\n\n");
1303 printf("fring_%i = %.16f\n", pWFHM->modeTag, pWFHM->fRING);
1304 printf("fdamp_%i = %.16f\n", pWFHM->modeTag, pWFHM->fDAMP);
1305 printf("fcutInsp_%i = %.16f \n", pWFHM->modeTag, pAmp->fAmpMatchIN);
1306 printf("fcutRD_%i = %.16f \n", pWFHM->modeTag, pAmp->fAmpMatchIM);
1307 #endif
1308
1309 /* Compute the frequencies for Intermediate collocation points. */
1310 /* It is done now because with the inspiral veto fAmpMatchIN will change, and we need the original here. */
1311 double df = pAmp->fAmpMatchIM - pAmp->fAmpMatchIN;
1312 pAmp->CollocationPointsFreqsAmplitudeInter[0] = pAmp->fAmpMatchIN + df/3. ;
1313 pAmp->CollocationPointsFreqsAmplitudeInter[1] = pAmp->fAmpMatchIN + df*2./3.;
1314
1315 int nCollocPtsInspAmp = pWFHM->nCollocPtsInspAmp;
1316 int nCollocPtsInterAmp = pWFHM->nCollocPtsInterAmp;
1317
1318 int modeint = pWFHM->modeInt;
1319
1320 #if DEBUG == 1
1321 printf("nCollocPtsInspAmp = %i \n",nCollocPtsInspAmp);
1322 printf("nCollocPtsInterAmp = %i \n",nCollocPtsInterAmp);
1323 #endif
1324
1325 /*****************/
1326 /* INSPIRAL */
1327 /*****************/
1328
1329 #if DEBUG == 1
1330 printf("\n**** INSPIRAL ****\n\n");
1331 printf("IMRPhenomXHMInspiralAmpVersion = %i\r\n",pWFHM->IMRPhenomXHMInspiralAmpVersion);
1332 #endif
1333
1334 /* Take Frequency Domain Post-Newtonian Amplitude Coefficients */
1335 IMRPhenomXHM_GetPNAmplitudeCoefficients(pAmp, pWFHM, pWF22);
1336
1337 #if DEBUG == 1
1338 printf("PN coeff %.16f %.16f\n", creal(pAmp->pnInitial),cimag(pAmp->pnInitial));
1339 printf("PN coeff %.16f %.16f\n", creal(pAmp->pnOneThird),cimag(pAmp->pnOneThird));
1340 printf("PN coeff %.16f %.16f\n", creal(pAmp->pnTwoThirds),cimag(pAmp->pnTwoThirds));
1341 printf("PN coeff %.16f %.16f\n", creal(pAmp->pnThreeThirds),cimag(pAmp->pnThreeThirds));
1342 printf("PN coeff %.16f %.16f\n", creal(pAmp->pnFourThirds),cimag(pAmp->pnFourThirds));
1343 printf("PN coeff %.16f %.16f\n", creal(pAmp->pnFiveThirds),cimag(pAmp->pnFiveThirds));
1344 printf("PN coeff %.16f %.16f\n", creal(pAmp->pnSixThirds),cimag(pAmp->pnSixThirds));
1345 #endif
1346
1347 // Initialize frequencies of colloc points !!!! ASSUMING YOU HAVE 3 COLLOC POINTS
1348 pAmp->CollocationPointsFreqsAmplitudeInsp[0] = 1.0 * pAmp->fAmpMatchIN;
1349 pAmp->CollocationPointsFreqsAmplitudeInsp[1] = 0.75 * pAmp->fAmpMatchIN;
1350 pAmp->CollocationPointsFreqsAmplitudeInsp[2] = 0.5 * pAmp->fAmpMatchIN;
1351
1352 double fcutInsp, f1, f2, f3;
1353 fcutInsp = pAmp->fAmpMatchIN; // Matching frequency between inspiral and intermediate
1354 f1 = pAmp->CollocationPointsFreqsAmplitudeInsp[0]; // Frequency of colloc point 1 = 1.*fcutInsp
1355 f2 = pAmp->CollocationPointsFreqsAmplitudeInsp[1]; // Frequency of colloc point 2 = 0.75*fcutInsp
1356 f3 = pAmp->CollocationPointsFreqsAmplitudeInsp[2]; // Frequency of colloc point 3 = 0.5*fcutInsp
1357
1358 // Compute the useful powers of fcutInsp, f1, f2, f3. Remembers: f3 < f2 < f1 = fcutInsp.
1359 IMRPhenomX_UsefulPowers powers_of_fcutInsp, powers_of_f1, powers_of_f2, powers_of_f3;
1360 IMRPhenomX_Initialize_Powers(&powers_of_fcutInsp, fcutInsp); // fcutInsp and f1 are the same but we keep them separated if in the future we change the frequencies of the collocatio points.
1361 IMRPhenomX_Initialize_Powers(&powers_of_f1, f1);
1362 IMRPhenomX_Initialize_Powers(&powers_of_f2, f2);
1363 IMRPhenomX_Initialize_Powers(&powers_of_f3, f3);
1364
1365 // Compute the values of Post-Newtoninan ansatz (without the pseudo-PN terms) at the frequencies of the collocation points.
1366 double PNf1, PNf2, PNf3;
1367 PNf1 = IMRPhenomXHM_Inspiral_PNAmp_Ansatz(&powers_of_f1, pWFHM, pAmp);
1368 PNf2 = IMRPhenomXHM_Inspiral_PNAmp_Ansatz(&powers_of_f2, pWFHM, pAmp);
1369 PNf3 = IMRPhenomXHM_Inspiral_PNAmp_Ansatz(&powers_of_f3, pWFHM, pAmp);
1370
1371 pAmp->PNAmplitudeInsp[0] = PNf1;
1372 pAmp->PNAmplitudeInsp[1] = PNf2;
1373 pAmp->PNAmplitudeInsp[2] = PNf3;
1374
1375 // Initialize values of collocation points at the previous 3 frequencies. They are taken from the parameter space fits.
1376 for(int i = 0; i<nCollocPtsInspAmp; i++)
1377 {
1378 pAmp->CollocationPointsValuesAmplitudeInsp[i] = fabs(pAmp->InspiralAmpFits[modeint*nCollocPtsInspAmp+i](pWF22,pWFHM->IMRPhenomXHMInspiralAmpFitsVersion));
1379 }
1380
1381
1382
1383 // Values of the collocation point minus the Post-Newtonian value. This gives a "collocation point" for the pseudo-PN part.
1384 // This way is more convenient becuase the reconstuction does not depended on the PN ansatz used.
1385 REAL8 iv1, iv2, iv3;
1386 iv1 = pAmp->CollocationPointsValuesAmplitudeInsp[0] - PNf1;
1387 iv2 = pAmp->CollocationPointsValuesAmplitudeInsp[1] - PNf2;
1388 iv3 = pAmp->CollocationPointsValuesAmplitudeInsp[2] - PNf3;
1389
1390
1391 #if DEBUG == 1
1392 printf("\nAmplitude pseudo collocation points before veto\n");
1393 printf("fAmpMatchIN = %.16f\n",pAmp->fAmpMatchIN);
1394 printf("F1 = %.16f\n", f1);
1395 printf("F2 = %.16f\n", f2);
1396 printf("F3 = %.16f\n\n", f3);
1397 printf("I1 = %.16f\n", pAmp->CollocationPointsValuesAmplitudeInsp[0]);
1398 printf("I2 = %.16f\n", pAmp->CollocationPointsValuesAmplitudeInsp[1]);
1399 printf("I3 = %.16f\n\n", pAmp->CollocationPointsValuesAmplitudeInsp[2]);
1400 printf("PNf1 = %.16f\n", PNf1);
1401 printf("PNf2 = %.16f\n", PNf2);
1402 printf("PNf3 = %.16f\n\n", PNf3);
1403 REAL8 piv1, piv2, piv3;
1404 piv1 = pAmp->CollocationPointsValuesAmplitudeInsp[0]*pWF22->ampNorm*powers_of_f1.m_seven_sixths;
1405 piv2 = pAmp->CollocationPointsValuesAmplitudeInsp[1]*pWF22->ampNorm*powers_of_f2.m_seven_sixths;
1406 piv3 = pAmp->CollocationPointsValuesAmplitudeInsp[2]*pWF22->ampNorm*powers_of_f3.m_seven_sixths;
1407 printf("p1 = %.16f\n", piv1);
1408 printf("p2 = %.16f\n", piv2);
1409 printf("p3 = %.16f\n\n", piv3);
1410 printf("V1 = %.16f\n", iv1);
1411 printf("V2 = %.16f\n", iv2);
1412 printf("V3 = %.16f\n\n", iv3);
1413 #endif
1414
1415 /*
1416 VETO: Choose the collocations points to use.
1417 Outside of the callibration region the collocation points can have a wavy behaviour so we remove some of them.
1418 */
1419 if(pWFHM->InspiralAmpVeto == 1){
1420 IMRPhenomXHM_Inspiral_Amplitude_Veto(&iv1, &iv2, &iv3, &powers_of_f1, &powers_of_f2, &powers_of_f3, pAmp, pWFHM);
1421 }
1422 #if DEBUG == 1
1423 printf("\nInspiral Veto: AmpVersion = %i",pWFHM->IMRPhenomXHMInspiralAmpVersion);
1424 #endif
1425 /* The 32 mode in this corner of the parameter space always uses the PN only */
1426 if(pWFHM->modeTag==32 && pWF22->q>2.5 && pWF22->chi1L<-0.9 && pWF22->chi2L<-0.9) pWFHM->IMRPhenomXHMInspiralAmpVersion = 0;
1427 /* The 32 mode in this corner of the parameter space removes the last collocation point */
1428 if(pWFHM->modeTag==32 && pWF22->q>2.5 && pWF22->chi1L<-0.6 && pWF22->chi2L>0. && iv1!=0){
1430 iv1 = 0.;
1431 }
1432 /* The 33 mode in this corner of the parameter space removes the last collocation point */
1433 if(pWFHM->modeTag==33 && (1.2>pWF22->q && pWF22->q>1. && pWF22->chi1L<-0.1 && pWF22->chi2L>0. && iv1!=0)){//November
1435 iv1 = 0.;
1436 }
1437 #if DEBUG == 1
1438 printf("\nInspiral Veto: AmpVersion = %i", pWFHM->IMRPhenomXHMInspiralAmpVersion);
1439 #endif
1440 //********* Remember that f3 < f2 < f1 !!!!! *****************
1441 /* Check for wavy collocation points. Only when we have 3 coll points. */
1442 if(pWFHM->IMRPhenomXHMInspiralAmpVersion == 3 && pAmp->WavyInsp == 1){
1444 iv2 = 0;
1445 #if DEBUG == 1
1446 printf("\nWavy Inspiral\n");
1447 #endif
1449 }
1450 }
1451 #if DEBUG == 1
1452 printf("\nInspiral Veto: AmpVersion = %i",pWFHM->IMRPhenomXHMInspiralAmpVersion);
1453 printf("\niv1 iv2 iv3 %e %e %e\n", iv1, iv2, iv3);
1454 #endif
1455 // Rename collocation points and frequencies.
1456 if(iv2==0){
1457 #if DEBUG == 1
1458 printf("\niv2 = 0\n");
1459 #endif
1460 iv2 = iv3;
1461 iv3 = 0.;
1462 powers_of_f2 = powers_of_f3;
1463 }
1464 if(iv1==0){
1465 #if DEBUG == 1
1466 printf("\niv1 = 0\n");
1467 #endif
1468 iv1 = iv2;
1469 powers_of_f1 = powers_of_f2;
1470 powers_of_fcutInsp = powers_of_f1;
1471 iv2 = iv3;
1472 iv3 = 0.;
1473 powers_of_f2 = powers_of_f3;
1474 }
1475 if(pWFHM->IMRPhenomXHMInspiralAmpVersion == 0) // when using PN we take fcutInsp = fMECO_lm
1476 {
1477 IMRPhenomX_Initialize_Powers(&powers_of_fcutInsp, (pWFHM->fMECOlm));
1478 }
1479
1480 // Update the Inspiral cutting frequency to the last collocation point alive.
1481 pAmp->fAmpMatchIN = powers_of_fcutInsp.itself;
1482 // End of VETO
1483
1484
1485 #if DEBUG == 1
1486 printf("\nAmplitude pseudo collocation points after veto\n");
1487 printf("fAmpMatchIN = %.16f\n", pAmp->fAmpMatchIN);
1488 printf("F1 = %.16f\n", powers_of_f1.itself);
1489 printf("F2 = %.16f\n", powers_of_f2.itself);
1490 printf("F3 = %.16f\n", powers_of_f3.itself);
1491 printf("V1 = %.16f\n", iv1);
1492 printf("V2 = %.16f\n", iv2);
1493 printf("V3 = %.16f\n", iv3);
1494 #endif
1495
1496 /* Get the pseudo-PN coefficients. */
1497 /* The whole inspiral ansatz is given by PNAnsatz(f) + pseudo-PN(f), with pseudo-PN(f) = rho1 *(f/fcutInsp)^(7/3) + rho2(f/fcutInsp)^(8/3) + rho3(f/fcutInsp)^(9/3).
1498 The coefficients are computed demanding that the pseudo-PN part at the three collocation points frequencies returns the actual collocation points: iv1, iv2, iv3. */
1499 pAmp->rho1 = IMRPhenomXHM_Inspiral_Amp_rho1(iv1, iv2, iv3, &powers_of_fcutInsp, &powers_of_f1, &powers_of_f2, &powers_of_f3, pWFHM);
1500 pAmp->rho2 = IMRPhenomXHM_Inspiral_Amp_rho2(iv1, iv2, iv3, &powers_of_fcutInsp, &powers_of_f1, &powers_of_f2, &powers_of_f3, pWFHM);
1501 pAmp->rho3 = IMRPhenomXHM_Inspiral_Amp_rho3(iv1, iv2, iv3, &powers_of_fcutInsp, &powers_of_f1, &powers_of_f2, &powers_of_f3, pWFHM);
1502
1503 #if DEBUG == 1
1504 printf("\nAmplitude pseudo PN coeffcients after veto\n");
1505 printf("rho1 = %.16f\n",pAmp->rho1);
1506 printf("rho2 = %.16f\n",pAmp->rho2);
1507 printf("rho3 = %.16f\n",pAmp->rho3);
1508 #endif
1509
1510 // To avoid passing an extra argument to IMRPhenomXHM_Inspiral_Amp_Ansatz we store this powers in pAmp.
1511 pAmp->fcutInsp_seven_thirds = powers_of_fcutInsp.seven_thirds;
1512 pAmp->fcutInsp_eight_thirds = powers_of_fcutInsp.eight_thirds;
1513 pAmp->fcutInsp_three = powers_of_fcutInsp.three;
1514
1515
1516 /*****************/
1517 /* RINGDOWN */
1518 /*****************/
1519 #if DEBUG == 1
1520 printf("\n\n**** RINGDOWN ****\n\n");
1521 #endif
1522
1523 // We have three "fitted" coefficients across parameter space: alambda, lambda and sigma (RDCoefficients[0,1,2]). Sigma will be constat for all the modes except the 21.
1524 pAmp->RDCoefficient[0] = fabs(pAmp->RingdownAmpFits[modeint*3](pWF22,pWFHM->IMRPhenomXHMRingdownAmpFitsVersion));
1525 pAmp->RDCoefficient[1] = pAmp->RingdownAmpFits[modeint*3+1](pWF22,pWFHM->IMRPhenomXHMRingdownAmpFitsVersion);
1526 pAmp->RDCoefficient[2] = pAmp->RingdownAmpFits[modeint*3+2](pWF22,pWFHM->IMRPhenomXHMRingdownAmpFitsVersion);
1527 pAmp->RDCoefficient[3] = 1./12.;
1528
1529 /* (IMRPhenomXPNRUseTunedCoprec) Amplitude deviations e.g. for PNR's copreessing model*/
1530 if(pWFHM->modeTag==33){
1531 // pAmp->lambda = pAmp->lambda + ( pWFHM->PNR_DEV_PARAMETER * pWFHM->MU1 );
1532 pAmp->RDCoefficient[0] = pAmp->RDCoefficient[0] + ( pWFHM->PNR_DEV_PARAMETER * pWFHM->MU3 );
1533 pAmp->RDCoefficient[2] = pAmp->RDCoefficient[2] + ( pWFHM->PNR_DEV_PARAMETER * pWFHM->MU4 );
1534 }
1535
1536 #if DEBUG == 1
1537 printf("\nuseInspAnsatzRigndown = %i\n",pAmp->useInspAnsatzRingdown);
1538 printf("alambda = %.16f\r\n",pAmp->RDCoefficient[0]);
1539 #endif
1540
1541 // For some cases with extreme spins there is almost no merger and the transition to the ringdown is very sharp.
1542 // The coefficients of the ringdown do not work well here and we take an approximation of the inspiral.
1543 if(pAmp->useInspAnsatzRingdown==1){
1544 // The Ringdown amp at fAmpMatchIM is set to be 0.9 the amplitude in the last inspiral collocation point
1545 IMRPhenomX_UsefulPowers powers_of_fRD;
1546 IMRPhenomX_Initialize_Powers(&powers_of_fRD, pAmp->fAmpMatchIM);
1547 pAmp->RDCoefficient[0] = 0.9*fabs(IMRPhenomXHM_Inspiral_Amp_Ansatz(&powers_of_fcutInsp, pWFHM, pAmp)/(IMRPhenomXHM_RD_Amp_Ansatz(&powers_of_fRD, pWFHM, pAmp)/pAmp->RDCoefficient[0]));
1548 }
1549
1550 #if DEBUG == 1
1551 printf("alambda = %.16f\r\n",pAmp->RDCoefficient[0]);
1552 printf("lambda = %.16f\r\n",pAmp->RDCoefficient[1]);
1553 printf("sigma = %.16f\r\n",pAmp->RDCoefficient[2]);
1554 printf("lc = %.16f\r\n",pAmp->RDCoefficient[3]);
1555 #endif
1556
1557
1558 /*******************/
1559 /* INTERMEDIATE */
1560 /*******************/
1561
1562 #if DEBUG == 1
1563 printf("\n\n**** INTERMEDIATE ****\n\n");
1564 #endif
1565
1566 // Initialize values of collocation points. They are taken from the parameter space fits.
1567 for(int i = 0; i<nCollocPtsInterAmp; i++)
1568 {
1569 pAmp->CollocationPointsValuesAmplitudeInter[i] = fabs(pAmp->IntermediateAmpFits[modeint*nCollocPtsInterAmp+i](pWF22,pWFHM->IMRPhenomXHMIntermediateAmpFitsVersion));
1570 }
1571
1572 /* For reconstructing the intermediate region we need value and derivative at the beginning and at the end of the intermediate region, given by the inspiral and ringdown ansatzs.
1573 And also the two values at the two intermediate collocation points. This gives 6 parameters, that will determine the 6 coefficients of the inverse of the 5th order polynomial
1574 that we use for reconstruction.
1575 This is the default behaviour, however the veto conditions can reduce the order of the polynomial.
1576 */
1577
1578 double F1, F2, F3, F4; //Frequencies of the collocation points in the intermediate region
1579 F1 = powers_of_fcutInsp.itself;
1582 F4 = pAmp->fAmpMatchIM;
1583
1584 #if DEBUG == 1
1585 printf("F1 = %.16f\n",F1);
1586 printf("F2 = %.16f\n",F2);
1587 printf("F3 = %.16f\n",F3);
1588 printf("F4 = %.16f\n\n",F4);
1589 #endif
1590
1591 // Initialize useful powers.
1592 IMRPhenomX_UsefulPowers powers_of_F1;
1593 IMRPhenomX_Initialize_Powers(&powers_of_F1,F1);
1594 IMRPhenomX_UsefulPowers powers_of_F4;
1595 IMRPhenomX_Initialize_Powers(&powers_of_F4,F4);
1596
1597 // Compute values at the boundaries (rescaled ansatz with the leading order of the 22).
1598 double inspF1 = IMRPhenomXHM_Inspiral_Amp_Ansatz(&powers_of_F1, pWFHM, pAmp);
1599 double rdF4;
1600 if (pWFHM->MixingOn == 1){
1601 rdF4 = cabs(SpheroidalToSpherical(&powers_of_F4, pAmp22, pPhase22, pAmp, pPhase, pWFHM, pWF22));
1602 }else{
1603 rdF4 = IMRPhenomXHM_RD_Amp_Ansatz(&powers_of_F4, pWFHM, pAmp);
1604 }
1605
1606 // Compute derivatives at the boundaries (rescaled ansatz with the leading order of the 22).
1607 double d1, d4;
1608 d1 = IMRPhenomXHM_Inspiral_Amp_NDAnsatz(&powers_of_F1, pWFHM, pAmp);
1609 if(pWFHM->MixingOn==1){
1610 d4 = IMRPhenomXHM_RD_Amp_NDAnsatz(&powers_of_F4, pAmp, pPhase, pWFHM, pAmp22, pPhase22, pWF22);
1611 }else{
1612 d4 = IMRPhenomXHM_RD_Amp_DAnsatz(&powers_of_F4, pWFHM, pAmp);
1613 }
1614 // Next use of Ansatz will be for return the full strain, set correct rescalefactor
1615 pAmp->InspRescaleFactor = 0;
1616 pAmp->InterRescaleFactor = 0;
1617 pAmp->RDRescaleFactor = 0;
1618
1619 #if DEBUG == 1
1620 printf("d1 = %.16f\n",d1);
1621 printf("d4 = %.16f\n",d4);
1622 #endif
1623
1624 // Derivatives at the boundaries of the whole ansatz.
1625 d1 = ((7.0/6.0) * pow(F1,1.0/6.0) / inspF1) - ( pow(F1,7.0/6.0) * d1 / (inspF1*inspF1) );
1626 d4 = ((7.0/6.0) * pow(F4,1.0/6.0) / rdF4) - ( pow(F4,7.0/6.0) * d4 / (rdF4*rdF4) );
1627
1628 #if DEBUG == 1
1629 printf("d1 = %.16f\n",d1);
1630 printf("d4 = %.16f\n\n",d4);
1631 #endif
1632
1633 // These values will feed the reconstruction function of the intermediate for getting the coefficients of the polynomial
1634 double V1, V2, V3, V4;
1635
1636 pWFHM->IMRPhenomXHMIntermediateAmpVersion = 105; // by default use 5th order polynomial
1637
1638 V1 = powers_of_F1.m_seven_sixths * inspF1;
1641 V4 = powers_of_F4.m_seven_sixths * rdF4;
1642
1643 #if DEBUG == 1
1644 printf("Before intermediate veto \n");
1645 printf("V1 = %.16f\n",V1);
1646 printf("V2 = %.16f\n",V2);
1647 printf("V3 = %.16f\n",V3);
1648 printf("V4 = %.16f\n",V4);
1649 printf("rdF4 = %.16f\n",rdF4);
1650 printf("F4.m_seven_sixths = %.16f\n\n",powers_of_F4.m_seven_sixths);
1651 #endif
1652
1653 // VETO: outside of the calibration region some collocation points are bad behaved. We veto them and change the type of reconstruction now.
1654
1655 if(pAmp->useInspAnsatzRingdown==1){
1656 //For these extreme cases we do not use intermediate collocation points -> third order polynomial.
1657 V2 = 1.;
1658 V3 = 1.;
1659 pWFHM->IMRPhenomXHMIntermediateAmpVersion = 101; /*Linear reconstruction*/
1660 #if DEBUG == 1
1661 printf("VETO: useInspAnsatzRingdown\n");
1662 printf("V2 = %.16f\n",V2);
1663 printf("V3 = %.16f\n",V3);
1664 #endif
1665 }
1666
1667 // The reconstruction function is the inverse of a polynomial. That we demand pass through the points V1, V2, V3, V4 and has the derivatives d1, d4 at the boundaries.
1668 // For simplicity we coded the reconstruction of the polynomial (Update_Intermediate_Amplitude_Coefficients), so we have to feed it with inverse values.
1669 V1 = 1.0 / V1;
1670 V2 = 1.0 / V2;
1671 V3 = 1.0 / V3;
1672 V4 = 1.0 / V4;
1673
1674 /* Apply NR tuning for precessing cases (IMRPhenomXPNRUseTunedCoprec) */
1675 if(pWFHM->modeTag==33){
1676 V2 = V2 + ( pWFHM->PNR_DEV_PARAMETER * pWFHM->MU1 );
1677 V3 = V3 + ( pWFHM->PNR_DEV_PARAMETER * pWFHM->MU2 );
1678 }
1679
1680
1681 #if DEBUG == 1
1682 printf("\nAfter Veto and inverse \n");
1683 printf("V1 = %.16f\n",V1);
1684 printf("V2 = %.16f\n",V2);
1685 printf("V3 = %.16f\n",V3);
1686 printf("V4 = %.16f\n",V4);
1687 printf("IMRPhenomXHMIntermediateAmpVersion = %i \n\n", pWFHM->IMRPhenomXHMIntermediateAmpVersion);
1688 #endif
1689
1690
1691 /* If the case does not have extreme mass ratio it skips the next block */
1692
1693
1694 /****** JUST EMR cases, 2 INTERMEDIATE REGIONS ******/
1695
1696 /* For EMR cases the amplitude shows a more pronounced drop off at the end of the inspiral.
1697 To model this we split the intermediate region in two.
1698 We add an extra collocation point between the end of the inspiral and the first intermediate collocation point, F0.
1699 From fcutInsp to FO we use a 4th order polynomial, and from F0 to fcutRD we use the usual 5th that is computed after this block. */
1700
1701 if(pWFHM->AmpEMR==1){
1702 #if DEBUG == 1
1703 printf("*** TWO INTERMEDIATE REGIONS ***\n\n");
1704 #endif
1705
1706 // Here we compute the first intermediate region with the inverse of a 4th polynomial.
1707 // For this we use point and derivative at fcutInsp and F0 (4 coefficients) and the value at the first collocation point (F2), with the total of 5 coefficients.
1708
1709 double F0, V0, d0, F0_seven_sixths; //Frequency, value, derivative and useful power at the extra collocation point for EMR
1710
1711 F0 = F1 + (F2-F1)/3.;
1712 pAmp->fAmpMatchInt12 = F0;
1713
1714 // Take the value and derivative from the parameter space fits.
1715 V0 = pAmp->IntermediateAmpFits[modeint*nCollocPtsInterAmp+8](pWF22,pWFHM->IMRPhenomXHMIntermediateAmpFitsVersion);
1716 d0 = pAmp->IntermediateAmpFits[modeint*nCollocPtsInterAmp+9](pWF22,pWFHM->IMRPhenomXHMIntermediateAmpFitsVersion);
1717
1718 #if DEBUG == 1
1719 printf("F0 = %.16f\n",F0);
1720 printf("V0 = %.16f\n",V0);
1721 printf("d0 = %.16f\n",d0);
1722 #endif
1723
1724 F0_seven_sixths = pow(F0,7.0/6.0);
1725
1726 d0 = ((7.0/6.0) / (V0*F0)) - ( d0 / (V0*V0*F0_seven_sixths) );
1727 V0 = 1. / V0;
1728
1729 #if DEBUG == 1
1730 printf("1/V0 = %.16f\n",V0);
1731 printf("1/d0 = %.16f\n",d0);
1732 #endif
1733
1734 // Get the coefficients of the polynomial for the first intermediate region
1735 pAmp->alpha0 = IMRPhenomXHM_Intermediate_Amp_delta0(d1,d0,V1,V2,V3,V0,F1,F2,F3,F0,104); //V3 and F3 will not be used when calling with 104
1736 pAmp->alpha1 = IMRPhenomXHM_Intermediate_Amp_delta1(d1,d0,V1,V2,V3,V0,F1,F2,F3,F0,104);
1737 pAmp->alpha2 = IMRPhenomXHM_Intermediate_Amp_delta2(d1,d0,V1,V2,V3,V0,F1,F2,F3,F0,104);
1738 pAmp->alpha3 = IMRPhenomXHM_Intermediate_Amp_delta3(d1,d0,V1,V2,V3,V0,F1,F2,F3,F0,104);
1739 pAmp->alpha4 = IMRPhenomXHM_Intermediate_Amp_delta4(d1,d0,V1,V2,V3,V0,F1,F2,F3,F0,104);
1740
1741 #if DEBUG == 1
1742 printf("Intermediate 1: feed values \n");
1743 printf("d1 = %.16f\n", d1);
1744 printf("d0 = %.16f\n", d0);
1745 printf("d4 = %.16f\n", d4);
1746 printf("V1 = %.16f\n", V1);
1747 printf("V0 = %.16f\n", V0);
1748 printf("V2 = %.16f\n", V2);
1749 printf("V3 = %.16f\n", V3);
1750 printf("V4 = %.16f\n", V4);
1751 printf("F1 = %.16f\n", F1);
1752 printf("F0 = %.16f\n", F0);
1753 printf("F2 = %.16f\n", F2);
1754 printf("F3 = %.16f\n", F3);
1755 printf("F4 = %.16f\n", F4);
1756 #endif
1757
1758 #if DEBUG == 1
1759 printf("\nIntermediate 1: polynomial coeffcients \r\n");
1760 printf("alpha0 = %.16f\n", pAmp->alpha0);
1761 printf("alpha1 = %.16f\n", pAmp->alpha1);
1762 printf("alpha2 = %.16f\n", pAmp->alpha2);
1763 printf("alpha3 = %.16f\n", pAmp->alpha3);
1764 printf("alpha4 = %.16f\n", pAmp->alpha4);
1765 #endif
1766
1767 //Update left collocation point for the 2nd intermediate region
1768 F1 = F0;
1769 V1 = V0;
1770 d1 = d0;
1771
1772 /**** END of first Intermediate region ****/
1773 }
1774 else{ /** This part is used both when we have a single intermediate region and for the second intermediate region **/
1775 pAmp->fAmpMatchInt12 = 0;
1776 pAmp->alpha0 = 1;
1777 pAmp->alpha1 = 1;
1778 pAmp->alpha2 = 1;
1779 pAmp->alpha3 = 1;
1780 pAmp->alpha4 = 1;
1781
1782 /** More vetos ***/
1783 if(pWFHM->IntermediateAmpVeto == 1 && pWFHM->IMRPhenomXHMIntermediateAmpVersion==105){ // only 21 mode
1784 IMRPhenomXHM_Intermediate_Amplitude_Veto(&V2, &V3, pWFHM, pWF22); // this changes the order of the polynomial to 4 or 3
1785 #if DEBUG == 1
1786 printf("VETO: Intermediate Amp Veto\n");
1787 printf("V2 = %.16f\n",V2);
1788 printf("V3 = %.16f\n",V3);
1789 #endif
1790 }
1791
1792 if(pWFHM->RingdownAmpVeto == 1){ // only 21, 32 mode
1793 IMRPhenomXHM_Ringdown_Amplitude_Veto(&V2, &V3, V4, pWFHM, pWF22); // If satisfied, remove the 2 inter collocation points
1794 #if DEBUG == 1
1795 printf("VETO: Ringdown Amp Veto\n");
1796 printf("V2 = %.16f\n",V2);
1797 printf("V3 = %.16f\n",V3);
1798 #endif
1799 }
1800
1801 if(pWFHM->IMRPhenomXHMIntermediateAmpVersion==105 && pAmp->WavyInt==1){
1802 if(WavyPoints(V2,V3,V4)==1){
1803 V3 = V2;
1804 F3 = F2;
1805 V2 = 1.;
1807 #if DEBUG == 1
1808 printf("VETO: Wavy Inter colloc points\n");
1809 printf("V2 = %.16f\n",V2);
1810 printf("V3 = %.16f\n",V3);
1811 #endif
1812 }
1813 }
1814
1815 if((pWF22->q>40. && pWF22->chi1L>0.9 && V2!=1 && V3!=1) ||
1816 (pWFHM->modeTag==32 && (pWFHM->IMRPhenomXHMIntermediateAmpVersion != 101) && ((pWF22->q>2.5 && pWF22->chi1L<-0.6 && pWF22->chi2L>0) || (pWF22->chi1L<-0.9&&pWF22->chi2L<-0.9))) ||
1817 (pWFHM->modeTag==21 && pWF22->eta<0.23 && pWF22->chi1L>0.7 && pWF22->chi2L<-0.5)){
1818 V2 = 1.;
1819 V3 = 1.;
1821 #if DEBUG == 1
1822 printf("VETO: veto regions\n");
1823 printf("V2 = %.16f\n",V2);
1824 printf("V3 = %.16f\n",V3);
1825 #endif
1826 }
1827 /*** End of vetos **/
1828 }
1829
1830 // The reconstruction function (Update_Intermediate_Amplitude_Coefficients) assumes that F3 is the point with value.
1831 // If F3 was removed in the veto the we replace it with F2
1832 if(V3 == 1.){
1833 V3 = V2;
1834 F3 = F2;
1835 V2 = 1.;
1836 }
1837
1838 /*
1839 Reconstruct the phenomenological coefficients for the intermediate ansatz
1840 */
1841
1842 // Store the values for the reconstruction.
1843 pAmp->v1 = V1;
1844 pAmp->v2 = V2;
1845 pAmp->v3 = V3;
1846 pAmp->v4 = V4;
1847 pAmp->f1 = F1;
1848 pAmp->f2 = F2;
1849 pAmp->f3 = F3;
1850 pAmp->f4 = F4;
1851 pAmp->d1 = d1;
1852 pAmp->d4 = d4;
1853
1854 #if DEBUG == 1
1855 printf("\nIntermediate Amplitude Input \r\n");
1856 printf("\nIMRPhenomXHMIntermediateAmpVersion = %i \r\n", pWFHM->IMRPhenomXHMIntermediateAmpVersion);
1857 printf("V1 = %.16f\n", V1);
1858 printf("V2 = %.16f\n", V2);
1859 printf("V3 = %.16f\n", V3);
1860 printf("V4 = %.16f\n", V4);
1861 printf("F1 = %.16f\n", F1);
1862 printf("F2 = %.16f\n", F2);
1863 printf("F3 = %.16f\n", F3);
1864 printf("F4 = %.16f\n", F4);
1865 printf("d1 = %.16f\n", d1);
1866 printf("d4 = %.16f\n", d4);
1867 printf("fAmpMatchIn = %.16f\n", pAmp->fAmpMatchIN);
1868 #endif
1869
1870 // Compute the coefficients of the polynomial
1872
1873
1874 #if DEBUG == 1
1875 printf("\nIntermediate polynomial coeffcients Before ChoosePolOrder \r\n");
1876 printf("\nIMRPhenomXHMIntermediateAmpVersion = %i \r\n", pWFHM->IMRPhenomXHMIntermediateAmpVersion);
1877 printf("delta0 = %.16f\n", pAmp->delta0);
1878 printf("delta1 = %.16f\n", pAmp->delta1);
1879 printf("delta2 = %.16f\n", pAmp->delta2);
1880 printf("delta3 = %.16f\n", pAmp->delta3);
1881 printf("delta4 = %.16f\n", pAmp->delta4);
1882 printf("delta5 = %.16f\n", pAmp->delta5);
1883 printf("fAmpMatchIn = %.16f\n", pAmp->fAmpMatchIN);
1884 #endif
1885
1886 // Check that the polynomial does not cross zero, because the actual reconstructing function is the inverse of this polynomial
1887 // If it crosses zero, then remove one collocation and lower the order of the polynomial.
1888 ChoosePolOrder(pWFHM, pAmp);
1889
1890 #if DEBUG == 1
1891 printf("\nIMRPhenomXHMIntermediateAmpVersion = %i \r\n", pWFHM->IMRPhenomXHMIntermediateAmpVersion);
1892 printf("\nIntermediate polynomial coeffcients After ChoosePolOrder\r\n");
1893 printf("delta0 = %.16f\n", pAmp->delta0);
1894 printf("delta1 = %.16f\n", pAmp->delta1);
1895 printf("delta2 = %.16f\n", pAmp->delta2);
1896 printf("delta3 = %.16f\n", pAmp->delta3);
1897 printf("delta4 = %.16f\n", pAmp->delta4);
1898 printf("delta5 = %.16f\n", pAmp->delta5);
1899 printf("fAmpMatchIn = %.16f\n", pAmp->fAmpMatchIN);
1900 #endif
1901 }// END of 122018 version
1902
1903 }
1904
1906 double factor = 0.;
1907 switch(rescalefactor){
1908 case 0:{ // Strain
1909 factor = 1.;
1910 break;
1911 }
1912 case 1:{ // 22factor
1913 factor = pAmp->ampNorm * powers_of_Mf->m_seven_sixths;
1914 break;
1915 }
1916 case 2:{ // lmfactor
1917 if (pAmp->PNdominantlmpower == 1) factor = pAmp->PNdominant * powers_of_Mf->m_seven_sixths * pAmp->PNdominantlm * powers_of_Mf->one_third;
1918 if (pAmp->PNdominantlmpower == 2) factor = pAmp->PNdominant * powers_of_Mf->m_seven_sixths * pAmp->PNdominantlm * powers_of_Mf->two_thirds;
1919 if (pAmp->PNdominantlmpower == 3) factor = pAmp->PNdominant * powers_of_Mf->m_seven_sixths * pAmp->PNdominantlm * powers_of_Mf->itself;
1920 if (pAmp->PNdominantlmpower == 4) factor = pAmp->PNdominant * powers_of_Mf->m_seven_sixths * pAmp->PNdominantlm * powers_of_Mf->four_thirds;
1921 break;
1922 }
1923 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in RescaleFactor: version %i is not valid. Recommended version is 1.", rescalefactor);}
1924 }
1925 return factor;
1926}
1927
1928
1929/*******************************************/
1930/* */
1931/* COMPUTE PHASE COEFFICIENTS */
1932/* */
1933/*******************************************/
1934
1935/* This function is an expansion in the frequency of the real part of the 21-amplitude, following the convention of Eq.(2.2) in the paper */
1936/* This function occasionally changes sign, so we need to keep track of this to correct the WF's phase, since the amplitude of IMRPhenomXHM is positive by construction */
1938
1939 double eta,chi1,chi2;
1940 eta =wf22->eta;
1941 chi1=wf22->chi1L;
1942 chi2=wf22->chi2L;
1943 double delta=sqrt(1 - 4*eta);
1944
1945 double output=(-16*delta*eta*ff*pow(LAL_PI,1.5))/(3.*sqrt(5)) + (4*pow(2,0.3333333333333333)*(chi1 - chi2 + delta *(chi1+ chi2))*eta*pow(ff,1.3333333333333333)*pow(LAL_PI,1.8333333333333333))/sqrt(5) + (2*pow(2,0.6666666666666666)*eta*(306*delta - 360*delta*eta)*pow(ff,1.6666666666666667)*pow(LAL_PI,2.1666666666666665))/(189.*sqrt(5));
1946
1947 if(output>=0) return(1);
1948 else return(-1);
1949 }
1950
1951/* compute phase coefficients */
1953
1954 int ell=pWFHM->ell;
1955 int emm=pWFHM->emm;
1956
1957 /* Pre-initialize all phenomenological coefficients */
1958
1959 //inspiral
1960 // these will be computed by rescaling IMRPhenomX inspiral coefficients
1961 for(int i=0; i<NMAX_INSPIRAL_COEFFICIENTS; i++)
1962 {
1963 pPhase->phi[i]=0.;
1964 pPhase->phiL[i]=0.;
1965 }
1966
1967 //intermediate
1968 // will be determined by solving linear system at collocation points (see Eqs. (5.5)-(5.6) in the paper)
1969 pPhase->c0 = 0.0;
1970 pPhase->c1 = 0.0;
1971 pPhase->c2 = 0.0;
1972 pPhase->c3 = 0.0;
1973 pPhase->c4 = 0.0;
1974 pPhase->cL = 0.0;
1975 // ringdown spherical
1976 // these will be obtained by rescaling the 22 ringdown parameters (see Eq. (6.5) in the paper)
1977 pPhase->alpha0 = 0.0;
1978 pPhase->alpha2 = 0.0;
1979 pPhase->alphaL = 0.0;
1980
1981
1982
1983 // set number of collocation points used (depends on the mode)
1984 int nCollocationPts_inter=pWFHM->nCollocPtsInterPhase;
1985 // define mass-ratio that discriminates between EMR and rest of cases
1986 double eta_m1=1./pWF22->eta;
1987 // store mode tag
1988 int modeint=pWFHM->modeInt;
1989
1990 IMRPhenomX_UsefulPowers powers_of_f;
1991
1992 //initialize frequencies of colloc points in intermediate region
1993 IMRPhenomXHM_Intermediate_CollocPtsFreqs(pPhase,pWFHM,pWF22);
1994
1995
1996 //for each collocation point, call fit giving the value of the phase derivative at that point
1997 for(int i = 0; i<N_MAX_COEFFICIENTS_PHASE_INTER; i++)
1998 {
2000 // time-shift wf so that modes peak around t=0
2001 // our hybrids are built so that their inverse FT peaks ~500 M before the end of the time-interval. Our reconstructed WFs will have the same normalization, so we shift the phase here so that the modes peak around t=0
2003 }
2004
2005 double fcutRD=pPhase->fPhaseMatchIM;
2006 double fcutInsp=pPhase->fPhaseMatchIN;
2007
2008
2009
2010 /************** Inspiral: rescale PhenX and apply PN corrections *********/
2011
2012 //collect all the phenX inspiral coefficients
2013 double phenXnonLog[]={pPhase22->phi0,pPhase22->phi1,pPhase22->phi2,pPhase22->phi3,pPhase22->phi4,pPhase22->phi5,pPhase22->phi6,pPhase22->phi7,pPhase22->phi8,pPhase22->phi9,0.,0,0,0};
2014 double phenXLog[]={0.,0.,0.,0.,0.,pPhase22->phi5L,pPhase22->phi6L,0.,pPhase22->phi8L,pPhase22->phi9L,0,0,0,0};
2015 double pseudoPN[]={0.,0.,0.,0.,0.,0.,0.,0.,pPhase22->sigma1,pPhase22->sigma2,pPhase22->sigma3,pPhase22->sigma4,pPhase22->sigma5};
2016
2017
2018
2019 //rescale the coefficients of phenX by applying phi_lm(f)~m/2 *phi_22(2/m f)
2020 // for more details, see Appendix D of the paper
2021 double m_over_2=emm*0.5, two_over_m=1./m_over_2;
2022
2023 double fact=pPhase22->phiNorm/pWF22->eta;
2024
2025 for(int i=0; i< NMAX_INSPIRAL_COEFFICIENTS; i++){
2026
2027 phenXnonLog[i]=phenXnonLog[i]*fact;
2028 phenXLog[i]=phenXLog[i]*fact;
2029 pseudoPN[i]=pseudoPN[i]*fact;
2030
2031 // scaling the logarithmic terms introduces some extra contributions in the non-log terms
2032 pPhase->phi[i]=(phenXnonLog[i]+pseudoPN[i]-phenXLog[i]*log(m_over_2))*pow(m_over_2,(8-i)/3.);
2033 pPhase->phiL[i]=phenXLog[i]*pow(m_over_2,(8-i)/3.);
2034 }
2035
2036
2037 // if the mass-ratio is not extreme, use the complex PN amplitudes to correct the orbital phase at linear order in f
2038 // see Eq. (4.12) in the paper
2039 if(pWF22->eta>0.01){
2040 pPhase->LambdaPN=IMRPhenomXHM_Insp_Phase_LambdaPN(pWF22->eta, pWFHM->modeTag);
2041 #if DEBUG == 1
2042 printf("Add phase shift from PN complex amplitude\n");
2043 #endif
2044
2045 }
2046 // else use some phenomenological fits: the fits give the coefficient of a linear-in-f term to be added to the orbital phase
2047 else{
2048
2049 pPhase->LambdaPN= pPhase->InspiralPhaseFits[pWFHM->modeInt](pWF22,pWFHM->IMRPhenomXHMInspiralPhaseVersion);
2050 #if DEBUG == 1
2051 printf("\nLambdaPN = %.16f\n", pPhase->LambdaPN);
2052 #endif
2053 }
2054
2055 pPhase->phi[8]= pPhase->phi[8]+ pPhase->LambdaPN;
2056
2057
2058 /****************************************
2059 *********** Intermediate-ringdown region ********
2060 ****************************************/
2061
2062 // below we set up the linear system to be solved in the intermediate region
2063
2064
2065 /* GSL objects for solving system of equations via LU decomposition */
2066 gsl_vector *b, *x;
2067 gsl_matrix *A;
2068 gsl_permutation *p;
2069 int s;
2070
2071 p = gsl_permutation_alloc(nCollocationPts_inter);
2072 b = gsl_vector_alloc(nCollocationPts_inter);
2073 x = gsl_vector_alloc(nCollocationPts_inter);
2074 A = gsl_matrix_alloc(nCollocationPts_inter,nCollocationPts_inter);
2075
2076
2077 // for high-spin cases: avoid sharp transitions for 21 mode by extending the inspiral region into the intermediate one
2078
2079 if(pWFHM->modeTag==21&&pWF22->STotR>=0.8){
2080
2081 double insp_vals[3], FF, diff12, diff23;
2082 for(int i=0; i<3; i++){
2083 FF=two_over_m*pPhase->CollocationPointsFreqsPhaseInter[i];
2084 IMRPhenomX_Initialize_Powers(&powers_of_f,FF);
2085 insp_vals[i]=1./pWF22->eta*IMRPhenomX_dPhase_22(FF,&powers_of_f,pPhase22,pWF22);
2086 }
2087
2088 diff12=insp_vals[0]-insp_vals[1];
2089 diff23=insp_vals[1]-insp_vals[2];
2090
2093
2094 }
2095
2096
2097 // choose collocation points according to spin/mass ratio
2098 // current catalogue of simulations include some cases that create unphysical effects in the fits -> we need to use different subset of collocation points according to the parameters (we have to pick 5 out of 6 available fits)
2099 /* cpoints_indices is an array of integers labelling the collocation points chosen in each case, e.g.
2100 cpoints_indices={0,1,3,4,5} would mean that we are discarding the 3rd collocation points in the reconstructio */
2101
2102 int cpoints_indices[nCollocationPts_inter];
2103 cpoints_indices[0]=0;
2104 cpoints_indices[1]=1;
2105 cpoints_indices[4]=5;
2106
2107
2108 if((pWF22->eta<pWFHM->etaEMR)||(emm==ell&&pWF22->STotR>=0.8)||(pWFHM->modeTag==33&&pWF22->STotR<0))
2109 {
2110 cpoints_indices[2]=3;
2111 cpoints_indices[3]=4;
2112 }
2113 else if(pWF22->STotR>=0.8&&pWFHM->modeTag==21){
2114
2115 cpoints_indices[2]=2;
2116 cpoints_indices[3]=4;
2117 }
2118
2119 else{
2120 cpoints_indices[2]=2;
2121 cpoints_indices[3]=3;
2122 }
2123
2124
2125
2126
2127 switch(pWFHM->MixingOn){
2128 // mode-mixing off: mode does not show significant mode-mixing
2129 // the MixingOn flag is automatically switched off for the modes 21,33,44
2130 case 0:
2131 {
2132
2133 /************ Ringdown **************/
2134
2135 /* ansatz:
2136 alpha0 + ((fRDlm^2) alpha2)/(f^2) + alphaL*(fdamplm)/((fdamplm)^2 + (f - fRDlm)^2)*/
2137
2138 //compute alpha2 by rescaling the coefficient of the 22-reconstruction
2139 double wlm;
2140 if(pWFHM->ell==pWFHM->emm) wlm=2;
2141 else wlm=pWFHM->emm/3.;
2142
2143 pPhase->alpha2=1./pow(pWFHM->fRING,2)*wlm*IMRPhenomXHM_RD_Phase_22_alpha2(pWF22,pWFHM->IMRPhenomXHMRingdownPhaseVersion);
2144
2145
2146 // compute alphaL
2148
2149 #if DEBUG == 1
2150 printf("alpha2_fit=%.16e\n",IMRPhenomXHM_RD_Phase_22_alpha2(pWF22,pWFHM->IMRPhenomXHMRingdownPhaseVersion));
2151 printf("Ringdown parameters: \n alpha2=%.16e \n alphaL=%.16e\n",pPhase->alpha2,pPhase->alphaL);
2152 # endif
2153
2154 // compute spherical-harmonic phase and its first derivative at matching frequency --used to glue RD and intermediate regions--
2155 IMRPhenomX_Initialize_Powers(&powers_of_f,fcutRD);
2156 pPhase->phi0RD=IMRPhenomXHM_RD_Phase_AnsatzInt(fcutRD, &powers_of_f,pWFHM, pPhase);
2157 pPhase->dphi0RD=IMRPhenomXHM_RD_Phase_Ansatz(fcutRD, &powers_of_f,pWFHM, pPhase);
2158
2159
2160 /************ Intermediate **************/
2161
2162
2163 // set up the linear system by calling the fits of the intermediate-region collocation points
2164 for(int i=0; i<nCollocationPts_inter; i++){
2165
2166 int ind=cpoints_indices[i];
2167 gsl_vector_set(b,i,pPhase->CollocationPointsValuesPhaseInter[ind]);
2168 REAL8 ff=pPhase->CollocationPointsFreqsPhaseInter[ind], ffm1=1./ff, ffm2=ffm1*ffm1;
2169 REAL8 fpowers[]={1., (pWFHM->fDAMP)/(pow(pWFHM->fDAMP,2)+pow(ff-(pWFHM->fRING),2)),ffm1,ffm2,ffm2*ffm2, ffm1*ffm2};
2170 for(int j=0; j<nCollocationPts_inter; j++)
2171 gsl_matrix_set(A,i,j,fpowers[j]);
2172 }
2173
2174 break;
2175 }
2176
2177 // mode-mixing on (the MixingOn flag is automatically switched on when reconstructing the 32 mode)
2178 case 1:
2179 {
2180
2181 // for the 32 mode, the ringdown waveform is computed outside of the dedicated amplitude and phase routines, by calling SpheroidalToSpherical
2182 // compute dphi/df and d2phi/df2 at fcutRD starting from rotation of spheroidal ansatz: these values will then enter the linear system that determines the coefficients in the intermediate region
2183
2184 // we compute derivatives by applying finite-difference schemes. In total we will need the value of the phase at three points around fcutRD
2185
2186 // we first compute the full spherical-harmonic-basis waveforms at three points
2187 double complex SphericalWF[3];
2188 double fstep=0.0000001;
2189
2190 for(int i=0; i<3; i++){
2191
2192 double FF=fcutRD+(i-1)*fstep;
2193 IMRPhenomX_UsefulPowers powers_of_FF;
2194 IMRPhenomX_Initialize_Powers(&powers_of_FF,FF);
2195 SphericalWF[i]=SpheroidalToSpherical(&powers_of_FF, pAmp22, pPhase22, pAmp, pPhase, pWFHM, pWF22);
2196 }
2197
2198 long double phase_args[]={fmodl(carg(SphericalWF[0]),2.*LAL_PI),fmodl(carg(SphericalWF[1]),2.*LAL_PI),fmodl(carg(SphericalWF[2]),2.*LAL_PI)};
2199
2200 // make sure that all the three points belong to the same branch of mod
2201 for(int i=0; i<3; i++){
2202
2203 if(phase_args[i]>0)
2204 phase_args[i]-=2.*LAL_PI;
2205
2206 }
2207
2208 // store spherical-harmonic phase at fcutRD
2209 pPhase->phi0RD=phase_args[1];
2210 long double fstep_m1= 1./fstep;
2211 // we apply the FD schemes to get first and second phase derivatives at fcutRD
2212 pPhase->dphi0RD=0.5*fstep_m1*(phase_args[2]-phase_args[0]);
2213 double d2phi0RD=fstep_m1*fstep_m1*(phase_args[2]-2.*phase_args[1]+phase_args[0]);
2214
2215 #if DEBUG == 1
2216 printf("dphi0ref=%.16e\t d2phi0ref=%.16e\n",pPhase->dphi0RD,d2phi0RD);
2217 # endif
2218
2219 // To achieve a smoother transition with the intermediate region (IR), we feed into the intermediate-region reconstruction the first and second derivative of the reconstructed ringdown phase, see Eq. (5.6)
2220
2221
2222 // first derivative
2223 // if the mass-ratio is not extreme, we use the derivative computed using the FD scheme above, else we keep the original value of the fit
2224 // in the second case, we will therefore need to glue intermediate and ringdown region later
2225 if(pWF22->eta>pWFHM->etaEMR){
2226 pPhase->CollocationPointsFreqsPhaseInter[nCollocationPts_inter-2]=fcutRD;
2227 pPhase->CollocationPointsValuesPhaseInter[nCollocationPts_inter-2]=pPhase->dphi0RD;
2228 }
2229 // second derivative
2230 pPhase->CollocationPointsFreqsPhaseInter[nCollocationPts_inter-1]=fcutRD;
2231 pPhase->CollocationPointsValuesPhaseInter[nCollocationPts_inter-1]=d2phi0RD;
2232
2233
2234 // set up the linear system to determine the intermediate region coefficients
2235 for(int i=0; i<nCollocationPts_inter; i++){
2236
2237 gsl_vector_set(b,i,pPhase->CollocationPointsValuesPhaseInter[i]);
2238 REAL8 ff=pPhase->CollocationPointsFreqsPhaseInter[i], ffm1=1./ff, ffm2=ffm1*ffm1;
2239 REAL8 fpowers[]={1., (pWFHM->fDAMP)/(pow(pWFHM->fDAMP,2)+pow(ff-(pWFHM->fRING),2)),ffm1,ffm2,ffm2*ffm2, ffm1*ffm2};
2240 for(int j=0; j<nCollocationPts_inter; j++)
2241 gsl_matrix_set(A,i,j,fpowers[j]);
2242
2243 }
2244
2245 // set the last point to equal the second derivative of the rotated spheroidal ansatz
2246 int cpoint_ind=nCollocationPts_inter-1;
2247 gsl_vector_set(b,cpoint_ind,pPhase->CollocationPointsValuesPhaseInter[cpoint_ind]);
2248 REAL8 ff=pPhase->CollocationPointsFreqsPhaseInter[cpoint_ind];
2249 REAL8 ffm1=1./ff, ffm2=ffm1*ffm1, ffm3=ffm2*ffm1, ffm4=ffm2*ffm2, ffm5=ffm3*ffm2;
2250
2251 REAL8 fpowers[]={0.,-2*(pWFHM->fDAMP)*(ff-(pWFHM->fRING))/pow((pow(pWFHM->fDAMP,2)+pow(ff-(pWFHM->fRING),2)),2),-ffm2, -2.* ffm3, -4.*ffm5,-3* ffm4};
2252 for(int j=0; j<nCollocationPts_inter; j++){
2253 gsl_matrix_set(A,cpoint_ind,j,fpowers[j]);
2254 }
2255
2256
2257 #if DEBUG == 1
2258 printf("Collocation points for intermediate ansatz %i:\n",nCollocationPts_inter);
2259 for(int i=0;i<nCollocationPts_inter; i++)
2260 printf("p%d={%.13f,%.13f};\n",i,pPhase->CollocationPointsFreqsPhaseInter[i],pPhase->CollocationPointsValuesPhaseInter[i]);
2261 #endif
2262
2263 break;
2264
2265
2266
2267 }
2268
2269 default:{XLALPrintError("Error in IMRPhenomXHM_GetPhaseCoefficients:mode-mixing not properly initialized.\n");}
2270
2271 }
2272
2273
2274 // at this point we have initizialized all the collocation points in the intermediate region, so we can solve to get the coefficients of the ansatz
2275
2276 /* In the intermediate region , the most general ansatz is
2277 (c0 + c1 /f + c2 /(f)^2 + c4 /(f)^4 + c3 /f^3 +cL fdamp/((fdamp)^2 + (f - fRD )^2))*/
2278
2279 /* We now solve the system A x = b via an LU decomposition */
2280 gsl_linalg_LU_decomp(A,p,&s);
2281 gsl_linalg_LU_solve(A,p,b,x);
2282
2283 pPhase->c0 = gsl_vector_get(x,0); // x[0]; // c0
2284 pPhase->cL = gsl_vector_get(x,1); // x[1] // cL
2285 pPhase->c1 = gsl_vector_get(x,2); // x[2]; // c1
2286 pPhase->c2 = gsl_vector_get(x,3); // x[3]; // c2
2287 pPhase->c4 = gsl_vector_get(x,4); // x[4]; // c4
2288
2289 /* (IMRPhenomXPNRUseTunedCoprec) Add PNR or input deviations (e.g. tuning for coprecessing frame model) to cL, c1 and c4 in the same manner as is done for the l=m=2 moment */
2290 if ( pWFHM->modeTag == 33){
2291 pPhase->c0 = pPhase->c0 + ( pWFHM->PNR_DEV_PARAMETER * pWFHM->NU0 );
2292 pPhase->cL = pPhase->cL + ( pWFHM->PNR_DEV_PARAMETER * pWFHM->NU4 );
2293 pPhase->c1 = pPhase->c1 + ( pWFHM->PNR_DEV_PARAMETER * pWFHM->ZETA2 );
2294 pPhase->c4 = pPhase->c4 + ( pWFHM->PNR_DEV_PARAMETER * pWFHM->ZETA1 );
2295 }
2296
2297
2298 // currently the 32 mode is calibrated using one extra point
2299 if ((pWFHM->modeTag)== 32)
2300 { pPhase->c3 = gsl_vector_get(x,5); //c3
2301
2302 // if the mass-ratio is extreme, we need to glue intermediate and ringdown region, as explained above
2303 if(pWF22->eta<pWFHM->etaEMR)
2304 {
2305 IMRPhenomX_Initialize_Powers(&powers_of_f,fcutRD);
2306 pPhase->c0=pPhase->c0+pPhase->dphi0RD-IMRPhenomXHM_Inter_Phase_Ansatz(fcutRD, &powers_of_f,pWFHM, pPhase);
2307
2308 }
2309 }
2310
2311 #if DEBUG == 1
2312 printf("Intermediate region coefficients:\n");
2313 printf("c0=%.16e\n c1=%.16e\n c2=%.16e\n c3=%.16e\n c4=%.16e\n cL=%.16e\n", pPhase->c0,pPhase->c1,pPhase->c2,pPhase->c3,pPhase->c4,pPhase->cL);
2314 #endif
2315
2316 gsl_vector_free(b);
2317 gsl_vector_free(x);
2318 gsl_matrix_free(A);
2319 gsl_permutation_free(p);
2320
2321
2322 /****************** end of ringdown and inspiral reconstruction **************/
2323
2324 //glue inspiral and intermediate
2325 IMRPhenomX_Initialize_Powers(&powers_of_f,fcutInsp);
2326 pPhase->C1INSP=IMRPhenomXHM_Inter_Phase_Ansatz(fcutInsp, &powers_of_f,pWFHM, pPhase)-IMRPhenomXHM_Inspiral_Phase_Ansatz(fcutInsp, &powers_of_f, pPhase);
2327 pPhase->CINSP=-(pPhase->C1INSP*fcutInsp)+IMRPhenomXHM_Inter_Phase_AnsatzInt(fcutInsp, &powers_of_f,pWFHM, pPhase)-IMRPhenomXHM_Inspiral_Phase_AnsatzInt(fcutInsp, &powers_of_f, pPhase);
2328 //glue ringdown and intermediate regions
2329 IMRPhenomX_Initialize_Powers(&powers_of_f,fcutRD);
2330
2331 pPhase->C1RD=IMRPhenomXHM_Inter_Phase_Ansatz(fcutRD, &powers_of_f,pWFHM, pPhase)-pPhase->dphi0RD;
2332 pPhase->CRD=-pPhase->C1RD*fcutRD+IMRPhenomXHM_Inter_Phase_AnsatzInt(fcutRD, &powers_of_f,pWFHM, pPhase)-pPhase->phi0RD;
2333
2334
2335 // we now have a C1 reconstruction of the phase
2336 // below we align each mode so that at low-f its relative phase wrt the 22 agrees with PN
2337 double falign;
2338 //somehow arbitary cutoff to pick frequency for alignment: must be in the inspiral region
2339 if(pWF22->eta>pWFHM->etaEMR)
2340 falign=0.6*m_over_2*pWF22->fMECO;
2341 else
2342 falign=m_over_2*pWF22->fMECO;
2343 // printf("0>> falign = %f\n",falign);
2344
2345 // // (ltl) 0.6 is found to be too large a factor, so let's try a smaller one
2346 // if (pWF22->IMRPhenomXPNRUseTunedCoprec && (pWFHM->modeTag==33)) {
2347 // if(pWF22->eta>pWFHM->etaEMR){
2348 // falign = 0.06*m_over_2*pWF22->fMECO;
2349 // printf("&>> falign = %f\n",falign);
2350 // } else {
2351 // falign = m_over_2*pWF22->fMECO;
2352 // }
2353 // }
2354
2355 IMRPhenomX_UsefulPowers powers_of_falign;
2356 IMRPhenomX_Initialize_Powers(&powers_of_falign,falign);
2357 IMRPhenomX_Initialize_Powers(&powers_of_f,two_over_m*falign);
2358
2359 /* compute explicitly the phase normalization to be applied to IMRPhenomX, when mode-mixing is on this will have been already computed in GetSpheroidalCoefficients */
2360 if(pWFHM->MixingOn==0){
2361 IMRPhenomX_UsefulPowers powers_of_MfRef;
2362 IMRPhenomX_Initialize_Powers(&powers_of_MfRef,pWF22->MfRef);
2364 pWFHM->timeshift=IMRPhenomX_TimeShift_22(pPhase22, pWF22);
2365 pWFHM->phiref22 = -1./pWF22->eta*IMRPhenomX_Phase_22(pWF22->MfRef, &powers_of_MfRef, pPhase22, pWF22) - pWFHM->timeshift*pWF22->MfRef - pWFHM->phaseshift + 2.0*pWF22->phi0 + LAL_PI_4;
2366 }
2367
2368 // we determine the phase normalization of each mode by imposing Eq. (4.13), i.e. phi_lm(fref)-m/2 phi_22(2/m fref)~3.*LAL_PI_4*(1-m_over_2)
2369 double deltaphiLM = m_over_2*(1./pWF22->eta*IMRPhenomX_Phase_22(two_over_m*falign, &powers_of_f,pPhase22,pWF22)+pWFHM->phaseshift + pWFHM->phiref22)+pWFHM->timeshift*falign-3.*LAL_PI_4*(1-m_over_2)-(IMRPhenomXHM_Inspiral_Phase_AnsatzInt(falign, &powers_of_falign,pPhase)+pPhase->C1INSP*falign+pPhase->CINSP);
2370 pPhase->deltaphiLM=fmod(deltaphiLM,2.*LAL_PI);
2371
2372 // /* (ll) We have commented out the code above in order to rewrite the same code in a more carefully formatted way */
2373
2374 // // (ll) Define the relative phase parameter used in Eq. (4.13)
2375 // double relative_phase_at_zero_freq = 3.*LAL_PI_4*(1-m_over_2);
2376
2377 // // (ll) Define the (2,2) inspiral phase evaluated at (2/m)*falign, and then rescaled by m/2
2378 // double rescaled_phi22_at_falign = m_over_2*(1./pWF22->eta*IMRPhenomX_Phase_22(two_over_m*falign, &powers_of_f,pPhase22,pWF22)+pWFHM->phaseshift + pWFHM->phiref22)+pWFHM->timeshift*falign;
2379
2380 // //
2381 // printf("*>> falign = %f\n",falign);
2382 // printf("*>> relative_phase_at_zero_freq = %f\n",relative_phase_at_zero_freq);
2383 // printf("*>> m_over_2 = %f\n\n",m_over_2);
2384
2385 // // (ll) Define XHM inspiral phase at align frequency
2386 // double phiLM_inspiral_at_falign = (IMRPhenomXHM_Inspiral_Phase_AnsatzInt(falign, &powers_of_falign,pPhase)+pPhase->C1INSP*falign+pPhase->CINSP);
2387
2388 // // (ll)
2389 // double deltaphiLM = rescaled_phi22_at_falign - ( relative_phase_at_zero_freq + phiLM_inspiral_at_falign );
2390
2391 // // (ll) Get the representation closest to 2*pi
2392 // pPhase->deltaphiLM=fmod(deltaphiLM,2.*LAL_PI);
2393
2394 #if DEBUG == 1
2395 printf("\n****** Connection coefficients of 22 in %i******\n", pWFHM->modeTag);
2396 printf("%.6f %.6f %.6f %.6f %.6f\n", pPhase22->C1Int, pPhase22->C2Int, pPhase22->C1MRD, pPhase22->C1MRD, IMRPhenomX_Phase_22(two_over_m*falign, &powers_of_f,pPhase22,pWF22));
2397 #endif
2398
2399
2400 // for the 21, we need to make sure the sign of the PN amplitude is positive, else we'll need to flip its phase by Pi
2401 if ((pWFHM->modeTag)== 21)
2402 {
2403 int ampsign=IMRPhenomXHM_PN21AmpSign(0.008,pWF22);
2404 // the sign of the 21 amplitude changes across the parameter space, we add Pi to the phase to account for this - remember that the amplitude of the model is positive by construction
2405 if(ampsign>0) pPhase->deltaphiLM+=LAL_PI;
2406
2407 }
2408
2409 // all the coefficients have been now stored
2410 // print parameters to file if in debugging mode
2411 #if DEBUG == 1
2412 FILE *file;
2413 char fileSpec[40];
2414 sprintf(fileSpec, "PhaseParameters%i.dat", pWFHM->modeTag);
2415
2416 file = fopen(fileSpec,"w");
2417
2418 fprintf(file,"\n*** %i Mode ***\n", pWFHM->modeTag);
2419
2420 fprintf(file,"\n*** Intrinsic Parameters ***\n");
2421 fprintf(file,"eta = %.16f\n", pWF22->eta);
2422 fprintf(file,"chi1z = %.16f\n", pWF22->chi1L);
2423 fprintf(file,"chi2z = %.16f\n", pWF22->chi2L);
2424
2425 fprintf(file,"\n*** QNM Frequencies ***\n");
2426 fprintf(file,"fRINGlm = %.16f\n", pWFHM->fRING);
2427 fprintf(file,"fDAMPlm = %.16f\n", pWFHM->fDAMP);
2428
2429
2430 fprintf(file,"\n*** Phase Inspiral ***\n");
2431 fprintf(file,"Lambda = %.16f\n", pPhase->LambdaPN);
2432
2433 fprintf(file,"\n*** Phase Intermediate ***\n f_i\t dphi(f_i)\n");
2434
2435 for(int i = 0; i<N_MAX_COEFFICIENTS_PHASE_INTER; i++)
2437
2438
2439 if(pWFHM->modeTag==32){
2440 fprintf(file,"\n*** Mixing Coefficients ***\n");
2441 fprintf(file,"222 = %.16f %.16f *I\n", creal(pWFHM->mixingCoeffs[0]), cimag(pWFHM->mixingCoeffs[0]));
2442 fprintf(file,"223 = %.16f %.16f *I\n", creal(pWFHM->mixingCoeffs[1]), cimag(pWFHM->mixingCoeffs[1]));
2443 fprintf(file,"322 = %.16f %.16f *I\n", creal(pWFHM->mixingCoeffs[2]), cimag(pWFHM->mixingCoeffs[2]));
2444 fprintf(file,"323 = %.16f %.16f *I\n", creal(pWFHM->mixingCoeffs[3]), cimag(pWFHM->mixingCoeffs[3]));
2445
2446 fprintf(file, "\n*** Ringdown Phase Spheroidal ***\n f_i \t dphi(f_i)\n");
2447 for(int i=0; i<pWFHM->nCollocPtsRDPhase; i++)
2448 fprintf(file, "%.16f \t %.16f \n", pPhase->CollocationPointsFreqsPhaseRD[i],pPhase->RingdownPhaseFits[i](pWF22,pWFHM->IMRPhenomXHMRingdownPhaseVersion));
2450
2451 }
2452
2453 else{
2454 fprintf(file, "\n*** Ringdown phase ***\n");
2455 fprintf(file,"alpha2=%.16f\n alphaL=%.16f",pPhase->alpha2,pPhase->alphaL);
2456 }
2457
2458 fclose(file);
2459 #endif
2460
2461 }
2462
2463
2464
2465/* Reconstruct the coefficients of the ringdown phase in a spheroidal-harmonics basis --only for the 32 mode */
2467
2468
2469 int nCollocationPts_RD_Phase = pWFHM->nCollocPtsRDPhase;
2470 double CollocValuesPhaseRingdown[nCollocationPts_RD_Phase];
2471 double CollocFreqsPhaseRingdown[nCollocationPts_RD_Phase];
2472
2473 #if DEBUG == 1
2474 printf("Declare Colloc Values and Freqs:\n");
2475 #endif
2476 // define gsl variables
2477 gsl_vector *b, *x;
2478 gsl_matrix *A;
2479 gsl_permutation *p;
2480 int s;
2481
2482 p = gsl_permutation_alloc(nCollocationPts_RD_Phase);
2483 b = gsl_vector_alloc(nCollocationPts_RD_Phase);
2484 x = gsl_vector_alloc(nCollocationPts_RD_Phase);
2485 A = gsl_matrix_alloc(nCollocationPts_RD_Phase,nCollocationPts_RD_Phase);
2486
2487 IMRPhenomXHM_Ringdown_CollocPtsFreqs(pPhase,pWFHM, pWF22);
2488
2489 #if DEBUG == 1
2490 printf("Initialize RD CollocPtsFreqs:\n");
2491 #endif
2492
2493 if(pWFHM->IMRPhenomXHMRingdownPhaseVersion == 122019){
2494 // first fill-in the collocation points for the phase
2495 for(int i=0; i<nCollocationPts_RD_Phase; i++)
2496 {
2497
2498 CollocValuesPhaseRingdown[i] =pPhase->RingdownPhaseFits[i](pWF22,pWFHM->IMRPhenomXHMRingdownPhaseFitsVersion);
2499 CollocFreqsPhaseRingdown[i] =pPhase->CollocationPointsFreqsPhaseRD[i];
2500 gsl_vector_set(b,i,CollocValuesPhaseRingdown[i]);
2501 REAL8 ff=CollocFreqsPhaseRingdown[i], ffm1=1./ff, ffm2=ffm1*ffm1;
2502 REAL8 fpowers[]={1., (pWFHM->fDAMP)/(pow(pWFHM->fDAMP,2)+pow(ff-(pWFHM->fRING),2)),ffm2,ffm2*ffm2};
2503 for(int j=0; j<nCollocationPts_RD_Phase; j++)
2504 gsl_matrix_set(A,i,j,fpowers[j]);
2505
2506 }
2507 }
2508 else{
2509 for(int i = 0; i < nCollocationPts_RD_Phase; i++){
2510 CollocValuesPhaseRingdown[i] = pPhase->RingdownPhaseFits[i](pWF22,pWFHM->IMRPhenomXHMRingdownPhaseFitsVersion);
2511 CollocFreqsPhaseRingdown[i] = pPhase->CollocationPointsFreqsPhaseRD[i];
2512 gsl_vector_set(b,i,CollocValuesPhaseRingdown[i]);
2513 REAL8 ff = CollocFreqsPhaseRingdown[i], ffm1 = 1./ff, ffm2 = ffm1 * ffm1;
2514 REAL8 lorentzian = pWFHM->fDAMP / (pWFHM->fDAMP * pWFHM->fDAMP + (ff - pWFHM->fRING) * (ff - pWFHM->fRING));
2515 REAL8 fpowers[] = {1., ffm1, ffm2, ffm2 * ffm2, lorentzian};
2516 for(int j = 0; j < nCollocationPts_RD_Phase; j++){
2517 gsl_matrix_set(A, i, j, fpowers[j]); // Ansatz spheroidal
2518 }
2519 }
2520 }
2521
2522 #if DEBUG == 1
2523 printf("Collocation points in ringdown region:\n");
2524 for(int i=0; i<nCollocationPts_RD_Phase; i++){
2525 printf("p%d=(%.16e,%.16e)\n",i+1,CollocFreqsPhaseRingdown[i],CollocValuesPhaseRingdown[i]);
2526 }
2527 #endif
2528
2529 // solve the linear system of Eq. (6.8)
2530 gsl_linalg_LU_decomp(A,p,&s);
2531 gsl_linalg_LU_solve(A,p,b,x);
2532
2533 if(pWFHM->IMRPhenomXHMRingdownPhaseVersion == 122019){
2534 /* ansatz: alpha0 + (alpha2)/(f^2)+ (alpha4)/(f^4) + alphaL*(fdamplm)/((fdamplm)^2 + (f - fRDlm)^2)*/
2535 pPhase->alpha0_S = gsl_vector_get(x,0);
2536 pPhase->alphaL_S = gsl_vector_get(x,1);
2537 pPhase->alpha2_S = gsl_vector_get(x,2);
2538 pPhase->alpha4_S = gsl_vector_get(x,3);
2539 }
2540 else{
2541 //printf("New Spheroidal Version. Getting coefficients\n");
2542 for(int i=0; i<nCollocationPts_RD_Phase; i++){
2543 pPhase->RDCoefficient[i] = gsl_vector_get(x, i);
2544 }
2545 }
2546
2547 #if DEBUG == 1
2548 printf("**********\n alpha0_S=%.16e \n alphaL_S=%.16e \n alpha2_S=%.16e \n alpha4_S=%.16e \n \n ", pPhase->alpha0_S, pPhase->alphaL_S, pPhase->alpha2_S, pPhase->alpha4_S);
2549 #endif
2550
2551 gsl_vector_free(b);
2552 gsl_vector_free(x);
2553 gsl_matrix_free(A);
2554 gsl_permutation_free(p);
2555
2556 // we have reconstructed the "shape" of the spheroidal ringdown phase derivative, dphiS, but we need to adjust the relative time and phase shift of the final phiS wrt IMRPhenomX
2557
2558 IMRPhenomX_UsefulPowers powers_of_FREF;
2559 /*************** time-shift of spheroidal ansatz *******************/
2560 double frefRD=pWF22->fRING+pWF22->fDAMP;
2561 IMRPhenomX_Initialize_Powers(&powers_of_FREF,frefRD);
2562 // here we call a fit for dphiS(fref)-dphi22(fref)
2564
2565 // we compute dphi22(fref)
2567 pWFHM->timeshift=IMRPhenomX_TimeShift_22(pPhase22, pWF22);
2568
2569 pPhase->phi0_S = 0;
2570
2571 // we impose that dphiS(fref)-dphi22(fref) has the value given by our fit
2572 if(pWFHM->IMRPhenomXHMRingdownPhaseVersion == 122019){
2573 double dphi22ref = 1./pWF22->eta*IMRPhenomX_dPhase_22(frefRD, &powers_of_FREF,pPhase22,pWF22)+pWFHM->timeshift;
2574 pPhase->alpha0_S = pPhase->alpha0_S +dphi22ref+tshift-IMRPhenomXHM_RD_Phase_Ansatz(frefRD,&powers_of_FREF,pWFHM,pPhase);
2575 }
2576 else{
2577 // REAL8 psi4tostrain=XLALSimIMRPhenomXPsi4ToStrain(pWF22->eta, pWF22->STotR, pWF22->dchi);
2578 // pPhase->RDCoefficient[0] -= 2.*LAL_PI*(500+psi4tostrain);
2579 double dphi22ref = 1./pWF22->eta*IMRPhenomX_dPhase_22(frefRD, &powers_of_FREF,pPhase22,pWF22)+pWFHM->timeshift;
2580 double dphi32ref = IMRPhenomXHM_RD_Phase_Ansatz(frefRD,&powers_of_FREF,pWFHM,pPhase);
2581 pPhase->RDCoefficient[0] += dphi22ref+tshift-dphi32ref;
2582 //printf("dphi32ref before time shift = %.16e\n", dphi32ref);
2583 // printf("fittimeshift = %.16e\n", tshift);
2584 // printf("Time shift = %.16e\n", dphi22ref+tshift-dphi32ref);
2585 // printf("dphi22ref = %.16e\n", dphi22ref);
2586 // printf("dphi32ref = %.16e\n", IMRPhenomXHM_RD_Phase_Ansatz(frefRD,&powers_of_FREF,pWFHM,pPhase));
2587 //pPhase->RDCoefficient[0] += (dphi22ref + tshift - IMRPhenomXHM_RD_Phase_Ansatz(frefRD, &powers_of_FREF, pWFHM, pPhase));
2588 //pPhase->RDCoefficient[5] -= 2.*LAL_PI*(500+psi4tostrain);
2589 INT4 exponent = 5; // a + b / f^5
2590 REAL8 ff = pWFHM->fRING + 2 * pWFHM->fDAMP;
2591 IMRPhenomX_UsefulPowers powers_of_ff;
2592 IMRPhenomX_Initialize_Powers(&powers_of_ff, ff);
2593 pPhase->RDCoefficient[nCollocationPts_RD_Phase + 1] = -IMRPhenomXHM_RD_Phase_DerAnsatz(ff, &powers_of_ff, pWFHM, pPhase) * pow(ff, exponent + 1) / exponent;
2594 pPhase->RDCoefficient[nCollocationPts_RD_Phase] = IMRPhenomXHM_RD_Phase_Ansatz(ff, &powers_of_ff, pWFHM, pPhase) - pPhase->RDCoefficient[nCollocationPts_RD_Phase + 1] / pow(ff, exponent);
2595 pPhase->RDCoefficient[nCollocationPts_RD_Phase + 2] = IMRPhenomXHM_RD_Phase_AnsatzInt(ff, &powers_of_ff, pWFHM, pPhase) - (pPhase->RDCoefficient[5]*ff - 0.25 * pPhase->RDCoefficient[6] * powers_of_ff.m_four);
2596 pWFHM->fPhaseRDflat = ff;
2597 }
2598
2599 /*************** phase-shift of spheroidal ansatz *******************/
2600 //if(pWFHM->IMRPhenomXHMRingdownPhaseVersion == 122019 || pWFHM->IMRPhenomXHMRingdownPhaseVersion == 20221114){
2601 frefRD = pWF22->fRING;
2602 // }
2603 // else{
2604 // frefRD = pWF22->fRING - pWF22->fDAMP;
2605 // //pWFHM->phaseshift += LAL_PI;
2606 // }
2607 IMRPhenomX_Initialize_Powers(&powers_of_FREF,frefRD);
2608
2609 /* we compute phi22(fref) */
2610 IMRPhenomX_UsefulPowers powers_of_MfRef;
2611 IMRPhenomX_Initialize_Powers(&powers_of_MfRef,pWF22->MfRef);
2612 pWFHM->phiref22 = -1./pWF22->eta*IMRPhenomX_Phase_22(pWF22->MfRef, &powers_of_MfRef, pPhase22, pWF22) - pWFHM->timeshift*pWF22->MfRef - pWFHM->phaseshift + 2.0*pWF22->phi0 + LAL_PI_4;
2613 double phi22ref=1./pWF22->eta*IMRPhenomX_Phase_22(frefRD, &powers_of_FREF,pPhase22,pWF22) + pWFHM->timeshift*frefRD + pWFHM->phaseshift + pWFHM->phiref22;
2614 // printf("phi22ref = %.16e\n", phi22ref);
2615 // printf("phi32ref = %.16e\n", IMRPhenomXHM_RD_Phase_AnsatzInt(frefRD,&powers_of_FREF,pWFHM,pPhase));
2616 // printf("frefRD = %.16e\n", frefRD);
2617
2618 // we call a fit for Mod[phiS(fref)-phi22(fref),2*Pi]
2620 //if (pWFHM->IMRPhenomXHMRingdownPhaseVersion == 20221114) phishift += LAL_PI;
2621 //printf("fitphaseshift = %.6f\n", phishift);
2622
2623 //we adjust the relative phase of our reconstruction
2624 pPhase->phi0_S= phi22ref - IMRPhenomXHM_RD_Phase_AnsatzInt(frefRD,&powers_of_FREF,pWFHM,pPhase) + phishift;
2625 // printf("phaseshift = %.16e\n", pPhase->phi0_S);
2626 //
2627 // printf("phi32ref = %.16e\n", IMRPhenomXHM_RD_Phase_AnsatzInt(frefRD,&powers_of_FREF,pWFHM,pPhase));
2628
2629
2630 #if DEBUG == 1
2631 printf("**********\n alpha0_S=%.16e \n alphaL_S=%.16e \n alpha2_S=%.16e \n alpha4_S=%.16e \n \n ", pPhase->alpha0_S, pPhase->alphaL_S, pPhase->alpha2_S, pPhase->alpha4_S);
2632 printf("**********\n **Spheroidal reconstruction parameters:**\n \n phi0_S=%.16e \n alpha0_S=%.16e \n alphaL_S=%.16e \n alpha2_S=%.16e \n alpha4_S=%.16e \n \n ",pPhase->phi0_S, pPhase->alpha0_S, pPhase->alphaL_S, pPhase->alpha2_S, pPhase->alpha4_S);
2633 #endif
2634
2635
2636 }
2637
2638
2639/**************************************/
2640/* */
2641/* Spheroidal -> Spherical rotation */
2642/* */
2643/**************************************/
2644/* The rotation consists of a linear transformation using the mixing coefficients given by Berti [10.1103/PhysRevD.90.064012].
2645
2646 h32_spherical = a1 * h22_spheroidal + a2 * h32_spheroidal, where a1, a2 are the mixing coefficients.
2647
2648 Since the 22 is the most dominant mode, it will not show a significant mixing with any other mode,
2649 so we can assume that h22_spheroidal = h22_spherical
2650*/
2651
2652 // In principle this could be generalized to the 43 mode: for the time being, assume the mode solved for is only the 32.
2654 {
2655 REAL8 Mf = powers_of_Mf->itself;
2656 // Compute the 22 mode using PhenomX functions. This gives the 22 mode rescaled with the leading order. This is because the 32 is also rescaled.
2657 REAL8 amp22=XLALSimIMRPhenomXRingdownAmplitude22AnsatzAnalytical(Mf, pWF22->fRING, pWF22->fDAMP, pAmp22->gamma1, pAmp22->gamma2, pAmp22->gamma3);
2658 REAL8 phi22=1./pWF22->eta*IMRPhenomX_Phase_22(Mf, powers_of_Mf, pPhase22,pWF22) + pWFlm->timeshift*Mf + pWFlm->phaseshift + pWFlm->phiref22;
2659 COMPLEX16 wf22R = amp22 * cexp(I * phi22);
2660 if (pWFlm->IMRPhenomXHMRingdownAmpVersion != 0){
2661 wf22R *= pWFlm->ampNorm * powers_of_Mf->m_seven_sixths;
2662 }
2663 // Compute 32 mode in spheroidal.
2664 REAL8 amplm=IMRPhenomXHM_RD_Amp_Ansatz(powers_of_Mf, pWFlm, pAmplm);
2665 REAL8 philm=IMRPhenomXHM_RD_Phase_AnsatzInt(Mf, powers_of_Mf,pWFlm, pPhaselm);
2666 // Do the rotation.
2667 COMPLEX16 sphericalWF_32=conj(pWFlm->mixingCoeffs[2]) * wf22R + conj(pWFlm->mixingCoeffs[3])*amplm*cexp(I*philm);
2668 return sphericalWF_32;
2669
2670 }
2671
2672 // If the 22 mode has been previously computed, we use it here for the rotation.
2674 {
2675 // The input 22 in the whole 22, and for the rotation we have to rescaled with the leading order. This is because the 32 is also rescaled.
2676 //complex double wf22R = wf22/(powers_of_f->m_seven_sixths * pWF22->amp0);
2677 COMPLEX16 wf22R = wf22 / pWFlm->Amp0;
2678 if (pWFlm->IMRPhenomXHMRingdownAmpVersion == 0){
2679 wf22R /= (pWFlm->ampNorm * powers_of_Mf->m_seven_sixths);
2680 }
2681 // Compute 32 mode in spheroidal.
2682 REAL8 amplm=IMRPhenomXHM_RD_Amp_Ansatz(powers_of_Mf, pWFlm, pAmplm);
2683 REAL8 philm=IMRPhenomXHM_RD_Phase_AnsatzInt(powers_of_Mf->itself, powers_of_Mf, pWFlm, pPhaselm);
2684 // Do the rotation
2685 COMPLEX16 sphericalWF_32 = conj(pWFlm->mixingCoeffs[2])*wf22R + conj(pWFlm->mixingCoeffs[3])*amplm*cexp(I*philm);
2686 return sphericalWF_32;
2687
2688 }
2689
2690
2691 /**************************************/
2692 /* */
2693 /* RECONSTRUCTION FUNCTIONS THROUGH */
2694 /* INSPIRAL, MERGER & RINGDOWN */
2695 /* */
2696 /**************************************/
2697
2698 /* These functions return the amplitude/phase for a given input frequency.
2699 They are piece-wise functions that distiguish between inspiral, intermediate and ringdown
2700 and call the corresponding reconstruction function.
2701 */
2702
2703 /*******************************************************************************/
2704 /* Compute IMRPhenomXHM AMPLITUDE given an amplitude coefficients struct pAmp */
2705 /*******************************************************************************/
2706
2707 // WITHOUT mode mixing. It returns the whole amplitude (in NR units) without the normalization factor of the 22: sqrt[2 * eta / (3 * pi^(1/3))]
2709 // If it is an odd mode and equal black holes case this mode is zero.
2710 if(pWF->Ampzero==1){
2711 return 0.;
2712 }
2713
2714 // Use step function to only calculate IMR regions in approrpiate frequency regime
2715 REAL8 Amp, Mf = powers_of_Mf->itself;
2716 // Inspiral range
2717 if (!IMRPhenomX_StepFuncBool(Mf, pAmp->fAmpMatchIN)){
2718 Amp = IMRPhenomXHM_Inspiral_Amp_Ansatz(powers_of_Mf, pWF, pAmp);
2719 }
2720 // MRD range
2721 else if (IMRPhenomX_StepFuncBool(Mf, pAmp->fAmpMatchIM)){
2722 // if (pWFHM->fAmpRDfalloff > 0 && IMRPhenomX_StepFuncBool(Mf, pWFHM->fAmpRDfalloff)){
2723 // Amp = pAmp->RDCoefficient[3] * exp(- pAmp->RDCoefficient[4] * (Mf - pWFHM->fAmpRDfalloff));
2724 // }
2725 //else{
2726 Amp = IMRPhenomXHM_RD_Amp_Ansatz(powers_of_Mf, pWF, pAmp);
2727 if (pWF->IMRPhenomXHMRingdownAmpVersion == 0) Amp *= powers_of_Mf->m_seven_sixths * pWF->ampNorm;
2728 //}
2729 }
2730 /* Intermediate range */
2731 // First intermediate region.
2732 else if ((pWF->AmpEMR==1) && !IMRPhenomX_StepFuncBool(Mf, pAmp->fAmpMatchInt12)){
2733 INT4 tmp = pAmp->InterAmpPolOrder;
2734 pAmp->InterAmpPolOrder = 1042;
2735 Amp = IMRPhenomXHM_Intermediate_Amp_Ansatz(powers_of_Mf, pWF, pAmp);
2736 pAmp->InterAmpPolOrder = tmp;
2737 }
2738 else{ //Second intermediate region
2739 Amp = IMRPhenomXHM_Intermediate_Amp_Ansatz(powers_of_Mf, pWF, pAmp);
2740 }
2741 if (Amp < 0 && pWF->IMRPhenomXHMReleaseVersion != 122019) Amp = FALSE_ZERO;
2742 return Amp;
2743 }
2744
2745 // WITH mode mixing. It returns the whole amplitude (in NR units) without the normalization factor of the 22: sqrt[2 * eta / (3 * pi^(1/3))].
2747 // Use step function to only calculate IMR regions in approrpiate frequency regime
2748 REAL8 Amp, Mf = powers_of_Mf->itself;
2749 // Inspiral range
2750 if (!IMRPhenomX_StepFuncBool(Mf, pAmp->fAmpMatchIN)){
2751 Amp = IMRPhenomXHM_Inspiral_Amp_Ansatz(powers_of_Mf, pWF, pAmp);
2752 }
2753 // MRD range
2754 else if (IMRPhenomX_StepFuncBool(Mf, pAmp->fAmpMatchIM)){
2755 Amp = cabs(SpheroidalToSpherical(powers_of_Mf, pAmp22, pPhase22, pAmp, pPhase, pWF, pWF22));
2756 if (pWF->IMRPhenomXHMRingdownAmpVersion == 0) Amp *= powers_of_Mf->m_seven_sixths * pWF->ampNorm;
2757 }
2758 /* Intermediate range */
2759 // First intermediate region
2760 else if ((pWF->AmpEMR==1) && !IMRPhenomX_StepFuncBool(Mf, pAmp->fAmpMatchInt12)){
2761 INT4 tmp = pAmp->InterAmpPolOrder;
2762 pAmp->InterAmpPolOrder = 1042;
2763 Amp = IMRPhenomXHM_Intermediate_Amp_Ansatz(powers_of_Mf, pWF, pAmp);
2764 pAmp->InterAmpPolOrder = tmp;
2765 }
2766 //Second intermediate region
2767 else{
2768 Amp = IMRPhenomXHM_Intermediate_Amp_Ansatz(powers_of_Mf, pWF, pAmp);
2769 }
2770 if (Amp < 0 && pWF->IMRPhenomXHMReleaseVersion != 122019) Amp = FALSE_ZERO;
2771 return Amp;
2772 }
2773
2774 // WITH mode mixing and recycling the previously computed 22 mode. It returns the whole amplitude (in NR units) without the normalization factor of the 22: sqrt[2 * eta / (3 * pi^(1/3))].
2776 // Use step function to only calculate IMR regions in approrpiate frequency regime
2777 REAL8 Amp, Mf = powers_of_Mf->itself;
2778 // Inspiral range
2779 if (!IMRPhenomX_StepFuncBool(Mf, pAmp->fAmpMatchIN)){
2780 Amp = IMRPhenomXHM_Inspiral_Amp_Ansatz(powers_of_Mf, pWF, pAmp);
2781 }
2782 // MRD range
2783 else if (IMRPhenomX_StepFuncBool(Mf, pAmp->fAmpMatchIM)){
2784 Amp = cabs(SpheroidalToSphericalRecycle(powers_of_Mf, wf22, pAmp, pPhase, pWF));
2785 if (pWF->IMRPhenomXHMRingdownAmpVersion == 0) Amp *= powers_of_Mf->m_seven_sixths * pWF->ampNorm;
2786 }
2787 /* Intermediate range */
2788 // First intermediate region
2789 else if ((pWF->AmpEMR==1) && !IMRPhenomX_StepFuncBool(Mf, pAmp->fAmpMatchInt12)){ // FIXME: why was the second condition commented?
2790 INT4 tmp = pAmp->InterAmpPolOrder;
2791 pAmp->InterAmpPolOrder = 1042;
2792 Amp = IMRPhenomXHM_Intermediate_Amp_Ansatz(powers_of_Mf, pWF, pAmp);
2793 pAmp->InterAmpPolOrder = tmp;
2794 }
2795 //Second intermediate region
2796 else{
2797 Amp = IMRPhenomXHM_Intermediate_Amp_Ansatz(powers_of_Mf, pWF, pAmp);
2798 }
2799 if (Amp < 0 && pWF->IMRPhenomXHMReleaseVersion != 122019) Amp = FALSE_ZERO;
2800 return Amp;
2801 }
2802
2803
2804 /*******************************************************************************/
2805 /* Compute IMRPhenomXHM PHASE given a phase coefficients struct pPhase */
2806 /*******************************************************************************/
2807
2808// WITHOUT mode mixing.
2810{
2811 REAL8 Mf = powers_of_Mf->itself;
2812 // Inspiral range, f < fPhaseInsMax
2813 if (!IMRPhenomX_StepFuncBool(Mf, pPhase->fPhaseMatchIN))
2814 {
2815 REAL8 PhiIns = IMRPhenomXHM_Inspiral_Phase_AnsatzInt(Mf, powers_of_Mf, pPhase);
2816 return PhiIns + pPhase->C1INSP*Mf + pPhase->CINSP + pPhase->deltaphiLM;
2817 }
2818 // MRD range, f > fPhaseIntMax
2819 if (IMRPhenomX_StepFuncBool(Mf, pPhase->fPhaseMatchIM))
2820 {
2821 REAL8 PhiMRD = IMRPhenomXHM_RD_Phase_AnsatzInt(Mf, powers_of_Mf, pWF, pPhase);
2822 return PhiMRD + pPhase->C1RD*Mf + pPhase->CRD + pPhase->deltaphiLM;
2823 }
2824 //Intermediate range, fPhaseInsMax < f < fPhaseIntMax
2825 REAL8 PhiInt = IMRPhenomXHM_Inter_Phase_AnsatzInt(Mf, powers_of_Mf, pWF, pPhase);
2826 return PhiInt + pPhase->deltaphiLM;
2827}
2828
2829// WITH mode mixing.
2831{
2832 REAL8 Mf = powers_of_Mf->itself;
2833 // Inspiral range, f < fPhaseInsMax
2834 if (!IMRPhenomX_StepFuncBool(Mf, pPhase->fPhaseMatchIN))
2835 {
2836 REAL8 PhiIns = IMRPhenomXHM_Inspiral_Phase_AnsatzInt(Mf, powers_of_Mf, pPhase);
2837 return PhiIns + pPhase->C1INSP*Mf + pPhase->CINSP + pPhase->deltaphiLM;
2838 }
2839 // MRD range, f > fPhaseIntMax
2840 if (IMRPhenomX_StepFuncBool(Mf, pPhase->fPhaseMatchIM))
2841 {
2842 REAL8 PhiMRD = carg(SpheroidalToSpherical(powers_of_Mf, pAmp22, pPhase22, pAmp, pPhase, pWF, pWF22));
2843 return PhiMRD + pPhase->C1RD*Mf + pPhase->CRD + pPhase->deltaphiLM;
2844 }
2845 //Intermediate range, fPhaseInsMax < f < fPhaseIntMax
2846 REAL8 PhiInt = IMRPhenomXHM_Inter_Phase_AnsatzInt(Mf, powers_of_Mf, pWF, pPhase);
2847 return PhiInt + pPhase->deltaphiLM;
2848}
2849
2850// WITH mode mixing and recycling the previously computed 22 mode.
2852{
2853 REAL8 Mf = powers_of_Mf->itself;
2854 // Inspiral range, f < fPhaseInsMax
2855 if (!IMRPhenomX_StepFuncBool(Mf, pPhase->fPhaseMatchIN))
2856 {
2857 REAL8 PhiIns = IMRPhenomXHM_Inspiral_Phase_AnsatzInt(Mf, powers_of_Mf, pPhase);
2858 return PhiIns + pPhase->C1INSP*Mf + pPhase->CINSP + pPhase->deltaphiLM;
2859 }
2860 // MRD range, f > fPhaseIntMax
2861 if (IMRPhenomX_StepFuncBool(Mf, pPhase->fPhaseMatchIM))
2862 {
2863 REAL8 PhiMRD = carg(SpheroidalToSphericalRecycle(powers_of_Mf, wf22, pAmp, pPhase, pWF));
2864 return PhiMRD + pPhase->C1RD*Mf + pPhase->CRD + pPhase->deltaphiLM;
2865 }
2866 //Intermediate range, fPhaseInsMax < f < fPhaseIntMax
2867 REAL8 PhiInt = IMRPhenomXHM_Inter_Phase_AnsatzInt(Mf, powers_of_Mf, pWF, pPhase);
2868 return PhiInt + pPhase->deltaphiLM;
2869}
2870
2871
2872 // WITHOUT mode mixing.
2874 {
2875 // Inspiral range, f < fPhaseInsMax
2876 if (!IMRPhenomX_StepFuncBool(f, pPhase->fPhaseMatchIN))
2877 {
2878 double dPhiIns = IMRPhenomXHM_Inspiral_Phase_Ansatz(f, powers_of_f, pPhase);
2879 return dPhiIns + pPhase->C1INSP;
2880 }
2881 // MRD range, f > fPhaseIntMax
2882 if (IMRPhenomX_StepFuncBool(f, pPhase->fPhaseMatchIM))
2883 {
2884 double dPhiMRD = IMRPhenomXHM_RD_Phase_Ansatz(f, powers_of_f, pWF, pPhase);
2885 return dPhiMRD + pPhase->C1RD;
2886 }
2887 //Intermediate range, fPhaseInsMax < f < fPhaseIntMax
2888 double dPhiInt = IMRPhenomXHM_Inter_Phase_Ansatz(f, powers_of_f, pWF, pPhase);
2889 return dPhiInt;
2890 }
2891
2892 // WITH mode mixing.
2894 {
2895 // Inspiral range, f < fPhaseInsMax
2896 if (!IMRPhenomX_StepFuncBool(f, pPhase->fPhaseMatchIN))
2897 {
2898 double dPhiIns = IMRPhenomXHM_Inspiral_Phase_Ansatz(f, powers_of_f, pPhase);
2899 return dPhiIns + pPhase->C1INSP;
2900 }
2901 // MRD range, f > fPhaseIntMax
2902 if (IMRPhenomX_StepFuncBool(f, pPhase->fPhaseMatchIM))
2903 {
2904 double dPhiMRD = carg(SpheroidalToSpherical(powers_of_f, pAmp22, pPhase22, pAmp, pPhase, pWF, pWF22));
2905 return dPhiMRD + pPhase->C1RD;
2906 }
2907 //Intermediate range, fPhaseInsMax < f < fPhaseIntMax
2908 double dPhiInt = IMRPhenomXHM_Inter_Phase_Ansatz(f, powers_of_f, pWF, pPhase);
2909 return dPhiInt;
2910 }
2911
2912
2913
2914 /*****************/
2915 /* DEBUGGING */
2916 /*****************/
2917
2918 // This is just for debugging. It prints some interesting parameters to a file.
2920 IMRPhenomXWaveformStruct *pWF, /**< Wf structure for the 22 mode*/
2921 IMRPhenomXHMWaveformStruct *pWFHM, /**< Wf structure for the lm mode*/
2922 IMRPhenomXHMAmpCoefficients *pAmp, /**< Coefficients struct of the lm Amplitude */
2923 UNUSED IMRPhenomXHMPhaseCoefficients *pPhase /**< Coefficients struct of the lm Phase */
2924 )
2925 {
2926
2927 FILE *file;
2928 char fileSpec[40];
2929 sprintf(fileSpec, "Parameters%i.dat", pWFHM->modeTag);
2930 printf("\nOutput Parameter file: %s\r\n",fileSpec);
2931 file = fopen(fileSpec,"w");
2932
2933 fprintf(file,"\n*** %i Mode ***\n", pWFHM->modeTag);
2934
2935 fprintf(file,"\n*** Intrinsic Parameters ***\n");
2936 fprintf(file,"eta = %.16f\n", pWF->eta);
2937 fprintf(file,"chi1z = %.16f\n", pWF->chi1L);
2938 fprintf(file,"chi2z = %.16f\n", pWF->chi2L);
2939
2940 fprintf(file,"\n*** Extrinsic Parameters ***\n");
2941 fprintf(file,"distance = %.16f\n", pWF->distance);
2942 fprintf(file,"inclination = %.16f\n", pWF->inclination);
2943 fprintf(file,"phiRef = %.16f\n", pWF->phifRef);
2944 fprintf(file,"fRef = %.16f\n", pWF->fRef);
2945 fprintf(file,"m1_SI = %.16f\n", pWF->m1_SI);
2946 fprintf(file,"m2_SI = %.16f\n", pWF->m2_SI);
2947
2948 fprintf(file,"\n*** Frequency Grid ***\n");
2949 fprintf(file,"deltaF = %.16f\n", pWF->deltaF);
2950 fprintf(file,"fMin = %.16f\n", pWF->fMin);
2951 fprintf(file,"fMax = %.16f\n", pWF->fMax);
2952
2953 fprintf(file,"\n*** QNM Frequencies ***\n");
2954 fprintf(file,"fRING22 = %.16f\n", pWF->fRING);
2955 fprintf(file,"fDAMP22 = %.16f\n", pWF->fDAMP);
2956 fprintf(file,"fRINGlm = %.16f\n", pWFHM->fRING);
2957 fprintf(file,"fDAMPlm = %.16f\n", pWFHM->fDAMP);
2958
2959 fprintf(file,"\n*** Amplitude Frequencies ***\n");
2960 fprintf(file,"fInsp1 = %.16f\n", pAmp->CollocationPointsFreqsAmplitudeInsp[0]);
2961 fprintf(file,"fInsp2 = %.16f\n", pAmp->CollocationPointsFreqsAmplitudeInsp[1]);
2962 fprintf(file,"fInsp3 = %.16f\n", pAmp->CollocationPointsFreqsAmplitudeInsp[2]);
2963 fprintf(file,"fAmpMatchIN = %.16f\n", pAmp->fAmpMatchIN);
2964 fprintf(file,"fAmpMatchInt12 = %.16f\n", pAmp->fAmpMatchInt12);
2965 fprintf(file,"fInt1 = %.16f\n", pAmp->CollocationPointsFreqsAmplitudeInter[0]);
2966 fprintf(file,"fInt2 = %.16f\n", pAmp->CollocationPointsFreqsAmplitudeInter[1]);
2967 fprintf(file,"fAmpMatchIM = %.16f\n", pAmp->fAmpMatchIM);
2968
2969 fprintf(file,"\n*** Amplitude Collocation Points ***\n");
2970 fprintf(file,"A_fInsp1 = %.16f\n", pAmp->CollocationPointsValuesAmplitudeInsp[0]);
2971 fprintf(file,"A_fInsp2 = %.16f\n", pAmp->CollocationPointsValuesAmplitudeInsp[1]);
2972 fprintf(file,"A_fInsp3 = %.16f\n", pAmp->CollocationPointsValuesAmplitudeInsp[2]);
2973 fprintf(file,"A_fInt1 = %.16f\n", pAmp->CollocationPointsValuesAmplitudeInter[0]);
2974 fprintf(file,"A_fInt2 = %.16f\n", pAmp->CollocationPointsValuesAmplitudeInter[1]);
2975
2976 fprintf(file,"\n*** Amplitude Coefficients ***\n");
2977 fprintf(file,"rho1 = %.16f\n", pAmp->rho1); //Inspiral coefficients
2978 fprintf(file,"rho2 = %.16f\n", pAmp->rho2);
2979 fprintf(file,"rho3 = %.16f\n", pAmp->rho3);
2980 fprintf(file,"delta0 = %.16f\n", pAmp->delta0); //Intermediate region
2981 fprintf(file,"delta1 = %.16f\n", pAmp->delta1);
2982 fprintf(file,"delta2 = %.16f\n", pAmp->delta2);
2983 fprintf(file,"delta3 = %.16f\n", pAmp->delta3);
2984 fprintf(file,"delta4 = %.16f\n", pAmp->delta4);
2985 fprintf(file,"delta5 = %.16f\n", pAmp->delta5);
2986 fprintf(file,"alambda = %.20f\n", pAmp->RDCoefficient[0]); //Ringdown region
2987 fprintf(file,"lambda = %.16f\n", pAmp->RDCoefficient[1]);
2988 fprintf(file,"sigma = %.16f\n", pAmp->RDCoefficient[2]);
2989 fprintf(file,"lc = %.16f\n", pAmp->RDCoefficient[3]);
2990 fprintf(file,"alpha0 = %.16f\n", pAmp->alpha0); //First Intermediate region (only for EMR)
2991 fprintf(file,"alpha1 = %.16f\n", pAmp->alpha1);
2992 fprintf(file,"alpha2 = %.16f\n", pAmp->alpha2);
2993 fprintf(file,"alpha3 = %.16f\n", pAmp->alpha3);
2994 fprintf(file,"alpha4 = %.16f\n", pAmp->alpha4);
2995
2996 fprintf(file,"\n*** Amplitude Intermediate Input ***\n");
2997 fprintf(file,"f1 = %.16f\n", pAmp->f1);
2998 fprintf(file,"f2 = %.16f\n", pAmp->f2);
2999 fprintf(file,"f3 = %.16f\n", pAmp->f3);
3000 fprintf(file,"f4 = %.16f\n", pAmp->f4);
3001 fprintf(file,"v1 = %.16f\n", pAmp->v1);
3002 fprintf(file,"v2 = %.16f\n", pAmp->v2);
3003 fprintf(file,"v3 = %.16f\n", pAmp->v3);
3004 fprintf(file,"v4 = %.16f\n", pAmp->v4);
3005 fprintf(file,"d1 = %.16f\n", pAmp->d1);
3006 fprintf(file,"d4 = %.16f\n", pAmp->d4);
3007
3008 fprintf(file,"\n*** PN Amplitude Inspiral ***\n");
3009 fprintf(file,"PN_f1 = %.16f\n", pAmp->PNAmplitudeInsp[0]);
3010 fprintf(file,"PN_f2 = %.16f\n", pAmp->PNAmplitudeInsp[1]);
3011 fprintf(file,"PN_f3 = %.16f\n", pAmp->PNAmplitudeInsp[2]);
3012
3013 fprintf(file,"\n*** Amplitude Versions ***\n");
3014 fprintf(file,"InspiralFits = %i\n", pWFHM->IMRPhenomXHMInspiralAmpFitsVersion);
3015 fprintf(file,"IntermediateFits = %i\n", pWFHM->IMRPhenomXHMIntermediateAmpFitsVersion);
3016 fprintf(file,"RDFits = %i\n", pWFHM->IMRPhenomXHMRingdownAmpFitsVersion);
3017 fprintf(file,"InspiralRecon = %i\n", pWFHM->IMRPhenomXHMInspiralAmpVersion);
3018 fprintf(file,"IntermediateRecon = %i\n", pWFHM->IMRPhenomXHMIntermediateAmpVersion);
3019 fprintf(file,"RDRecon = %i\n", pWFHM->IMRPhenomXHMRingdownAmpVersion);
3020
3021 fprintf(file,"\n*** Amplitude Extra ***\n");
3022 fprintf(file,"AmpEMR = %i\n", pWFHM->AmpEMR);
3023 fprintf(file,"Ampzero = %i\n", pWFHM->Ampzero);
3024
3025 if(pWFHM->modeTag==32){
3026 fprintf(file,"\n*** Mixing Coefficients ***\n");
3027 fprintf(file,"222 = %.16f %.16f\n", creal(pWFHM->mixingCoeffs[0]), cimag(pWFHM->mixingCoeffs[0]));
3028 fprintf(file,"223 = %.16f %.16f\n", creal(pWFHM->mixingCoeffs[1]), cimag(pWFHM->mixingCoeffs[1]));
3029 fprintf(file,"322 = %.16f %.16f\n", creal(pWFHM->mixingCoeffs[2]), cimag(pWFHM->mixingCoeffs[2]));
3030 fprintf(file,"323 = %.16f %.16f\n", creal(pWFHM->mixingCoeffs[3]), cimag(pWFHM->mixingCoeffs[3]));
3031 }
3032
3033 fclose(file);
3034
3035 return 0;
3036 }
double XLALSimIMRPhenomXCP_ZETA1_l3m3(double theta, double eta, double a1)
double XLALSimIMRPhenomXCP_NU5_l3m3(double theta, double eta, double a1)
double XLALSimIMRPhenomXCP_MU1_l3m3(double theta, double eta, double a1)
double XLALSimIMRPhenomXCP_MU3_l3m3(double theta, double eta, double a1)
double XLALSimIMRPhenomXCP_MU4_l3m3(double theta, double eta, double a1)
double XLALSimIMRPhenomXCP_NU6_l3m3(double theta, double eta, double a1)
double XLALSimIMRPhenomXCP_ZETA2_l3m3(double theta, double eta, double a1)
double XLALSimIMRPhenomXCP_NU4_l3m3(double theta, double eta, double a1)
double XLALSimIMRPhenomXCP_MU2_l3m3(double theta, double eta, double a1)
static double fdamp(double eta, double chi1, double chi2, double finspin)
fdamp is the complex part of the ringdown frequency 1508.07250 figure 9
static double fring(double eta, double chi1, double chi2, double finspin)
fring is the real part of the ringdown frequency 1508.07250 figure 9
static double evaluate_QNMfit_fring33(double finalDimlessSpin)
static double evaluate_QNMfit_fdamp33(double finalDimlessSpin)
static double evaluate_QNMfit_fdamp44(double finalDimlessSpin)
static double evaluate_QNMfit_fdamp21(double finalDimlessSpin)
static double evaluate_QNMfit_fring21(double finalDimlessSpin)
static double evaluate_QNMfit_fring44(double finalDimlessSpin)
static double double delta
double IMRPhenomX_TimeShift_22(IMRPhenomXPhaseCoefficients *pPhase, IMRPhenomXWaveformStruct *pWF)
void IMRPhenomX_Phase_22_ConnectionCoefficients(IMRPhenomXWaveformStruct *pWF, IMRPhenomXPhaseCoefficients *pPhase)
int IMRPhenomX_Initialize_Powers(IMRPhenomX_UsefulPowers *p, REAL8 number)
double IMRPhenomX_Phase_22(double ff, IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXPhaseCoefficients *pPhase, IMRPhenomXWaveformStruct *pWF)
double IMRPhenomX_dPhase_22(double ff, IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXPhaseCoefficients *pPhase, IMRPhenomXWaveformStruct *pWF)
IMRPhenomX_UsefulPowers powers_of_lalpiHM
void IMRPhenomXHM_Get_Inspiral_Amp_Coefficients(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
static double IMRPhenomXHM_Inspiral_Phase_AnsatzInt(double Mf, IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMPhaseCoefficients *pPhase)
static double IMRPhenomXHM_Insp_Amp_21_iv2(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Inspiral_Amp_rho3(double v1, double v2, double v3, IMRPhenomX_UsefulPowers *powers_of_fcutInsp, IMRPhenomX_UsefulPowers *powers_of_f1, IMRPhenomX_UsefulPowers *powers_of_f2, IMRPhenomX_UsefulPowers *powers_of_f3, IMRPhenomXHMWaveformStruct *pWFHM)
static double IMRPhenomXHM_Inspiral_Amp_NDAnsatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
static double IMRPhenomXHM_Inspiral_Amp_Ansatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
static double IMRPhenomXHM_Insp_Amp_44_iv2(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Inspiral_Amp_rho2(double v1, double v2, double v3, IMRPhenomX_UsefulPowers *powers_of_fcutInsp, IMRPhenomX_UsefulPowers *powers_of_f1, IMRPhenomX_UsefulPowers *powers_of_f2, IMRPhenomX_UsefulPowers *powers_of_f3, IMRPhenomXHMWaveformStruct *pWFHM)
static double IMRPhenomXHM_Insp_Amp_33_iv1(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Insp_Phase_33_lambda(IMRPhenomXWaveformStruct *pWF, int InspPhaseFlag)
static double IMRPhenomXHM_Insp_Amp_33_iv2(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Insp_Amp_32_iv3(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Insp_Phase_21_lambda(IMRPhenomXWaveformStruct *pWF, int InspPhaseFlag)
static double IMRPhenomXHM_Insp_Phase_44_lambda(IMRPhenomXWaveformStruct *pWF, int InspPhaseFlag)
static double IMRPhenomXHM_Insp_Amp_32_iv1(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Insp_Phase_32_lambda(IMRPhenomXWaveformStruct *pWF, int InspPhaseFlag)
static double IMRPhenomXHM_Insp_Amp_33_iv3(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Insp_Amp_44_iv1(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
void IMRPhenomXHM_Inspiral_Amplitude_Veto(double *iv1, double *iv2, double *iv3, IMRPhenomX_UsefulPowers *powers_of_f1, IMRPhenomX_UsefulPowers *powers_of_f2, IMRPhenomX_UsefulPowers *powers_of_f3, IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMWaveformStruct *pWFHM)
static double IMRPhenomXHM_Insp_Amp_21_iv3(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
int WavyPoints(double p1, double p2, double p3)
static double IMRPhenomXHM_Insp_Amp_44_iv3(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Inspiral_PNAmp_Ansatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
static double IMRPhenomXHM_Insp_Amp_21_iv1(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Inspiral_Phase_Ansatz(double Mf, IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMPhaseCoefficients *pPhase)
static double IMRPhenomXHM_Inspiral_Amp_rho1(double v1, double v2, double v3, IMRPhenomX_UsefulPowers *powers_of_fcutInsp, IMRPhenomX_UsefulPowers *powers_of_f1, IMRPhenomX_UsefulPowers *powers_of_f2, IMRPhenomX_UsefulPowers *powers_of_f3, IMRPhenomXHMWaveformStruct *pWFHM)
static double IMRPhenomXHM_Insp_Phase_LambdaPN(double eta, int modeInt)
static double IMRPhenomXHM_Insp_Amp_32_iv2(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Intermediate_Amp_delta0(double d1, double d4, double v1, double v2, double v3, double v4, double f1, double f2, double f3, double f4, int IntAmpFlag)
static double IMRPhenomXHM_Inter_Phase_21_p3(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Amp_21_int4(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
void ChoosePolOrder(IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
static double IMRPhenomXHM_Inter_Amp_44_int2(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Amp_33_int1(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Amp_33_int2(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Amp_32_int2(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Amp_21_dint0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Phase_44_p5(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Amp_32_int0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Phase_44_p6(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Phase_32_p6(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Phase_21_p6(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Amp_33_int3(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Amp_44_int4(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Phase_Ansatz(double ff, IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMPhaseCoefficients *pPhase)
void IMRPhenomXHM_Intermediate_Amplitude_Veto(double *int1, double *int2, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
static double IMRPhenomXHM_Inter_Amp_33_int4(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Phase_33_p3(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Amp_21_int1(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Amp_32_dint0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
void Update_Intermediate_Amplitude_Coefficients(IMRPhenomXHMAmpCoefficients *pAmp, int IntAmpFlag)
static double IMRPhenomXHM_Inter_Phase_44_p1(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Phase_33_p4(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Phase_32_p4(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Amp_32_int3(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Phase_44_p2(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Phase_44_p4(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Phase_21_p5(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Amp_21_int0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Phase_32_p1(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Intermediate_Amp_delta3(double d1, double d4, double v1, double v2, double v3, double v4, double f1, double f2, double f3, double f4, int IntAmpFlag)
static double IMRPhenomXHM_Inter_Amp_21_int2(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Amp_44_int0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Intermediate_Amp_delta2(double d1, double d4, double v1, double v2, double v3, double v4, double f1, double f2, double f3, double f4, int IntAmpFlag)
static double IMRPhenomXHM_Inter_Amp_44_int3(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Phase_AnsatzInt(double ff, IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMPhaseCoefficients *pPhase)
static double IMRPhenomXHM_Inter_Amp_32_int1(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Amp_32_int4(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Intermediate_Amp_delta4(double d1, double d4, double v1, double v2, double v3, double v4, double f1, double f2, double f3, double f4, int IntAmpFlag)
static double IMRPhenomXHM_Inter_Phase_44_p3(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Phase_33_p1(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Amp_44_dint0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Intermediate_Amp_Ansatz(IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
static double IMRPhenomXHM_Inter_Phase_32_p5(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Amp_21_int3(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Phase_33_p5(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Amp_33_dint0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Phase_33_p2(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Phase_21_p2(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Phase_33_p6(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Amp_33_int0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Amp_44_int1(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
static double IMRPhenomXHM_Inter_Phase_32_p3(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Phase_21_p4(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Inter_Phase_32_p2(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
static double IMRPhenomXHM_Intermediate_Amp_delta1(double d1, double d4, double v1, double v2, double v3, double v4, double f1, double f2, double f3, double f4, int IntAmpFlag)
void IMRPhenomXHM_Intermediate_Amp_Coefficients(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22)
static double IMRPhenomXHM_Inter_Phase_21_p1(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
int IMRPhenomXHM_PN21AmpSign(double ff, IMRPhenomXWaveformStruct *wf22)
REAL8 IMRPhenomXHM_Amplitude_ModeMixing(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXHMWaveformStruct *pWF, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22, IMRPhenomXWaveformStruct *pWF22)
void IMRPhenomXHM_SetHMWaveformVariables(int ell, int emm, IMRPhenomXHMWaveformStruct *wf, IMRPhenomXWaveformStruct *wf22, QNMFits *qnms, LALDict *LALParams)
double IMRPhenomXHM_Amplitude_fcutInsp(IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
int ParametersToFile(IMRPhenomXWaveformStruct *pWF, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp, UNUSED IMRPhenomXHMPhaseCoefficients *pPhase)
void IMRPhenomXHM_GetPNAmplitudeCoefficients(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
void IMRPhenomXHM_Ringdown_CollocPtsFreqs(IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
void Get21PNAmplitudeCoefficients(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXWaveformStruct *pWF22)
void IMRPhenomXHM_Initialize_MixingCoeffs(IMRPhenomXHMWaveformStruct *wf, IMRPhenomXWaveformStruct *wf22)
double IMRPhenomXHM_dPhase_ModeMixing(double f, IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXHMWaveformStruct *pWF, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22, IMRPhenomXWaveformStruct *pWF22)
REAL8 IMRPhenomXHM_Phase_ModeMixingRecycle(IMRPhenomX_UsefulPowers *powers_of_Mf, COMPLEX16 wf22, IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXHMWaveformStruct *pWF)
void IMRPhenomXHM_GetAmplitudeCoefficients(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
void IMRPhenomXHM_GetPhaseCoefficients(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22, UNUSED LALDict *lalParams)
COMPLEX16 SpheroidalToSpherical(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22, IMRPhenomXHMAmpCoefficients *pAmplm, IMRPhenomXHMPhaseCoefficients *pPhaselm, IMRPhenomXHMWaveformStruct *pWFlm, IMRPhenomXWaveformStruct *pWF22)
REAL8 IMRPhenomXHM_Amplitude_ModeMixingRecycle(IMRPhenomX_UsefulPowers *powers_of_Mf, COMPLEX16 wf22, IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXHMWaveformStruct *pWF)
REAL8 IMRPhenomXHM_Amplitude_noModeMixing(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMWaveformStruct *pWF)
void IMRPhenomXHM_Intermediate_CollocPtsFreqs(IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
double GetfcutInsp(IMRPhenomXWaveformStruct *pWF22, IMRPhenomXHMWaveformStruct *pWFHM)
double IMRPhenomXHM_Amplitude_fcutRD(IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
void GetSpheroidalCoefficients(IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXPhaseCoefficients *pPhase22, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
void IMRPhenomXHM_FillAmpFitsArray(IMRPhenomXHMAmpCoefficients *pAmp)
void IMRPhenomXHM_FillPhaseFitsArray(IMRPhenomXHMPhaseCoefficients *pPhase)
REAL8 IMRPhenomXHM_Phase_noModeMixing(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXHMWaveformStruct *pWF, UNUSED IMRPhenomXWaveformStruct *pWF22)
double IMRPhenomXHM_dPhase_noModeMixing(double f, IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXHMWaveformStruct *pWF, UNUSED IMRPhenomXWaveformStruct *pWF22)
REAL8 IMRPhenomXHM_Phase_ModeMixing(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXHMWaveformStruct *pWF, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22, IMRPhenomXWaveformStruct *pWF22)
REAL8 IMRPhenomXHM_GenerateRingdownFrequency(UINT4 ell, UINT4 emm, IMRPhenomXWaveformStruct *wf22)
COMPLEX16 SpheroidalToSphericalRecycle(IMRPhenomX_UsefulPowers *powers_of_Mf, COMPLEX16 wf22, IMRPhenomXHMAmpCoefficients *pAmplm, IMRPhenomXHMPhaseCoefficients *pPhaselm, IMRPhenomXHMWaveformStruct *pWFlm)
double RescaleFactor(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMAmpCoefficients *pAmp, UINT2 rescalefactor)
void IMRPhenomXHM_Initialize_QNMs(QNMFits *qnms)
#define NMAX_INSPIRAL_COEFFICIENTS
#define FALSE_ZERO
double evaluate_QNMfit_re_l2m2lp3(double finalDimlessSpin)
double evaluate_QNMfit_re_l2m2lp2(double finalDimlessSpin)
double evaluate_QNMfit_im_l2m2lp3(double finalDimlessSpin)
double evaluate_QNMfit_fdamp32(double finalDimlessSpin)
double evaluate_QNMfit_im_l3m2lp3(double finalDimlessSpin)
double evaluate_QNMfit_im_l2m2lp2(double finalDimlessSpin)
double evaluate_QNMfit_re_l3m2lp3(double finalDimlessSpin)
double evaluate_QNMfit_fring32(double finalDimlessSpin)
double evaluate_QNMfit_im_l3m2lp2(double finalDimlessSpin)
double evaluate_QNMfit_re_l3m2lp2(double finalDimlessSpin)
static double IMRPhenomXHM_RD_Amp_44_rdcp1(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Phase_32_SpheroidalTimeShift(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag)
static double IMRPhenomXHM_RD_Phase_DerAnsatz(double ff, IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMPhaseCoefficients *pPhase)
static double IMRPhenomXHM_RD_Amp_33_rdcp1(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Phase_22_alpha2(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag)
static double IMRPhenomXHM_RD_Amp_DAnsatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWF, IMRPhenomXHMAmpCoefficients *pAmp)
static void IMRPhenomXHM_RD_Amp_Coefficients(IMRPhenomXWaveformStruct *pWF22, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
static double IMRPhenomXHM_RD_Amp_21_rdcp1(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_44_lambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
void IMRPhenomXHM_Ringdown_Amplitude_Veto(double *V2, double *V3, double V4, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
static double IMRPhenomXHM_RD_Amp_32_lambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_21_alambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Phase_32_p5(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag)
static double IMRPhenomXHM_RD_Phase_32_SpheroidalPhaseShift(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag)
static double IMRPhenomXHM_RD_Amp_32_rdcp1(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Phase_AnsatzInt(double ff, IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMPhaseCoefficients *pPhase)
static double IMRPhenomXHM_RD_Amp_33_rdcp3(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_21_rdcp2(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_33_lambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_Ansatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
static double IMRPhenomXHM_RD_Phase_32_p1(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag)
static double IMRPhenomXHM_RD_Amp_32_rdaux1(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_21_rdcp3(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_32_rdaux2(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_32_rdcp3(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_32_rdcp2(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Phase_32_p3(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag)
static double IMRPhenomXHM_RD_Amp_33_alambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_33_sigma(UNUSED IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Phase_32_p2(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag)
static double IMRPhenomXHM_RD_Amp_44_rdcp2(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_32_sigma(UNUSED IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_33_rdcp2(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Phase_22_alphaL(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag)
static double IMRPhenomXHM_RD_Amp_NDAnsatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22, IMRPhenomXWaveformStruct *pWF22)
static double IMRPhenomXHM_RD_Phase_Ansatz(double ff, IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMPhaseCoefficients *pPhase)
static double IMRPhenomXHM_RD_Phase_32_p4(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag)
static double IMRPhenomXHM_RD_Amp_44_sigma(UNUSED IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_32_alambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_44_rdcp3(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_21_sigma(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_21_lambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
static double IMRPhenomXHM_RD_Amp_44_alambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag)
#define N_MAX_COEFFICIENTS_PHASE_INTER
bool IMRPhenomX_StepFuncBool(const double t, const double t1)
REAL8 XLALSimIMRPhenomXRingdownAmplitude22AnsatzAnalytical(REAL8 ff, REAL8 fRD, REAL8 fDA, REAL8 gamma1, REAL8 gamma2, REAL8 gamma3)
"Analytical" phenomenological ringdown ansatz for amplitude.
REAL8 XLALSimIMRPhenomXPsi4ToStrain(double eta, double S, double dchi)
This is a fit of the time-difference between t_peak of strain and t_peak of psi4 needed to align in t...
#define c
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMInspiralAmpVersion(LALDict *params)
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMInspiralPhaseVersion(LALDict *params)
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMInspiralAmpFitsVersion(LALDict *params)
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMIntermediateAmpFreqsVersion(LALDict *params)
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMReleaseVersion(LALDict *params)
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMIntermediatePhaseVersion(LALDict *params)
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMRingdownPhaseVersion(LALDict *params)
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMRingdownAmpFreqsVersion(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupPhenomXCPMU1l3m3(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupPhenomXCPNU5l3m3(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupPhenomXCPZETA2l3m3(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupPhenomXCPNU6l3m3(LALDict *params)
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMIntermediateAmpFitsVersion(LALDict *params)
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMIntermediateAmpVersion(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupPhenomXCPZETA1l3m3(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupPhenomXCPNU0l3m3(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupPhenomXCPMU3l3m3(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupPhenomXCPNU4l3m3(LALDict *params)
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMRingdownAmpFitsVersion(LALDict *params)
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMRingdownAmpVersion(LALDict *params)
INT4 XLALSimInspiralWaveformParamsLookupPhenomXHMInspiralAmpFreqsVersion(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupPhenomXCPMU2l3m3(LALDict *params)
REAL8 XLALSimInspiralWaveformParamsLookupPhenomXCPMU4l3m3(LALDict *params)
#define fprintf
int s
Definition: bh_qnmode.c:137
double i
Definition: bh_ringdown.c:118
#define LAL_PI
#define LAL_PI_4
double complex COMPLEX16
double REAL8
uint16_t UINT2
uint32_t UINT4
int32_t INT4
static const INT4 r
#define XLAL_ERROR_VOID(...)
#define XLAL_ERROR_REAL8(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_EFUNC
XLAL_EDOM
XLAL_EINVAL
string version
p
x
REAL8 RDCoefficient[N_MAX_COEFFICIENTS_AMPLITUDE_RING]
ParameterSpaceFit RingdownAmpFits[N_HIGHERMODES_IMPLEMENTED *N_MAX_COEFFICIENTS_AMPLITUDE_RING+N_MAX_COEFFICIENTS_AMPLITUDE_RDAUX]
REAL8 CollocationPointsFreqsAmplitudeInter[N_MAX_COEFFICIENTS_AMPLITUDE_INTER]
REAL8 CollocationPointsFreqsAmplitudeInsp[N_MAX_COEFFICIENTS_AMPLITUDE_INS]
UINT2 VersionCollocPtsInter[N_MAX_COEFFICIENTS_AMPLITUDE_INTER]
ParameterSpaceFit IntermediateAmpFits[N_HIGHERMODES_IMPLEMENTED *N_MAX_COEFFICIENTS_AMPLITUDE_INTER]
REAL8 PNAmplitudeInsp[N_MAX_COEFFICIENTS_AMPLITUDE_INS]
REAL8 CollocationPointsValuesAmplitudeInter[N_MAX_COEFFICIENTS_AMPLITUDE_INTER]
REAL8 CollocationPointsValuesAmplitudeInsp[N_MAX_COEFFICIENTS_AMPLITUDE_INS]
ParameterSpaceFit InspiralAmpFits[N_HIGHERMODES_IMPLEMENTED *N_MAX_COEFFICIENTS_AMPLITUDE_INS]
REAL8 phiL[N_MAX_COEFFICIENTS_PHASE_INS]
ParameterSpaceFit RingdownPhaseFits[N_MAX_COEFFICIENTS_PHASE_RING]
REAL8 CollocationPointsFreqsPhaseRD[N_MAX_COEFFICIENTS_PHASE_RING]
REAL8 CollocationPointsFreqsPhaseInter[N_MAX_COEFFICIENTS_PHASE_INTER]
REAL8 phi[N_MAX_COEFFICIENTS_PHASE_INS]
ParameterSpaceFit IntermediatePhaseFits[N_HIGHERMODES_IMPLEMENTED *N_MAX_COEFFICIENTS_PHASE_INTER]
ParameterSpaceFit InspiralPhaseFits[N_HIGHERMODES_IMPLEMENTED]
REAL8 CollocationPointsValuesPhaseInter[N_MAX_COEFFICIENTS_PHASE_INTER]
REAL8 RDCoefficient[N_MAX_COEFFICIENTS_PHASE_RING+3]
fitQNM_fring fring_lm[N_HIGHERMODES_IMPLEMENTED]
fitQNM_fdamp fdamp_lm[N_HIGHERMODES_IMPLEMENTED]
char output[FILENAME_MAX]
Definition: unicorn.c:26