LALSimulation  5.4.0.1-fe68b98
LALSimIMRPhenomXHM_ringdown.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 
22 
23 /* Include fits of the Ringdown related quantities and the ansatz for amplitude and phase */
24 
25 /* These parameter-space fits are documented in the supplementary material of https://dcc.ligo.org/P2000011-v2.
26  There are 2 Mathematica notebooks (one for amplitude and one for phase) that read the fits data and automatically generate the C-code below.
27  For more information read https://git.ligo.org/waveforms/reviews/imrphenomx/blob/master/documentation/ParspaceFits/README and the documentation in the notebooks. */
28 
29 
30 /****************************************/
31 /* */
32 /* AMPLITUDE */
33 /* */
34 /****************************************/
35 
36 /* Fits over parameter space for the ringdown amplitude coefficients (alambda, lambda, sigma) for each mode. */
37 
38 // The spin parameter S = (m1^2*chi1 + m2^2*chi2)/(m1^2 + m2^2)
39 
40 
41 // alambda, lambda and sigma are the coefficients of the ringdown ansatz, see Eq. (6.2)
42 
43 /* Start of Amp Parameter Space Fits */
44 
45 
46 static double IMRPhenomXHM_RD_Amp_21_alambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag) {
47  double total=0;
48  switch (RDAmpFlag){
49  case 122018:{
50  double eta = pWF->eta;
51  double S = pWF->STotR;
52  double delta=sqrt(1.-4.*eta),eta2,S2;
53  eta2 = pow(eta,2);
54  S2 = pow(S,2);
55  double noSpin = sqrt(eta - 4.*eta2)*(0.00734983387668636 - 0.0012619735607202085*eta + 0.01042318959002753*eta2);
56  double eqSpin = sqrt(eta - 4.*eta2)*S*(-0.004839645742570202 - 0.0013927779195756036*S + eta2*(-0.054621206928483663 + 0.025956604949552205*S + 0.020360826886107204*S2));
57  double uneqSpin = -0.018115657394753674*pWF->dchi*eta2*(1. - 10.539795474715346*eta2*delta);
58  total = noSpin + eqSpin + uneqSpin;
59  break;
60  }
61  case 122022:{
62  double eta = pWF->eta;
63  double delta = pWF->delta;
64  double S = pWF->chiPNHat;
65  double chidiff = pWF->dchi_half;
66  double eta1 = eta;
67  double eta2 = eta1 * eta1;
68  double eta3 = eta1 * eta2;
69  double eta4 = eta1 * eta3;
70  double eta5 = eta1 * eta4;
71  double S1 = S;
72  double S2 = S1 * S1;
73  double chidiff1 = chidiff;
74  total = fabs(delta*(0.24548180919287976 - 0.25565119457386487*eta1)*eta1 + chidiff1*delta*eta1*(0.5670798742968471*eta1 - 14.276514548218454*eta2 + 45.014547333879136*eta3) + chidiff1*delta*eta1*(0.4580805242442763*eta1 - 4.859294663135058*eta2 + 14.995447609839573*eta3)*S1 + chidiff1*eta5*(-27.031582936285528 + 6.468164760468401*S1 + 0.34222101136488015*S2) + delta*eta1*S1*(-0.2204878224611389*(1.0730799832007898 - 3.44643820338605*eta1 + 32.520429274459836*eta2 - 83.21097158567372*eta3) + 0.008901444811471891*(-5.876973170072921 + 120.70115519895002*eta1 - 916.5281661566283*eta2 + 2306.8425350489847*eta3)*S1 + 0.015541783867953005*(2.4780170686140455 + 17.377013149762398*eta1 - 380.91157168170236*eta2 + 1227.5332509075172*eta3)*S2));
75  break;
76  }
77  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_21_alambda: version %i is not valid.", RDAmpFlag);}
78  }
79  return total;
80 }
81 
82 static double IMRPhenomXHM_RD_Amp_21_lambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag) {
83  double total=0;
84  switch (RDAmpFlag){
85  case 122018:{
86  double eta = pWF->eta;
87  double S = pWF->STotR;
88  double delta=sqrt(1.-4.*eta), eta2;
89  eta2 = pow(eta,2);
90  double noSpin = 0.5566284518926176 + 0.12651770333481904*eta + 1.8084545267208734*eta2;
91  double eqSpin = (0.29074922226651545 + eta2*(-2.101111399437034 - 3.4969956644617946*S) + eta*(0.059317243606471406 - 0.31924748117518226*S) + 0.27420263462336675*S)*S;
92  double uneqSpin = 1.0122975748481835*pWF->dchi*eta2*delta;
93  total = noSpin + eqSpin + uneqSpin;
94  break;
95  }
96  case 122022:{
97  double eta = pWF->eta;
98  double S = pWF->chiPNHat;
99  double chidiff = pWF->dchi_half;
100  double eta1 = eta;
101  double eta2 = eta1 * eta1;
102  double eta3 = eta1 * eta2;
103  double eta4 = eta1 * eta3;
104  double S1 = S;
105  double chidiff1 = chidiff;
106  total = fabs(1.0092933052569775 - 0.2791855444800297*eta1 + 1.7110615047319937*eta2 + chidiff1*(-0.1054835719277311*eta1 + 7.506083919925026*eta2 - 30.1595680078279*eta3) + chidiff1*(2.078267611384239*eta1 - 10.166026002515457*eta2 - 1.2091616330625208*eta3)*S1 + S1*(0.17250873250247642*(1.0170226856985174 + 1.0395650952176598*eta1 - 35.73623734051525*eta2 + 403.68074286921444*eta3 - 1194.6152711219886*eta4) + 0.06850746964805364*(1.507796537056924 + 37.81075363806507*eta1 - 863.117144661059*eta2 + 6429.543634627373*eta3 - 15108.557419182316*eta4)*S1));
107  break;
108  }
109  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_21_lambda: version %i is not valid.", RDAmpFlag);}
110  }
111  return total;
112 }
113 
114 static double IMRPhenomXHM_RD_Amp_33_alambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag) {
115  double total=0;
116  switch (RDAmpFlag){
117  case 122018:{
118  double eta = pWF->eta;
119  double S = pWF->STotR;
120  double eta2,eta4,delta=sqrt(1-4*eta);
121  eta2 = pow(eta,2);
122  eta4 = pow(eta,4);
123  double noSpin = sqrt(eta - 4.*eta2)*(0.013700854227665184 + 0.01202732427321774*eta + 0.0898095508889557*eta2);
124  double eqSpin = sqrt(eta - 4.*eta2)*(0.0075858980586079065 + eta*(-0.013132320758494439 - 0.018186317026076343*S) + 0.0035617441651710473*S)*S;
125  double uneqSpin = eta4*(pWF->chi2L*(-0.09802218411554885 - 0.05745949361626237*S) + pWF->chi1L*(0.09802218411554885 + 0.05745949361626237*S) + eta2*(pWF->chi1L*(-4.2679864481479886 - 11.877399902871485*S) + pWF->chi2L*(4.2679864481479886 + 11.877399902871485*S))*delta);
126  total = noSpin + eqSpin + uneqSpin;
127  break;
128  }
129  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_33_alambda: version %i is not valid.", RDAmpFlag);}
130  }
131  return total;
132 }
133 
134 static double IMRPhenomXHM_RD_Amp_33_lambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag) {
135  double total=0;
136  switch (RDAmpFlag){
137  case 122018:{
138  double eta = pWF->eta;
139  double S = pWF->STotR;
140  double eta2,S2,delta=sqrt(1-4*eta);
141  eta2 = pow(eta,2);
142  S2 = pow(S,2);
143  double noSpin = 0.7435306475478924 - 0.06688558533374556*eta + 1.471989765837694*eta2;
144  double eqSpin = S*(0.19457194111990656 + 0.07564220573555203*S + eta*(-0.4809350398289311 + 0.17261430318577403*S - 0.1988991467974821*S2));
145  double uneqSpin = 1.8881959341735146*pWF->dchi*eta2*delta;
146  total = noSpin + eqSpin + uneqSpin;
147  break;
148  }
149  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_33_lambda: version %i is not valid.", RDAmpFlag);}
150  }
151  return total;
152 }
153 
154 static double IMRPhenomXHM_RD_Amp_32_alambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag) {
155  double total=0;
156  switch (RDAmpFlag){
157  case 122018:{
158  double eta = pWF->eta;
159  double S = pWF->STotR;
160  double eta2,eta3,eta4,eta5,eta6;
161  eta2 = pow(eta,2);
162  eta3 = pow(eta,3);
163  eta4 = pow(eta,4);
164  eta5 = pow(eta,5);
165  eta6 = pow(eta,6);
166  double noSpin = 0.00012587900257140724 + 0.03927886286971654*eta - 0.8109309606583066*eta2 + 8.820604907164254*eta3 - 51.43344812454074*eta4 + 141.81940900657446*eta5 - 140.0426973304466*eta6;
167  double eqSpin = S*(-0.00006001471234796344 + eta4*(-0.7849112300598181 - 2.09188976953315*S) + eta2*(0.08311497969032984 - 0.15569578955822236*S) + eta*(-0.01083175709906557 + 0.00568899459837252*S) - 0.00009363591928190229*S + 1.0670798489407887*eta3*S);
168  double uneqSpin = -0.04537308968659669*pow(pWF->dchi,2)*eta2*(1. - 8.711096029480697*eta + 18.362371966229926*eta2) + pWF->dchi*(-297.36978685672733 + 3103.2516759087644*eta - 10001.774055779177*eta2 + 9386.734883473799*eta3)*eta6;
169  total = noSpin + eqSpin + uneqSpin;
170  break;
171  }
172  case 122022:{
173  double eta = pWF->eta;
174  double delta = pWF->delta;
175  double S = pWF->STotR;
176  double chidiff = pWF->dchi_half;
177  double eta1 = eta;
178  double eta2 = eta1 * eta1;
179  double eta3 = eta1 * eta2;
180  double eta4 = eta1 * eta3;
181  double eta5 = eta1 * eta4;
182  double eta6 = eta1 * eta5;
183  double S1 = S;
184  double S2 = S1 * S1;
185  double chidiff1 = chidiff;
186  double chidiff2 = chidiff1 * chidiff1;
187  total = chidiff2*(-3.4614418482110163*eta3 + 35.464117772624164*eta4 - 85.19723511005235*eta5) + chidiff1*delta*(2.0328561081997463*eta3 - 46.18751757691501*eta4 + 170.9266105597438*eta5) + chidiff2*(-0.4600401291210382*eta3 + 12.23450117663151*eta4 - 42.74689906831975*eta5)*S1 + chidiff1*delta*(5.786292428422767*eta3 - 53.60467819078566*eta4 + 117.66195692191727*eta5)*S1 + S1*(-0.0013330716557843666*(56.35538385647113*eta1 - 1218.1550992423377*eta2 + 16509.69605686402*eta3 - 102969.88022112886*eta4 + 252228.94931931415*eta5 - 150504.2927996263*eta6) + 0.0010126460331462495*(-33.87083889060834*eta1 + 502.6221651850776*eta2 - 1304.9210590188136*eta3 - 36980.079328277505*eta4 + 295469.28617550555*eta5 - 597155.7619486618*eta6)*S1 - 0.00043088431510840695*(-30.014415072587354*eta1 - 1900.5495690280086*eta2 + 76517.21042363928*eta3 - 870035.1394696251*eta4 + 3.9072674134789007e6*eta5 - 6.094089675611567e6*eta6)*S2) + (0.08408469319155859*eta1 - 1.223794846617597*eta2 + 6.5972460654253515*eta3 - 15.707327897569396*eta4 + 14.163264397061505*eta5)*pow(1 - 8.612447115134758*eta1 + 18.93655612952139*eta2,-1);
188  break;
189  }
190  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_32_alambda: version %i is not valid.", RDAmpFlag);}
191  }
192  return total;
193 }
194 
195 static double IMRPhenomXHM_RD_Amp_32_lambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag) {
196  double total=0;
197  switch (RDAmpFlag){
198  case 122018:{
199  double eta = pWF->eta;
200  double S = pWF->STotR;
201  double eta2,eta3,delta=sqrt(1.-4*eta);
202  eta2 = pow(eta,2);
203  eta3 = pow(eta,3);
204  double noSpin = (sqrt(1. - 3.*eta)*(0.0341611244787871 - 0.3197209728114808*eta + 0.7689553234961991*eta2))/(0.048429644168112324 - 0.43758296068790314*eta + eta2);
205  double eqSpin = sqrt(1. - 3.*eta)*S*(0.11057199932233873 + eta2*(25.536336676250748 - 71.18182757443142*S) + 9.790509295728649*eta*S + eta3*(-56.96407763839491 + 175.47259563543165*S));
206  double uneqSpin = -5.002106168893265*pow(pWF->dchi,2)*eta2*delta;
207  total = noSpin + eqSpin + uneqSpin;
208  break;
209  }
210  case 122022:{
211  double eta = pWF->eta;
212  double delta = pWF->delta;
213  double S = pWF->STotR;
214  double chidiff = pWF->dchi_half;
215  double eta1 = eta;
216  double eta2 = eta1 * eta1;
217  double eta3 = eta1 * eta2;
218  double eta4 = eta1 * eta3;
219  double eta5 = eta1 * eta4;
220  double S1 = S;
221  double chidiff1 = chidiff;
222  double chidiff2 = chidiff1 * chidiff1;
223  total = 0.978510781593996 + 0.36457571743142897*eta1 - 12.259851752618998*eta2 + 49.19719473681921*eta3 + chidiff1*delta*(-188.37119473865533*eta3 + 2151.8731700399308*eta4 - 6328.182823770599*eta5) + chidiff2*(115.3689949926392*eta3 - 1159.8596972989067*eta4 + 2657.6998831179444*eta5) + S1*(0.22358643406992756*(0.48943645614341924 - 32.06682257944444*eta1 + 365.2485484044132*eta2 - 915.2489655397206*eta3) + 0.0792473022309144*(1.877251717679991 - 103.65639889587327*eta1 + 1202.174780792418*eta2 - 3206.340850767219*eta3)*S1);
224  break;
225  }
226  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_32_lambda: version %i is not valid.", RDAmpFlag);}
227  }
228  return total;
229 }
230 
231 static double IMRPhenomXHM_RD_Amp_44_alambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag) {
232  double total=0;
233  switch (RDAmpFlag){
234  case 122018:{
235  double eta = pWF->eta;
236  double S = pWF->STotR;
237  double delta=sqrt(1.-4.*eta),eta2,eta3,eta4,eta5,S2;
238  eta2 = pow(eta,2);
239  eta3 = pow(eta,3);
240  eta4 = pow(eta,4);
241  eta5 = pow(eta,5);
242  S2 = pow(S,2);
243  double noSpin = sqrt(eta - 3.*eta2)*(0.007904587819112173 + 0.09558474985614368*eta - 2.663803397359775*eta2 + 28.298192768381554*eta3 - 136.10446022757958*eta4 + 233.23167528016833*eta5);
244  double eqSpin = sqrt(eta - 3.*eta2)*S*(0.0049703757209330025 + 0.004122811292229324*S + eta*(-0.06166686913913691 + 0.014107365722576927*S)*S + eta2*(-0.2945455034809188 + 0.4139026619690879*S - 0.1389170612199015*S2) + eta3*(0.9225758392294605 - 0.9656098473922222*S + 0.19708289555425246*S2) + 0.000657528128497184*S2);
245  double uneqSpin = 0.00659873279539475*pWF->dchi*eta2*delta;
246  total = noSpin + eqSpin + uneqSpin;
247  break;
248  }
249  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_44_a: version %i is not valid.", RDAmpFlag);}
250  }
251  return total;
252 }
253 
254 static double IMRPhenomXHM_RD_Amp_44_lambda(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag) {
255  double total=0;
256  switch (RDAmpFlag){
257  case 122018:{
258  double eta = pWF->eta;
259  double S = pWF->STotR;
260  double delta=sqrt(1.-4.*eta),eta2,eta3,eta4,eta5,eta6,eta7;
261  eta2 = pow(eta,2);
262  eta3 = pow(eta,3);
263  eta4 = pow(eta,4);
264  eta5 = pow(eta,5);
265  eta6 = pow(eta,6);
266  eta7 = pow(eta,7);
267  double noSpin = 0.7702864948772887 + 32.81532373698395*eta - 1625.1795901450212*eta2 + 31305.458876573215*eta3 - 297375.5347399236*eta4 + 1.4726521941846698e6*eta5 - 3.616582470072637e6*eta6 + 3.4585865843680725e6*eta7;
268  double eqSpin = (-0.03011582009308575*S + 0.09034746027925727*eta*S + 1.8738784391649446*eta2*S - 5.621635317494836*eta3*S)/(-1.1340218677260014 + S);
269  double uneqSpin = 0.959943270591552*pow(pWF->dchi,2)*eta2 + 0.853573071529436*pWF->dchi*eta2*delta;
270  total = noSpin + eqSpin + uneqSpin;
271  break;
272  }
273  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_44_lambda: version %i is not valid.", RDAmpFlag);}
274  }
275  return total;
276 }
277 
278 static double IMRPhenomXHM_RD_Amp_21_sigma(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag) {
279  double total=0;
280  switch (RDAmpFlag){
281  case 122018:{
282  double eta = pWF->eta;
283  double S = pWF->STotR;
284  double delta=sqrt(1.-4.*eta),eta2;
285  eta2 = pow(eta,2);
286  double noSpin = 1.2922261617161441 + 0.0019318405961363861*eta;
287  double eqSpin = (0.04927982551108649 - 0.6703778360948937*eta + 2.6625014134659772*eta2)*S;
288  double uneqSpin = 1.2001101665670462*(pWF->chi1L + pow(pWF->chi1L,2) - 2.*pWF->chi1L*pWF->chi2L + (-1. + pWF->chi2L)*pWF->chi2L)*eta2*delta;
289  total = noSpin + eqSpin + uneqSpin;
290  break;
291  }
292  case 122022:{
293  double eta = pWF->eta;
294  double S = pWF->STotR;
295  double chidiff = pWF->dchi_half;
296  double eta1 = eta;
297  double eta2 = eta1 * eta1;
298  double eta3 = eta1 * eta2;
299  double S1 = S;
300  double chidiff1 = chidiff;
301  total = fabs(1.374451177213076 - 0.1147381625630186*eta1 + chidiff1*(0.6646459256372743*eta1 - 5.020585319906719*eta2 + 9.817281653770431*eta3) + chidiff1*(3.8734254747587973*eta1 - 39.880716190740465*eta2 + 99.05511583518896*eta3)*S1 + S1*(0.013272603498067647*(1.809972721953344 - 12.560287006325837*eta1 - 134.597005438578*eta2 + 786.2235720637008*eta3) + 0.006850483944311038*(-6.478737679813189 - 200.29813775611166*eta1 + 2744.3629484255357*eta2 - 7612.096007280672*eta3)*S1));
302  break;
303  }
304  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_21_sigma: version %i is not valid.", RDAmpFlag);}
305  }
306  return total;
307 }
308 
309 static double IMRPhenomXHM_RD_Amp_33_sigma(UNUSED IMRPhenomXWaveformStruct *pWF, int RDAmpFlag) {
310  double total=0;
311  switch (RDAmpFlag){
312  case 122018:{
313  double noSpin = 1.3;
314  double eqSpin = 0.;
315  double uneqSpin = 0.;
316  total = noSpin + eqSpin + uneqSpin;
317  break;
318  }
319  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_33_sigma: version %i is not valid.", RDAmpFlag);}
320  }
321  return total;
322 }
323 
324 static double IMRPhenomXHM_RD_Amp_32_sigma(UNUSED IMRPhenomXWaveformStruct *pWF, int RDAmpFlag) {
325  double total=0;
326  switch (RDAmpFlag){
327  case 122018:{
328  double noSpin = 1.33;
329  double eqSpin = 0.;
330  double uneqSpin = 0.;
331  total = noSpin + eqSpin + uneqSpin;
332  break;
333  }
334  case 122022:{
335  double eta = pWF->eta;
336  double delta = pWF->delta;
337  double S = pWF->STotR;
338  double chidiff = pWF->dchi_half;
339  double eta1 = eta;
340  double eta2 = eta1 * eta1;
341  double eta3 = eta1 * eta2;
342  double eta4 = eta1 * eta3;
343  double eta5 = eta1 * eta4;
344  double S1 = S;
345  double chidiff1 = chidiff;
346  double chidiff2 = chidiff1 * chidiff1;
347  total = 1.3353917551819414 + 0.13401718687342024*eta1 + chidiff1*delta*(144.37065005786636*eta3 - 754.4085447486738*eta4 + 123.86194078913776*eta5) + chidiff2*(209.09202210427972*eta3 - 1769.4658099037918*eta4 + 3592.287297392387*eta5) + S1*(-0.012086025709597246*(-6.230497473791485 + 600.5968613752918*eta1 - 6606.1009717965735*eta2 + 17277.60594350428*eta3) - 0.06066548829900489*(-0.9208054306316676 + 142.0346574366267*eta1 - 1567.249168668069*eta2 + 4119.373703246675*eta3)*S1);
348  break;
349  }
350  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_32_sigma: version %i is not valid.", RDAmpFlag);}
351  }
352  return total;
353 }
354 
355 static double IMRPhenomXHM_RD_Amp_44_sigma(UNUSED IMRPhenomXWaveformStruct *pWF, int RDAmpFlag) {
356  double total=0;
357  switch (RDAmpFlag){
358  case 122018:{
359  double noSpin = 1.33;
360  double eqSpin = 0.;
361  double uneqSpin = 0.;
362  total = noSpin + eqSpin + uneqSpin;
363  break;
364  }
365  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_44_sigma: version %i is not valid.", RDAmpFlag);}
366  }
367  return total;
368 }
369 
370 
371 static double IMRPhenomXHM_RD_Amp_21_rdcp1(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
372  double total=0;
373  switch (RDAmpFlag){
374  case 122022:{
375  double eta = pWF->eta;
376  double delta = pWF->delta;
377  double S = pWF->chiPNHat;
378  double chidiff = pWF->dchi_half;
379  double eta1 = eta;
380  double eta2 = eta1 * eta1;
381  double eta3 = eta1 * eta2;
382  double eta4 = eta1 * eta3;
383  double eta5 = eta1 * eta4;
384  double eta6 = eta1 * eta5;
385  double S1 = S;
386  double S2 = S1 * S1;
387  double chidiff1 = chidiff;
388  total = fabs(delta*eta1*(12.880905080761432 - 23.5291063016996*eta1 + 92.6090002736012*eta2 - 175.16681482428694*eta3) + chidiff1*delta*eta1*(26.89427230731867*eta1 - 710.8871223808559*eta2 + 2255.040486907459*eta3) + chidiff1*delta*eta1*(21.402708785047853*eta1 - 232.07306353130417*eta2 + 591.1097623278739*eta3)*S1 + delta*eta1*S1*(-10.090867481062709*(0.9580746052260011 + 5.388149112485179*eta1 - 107.22993216128548*eta2 + 801.3948756800821*eta3 - 2688.211889175019*eta4 + 3950.7894052628735*eta5 - 1992.9074348833092*eta6) - 0.42972412296628143*(1.9193131231064235 + 139.73149069609775*eta1 - 1616.9974609915555*eta2 - 3176.4950303461164*eta3 + 107980.65459735804*eta4 - 479649.75188253267*eta5 + 658866.0983367155*eta6)*S1) + chidiff1*eta5*(-1512.439342647443 + 175.59081294852444*S1 + 10.13490934572329*S2));
389  break;
390  }
391  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_21_rdcp1: version %i is not valid.", RDAmpFlag);}
392  }
393  return total;
394 }
395 
396 static double IMRPhenomXHM_RD_Amp_21_rdcp2(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
397  double total=0;
398  switch (RDAmpFlag){
399  case 122022:{
400  double eta = pWF->eta;
401  double delta = pWF->delta;
402  double S = pWF->chiPNHat;
403  double chidiff = pWF->dchi_half;
404  double eta1 = eta;
405  double eta2 = eta1 * eta1;
406  double eta3 = eta1 * eta2;
407  double eta4 = eta1 * eta3;
408  double eta5 = eta1 * eta4;
409  double S1 = S;
410  double S2 = S1 * S1;
411  double chidiff1 = chidiff;
412  total = fabs(delta*(9.112452928978168 - 7.5304766811877455*eta1)*eta1 + chidiff1*delta*eta1*(16.236533863306132*eta1 - 500.11964987628926*eta2 + 1618.0818430353293*eta3) + chidiff1*delta*eta1*(2.7866868976718226*eta1 - 0.4210629980868266*eta2 - 20.274691328125606*eta3)*S1 + chidiff1*eta5*(-1116.4039232324135 + 245.73200219767514*S1 + 21.159179960295855*S2) + delta*eta1*S1*(-8.236485576091717*(0.8917610178208336 + 5.1501231412520285*eta1 - 87.05136337926156*eta2 + 519.0146702141192*eta3 - 997.6961311502365*eta4) + 0.2836840678615208*(-0.19281297100324718 - 57.65586769647737*eta1 + 586.7942442434971*eta2 - 1882.2040277496196*eta3 + 2330.3534917059906*eta4)*S1 + 0.40226131643223145*(-3.834742668014861 + 190.42214703482531*eta1 - 2885.5110686004946*eta2 + 16087.433824017446*eta3 - 29331.524552164105*eta4)*S2));
413  break;
414  }
415  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_21_rdcp2: version %i is not valid.", RDAmpFlag);}
416  }
417  return total;
418 }
419 
420 static double IMRPhenomXHM_RD_Amp_21_rdcp3(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
421  double total=0;
422  switch (RDAmpFlag){
423  case 122022:{
424  double eta = pWF->eta;
425  double delta = pWF->delta;
426  double S = pWF->chiPNHat;
427  double chidiff = pWF->dchi_half;
428  double eta1 = eta;
429  double eta2 = eta1 * eta1;
430  double eta3 = eta1 * eta2;
431  double eta4 = eta1 * eta3;
432  double eta5 = eta1 * eta4;
433  double S1 = S;
434  double S2 = S1 * S1;
435  double chidiff1 = chidiff;
436  total = fabs(delta*(2.920930733198033 - 3.038523690239521*eta1)*eta1 + chidiff1*delta*eta1*(6.3472251472354975*eta1 - 171.23657247338042*eta2 + 544.1978232314333*eta3) + chidiff1*delta*eta1*(1.9701247529688362*eta1 - 2.8616711550845575*eta2 - 0.7347258030219584*eta3)*S1 + chidiff1*eta5*(-334.0969956136684 + 92.91301644484749*S1 - 5.353399481074393*S2) + delta*eta1*S1*(-2.7294297839371824*(1.148166706456899 - 4.384077347340523*eta1 + 36.120093043420326*eta2 - 87.26454353763077*eta3) + 0.23949142867803436*(-0.6931516433988293 + 33.33372867559165*eta1 - 307.3404155231787*eta2 + 862.3123076782916*eta3)*S1 + 0.1930861073906724*(3.7735099269174106 - 19.11543562444476*eta1 - 78.07256429516346*eta2 + 485.67801863289293*eta3)*S2));
437  break;
438  }
439  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_21_rdcp3: version %i is not valid.", RDAmpFlag);}
440  }
441  return total;
442 }
443 
444 static double IMRPhenomXHM_RD_Amp_33_rdcp1(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
445  double total=0;
446  switch (RDAmpFlag){
447  case 122022:{
448  double eta = pWF->eta;
449  double delta = pWF->delta;
450  double S = pWF->STotR;
451  double chidiff = pWF->dchi_half;
452  double eta1 = eta;
453  double eta2 = eta1 * eta1;
454  double eta3 = eta1 * eta2;
455  double eta4 = eta1 * eta3;
456  double eta5 = eta1 * eta4;
457  double S1 = S;
458  double S2 = S1 * S1;
459  double chidiff1 = chidiff;
460  double chidiff2 = chidiff1 * chidiff1;
461  total = delta*eta1*(12.439702602599235 - 4.436329538596615*eta1 + 22.780673360839497*eta2) + delta*eta1*(chidiff1*(-41.04442169938298*eta1 + 502.9246970179746*eta2 - 1524.2981907688634*eta3) + chidiff2*(32.23960072974939*eta1 - 365.1526474476759*eta2 + 1020.6734178547847*eta3)) + chidiff1*delta*eta1*(-52.85961155799673*eta1 + 577.6347407795782*eta2 - 1653.496174539196*eta3)*S1 + chidiff1*eta5*(257.33227387984863 - 34.5074027042393*chidiff2 - 21.836905132600755*S1 - 15.81624534976308*S2) + 13.499999999999998*delta*eta1*S1*(-0.13654149379906394*(2.719687834084113 + 29.023992126142304*eta1 - 742.1357702210267*eta2 + 4142.974510926698*eta3 - 6167.08766058184*eta4 - 3591.1757995710486*eta5) - 0.06248535354306988*(6.697567446351289 - 78.23231700361792*eta1 + 444.79350113344543*eta2 - 1907.008984765889*eta3 + 6601.918552659412*eta4 - 10056.98422430965*eta5)*S1)*pow(-3.9329308614837704 + S1,-1);
462  break;
463  }
464  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_33_rdcp1: version %i is not valid.", RDAmpFlag);}
465  }
466  return total;
467 }
468 
469 static double IMRPhenomXHM_RD_Amp_33_rdcp2(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
470  double total=0;
471  switch (RDAmpFlag){
472  case 122022:{
473  double eta = pWF->eta;
474  double delta = pWF->delta;
475  double S = pWF->STotR;
476  double chidiff = pWF->dchi_half;
477  double eta1 = eta;
478  double eta2 = eta1 * eta1;
479  double eta3 = eta1 * eta2;
480  double eta4 = eta1 * eta3;
481  double eta5 = eta1 * eta4;
482  double S1 = S;
483  double S2 = S1 * S1;
484  double chidiff1 = chidiff;
485  double chidiff2 = chidiff1 * chidiff1;
486  total = delta*eta1*(8.425057692276933 + 4.543696144846763*eta1) + chidiff1*delta*eta1*(-32.18860840414171*eta1 + 412.07321398189293*eta2 - 1293.422289802462*eta3) + chidiff1*delta*eta1*(-17.18006888428382*eta1 + 190.73514518113845*eta2 - 636.4802385540647*eta3)*S1 + delta*eta1*S1*(0.1206817303851239*(8.667503604073314 - 144.08062755162752*eta1 + 3188.189172446398*eta2 - 35378.156133055556*eta3 + 163644.2192178668*eta4 - 265581.70142471837*eta5) + 0.08028332044013944*(12.632478544060636 - 322.95832000179297*eta1 + 4777.45310151897*eta2 - 35625.58409457366*eta3 + 121293.97832549023*eta4 - 148782.33687815256*eta5)*S1) + chidiff1*eta5*(159.72371180117415 - 29.10412708633528*chidiff2 - 1.873799747678187*S1 + 41.321480132899524*S2);
487  break;
488  }
489  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_33_rdcp2: version %i is not valid.", RDAmpFlag);}
490  }
491  return total;
492 }
493 
494 static double IMRPhenomXHM_RD_Amp_33_rdcp3(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
495  double total=0;
496  switch (RDAmpFlag){
497  case 122022:{
498  double eta = pWF->eta;
499  double delta = pWF->delta;
500  double S = pWF->STotR;
501  double chidiff = pWF->dchi_half;
502  double eta1 = eta;
503  double eta2 = eta1 * eta1;
504  double eta3 = eta1 * eta2;
505  double eta4 = eta1 * eta3;
506  double eta5 = eta1 * eta4;
507  double S1 = S;
508  double S2 = S1 * S1;
509  double chidiff1 = chidiff;
510  double chidiff2 = chidiff1 * chidiff1;
511  total = delta*eta1*(2.485784720088995 + 2.321696430921996*eta1) + delta*eta1*(chidiff1*(-10.454376404653859*eta1 + 147.10344302665484*eta2 - 496.1564538739011*eta3) + chidiff2*(-5.9236399792925996*eta1 + 65.86115501723127*eta2 - 197.51205149250532*eta3)) + chidiff1*delta*eta1*(-10.27418232676514*eta1 + 136.5150165348149*eta2 - 473.30988537734174*eta3)*S1 + chidiff1*eta5*(32.07819766300362 - 3.071422453072518*chidiff2 + 35.09131921815571*S1 + 67.23189816732847*S2) + 13.499999999999998*delta*eta1*S1*(0.0011484326782460882*(4.1815722950796035 - 172.58816646768219*eta1 + 5709.239330076732*eta2 - 67368.27397765424*eta3 + 316864.0589150127*eta4 - 517034.11171277676*eta5) - 0.009496797093329243*(0.9233282181397624 - 118.35865186626413*eta1 + 2628.6024206791726*eta2 - 23464.64953722729*eta3 + 94309.57566199072*eta4 - 140089.40725211444*eta5)*S1)*pow(0.09549360183532198 - 0.41099904730526465*S1 + S2,-1);
512  break;
513  }
514  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_33_rdcp3: version %i is not valid.", RDAmpFlag);}
515  }
516  return total;
517 }
518 
519 static double IMRPhenomXHM_RD_Amp_44_rdcp1(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
520  double total=0;
521  switch (RDAmpFlag){
522  case 122022:{
523  double eta = pWF->eta;
524  double delta = pWF->delta;
525  double S = pWF->STotR;
526  double chidiff = pWF->dchi_half;
527  double eta1 = eta;
528  double eta2 = eta1 * eta1;
529  double eta3 = eta1 * eta2;
530  double eta4 = eta1 * eta3;
531  double eta5 = eta1 * eta4;
532  double eta6 = eta1 * eta5;
533  double S1 = S;
534  double S2 = S1 * S1;
535  double chidiff1 = chidiff;
536  double chidiff2 = chidiff1 * chidiff1;
537  total = eta1*(chidiff1*delta*(-8.51952446214978*eta1 + 117.76530248141987*eta2 - 297.2592736781142*eta3) + chidiff2*(-0.2750098647982238*eta1 + 4.456900599347149*eta2 - 8.017569928870929*eta3)) + eta1*(5.635069974807398 - 33.67252878543393*eta1 + 287.9418482197136*eta2 - 3514.3385364216438*eta3 + 25108.811524802128*eta4 - 98374.18361532023*eta5 + 158292.58792484726*eta6) + eta1*S1*(-0.4360849737360132*(-0.9543114627170375 - 58.70494649755802*eta1 + 1729.1839588870455*eta2 - 16718.425586396803*eta3 + 71236.86532610047*eta4 - 111910.71267453219*eta5) - 0.024861802943501172*(-52.25045490410733 + 1585.462602954658*eta1 - 15866.093368857853*eta2 + 35332.328181283*eta3 + 168937.32229060197*eta4 - 581776.5303770923*eta5)*S1 + 0.005856387555754387*(186.39698091707513 - 9560.410655118145*eta1 + 156431.3764198244*eta2 - 1.0461268207440731e6*eta3 + 3.054333578686424e6*eta4 - 3.2369858387064277e6*eta5)*S2);
538  break;
539  }
540  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_44_rdcp1: version %i is not valid.", RDAmpFlag);}
541  }
542  return total;
543 }
544 
545 static double IMRPhenomXHM_RD_Amp_44_rdcp2(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
546  double total=0;
547  switch (RDAmpFlag){
548  case 122022:{
549  double eta = pWF->eta;
550  double delta = pWF->delta;
551  double S = pWF->STotR;
552  double chidiff = pWF->dchi_half;
553  double eta1 = eta;
554  double eta2 = eta1 * eta1;
555  double eta3 = eta1 * eta2;
556  double eta4 = eta1 * eta3;
557  double eta5 = eta1 * eta4;
558  double S1 = S;
559  double S2 = S1 * S1;
560  double chidiff1 = chidiff;
561  double chidiff2 = chidiff1 * chidiff1;
562  total = eta1*(chidiff1*delta*(-2.861653255976984*eta1 + 50.50227103211222*eta2 - 123.94152825700999*eta3) + chidiff2*(2.9415751419018865*eta1 - 28.79779545444817*eta2 + 72.40230240887851*eta3)) + eta1*(3.2461722686239307 + 25.15310593958783*eta1 - 792.0167314124681*eta2 + 7168.843978909433*eta3 - 30595.4993786313*eta4 + 49148.57065911245*eta5) + eta1*S1*(-0.23311779185707152*(-1.0795711755430002 - 20.12558747513885*eta1 + 1163.9107546486134*eta2 - 14672.23221502075*eta3 + 73397.72190288734*eta4 - 127148.27131388368*eta5) + 0.025805905356653*(11.929946153728276 + 350.93274421955806*eta1 - 14580.02701600596*eta2 + 174164.91607515427*eta3 - 819148.9390278616*eta4 + 1.3238624538095295e6*eta5)*S1 + 0.019740635678180102*(-7.046295936301379 + 1535.781942095697*eta1 - 27212.67022616794*eta2 + 201981.0743810629*eta3 - 696891.1349708183*eta4 + 910729.0219043035*eta5)*S2);
563  break;
564  }
565  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_44_rdcp2: version %i is not valid.", RDAmpFlag);}
566  }
567  return total;
568 }
569 
570 static double IMRPhenomXHM_RD_Amp_44_rdcp3(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
571  double total=0;
572  switch (RDAmpFlag){
573  case 122022:{
574  double eta = pWF->eta;
575  double delta = pWF->delta;
576  double S = pWF->STotR;
577  double chidiff = pWF->dchi_half;
578  double eta1 = eta;
579  double eta2 = eta1 * eta1;
580  double eta3 = eta1 * eta2;
581  double eta4 = eta1 * eta3;
582  double eta5 = eta1 * eta4;
583  double S1 = S;
584  double chidiff1 = chidiff;
585  double chidiff2 = chidiff1 * chidiff1;
586  total = eta1*(chidiff1*delta*(2.4286414692113816*eta1 - 23.213332913737403*eta2 + 66.58241012629095*eta3) + chidiff2*(3.085167288859442*eta1 - 31.60440418701438*eta2 + 78.49621016381445*eta3)) + eta1*(0.861883217178703 + 13.695204704208976*eta1 - 337.70598252897696*eta2 + 2932.3415281149432*eta3 - 12028.786386004691*eta4 + 18536.937955014455*eta5) + eta1*S1*(-0.048465588779596405*(-0.34041762314288154 - 81.33156665674845*eta1 + 1744.329802302927*eta2 - 16522.343895064576*eta3 + 76620.18243090731*eta4 - 133340.93723954144*eta5) + 0.024804027856323612*(-8.666095805675418 + 711.8727878341302*eta1 - 13644.988225595187*eta2 + 112832.04975245205*eta3 - 422282.0368440555*eta4 + 584744.0406581408*eta5)*S1);
587  break;
588  }
589  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_44_rdcp3: version %i is not valid.", RDAmpFlag);}
590  }
591  return total;
592 }
593 
594 static double IMRPhenomXHM_RD_Amp_32_rdaux1(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
595  double total=0;
596  switch (RDAmpFlag){
597  case 122022:{
598  double eta = pWF->eta;
599  double delta = pWF->delta;
600  double S = pWF->STotR;
601  double chidiff = pWF->dchi_half;
602  double eta1 = eta;
603  double eta2 = eta1 * eta1;
604  double eta3 = eta1 * eta2;
605  double eta4 = eta1 * eta3;
606  double eta5 = eta1 * eta4;
607  double eta6 = eta1 * eta5;
608  double S1 = S;
609  double chidiff1 = chidiff;
610  double chidiff2 = chidiff1 * chidiff1;
611  total = chidiff2*(-4.188795724777721*eta2 + 53.39200466700963*eta3 - 131.19660856923554*eta4) + chidiff1*delta*(14.284921364132623*eta2 - 321.26423637658746*eta3 + 1242.865584938088*eta4) + S1*(-0.022968727462555794*(83.66854837403105*eta1 - 3330.6261333413177*eta2 + 77424.12614733395*eta3 - 710313.3016672594*eta4 + 2.6934917075009225e6*eta5 - 3.572465179268999e6*eta6) + 0.0014795114305436387*(-1672.7273629876313*eta1 + 90877.38260964208*eta2 - 1.6690169155105734e6*eta3 + 1.3705532554135624e7*eta4 - 5.116110998398143e7*eta5 + 7.06066766311127e7*eta6)*S1) + (4.45156488896258*eta1 - 77.39303992494544*eta2 + 522.5070635563092*eta3 - 1642.3057499049708*eta4 + 2048.333892310575*eta5)*pow(1 - 9.611489164758915*eta1 + 24.249594730050312*eta2,-1);
612  break;
613  }
614  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_32_rdaux1: version %i is not valid.", RDAmpFlag);}
615  }
616  return total;
617 }
618 
619 static double IMRPhenomXHM_RD_Amp_32_rdaux2(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
620  double total=0;
621  switch (RDAmpFlag){
622  case 122022:{
623  double eta = pWF->eta;
624  double delta = pWF->delta;
625  double S = pWF->chiPNHat;
626  double chidiff = pWF->dchi_half;
627  double eta1 = eta;
628  double eta2 = eta1 * eta1;
629  double eta3 = eta1 * eta2;
630  double eta4 = eta1 * eta3;
631  double eta5 = eta1 * eta4;
632  double eta6 = eta1 * eta5;
633  double eta7 = eta1 * eta6;
634  double S1 = S;
635  double chidiff1 = chidiff;
636  double chidiff2 = chidiff1 * chidiff1;
637  total = chidiff2*(-18.550171209458394*eta2 + 188.99161055445936*eta3 - 440.26516625611*eta4) + chidiff1*delta*(13.132625215315063*eta2 - 340.5204040505528*eta3 + 1327.1224176812448*eta4) + S1*(-0.16707403272774676*(6.678916447469937*eta1 + 1331.480396625797*eta2 - 41908.45179140144*eta3 + 520786.0225074669*eta4 - 3.1894624909922685e6*eta5 + 9.51553823212259e6*eta6 - 1.1006903622406831e7*eta7) + 0.015205286051218441*(108.10032279461095*eta1 - 16084.215590200103*eta2 + 462957.5593513407*eta3 - 5.635028227588545e6*eta4 + 3.379925277713386e7*eta5 - 9.865815275452062e7*eta6 + 1.1201307979786257e8*eta7)*S1) + (3.902154247490771*eta1 - 55.77521071924907*eta2 + 294.9496843041973*eta3 - 693.6803787318279*eta4 + 636.0141528226893*eta5)*pow(1 - 8.56699762573719*eta1 + 19.119341007236955*eta2,-1);
638  break;
639  }
640  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_32_rdaux2: version %i is not valid.", RDAmpFlag);}
641  }
642  return total;
643 }
644 
645 static double IMRPhenomXHM_RD_Amp_32_rdcp1(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
646  double total=0;
647  switch (RDAmpFlag){
648  case 122022:{
649  double eta = pWF->eta;
650  double delta = pWF->delta;
651  double S = pWF->STotR;
652  double chidiff = pWF->dchi_half;
653  double eta1 = eta;
654  double eta2 = eta1 * eta1;
655  double eta3 = eta1 * eta2;
656  double eta4 = eta1 * eta3;
657  double eta5 = eta1 * eta4;
658  double eta6 = eta1 * eta5;
659  double S1 = S;
660  double S2 = S1 * S1;
661  double chidiff1 = chidiff;
662  double chidiff2 = chidiff1 * chidiff1;
663  total = chidiff2*(-261.63903838092017*eta3 + 2482.4929818200458*eta4 - 5662.765952006266*eta5) + chidiff1*delta*(200.3023530582654*eta3 - 3383.07742098347*eta4 + 11417.842708417566*eta5) + chidiff2*(-177.2481070662751*eta3 + 1820.8637746828358*eta4 - 4448.151940319403*eta5)*S1 + chidiff1*delta*(412.749304734278*eta3 - 4156.641392955615*eta4 + 10116.974216563232*eta5)*S1 + S1*(-0.07383539239633188*(40.59996146686051*eta1 - 527.5322650311067*eta2 + 4167.108061823492*eta3 - 13288.883172763119*eta4 - 23800.671572828596*eta5 + 146181.8016013141*eta6) + 0.03576631753501686*(-13.96758180764024*eta1 - 797.1235306450683*eta2 + 18007.56663810595*eta3 - 151803.40642097822*eta4 + 593811.4596071478*eta5 - 878123.747877138*eta6)*S1 + 0.01007493097350273*(-27.77590078264459*eta1 + 4011.1960424049857*eta2 - 152384.01804465035*eta3 + 1.7595145936445233e6*eta4 - 7.889230647117076e6*eta5 + 1.2172078072446395e7*eta6)*S2) + (4.146029818148087*eta1 - 61.060972560568054*eta2 + 336.3725848841942*eta3 - 832.785332776221*eta4 + 802.5027431944313*eta5)*pow(1 - 8.662174796705683*eta1 + 19.288918757536685*eta2,-1);
664  break;
665  }
666  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_32_rdcp1: version is not valid. Recommended version is 122022.");}
667  }
668  return total;
669 }
670 
671 static double IMRPhenomXHM_RD_Amp_32_rdcp2(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
672  double total=0;
673  switch (RDAmpFlag){
674  case 122022:{
675  double eta = pWF->eta;
676  double delta = pWF->delta;
677  double S = pWF->STotR;
678  double chidiff = pWF->dchi_half;
679  double eta1 = eta;
680  double eta2 = eta1 * eta1;
681  double eta3 = eta1 * eta2;
682  double eta4 = eta1 * eta3;
683  double eta5 = eta1 * eta4;
684  double eta6 = eta1 * eta5;
685  double S1 = S;
686  double S2 = S1 * S1;
687  double chidiff1 = chidiff;
688  double chidiff2 = chidiff1 * chidiff1;
689  total = chidiff2*(-220.42133216774002*eta3 + 2082.031407555522*eta4 - 4739.292554291661*eta5) + chidiff1*delta*(179.07548162694007*eta3 - 2878.2078963030094*eta4 + 9497.998559135678*eta5) + chidiff2*(-128.07917402087625*eta3 + 1392.4598433465628*eta4 - 3546.2644951338134*eta5)*S1 + chidiff1*delta*(384.31792882093424*eta3 - 3816.5687272960417*eta4 + 9235.479593415908*eta5)*S1 + S1*(-0.06144774696295017*(35.72693522898656*eta1 - 168.08433700852038*eta2 - 3010.678442066521*eta3 + 45110.034521934074*eta4 - 231569.4154711447*eta5 + 414234.84895584086*eta6) + 0.03663881822701642*(-22.057692852225696*eta1 + 223.9912685075838*eta2 - 1028.5261783449762*eta3 - 12761.957255385*eta4 + 141784.13567610556*eta5 - 328718.5349981628*eta6)*S1 + 0.004849853669413881*(-90.35491669965123*eta1 + 19286.158446325957*eta2 - 528138.5557827373*eta3 + 5.175061086459432e6*eta4 - 2.1142182400264673e7*eta5 + 3.0737963347449116e7*eta6)*S2) + (3.133378729082171*eta1 - 45.83572706555282*eta2 + 250.23275606463622*eta3 - 612.0498767005383*eta4 + 580.3574091493459*eta5)*pow(1 - 8.698032720488515*eta1 + 19.38621948411302*eta2,-1);
690  break;
691  }
692  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_32_rdcp2: version is not valid. Recommended version is 122022.");}
693  }
694  return total;
695 }
696 
697 static double IMRPhenomXHM_RD_Amp_32_rdcp3(IMRPhenomXWaveformStruct *pWF, int RDAmpFlag){
698  double total=0;
699  switch (RDAmpFlag){
700  case 122022:{
701  double eta = pWF->eta;
702  double delta = pWF->delta;
703  double S = pWF->STotR;
704  double chidiff = pWF->dchi_half;
705  double eta1 = eta;
706  double eta2 = eta1 * eta1;
707  double eta3 = eta1 * eta2;
708  double eta4 = eta1 * eta3;
709  double eta5 = eta1 * eta4;
710  double eta6 = eta1 * eta5;
711  double S1 = S;
712  double chidiff1 = chidiff;
713  double chidiff2 = chidiff1 * chidiff1;
714  total = chidiff2*(-79.14146757219045*eta3 + 748.8207876524461*eta4 - 1712.3401586150026*eta5) + chidiff1*delta*(65.1786095079065*eta3 - 996.4553252426255*eta4 + 3206.5675278160684*eta5) + chidiff2*(-36.474455088940225*eta3 + 421.8842792746865*eta4 - 1117.0227933265749*eta5)*S1 + chidiff1*delta*(169.07368933925878*eta3 - 1675.2562326502878*eta4 + 4040.0077967763787*eta5)*S1 + S1*(-0.01992370601225598*(36.307098892574196*eta1 - 846.997262853445*eta2 + 16033.60939445582*eta3 - 138800.53021166887*eta4 + 507922.88946543116*eta5 - 647376.1499824544*eta6) + 0.014207919520826501*(-33.80287899746716*eta1 + 1662.2913368534057*eta2 - 31688.885017467597*eta3 + 242813.43893659746*eta4 - 793178.4767168422*eta5 + 929016.897093022*eta6)*S1) + (0.9641853854287679*eta1 - 13.801372413989519*eta2 + 72.80610853168994*eta3 - 168.65551450831953*eta4 + 147.2372582604103*eta5)*pow(1 - 8.65963828355163*eta1 + 19.112920222001367*eta2,-1);
715  break;
716  }
717  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Amp_32_rdcp3: version is not valid. Recommended version is 122022.");}
718  }
719  return total;
720 }
721 
722 /* End of Amp Parameter Space Fits */
723 
724 
725 /************** Ringdown coefficients from collocation points *************/
726 
728  switch (pWFHM->IMRPhenomXHMRingdownAmpVersion){
729  case 0:{
730  // We have three "fitted" coefficients across parameter space: alambda, lambda and sigma. Sigma will be constat for all the modes except the 21.
731  pAmp->RDCoefficient[0] = fabs(pAmp->RingdownAmpFits[pWFHM->modeInt*3](pWF22,pWFHM->IMRPhenomXHMRingdownAmpFitsVersion));
732  pAmp->RDCoefficient[1] = pAmp->RingdownAmpFits[pWFHM->modeInt*3+1](pWF22,pWFHM->IMRPhenomXHMRingdownAmpFitsVersion);
733  pAmp->RDCoefficient[2] = pAmp->RingdownAmpFits[pWFHM->modeInt*3+2](pWF22,pWFHM->IMRPhenomXHMRingdownAmpFitsVersion);
734  pAmp->RDCoefficient[3] = 1./12.;
735  break;
736  }
737  case 1:
738  case 2:{
739  if (pWFHM->IMRPhenomXHMRingdownAmpVersion == 1){
740  pAmp->RDCoefficient[0] = fabs(pAmp->RingdownAmpFits[pWFHM->modeInt*3](pWF22,pWFHM->IMRPhenomXHMRingdownAmpFitsVersion));
741  pAmp->RDCoefficient[1] = pAmp->RingdownAmpFits[pWFHM->modeInt*3+1](pWF22,pWFHM->IMRPhenomXHMRingdownAmpFitsVersion);
742  pAmp->RDCoefficient[2] = pAmp->RingdownAmpFits[pWFHM->modeInt*3+2](pWF22,pWFHM->IMRPhenomXHMRingdownAmpFitsVersion);
743  }
744  else if (pWFHM->IMRPhenomXHMRingdownAmpVersion == 2){
745  double rdcp1 = fabs(pAmp->RingdownAmpFits[12 + pWFHM->modeInt * 3](pWF22, pWFHM->IMRPhenomXHMRingdownAmpFitsVersion));
746  double rdcp2 = fabs(pAmp->RingdownAmpFits[13 + pWFHM->modeInt * 3](pWF22, pWFHM->IMRPhenomXHMRingdownAmpFitsVersion));
747  double rdcp3 = fabs(pAmp->RingdownAmpFits[14 + pWFHM->modeInt * 3](pWF22, pWFHM->IMRPhenomXHMRingdownAmpFitsVersion));
748 
749  pAmp->CollocationPointsFreqsAmplitudeRD[0] = pWFHM->fRING - pWFHM->fDAMP;
750  pAmp->CollocationPointsFreqsAmplitudeRD[1] = pWFHM->fRING;
751  pAmp->CollocationPointsFreqsAmplitudeRD[2] = pWFHM->fRING + pWFHM->fDAMP;
752  /* Apply vetos to RDCP. Assuming they are strain */
753  // float rdveto = 0.01; // Applied over the strain / RescaleFactor_lm
754  // IMRPhenomX_UsefulPowers powers_of_RDCP1, powers_of_RDCP2, powers_of_RDCP3; // PN power v^{1/3} = (2pif/m)
755  // IMRPhenomX_Initialize_Powers(&powers_of_RDCP1, pAmp->CollocationPointsFreqsAmplitudeRD[0]);
756  // IMRPhenomX_Initialize_Powers(&powers_of_RDCP2, pAmp->CollocationPointsFreqsAmplitudeRD[1]);
757  // IMRPhenomX_Initialize_Powers(&powers_of_RDCP3, pAmp->CollocationPointsFreqsAmplitudeRD[2]);
758  // double rescale_factor_lm;
759  // rescale_factor_lm = RescaleFactor(&powers_of_RDCP1, pAmp, 2);
760  // if ( rdcp1 / rescale_factor_lm < rdveto ){
761  // rdcp1 = 0.9 * rdcp2;
762  // }
763  // rescale_factor_lm = RescaleFactor(&powers_of_RDCP2, pAmp, 2);
764  // if ( rdcp2 / rescale_factor_lm < rdveto ){
765  // rdcp2 = 0.9 * rdcp1;
766  // }
767  // rescale_factor_lm = RescaleFactor(&powers_of_RDCP3, pAmp, 2);
768  // if ( rdcp3 / rescale_factor_lm < rdveto ){
769  // rdcp3 = 0.9 * rdcp2;
770  // }
771  if ( rdcp3 >= rdcp2 * rdcp2 / rdcp1 ){
772  rdcp3 = 0.5 * rdcp2 * rdcp2 / rdcp1;
773  }
774  if ( rdcp3 > rdcp2 ){
775  rdcp3 = 0.5 * rdcp2;
776  }
777  if (rdcp1 < rdcp2 && rdcp3 > rdcp1){
778  rdcp3 = rdcp1;
779  }
780  /* End of vetos */
781  pAmp->CollocationPointsValuesAmplitudeRD[0] = rdcp1;
782  pAmp->CollocationPointsValuesAmplitudeRD[1] = rdcp2;
783  pAmp->CollocationPointsValuesAmplitudeRD[2] = rdcp3;
784  REAL8 deno = (sqrt(rdcp1 / rdcp3) - (rdcp1 / rdcp2));
785  if (deno <= 0){
786  deno = 1e-16;
787  }
788  pAmp->RDCoefficient[0] = rdcp1 * pWFHM->fDAMP / deno;
789  pAmp->RDCoefficient[2] = sqrt(pAmp->RDCoefficient[0] / (rdcp2 * pWFHM->fDAMP));
790  pAmp->RDCoefficient[1] = 0.5 * pAmp->RDCoefficient[2] * log(rdcp1 / rdcp3);
791  }
792 
793  if (pWFHM->RingdownAmpVeto == 1){
794  pAmp->RDCoefficient[1] = 1;
795  pAmp->RDCoefficient[2] = 1.35;
796  pAmp->RDCoefficient[0] = pAmp->CollocationPointsValuesAmplitudeRD[2] * pWFHM->fDAMP * pAmp->RDCoefficient[2] * pAmp->RDCoefficient[2];
797  }
798 
799  if(pWFHM->fAmpRDfalloff > 0){
800  IMRPhenomX_UsefulPowers powers_of_RDfalloff;
801  IMRPhenomX_Initialize_Powers(&powers_of_RDfalloff, pWFHM->fAmpRDfalloff);
802  REAL8 tmp = pWFHM->fAmpRDfalloff;
803  pWFHM->fAmpRDfalloff = 0;
804  pAmp->RDCoefficient[3] = IMRPhenomXHM_RD_Amp_Ansatz(&powers_of_RDfalloff, pWFHM, pAmp);
805  pAmp->RDCoefficient[4] = -1. * IMRPhenomXHM_RD_Amp_DAnsatz(&powers_of_RDfalloff, pWFHM, pAmp) / pAmp->RDCoefficient[3];
806  pWFHM->fAmpRDfalloff = tmp;
807  }
808  if(pAmp->nCoefficientsRDAux > 0)
809  IMRPhenomXHM_RDAux_Amp_Coefficients(pWF22, pWFHM, pAmp);
810  break;
811  }
812  default:{
813  XLAL_ERROR_VOID(XLAL_EINVAL, "Error in IMRPhenomXHM_RD_Amp_Coefficients: IMRPhenomXHMRingdownAmpVersion is not valid.\n");
814  }
815  }
816 }
817 
819 
820  for(UINT2 i = 0; i < pAmp->nCollocPtsRDAux; i++)
822  IMRPhenomX_UsefulPowers powers_of_fRDAux;
823  IMRPhenomX_Initialize_Powers(&powers_of_fRDAux, pAmp->fRDAux);
824  pAmp->CollocationPointsValuesAmplitudeRDAux[pAmp->nCollocPtsRDAux] = IMRPhenomXHM_RD_Amp_Ansatz(&powers_of_fRDAux, pWFHM, pAmp); //pAmp->CollocationPointsValuesAmplitudeRD[0];
825  pAmp->CollocationPointsValuesAmplitudeRDAux[pAmp->nCollocPtsRDAux + 1] = IMRPhenomXHM_RD_Amp_DAnsatz(&powers_of_fRDAux, pWFHM, pAmp);
828  }
829 
831  pAmp->CollocationPointsFreqsAmplitudeRDAux[1] = 0.5 * (pAmp->fAmpMatchIM + pAmp->fRDAux); // First Chebyshev node
834 
835 
836  /* GSL objects for solving system of equations via LU decomposition */
837  gsl_vector *b, *x;
838  gsl_matrix *A;
839  gsl_permutation *p;
840  int signum; // No need to set, used internally by gsl_linalg_LU_decomp
841 
842  p = gsl_permutation_alloc(pAmp->nCoefficientsRDAux);
843  b = gsl_vector_alloc(pAmp->nCoefficientsRDAux);
844  x = gsl_vector_alloc(pAmp->nCoefficientsRDAux);
845  A = gsl_matrix_alloc(pAmp->nCoefficientsRDAux, pAmp->nCoefficientsRDAux);
846 
847  /* Define linear system of equations */
848 
849  //FIXME: be careful with the indexing, where are the RDaux CollocPoints in CollocationPointsValuesAmplitudeRD?
850  // Should be at the end, although this region goes before than the the "normal RD region".
851  for(INT4 i = 0; i < pAmp->nCoefficientsRDAux; i++){
852  // b is the vector with the values of collocation points
853  gsl_vector_set(b, i, pAmp->CollocationPointsValuesAmplitudeRDAux[i]);
854  //FIXME: distinguish InterAmp ansatzaes versions
855  // Set system matrix: Polynomial at the collocation points frequencies + derivative at the right boundary
856  /* A = (1, f1, f1^2, f1^3, f1^4)
857  (1, f2, f2^2, f2^3, f2^4)
858  (1, f3, f3^2, f3^3, f3^4)
859  (0, 1, f3, f3^2, f3^3)
860  Until number of collocation points
861  */
862  REAL8 fcollpoint = pAmp->CollocationPointsFreqsAmplitudeRDAux[i];
863  REAL8 fpower = 1.; // 1, f, f^2, f^3, f^4, ...
864  if (i < pAmp->nCoefficientsRDAux - 1){
865  for(INT4 j = 0; j < pAmp->nCoefficientsRDAux; j++){
866  gsl_matrix_set(A, i, j, fpower);
867  fpower *= fcollpoint;
868  }
869  }
870  else{ // Last row of the matrix for the derivative
871  fpower = 1.;
872  gsl_matrix_set(A, i, 0, 0.);
873  for(INT4 j = 1; j < pAmp->nCoefficientsRDAux; j++){
874  gsl_matrix_set(A, i, j, j * fpower);
875  fpower *= fcollpoint;
876  }
877  }
878  }
879 
880  /* We now solve the system A x = b via an LU decomposition. x is the solution vector */
881  gsl_linalg_LU_decomp(A, p, &signum);
882  gsl_linalg_LU_solve(A, p, b, x);
883 
884  for (INT4 i = 0; i < pAmp->nCoefficientsRDAux; i++){
885  pAmp->RDAuxCoefficient[i] = gsl_vector_get(x, i);
886  }
887 
888  gsl_vector_free(b);
889  gsl_vector_free(x);
890  gsl_matrix_free(A);
891  gsl_permutation_free(p);
892 }
893 
894 /************** Amplitude Ringdown Ansatz *************/
895 
896 // For the modes with mixing this is the ansatz of the spheroidal part.
898 
899  double ff = powers_of_Mf->itself;
900  int RDAmpFlag = pWFHM->IMRPhenomXHMRingdownAmpVersion;
901  double frd = pWFHM->fRING;
902  double fda = pWFHM->fDAMP;
903  double dfr = ff - frd;
904  double ampRD = 0.;
905 
906  switch ( RDAmpFlag )
907  {
908  case 0: /* Canonical, 3 fitted coefficients + fring, fdamp, lc that are fixed. sigma is also fixed except for the 21 mode. */
909  { // Only for the 122018 release.
910  double dfd = fda * pAmp->RDCoefficient[2];
911  double lc = pAmp->RDCoefficient[3];
912  ampRD = (fda *fabs(pAmp->RDCoefficient[0]) * pAmp->RDCoefficient[2])*exp(- dfr * pAmp->RDCoefficient[1] / dfd )/ (dfr*dfr + dfd*dfd)*pow(ff,-lc);
913  // The line below returns the strain amplitude
914  // if (pAmp->RDRescaleFactor == 0){
915  // ampRD *= (pWFHM->ampNorm * powers_of_Mf->m_seven_sixths);
916  // //printf("%.10f %.16e\n", ff, ampRD);
917  // }
918  break;
919  }
920  case 1:
921  case 2:
922  {
923  if(pAmp->nCoefficientsRDAux > 0 && !IMRPhenomX_StepFuncBool(ff, pAmp->fRDAux)){
924  /* Polynomial */
925  double fpower = 1.;
926  for (UINT2 i = 0; i < pAmp->nCoefficientsRDAux; i++){
927  ampRD += fpower * pAmp->RDAuxCoefficient[i];
928  fpower *= ff;
929  }
930  }
931  else if (pWFHM->fAmpRDfalloff > 0 && IMRPhenomX_StepFuncBool(ff, pWFHM->fAmpRDfalloff)){
932  ampRD = pAmp->RDCoefficient[3] * exp(- pAmp->RDCoefficient[4] * (ff - pWFHM->fAmpRDfalloff));
933  }
934  else{ /* Lorentzian with exponential falloff */
935  double dfd = fda * pAmp->RDCoefficient[2];
936  ampRD = pAmp->RDCoefficient[0] * fda / ( exp(pAmp->RDCoefficient[1] / dfd * dfr) * (dfr * dfr + dfd * dfd)); // * pWF->ampNorm * factor;
937  }
938  if (pAmp->RDRescaleFactor !=0){
939  ampRD /= RescaleFactor(powers_of_Mf, pAmp, pAmp->RDRescaleFactor);
940  }
941  break;
942  }
943  default:
944  {
945  XLAL_ERROR_REAL8(XLAL_EINVAL, "Error in IMRPhenomXHM_RD_Amp_Ansatz: IMRPhenomXHMRingdownAmpVersion = %i is not valid. \n", RDAmpFlag);
946  }
947  }
948 
949  return ampRD;
950 }
951 
952 
953 /*** Derivative of the RD Ansatz for modes without mixing ***/
954 
956 
957  double ff = powers_of_Mf->itself;
958  int RDAmpFlag = pWF->IMRPhenomXHMRingdownAmpVersion;
959  double frd = pWF->fRING;
960  double fda = pWF->fDAMP;
961  double DampRD;
962  double numerator,denom;
963 
964  switch ( RDAmpFlag )
965  {
966  case 0: /* Canonical, 3 fitted coefficients + fring, fdamp, lc that are fixed. sigma is also fixed except for the 21 mode. */
967  {
968  double dfd = fda * pAmp->RDCoefficient[2];
969  double lc = pAmp->RDCoefficient[3];
970  double lambda = pAmp->RDCoefficient[1], expon;
971  numerator = fabs(pAmp->RDCoefficient[0])*pow(ff,-1.-lc)*( dfd*lc*(frd*frd + dfd*dfd)
972  + ff*( frd*frd*lambda - 2.*dfd*frd*(1.+lc) + dfd*dfd*lambda)
973  + ff*ff*( -2.*lambda*frd + dfd*(2.+lc) )
974  + ff*ff*ff*( lambda )
975  );
976  denom = frd*frd + dfd*dfd - ff*2.*frd + ff*ff;
977  expon = exp(-(ff-frd)*lambda/dfd);
978  DampRD = -numerator*expon/(denom*denom);
979  break;
980  }
981  case 1:
982  case 2:
983  {
984  double dfr = ff - frd;
985  numerator = pAmp->RDCoefficient[0] * (dfr * dfr * pAmp->RDCoefficient[1] + 2 * fda * dfr * pAmp->RDCoefficient[2] + fda * fda * pAmp->RDCoefficient[1] * pAmp->RDCoefficient[2] * pAmp->RDCoefficient[2]);
986  denom = (dfr * dfr + fda * fda * pAmp->RDCoefficient[2] * pAmp->RDCoefficient[2]);
987  denom = pAmp->RDCoefficient[2] * denom * denom * exp(dfr * pAmp->RDCoefficient[1] / (fda * pAmp->RDCoefficient[2]));
988  DampRD = - numerator / denom;
989  break;
990  }
991  default:
992  {
993  XLAL_ERROR_REAL8(XLAL_EINVAL, "Error in IMRPhenomXHM_RD_Amp_DAnsatz: IMRPhenomXHMRingdownAmpVersion = %i is not valid. \n", RDAmpFlag);
994  }
995  }
996 
997  return DampRD;
998 }
999 
1000 /*** Derivative of the RD Ansatz for modes with mixing ***/
1001 // It can not be obtained analytically, so we use finite difference of 4th order
1003 
1004  double ff = powers_of_Mf->itself;
1005  double df = 10e-10;
1006  double Nder;
1007  double fun2R = ff + 2*df;
1008  double funR = ff + df;
1009  double funL = ff - df;
1010  double fun2L = ff - 2*df;
1011  IMRPhenomX_UsefulPowers powers_of_fun;
1012 
1013  IMRPhenomX_Initialize_Powers(&powers_of_fun,fun2R);
1014  fun2R = cabs(SpheroidalToSpherical(&powers_of_fun, pAmp22, pPhase22, pAmp, pPhase, pWFHM, pWF22));
1015 
1016  IMRPhenomX_Initialize_Powers(&powers_of_fun,funR);
1017  funR = cabs(SpheroidalToSpherical(&powers_of_fun, pAmp22, pPhase22, pAmp, pPhase, pWFHM, pWF22));
1018 
1019  IMRPhenomX_Initialize_Powers(&powers_of_fun,funL);
1020  funL = cabs(SpheroidalToSpherical(&powers_of_fun, pAmp22, pPhase22, pAmp, pPhase, pWFHM, pWF22));
1021 
1022  IMRPhenomX_Initialize_Powers(&powers_of_fun,fun2L);
1023  fun2L = cabs(SpheroidalToSpherical(&powers_of_fun, pAmp22, pPhase22, pAmp, pPhase, pWFHM, pWF22));
1024 
1025  Nder = (-fun2R + 8*funR - 8*funL + fun2L )/(12*df);
1026 
1027  return Nder;
1028 
1029 }
1030 
1031 /* There are cases for some modes where the ringdown amplitude is very low compared to the inspiral and the intermediate reconstruction fails to connect them.
1032  For those cases we remove the two collocation points and reconstruct with a third order polynomial or with a linear one for the 32 mode. */
1033 void IMRPhenomXHM_Ringdown_Amplitude_Veto(double *V2, double *V3, double V4, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22){
1034 
1035  double threshold;
1036  if(pWFHM->modeTag == 32){
1037  threshold = 0.1/(pWF22->ampNorm);
1038  if(1./V4 < threshold){
1039  *V2 = 1.;
1040  *V3 = 1.;
1042  }
1043  }
1044  else{
1045  threshold = 0.01/(pWF22->ampNorm);
1046  if(1./V4 < threshold){
1047  *V2 = 1.;
1048  *V3 = 1.;
1049  pWFHM->IMRPhenomXHMIntermediateAmpVersion = 1032;
1050  }
1051  }
1052 }
1053 
1054 
1055 
1056 /****************************************/
1057 /* */
1058 /* PHASE */
1059 /* */
1060 /****************************************/
1061 
1062 /* Fits over parameter space for the ringdown phase quantities. */
1063 
1064 // The spin parameter S = (m1^2*chi1 + m2^2*chi2)/(m1^2 + m2^2)
1065 
1066 /* Start of Phase Parameter Space Fits */
1067 
1068 static double IMRPhenomXHM_RD_Phase_22_alpha2(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag) {
1069  double total=0;
1070  switch (RDPhaseFlag){
1071  case 122019:{
1072  double eta = pWF->eta;
1073  double S = pWF->STotR;
1074  double eta2,eta3,eta4,delta=sqrt(1-4*eta);
1075  eta2 = pow(eta,2);
1076  eta3 = pow(eta,3);
1077  eta4 = pow(eta,4);
1078  double noSpin = 0.2088669311744758 - 0.37138987533788487*eta + 6.510807976353186*eta2 - 31.330215053905395*eta3 + 55.45508989446867*eta4;
1079  double eqSpin = ((0.2393965714370633 + 1.6966740823756759*eta - 16.874355161681766*eta2 + 38.61300158832203*eta3)*S)/(1. - 0.633218538432246*S);
1080  double uneqSpin = pWF->dchi*(0.9088578269496244*pow(eta,2.5) + 15.619592332008951*pWF->dchi*pow(eta,3.5))*delta;
1081  total = noSpin + eqSpin + uneqSpin;
1082  break;
1083  }
1084  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Phase_22_alpha2: version % is not valid.", RDPhaseFlag);}
1085  }
1086  return total;
1087 }
1088 
1089 static double IMRPhenomXHM_RD_Phase_22_alphaL(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag) {
1090  double total=0;
1091  switch (RDPhaseFlag){
1092  case 122019:{
1093  double eta = pWF->eta;
1094  double S = pWF->STotR;
1095  double delta=sqrt(1.- 4.*eta),eta2,eta3,eta4,S2;
1096  eta2 = pow(eta,2);
1097  eta3 = pow(eta,3);
1098  eta4 = pow(eta,4);
1099  S2 = pow(S,2);
1100  double noSpin = eta*(-1.1926122248825484 + 2.5400257699690143*eta - 16.504334734464244*eta2 + 27.623649807617376*eta3);
1101  double eqSpin = eta3*S*(35.803988443700824 + 9.700178927988006*S - 77.2346297158916*S2) + eta*S*(0.1034526554654983 - 0.21477847929548569*S - 0.06417449517826644*S2) + eta2*S*(-4.7282481007397825 + 0.8743576195364632*S + 8.170616575493503*S2) + eta4*S*(-72.50310678862684 - 39.83460092417137*S + 180.8345521274853*S2);
1102  double uneqSpin = (-0.7428134042821221*pWF->chi1L*pow(eta,3.5) + 0.7428134042821221*pWF->chi2L*pow(eta,3.5) + 17.588573345324154*pow(pWF->chi1L,2)*pow(eta,4.5) - 35.17714669064831*pWF->chi1L*pWF->chi2L*pow(eta,4.5) + 17.588573345324154*pow(pWF->chi2L,2)*pow(eta,4.5))*delta;
1103  total = noSpin + eqSpin + uneqSpin;
1104  break;
1105  }
1106  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Phase_22_alphaL: version % is not valid.", RDPhaseFlag);}
1107  }
1108  return total;
1109 }
1110 
1111 /**************** 32 specific fits ***************/
1113  double total;
1114  switch (RDPhaseFlag){
1115  case 122019:{
1116  double eta = pWF->eta;
1117  double S = pWF->STotR;
1118  double eta2,eta3,eta4,eta5,S2,S3,S4;
1119  eta2 = pow(eta,2);
1120  eta3 = pow(eta,3);
1121  eta4 = pow(eta,4);
1122  eta5 = pow(eta,5);
1123  S2 = pow(S,2);
1124  S3 = pow(S,3);
1125  S4 = pow(S,4);
1126  double noSpin = 11.851438981981772 + 167.95086712701223*eta - 4565.033758777737*eta2 + 61559.132976189896*eta3 - 364129.24735853914*eta4 + 739270.8814129328*eta5;
1127  double eqSpin = (9.506768471271634 + 434.31707030999445*eta - 8046.364492927503*eta2 + 26929.677144312944*eta3)*S + (-5.949655484033632 - 307.67253970367034*eta + 1334.1062451631644*eta2 + 3575.347142399199*eta3)*S2 + (3.4881615575084797 - 2244.4613237912527*eta + 24145.932943269272*eta2 - 60929.87465551446*eta3)*S3 + (15.585154698977842 - 2292.778112523392*eta + 24793.809334683185*eta2 - 65993.84497923202*eta3)*S4;
1128  double uneqSpin = 465.7904934097202*pWF->dchi*sqrt(1.-4.*eta)*eta2;
1129  total = noSpin + eqSpin + uneqSpin;
1130  break;
1131  }
1132  case 122022:{
1133  double eta = pWF->eta;
1134  double delta = pWF->delta;
1135  double S = pWF->STotR;
1136  double chidiff = pWF->dchi_half;
1137  double eta1 = eta;
1138  double eta2 = eta1 * eta1;
1139  double eta3 = eta1 * eta2;
1140  double eta4 = eta1 * eta3;
1141  double S1 = S;
1142  double chidiff1 = chidiff;
1143  total = chidiff1*delta*(-3437.4083682807154*eta2 + 29349.322789666763*eta3 - 46822.26853496922*eta4) + S1*(18.280316689743625*(46.16185108285708*eta1 - 515.3588338752849*eta2 + 1475.462268010926*eta3) + 2.1970246269696427*(370.70040479593024*eta1 - 4329.153306191607*eta2 + 12065.631680744644*eta3)*S1) + (12.080898026205173 - 79.10914761468462*eta1 - 89.82426456495799*eta2 + 915.6864093792078*eta3)*pow(1 - 9.443713298364061*eta1 + 22.353898754970686*eta2,-1);
1144  break;
1145  }
1146  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Phase_32_SpheroidalTimeShift: version %i is not valid.", RDPhaseFlag);}
1147  }
1148  return total;
1149 }
1150 
1152  double total;
1153  switch (RDPhaseFlag){
1154  case 122019:
1155  {
1156  double eta = pWF->eta;
1157  double S = pWF->STotR;
1158  double eta2,eta3,eta4,eta5,eta6,eta7,S2,S3,S4;
1159  eta2 = pow(eta,2);
1160  eta3 = pow(eta,3);
1161  eta4 = pow(eta,4);
1162  eta5 = pow(eta,5);
1163  eta6 = pow(eta,6);
1164  eta7 = pow(eta,7);
1165  S2 = pow(S,2);
1166  S3 = pow(S,3);
1167  S4 = pow(S,4);
1168  double noSpin = -1.3328895897490733 - 22.209549522908667*eta + 1056.2426481245027*eta2 - 21256.376324666326*eta3 + 246313.12887984765*eta4 - 1.6312968467540336e6*eta5 + 5.614617173188322e6*eta6 - 7.612233821752137e6*eta7;
1169  double eqSpin = (S*(-1.622727240110213 + 0.9960210841611344*S - 1.1239505323267036*S2 - 1.9586085340429995*S3 + eta2*(196.7055281997748 + 135.25216875394943*S + 1086.7504825459278*S2 + 546.6246807461155*S3 - 312.1010566468068*S4) + 0.7638287749489343*S4 + eta*(-47.475568056234245 - 35.074072557604445*S - 97.16014978329918*S2 - 34.498125910065156*S3 + 24.02858084544326*S4) + eta3*(62.632493533037625 - 22.59781899512552*S - 2683.947280170815*S2 - 1493.177074873678*S3 + 805.0266029288334*S4)))/(-2.950271397057221 + 1.*S);
1170  double uneqSpin = (sqrt(1.-4.*eta)*(pWF->chi2L*pow(eta,2.5)*(88.56162028006072 - 30.01812659282717*S) + pWF->chi2L*eta2*(43.126266433486435 - 14.617728550838805*S) + pWF->chi1L*eta2*(-43.126266433486435 + 14.617728550838805*S) + pWF->chi1L*pow(eta,2.5)*(-88.56162028006072 + 30.01812659282717*S)))/(-2.950271397057221 + 1.*S);
1171  total = noSpin + eqSpin + uneqSpin;
1172  break;
1173  }
1174  case 122022:{
1175  double eta = pWF->eta;
1176  double delta = pWF->delta;
1177  double S = pWF->chiPNHat;
1178  double chidiff = pWF->dchi_half;
1179  double eta1 = eta;
1180  double eta2 = eta1 * eta1;
1181  double eta3 = eta1 * eta2;
1182  double eta4 = eta1 * eta3;
1183  double eta5 = eta1 * eta4;
1184  double S1 = S;
1185  double S2 = S1 * S1;
1186  double chidiff1 = chidiff;
1187  double chidiff2 = chidiff1 * chidiff1;
1188  total = chidiff1*delta*(-4055.661463620154*eta3 + 48053.79215999518*eta4 - 126644.96534635697*eta5) + chidiff2*(1489.2939046368886*eta3 - 11333.726790227513*eta4 + 21470.681301598026*eta5) + S1*(-0.09759112086805133*(-3.6379140102351064 + 510.0158912180661*eta1 - 12528.715040030444*eta2 + 82416.24893999398*eta3 - 161740.0041427807*eta4) - 0.20117612026208484*(-5.676646590653427 + 183.78422258983136*eta1 - 3807.617101722895*eta2 + 24219.360141326677*eta3 - 45892.369216999985*eta4)*S1 + 0.06537048555519695*(21.388069487574892 + 408.2245599781871*eta1 - 9652.666650065075*eta2 + 57782.086859487965*eta3 - 110112.73697613904*eta4)*S2) + (-1.3903824533899325 + 13.761709564309667*eta1 - 10.633427224975128*eta2 - 263.01839936998965*eta3 + 736.1912896690361*eta4)*pow(1 - 9.458921853648649*eta1 + 22.932673653997934*eta2,-1);
1189  break;
1190  }
1191  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Phase_32_SpheroidalPhaseShift: version %i is not valid.", RDPhaseFlag);}
1192  }
1193  return total;
1194 }
1195 
1196 // collocation points
1197 static double IMRPhenomXHM_RD_Phase_32_p1(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag) {
1198  double total;
1199  switch (RDPhaseFlag){
1200  case 122019:{
1201  double eta = pWF->eta;
1202  double S = pWF->STotR;
1203  double eta2,eta3,eta4,eta5,S2,S3,S4,S5;
1204  eta2 = pow(eta,2);
1205  eta3 = pow(eta,3);
1206  eta4 = pow(eta,4);
1207  eta5 = pow(eta,5);
1208  S2 = pow(S,2);
1209  S3 = pow(S,3);
1210  S4 = pow(S,4);
1211  S5 = pow(S,5);
1212  double noSpin = 3169.372056189274 + 426.8372805022653*eta - 12569.748101922158*eta2 + 149846.7281073725*eta3 - 817182.2896823225*eta4 + 1.5674053633767858e6*eta5;
1213  double eqSpin = (19.23408352151287 - 1762.6573670619173*eta + 7855.316419853637*eta2 - 3785.49764771212*eta3)*S + (-42.88446003698396 + 336.8340966473415*eta - 5615.908682338113*eta2 + 20497.5021807654*eta3)*S2 + (13.918237996338371 + 10145.53174542332*eta - 91664.12621864353*eta2 + 201204.5096556517*eta3)*S3 + (-24.72321125342808 - 4901.068176970293*eta + 53893.9479532688*eta2 - 139322.02687945773*eta3)*S4 + (-61.01931672442576 - 16556.65370439302*eta + 162941.8009556697*eta2 - 384336.57477596396*eta3)*S5;
1214  double uneqSpin = pWF->dchi*sqrt(1.-4.*eta)*eta2*(641.2473192044652 - 1600.240100295189*pWF->chi1L*eta + 1600.240100295189*pWF->chi2L*eta + 13275.623692212472*eta*S);
1215  total = noSpin + eqSpin + uneqSpin;
1216  break;
1217  }
1218  case 122022:{
1219  double eta = pWF->eta;
1220  double delta = pWF->delta;
1221  double S = pWF->STotR;
1222  double chidiff = pWF->dchi_half;
1223  double eta1 = eta;
1224  double eta2 = eta1 * eta1;
1225  double eta3 = eta1 * eta2;
1226  double eta4 = eta1 * eta3;
1227  double S1 = S;
1228  double S2 = S1 * S1;
1229  double chidiff1 = chidiff;
1230  total = 3221.932636435056 - 134.59521937837278*eta1 + chidiff1*delta*(319.06177738567885*eta2 - 20578.40007594641*eta3 + 101859.1659970414*eta4) + S1*(24.059115241154707*(436.41245673494626*eta2 - 2938.503437122471*eta3 + 5027.414440730744*eta4) + 21.848741292340005*(1251.706577839354*eta2 - 14171.490147583942*eta3 + 36914.6553449061*eta4)*S1 + 15.901300902033508*(-149.2789474539545*eta2 + 3483.608736789833*eta3 - 11289.97178789606*eta4)*S2)*pow(-1.803337035190313 + S1,-1) - 0.012868288384497346*pow(0.000187943994873322 + pow(-0.1923690355128322 + eta1,2),-1);
1231  break;
1232  }
1233  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Phase_32_p1: version % is not valid.", RDPhaseFlag);}
1234  }
1235  return total;
1236 }
1237 // collocation points
1238 static double IMRPhenomXHM_RD_Phase_32_p2(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag) {
1239  double total;
1240  switch (RDPhaseFlag){
1241  case 122019:{
1242  double eta = pWF->eta;
1243  double S = pWF->STotR;
1244  double eta2,eta3,S2,S3,S4;
1245  eta2 = pow(eta,2);
1246  eta3 = pow(eta,3);
1247  S2 = pow(S,2);
1248  S3 = pow(S,3);
1249  S4 = pow(S,4);
1250  double noSpin = 3131.0260952676376 + 206.09687819102305*eta - 2636.4344627081873*eta2 + 7475.062269742079*eta3;
1251  double eqSpin = (49.90874152040307 - 691.9815135740145*eta - 434.60154548208334*eta2 + 10514.68111669422*eta3)*S + (97.3078084654917 - 3458.2579971189534*eta + 26748.805404989867*eta2 - 56142.13736008524*eta3)*S2 + (-132.49105074500454 + 429.0787542102207*eta + 7269.262546204149*eta2 - 27654.067482558712*eta3)*S3 + (-227.8023564332453 + 5119.138772157134*eta - 34444.2579678986*eta2 + 69666.01833764123*eta3)*S4;
1252  double uneqSpin = 477.51566939885424*pWF->dchi*sqrt(1.-4.*eta)*eta2;
1253  total = noSpin + eqSpin + uneqSpin;
1254  break;
1255  }
1256  case 122022:{
1257  double eta = pWF->eta;
1258  double delta = pWF->delta;
1259  double S = pWF->STotR;
1260  double chidiff = pWF->dchi_half;
1261  double eta1 = eta;
1262  double eta2 = eta1 * eta1;
1263  double eta3 = eta1 * eta2;
1264  double eta4 = eta1 * eta3;
1265  double S1 = S;
1266  double chidiff1 = chidiff;
1267  total = 3169.5772463611165 + 56.56534589293562*eta1 - 863.5731390762933*eta2 + 2693.8619211321557*eta3 + chidiff1*delta*(-2818.2944800258847*eta2 + 15684.658457287562*eta3 + 14379.128341035908*eta4) + S1*(-0.16388886708177886*(334.30009385854424*eta2 - 154749.10305716714*eta3 + 613903.6107269318*eta4) + 11.950465013745157*(1079.481585746054*eta2 - 11981.85336876442*eta3 + 30911.708103120814*eta4)*S1)*pow(-1.169876031327984 + S1,-1) - 0.009425837438775205*pow(0.00016960223009674388 + pow(-0.20083535429185695 + eta1,2),-1);
1268  break;
1269  }
1270  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Phase_32_p2: version % is not valid.", RDPhaseFlag);}
1271  }
1272  return total;
1273 }
1274 // collocation points
1275 static double IMRPhenomXHM_RD_Phase_32_p3(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag) {
1276  double total;
1277  switch (RDPhaseFlag){
1278  case 122019:{
1279  double eta = pWF->eta;
1280  double S = pWF->STotR;
1281  double eta2,eta3,S2,S3,S4;
1282  eta2 = pow(eta,2);
1283  eta3 = pow(eta,3);
1284  S2 = pow(S,2);
1285  S3 = pow(S,3);
1286  S4 = pow(S,4);
1287  double noSpin = 3082.803556599222 + 76.94679795837645*eta - 586.2469821978381*eta2 + 977.6115755788503*eta3;
1288  double eqSpin = (45.08944710349874 - 807.7353772747749*eta + 1775.4343704616288*eta2 + 2472.6476419567534*eta3)*S + (95.57355060136699 - 2224.9613131172046*eta + 13821.251641893134*eta2 - 25583.314298758105*eta3)*S2 + (-144.96370424517866 + 2268.4693587493093*eta - 10971.864789147161*eta2 + 16259.911572457446*eta3)*S3 + (-227.8023564332453 + 5119.138772157134*eta - 34444.2579678986*eta2 + 69666.01833764123*eta3)*S4;
1289  double uneqSpin = 378.2359918274837*pWF->dchi*sqrt(1.-4.*eta)*eta2;
1290  total = noSpin + eqSpin + uneqSpin;
1291  break;
1292  }
1293  case 122022:{
1294  double eta = pWF->eta;
1295  double delta = pWF->delta;
1296  double S = pWF->STotR;
1297  double chidiff = pWF->dchi_half;
1298  double eta1 = eta;
1299  double eta2 = eta1 * eta1;
1300  double eta3 = eta1 * eta2;
1301  double eta4 = eta1 * eta3;
1302  double eta5 = eta1 * eta4;
1303  double eta6 = eta1 * eta5;
1304  double S1 = S;
1305  double chidiff1 = chidiff;
1306  total = 3119.6603946770488 + 168.61554447853712*eta1 - 1777.6654596491376*eta2 + 5037.407962552042*eta3 + chidiff1*delta*(45693.135566736484*eta2 - 789332.4959926775*eta3 + 4.460496312695218e6*eta4 - 8.176309211912101e6*eta5) + chidiff1*delta*(7840.121424232572*eta2 - 47166.09840761356*eta3 + 66597.52917033392*eta4)*S1 + S1*(-6.019579546899472*(14538.163921822728*eta2 - 318911.4362763759*eta3 + 2.6041867832020866e6*eta4 - 9.288489508236282e6*eta5 + 1.2170972980338342e7*eta6) - 7.739304888898913*(11114.219992659659*eta2 - 231541.9569739445*eta3 + 1.8069370995120746e6*eta4 - 6.203273456127891e6*eta5 + 7.874294046591697e6*eta6)*S1) - 0.0003538235766427527*pow(8.2810517855422e-6 + pow(-0.2075868299995718 + eta1,2),-1);
1307  break;
1308  }
1309  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Phase_32_p3: version % is not valid.", RDPhaseFlag);}
1310  }
1311  return total;
1312 }
1313 // collocation points
1314 static double IMRPhenomXHM_RD_Phase_32_p4(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag) {
1315  double total;
1316  switch (RDPhaseFlag){
1317  case 122019:{
1318  double eta = pWF->eta;
1319  double S = pWF->STotR;
1320  double eta2,S2,S3,S4;
1321  eta2 = pow(eta,2);
1322  S2 = pow(S,2);
1323  S3 = pow(S,3);
1324  S4 = pow(S,4);
1325  double noSpin = 3077.0657367004565 + 64.99844502520415*eta - 357.38692756785395*eta2;
1326  double eqSpin = (34.793450080444714 - 986.7751755509875*eta - 9490.641676924794*pow(eta,3) + 5700.682624203565*eta2)*S + (57.38106384558743 - 1644.6690499868596*eta - 19906.416384606226*pow(eta,3) + 11008.881935880598*eta2)*S2 + (-126.02362949830213 + 3169.3397351803583*eta + 62863.79877094988*pow(eta,3) - 26766.730897942085*eta2)*S3 + (-169.30909412804587 + 4900.706039920717*eta + 95314.99988114933*pow(eta,3) - 41414.05689348732*eta2)*S4;
1327  double uneqSpin = 390.5443469721231*pWF->dchi*sqrt(1.-4.*eta)*eta2;
1328  total = noSpin + eqSpin + uneqSpin;
1329  break;
1330  }
1331  case 122022:{
1332  double eta = pWF->eta;
1333  double delta = pWF->delta;
1334  double S = pWF->STotR;
1335  double chidiff = pWF->dchi_half;
1336  double eta1 = eta;
1337  double eta2 = eta1 * eta1;
1338  double eta3 = eta1 * eta2;
1339  double eta4 = eta1 * eta3;
1340  double eta5 = eta1 * eta4;
1341  double eta6 = eta1 * eta5;
1342  double eta7 = eta1 * eta6;
1343  double S1 = S;
1344  double chidiff1 = chidiff;
1345  total = 3063.2702533356364 + 74.20321762511647*eta1 - 346.6653326379183*eta2 + chidiff1*delta*(2604.6711121030685*eta2 - 25322.83641432119*eta3 + 64521.907625802785*eta4) + 17.748987975469845*(6659.676835974436*eta2 - 207712.54687648916*eta3 + 2.5247192995989644e6*eta4 - 1.4825576629165621e7*eta5 + 4.215660954626601e7*eta6 - 4.662796037240443e7*eta7)*S1*pow(-1.2617538728082525 + S1,-1) - 2.211595095940034e-6*pow(-0.000010240443266599763 + pow(-0.17594819839300638 + eta1,2),-1);
1346  break;
1347  }
1348  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Phase_32_p4: version % is not valid.", RDPhaseFlag);}
1349  }
1350  return total;
1351 }
1352 
1353 
1354 static double IMRPhenomXHM_RD_Phase_32_p5(IMRPhenomXWaveformStruct *pWF, int RDPhaseFlag){
1355  double total=0;
1356  switch (RDPhaseFlag){
1357  case 122022:{
1358  double eta = pWF->eta;
1359  double delta = pWF->delta;
1360  double S = pWF->STotR;
1361  double chidiff = pWF->dchi_half;
1362  double eta1 = eta;
1363  double eta2 = eta1 * eta1;
1364  double eta3 = eta1 * eta2;
1365  double eta4 = eta1 * eta3;
1366  double eta5 = eta1 * eta4;
1367  double eta6 = eta1 * eta5;
1368  double eta7 = eta1 * eta6;
1369  double S1 = S;
1370  double chidiff1 = chidiff;
1371  total = 3099.009925231132 - 5.823302919769178*eta1 + chidiff1*delta*(1936.2628510750844*eta2 - 32280.62038532877*eta3 + 97943.49145078743*eta4) + 2.8281136508769498*(6856.762845446495*eta2 - 201653.94475043533*eta3 + 2.4217205751584964e6*eta4 - 1.4582166806161262e7*eta5 + 4.377208838391116e7*eta6 - 5.214533274483399e7*eta7)*S1*pow(-1.0624647822393556 + S1,-1);
1372  break;
1373  }
1374  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_RD_Phase_32_p5: version % is not valid.", RDPhaseFlag);}
1375  }
1376  return total;
1377 }
1378 
1379 /* End of Phase Parameter Space Fits */
1380 
1381 /************** ANSATZ PHASE DERIVATIVE **************/
1382 
1384 
1385  double frd = pWFHM->fRING;
1386  double fda = pWFHM->fDAMP;
1387  double dphaseRD;
1388 
1389  switch ( pWFHM->MixingOn )
1390  {
1391  case 0:
1392  {
1393  /* rescaling of the 22 ansatz -- used for (21),(33),(44) */
1394  /*ansatz:
1395  alpha0 + ((fRDlm^2) alpha2)/(f^2) + alphaL*(fdamplm)/((fdamplm)^2 + (f - fRDlm)^2)*/
1396  dphaseRD = ( pPhase->alpha0 + frd*frd*(pPhase->alpha2)*powers_of_f->m_two + ( (pPhase->alphaL)* fda/(fda*fda +(ff - frd)*(ff - frd)) ) );
1397  break;
1398 
1399  }
1400  case 1:
1401  {
1402  /* calibration of spheroidal ringdown waveform for (32) */
1403 
1404  if(pWFHM->IMRPhenomXHMRingdownPhaseVersion == 122019){
1405  /* ansatz: alpha0 + (alpha2)/(f^2)+ (alpha4)/(f^4) + alphaL*(fdamplm)/((fdamplm)^2 + (f - fRDlm)^2)*/
1406  dphaseRD = ( pPhase->alpha0_S + (pPhase->alpha2_S)*powers_of_f->m_two + (pPhase->alpha4_S)*powers_of_f->m_four +( (pPhase->alphaL_S)* fda/(fda*fda +(ff - frd)*(ff - frd)) ) );
1407  }
1408  else{ // FIXME: 1/eta???
1409  if(pWFHM->fPhaseRDflat > 0 && IMRPhenomX_StepFuncBool(ff, pWFHM->fPhaseRDflat)){
1410  dphaseRD = pPhase->RDCoefficient[5] + pPhase->RDCoefficient[6] * powers_of_f->m_five;
1411  }
1412  else{
1413  /* ansatz: a0 + a1/f + a2/f^2 + a3/f^4 + a4*fdamplm/(fdamplm^2 + (f - fRDlm)^2) */
1414  dphaseRD = ( pPhase->RDCoefficient[0] + pPhase->RDCoefficient[1] * powers_of_f->m_one + pPhase->RDCoefficient[2] * powers_of_f->m_two + pPhase->RDCoefficient[3] * powers_of_f->m_four + pPhase->RDCoefficient[4]*fda / (fda*fda + (ff - frd)*(ff - frd)) );
1415  }
1416  }
1417  break;
1418  }
1419  default:
1420  {XLAL_ERROR(XLAL_EDOM, "Error in IMRPhenomXHM_RD_Phase_Ansatz: version is not valid. Use version 0 for modes (2,1),(3,3),(4,4) and 1 for (3,2).\n");}
1421  }
1422  return dphaseRD;
1423 }
1424 
1425 /************** ANSATZ INTEGRATED PHASE **************/
1426 
1428 
1429  double invf = powers_of_f->m_one;
1430  double phaseRD;
1431  double frd=pWFHM->fRING;
1432  double fda=pWFHM->fDAMP;
1433 
1434  switch ( pWFHM->MixingOn )
1435  {
1436  case 0:
1437  {
1438  /* rescaling of the 22 ansatz -- used for (21),(33),(44) */
1439  /*ansatz:
1440  alpha0 f - fRDlm^2*alpha2)/f + alphaL*ArcTan[(f - fRDlm)/fdamplm]*/
1441  phaseRD = pPhase->alpha0*ff -frd*frd *(pPhase->alpha2)*invf + (pPhase->alphaL)* atan((ff-frd)/fda);
1442  break;
1443  }
1444  case 1:
1445  {
1446  /* calibration of spheroidal ringdown waveform for (32) */
1447  double invf3 = powers_of_f->m_three;
1448 
1449  if(pWFHM->IMRPhenomXHMRingdownPhaseVersion == 122019){
1450  /* ansatz: f alpha0 - (alpha4)/(3 f^3) - (alpha2)/f + alphaL ArcTan[(f - fRDlm)/fdamplm]*/
1451  phaseRD = pPhase->phi0_S+pPhase->alpha0_S*ff -(pPhase->alpha2_S)*invf -1./3.*(pPhase->alpha4_S)*invf3 +(pPhase->alphaL_S)* atan((ff-frd)/fda);
1452  }
1453  else{
1454  if(pWFHM->fPhaseRDflat > 0 && IMRPhenomX_StepFuncBool(ff, pWFHM->fPhaseRDflat)){
1455  // a + b / f^5
1456  phaseRD = pPhase->phi0_S + pPhase->RDCoefficient[7] + pPhase->RDCoefficient[5]*ff - 0.25 * pPhase->RDCoefficient[6] * powers_of_f->m_four;
1457  }
1458  else{
1459  /* ansatz: f a0 + a1 Log(f) - a2/f - a3/(3f^3) + a4*ArcTan[(f - fRDlm)/fdamplm] */
1460  phaseRD = pPhase->phi0_S + pPhase->RDCoefficient[0]*ff + pPhase->RDCoefficient[1]*powers_of_f->log - pPhase->RDCoefficient[2]*invf - 1/3.*pPhase->RDCoefficient[3]*invf3 + pPhase->RDCoefficient[4]*atan( (ff-frd)/fda );
1461  }
1462  }
1463  break;
1464  }
1465  default:
1466  {XLAL_ERROR(XLAL_EDOM, "Error in IMRPhenomXHM_RD_Phase_AnsatzInt: version is not valid. Use version 0 for modes (2,1),(3,3),(4,4) and 1 for (3,2).\n");}
1467  }
1468  return phaseRD;
1469 }
1470 
1472 
1473  double frd = pWFHM->fRING;
1474  double fda = pWFHM->fDAMP;
1475  double ddphaseRD;
1476 
1477 
1478  switch ( pWFHM->MixingOn )
1479  {
1480  case 0:
1481  {
1482  /* rescaling of the 22 ansatz -- used for (21),(33),(44) */
1483  /*ansatz:
1484  alpha0 + ((fRDlm^2) alpha2)/(f^2) + alphaL*(fdamplm)/((fdamplm)^2 + (f - fRDlm)^2)*/
1485  ddphaseRD = -2*frd*frd*(pPhase->alpha2)*powers_of_f->m_three + -2*pPhase->alphaL* fda *(ff-frd)/pow(fda*fda +(ff - frd)*(ff - frd), 2) ;
1486  break;
1487  }
1488  case 1:
1489  {
1490  /* calibration of spheroidal ringdown waveform for (32) */
1491  if(pWFHM->IMRPhenomXHMRingdownPhaseVersion == 122019){
1492  /* ansatz: alpha0 + (alpha2)/(f^2)+ (alpha4)/(f^4) + alphaL*(fdamplm)/((fdamplm)^2 + (f - fRDlm)^2)*/
1493  ddphaseRD = -2*pPhase->alpha2_S*powers_of_f->m_three - 4*(pPhase->alpha4_S)*powers_of_f->m_five - 2*( pPhase->alphaL_S* fda * (ff-frd)/pow(fda*fda +(ff - frd)*(ff - frd), 2) ) ;
1494  }
1495  else{
1496  /* ansatz: a0 + a1/f + a2/f^2 + a3/f^4 + a4*fdamplm/(fdamplm^2 + (f - fRDlm)^2) */
1497  ddphaseRD = ( - pPhase->RDCoefficient[1] * powers_of_f->m_two - 2 * pPhase->RDCoefficient[2] * powers_of_f->m_three - 4 * pPhase->RDCoefficient[3] * powers_of_f->m_five - 2 * pPhase->RDCoefficient[4]*fda*(ff-frd) / pow(fda*fda + (ff - frd)*(ff - frd), 2) );
1498  }
1499  break;
1500  }
1501  default:
1502  {XLAL_ERROR(XLAL_EDOM, "Error in IMRPhenomXHM_RD_Phase_Ansatz: version is not valid. Use version 0 for modes (2,1),(3,3),(4,4) and 1 for (3,2).\n");}
1503  }
1504  return ddphaseRD;
1505 }
static double double delta
int IMRPhenomX_Initialize_Powers(IMRPhenomX_UsefulPowers *p, REAL8 number)
COMPLEX16 SpheroidalToSpherical(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22, IMRPhenomXHMAmpCoefficients *pAmplm, IMRPhenomXHMPhaseCoefficients *pPhaselm, IMRPhenomXHMWaveformStruct *pWFlm, IMRPhenomXWaveformStruct *pWF22)
double RescaleFactor(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMAmpCoefficients *pAmp, UINT2 rescalefactor)
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 void IMRPhenomXHM_RDAux_Amp_Coefficients(IMRPhenomXWaveformStruct *pWF22, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
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)
bool IMRPhenomX_StepFuncBool(const double t, const double t1)
double i
Definition: bh_ringdown.c:118
double e
Definition: bh_ringdown.c:117
double REAL8
uint16_t UINT2
int32_t INT4
#define XLAL_ERROR_VOID(...)
#define XLAL_ERROR_REAL8(...)
#define XLAL_ERROR(...)
XLAL_EDOM
XLAL_EINVAL
list x
list p
REAL8 RDCoefficient[N_MAX_COEFFICIENTS_AMPLITUDE_RING]
REAL8 RDAuxCoefficient[N_MAX_COEFFICIENTS_AMPLITUDE_RDAUX]
ParameterSpaceFit RingdownAmpFits[N_HIGHERMODES_IMPLEMENTED *N_MAX_COEFFICIENTS_AMPLITUDE_RING+N_MAX_COEFFICIENTS_AMPLITUDE_RDAUX]
REAL8 CollocationPointsValuesAmplitudeRD[N_MAX_COEFFICIENTS_AMPLITUDE_RING]
REAL8 CollocationPointsFreqsAmplitudeRDAux[N_MAX_COEFFICIENTS_AMPLITUDE_RDAUX]
REAL8 CollocationPointsFreqsAmplitudeRD[N_MAX_COEFFICIENTS_AMPLITUDE_RING]
REAL8 CollocationPointsValuesAmplitudeRDAux[N_MAX_COEFFICIENTS_AMPLITUDE_RDAUX]
REAL8 RDCoefficient[N_MAX_COEFFICIENTS_PHASE_RING+3]