LALSimulation  5.4.0.1-fe68b98
LALSimIMRPhenomXHM_intermediate.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 //
21 // LALSimIMRPhenomXHM_Intermediate.c
22 //
23 
24 #include <gsl/gsl_poly.h>
25 
26 #include "LALSimIMRPhenomXHM_intermediate.h"
27 
28 /***********************************************/
29 /* */
30 /* AMPLITUDE */
31 /* */
32 /***********************************************/
33 
34 
35 // Fits of the intermediate collocation points over parameter space.
36 
37 /* These parameter-space fits are documented in the supplementary material of https://dcc.ligo.org/P2000011-v2.
38  There are 2 Mathematica notebooks (one for amplitude and one for phase) that read the fits data and automatically generate the C-code below.
39  For more information read https://git.ligo.org/waveforms/reviews/imrphenomx/blob/master/documentation/ParspaceFits/README and the documentation in the notebooks. */
40 
41 // IMRPhenomXHM_Inter_Amp_lm_intX returns the value of the amplitude at the collocation point intX
42 // IMRPhenomXHM_Inter_Amp_lm_dintX returns the value of the derivative of the amplitude at the collocation point intX
43 // for a definition of the frequencies corresponding to int0/1/2, see Sec V.A.
44 
45 // The dominant spin parameter is S = (m1^2*chi1 + m2^2*chi2)/(m1^2 + m2^2)
46 
47 /* Start of Amp Parameter Space Fits */
48 
49 static double IMRPhenomXHM_Inter_Amp_21_int1(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
50  double total=0;
51  switch (InterAmpFlag){
52  case 122018:{
53  double eta = pWF->eta;
54  double S = pWF->STotR;
55  double eta2,S2;
56  eta2 = pow(eta,2);
57  S2 = pow(S,2);
58  double noSpin = sqrt(eta - 4.*eta2)*(21.256776327599113 - 25.594352690383847*eta + 30.14761650482866*eta2);
59  double eqSpin = sqrt(eta - 4.*eta2)*S*(-11.262044985632757 - 1.8167045597937677*S + eta*(-1.1798437990445079 + 6.344825546437461*S - 4.881427482271166*S2));
60  double uneqSpin = -3.6366100759176696*pow(pWF->dchi,2)*(1. - 4.*eta)*eta - 31.60048733143782*(pWF->dchi)*eta2*(1. + 2.1502870640831855*eta2);
61  total = noSpin + eqSpin + uneqSpin;
62  break;
63  }
64  case 122022:{
65  double eta = pWF->eta;
66  double delta = pWF->delta;
67  double S = pWF->STotR;
68  double chidiff = pWF->dchi_half;
69  double eta1 = eta;
70  double eta2 = eta1 * eta1;
71  double eta3 = eta1 * eta2;
72  double eta4 = eta1 * eta3;
73  double eta5 = eta1 * eta4;
74  double S1 = S;
75  double S2 = S1 * S1;
76  double chidiff1 = chidiff;
77  double chidiff2 = chidiff1 * chidiff1;
78  total = fabs(delta*eta1*(chidiff2*(5.159755997682368*eta1 - 30.293198248154948*eta2 + 63.70715919820867*eta3) + chidiff1*(8.262642080222694*eta1 - 415.88826990259116*eta2 + 1427.5951158851076*eta3)) + delta*eta1*(18.55363583212328 - 66.46950491124205*eta1 + 447.2214642597892*eta2 - 1614.178472020212*eta3 + 2199.614895727586*eta4) + chidiff1*eta5*(-1698.841763891122 - 195.27885562092342*S1 - 1.3098861736238572*S2) + delta*eta1*(chidiff1*(34.17829404207186*eta1 - 386.34587928670015*eta2 + 1022.8553774274128*eta3)*S1 + chidiff1*(56.76554600963724*eta1 - 491.4593694689354*eta2 + 1016.6019654342113*eta3)*S2) + delta*eta1*S1*(-8.276366844994188*(1.0677538075697492 - 24.12941323757896*eta1 + 516.7886322104276*eta2 - 4389.799658723288*eta3 + 16770.447637953577*eta4 - 23896.392706809565*eta5) - 1.6908277400304084*(3.4799140066657928 - 29.00026389706585*eta1 + 114.8330693231833*eta2 - 184.13091281984674*eta3 + 592.300353344717*eta4 - 2085.0821513466053*eta5)*S1 - 0.46006975902558517*(-2.1663474937625975 + 826.026625945615*eta1 - 17333.549622759732*eta2 + 142904.08962903373*eta3 - 528521.6231015554*eta4 + 731179.456702448*eta5)*S2));
79  break;
80  }
81  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_21_int1: version %i is not valid.", InterAmpFlag);}
82  }
83  return total;
84 }
85 
86 static double IMRPhenomXHM_Inter_Amp_21_int2(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
87  double total=0;
88  switch (InterAmpFlag){
89  case 122018:{
90  double eta = pWF->eta;
91  double S = pWF->STotR;
92  double eta2 = pow(eta,2);
93  double noSpin = sqrt(eta - 4.*eta2)*(19.15445065708005 - 21.13596229438309*eta + 29.742565944285772*eta2);
94  double eqSpin = sqrt(eta - 4.*eta2)*S*(-12.766814596085734 - 2.123816950673979*S + eta*(-2.913184982025043 + 6.006571549661901*S));
95  double uneqSpin = -25.856046423804255*(pWF->dchi)*eta2*(1. + 5.7871199275552*eta2);
96  total = noSpin + eqSpin + uneqSpin;
97  break;
98  }
99  case 122022:{
100  double eta = pWF->eta;
101  double delta = pWF->delta;
102  double S = pWF->STotR;
103  double chidiff = pWF->dchi_half;
104  double eta1 = eta;
105  double eta2 = eta1 * eta1;
106  double eta3 = eta1 * eta2;
107  double eta4 = eta1 * eta3;
108  double eta5 = eta1 * eta4;
109  double S1 = S;
110  double S2 = S1 * S1;
111  double chidiff1 = chidiff;
112  total = fabs(delta*eta1*(13.757856231617446 - 12.783698329428516*eta1 + 12.048194546899204*eta2) + chidiff1*delta*eta1*(15.107530092096438*eta1 - 416.811753638553*eta2 + 1333.6181181686939*eta3) + chidiff1*eta5*(-1549.6199518612063 - 102.34716990474509*S1 - 3.3637011939285015*S2) + delta*eta1*(chidiff1*(36.358142200869295*eta1 - 384.2123173145321*eta2 + 984.6826660818275*eta3)*S1 + chidiff1*(4.159271594881928*eta1 + 105.10911749116399*eta2 - 639.190132707115*eta3)*S2) + delta*eta1*S1*(-8.097876227116853*(0.6569459700232806 + 9.861355377849485*eta1 - 116.88834714736281*eta2 + 593.8035334117192*eta3 - 1063.0692862578455*eta4) - 1.0546375154878165*(0.745557030602097 + 65.25215540635162*eta1 - 902.5751736558435*eta2 + 4350.442990924205*eta3 - 7141.611333893155*eta4)*S1 - 0.5006664599166409*(10.289020582277626 - 212.00728173197498*eta1 + 2334.0029399672358*eta2 - 11939.621138801092*eta3 + 21974.8201355744*eta4)*S2));
113  break;
114  }
115  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_21_int2: version %i is not valid.", InterAmpFlag);}
116  }
117  return total;
118 }
119 
120 static double IMRPhenomXHM_Inter_Amp_33_int1(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
121  double total=0;
122  switch (InterAmpFlag){
123  case 122018:{
124  double eta = pWF->eta;
125  double S = pWF->STotR;
126  double eta2,eta3,eta4,eta6;
127  eta2 = pow(eta,2);
128  eta3 = pow(eta,3);
129  eta4 = pow(eta,4);
130  eta6 = pow(eta,6);
131  double noSpin = sqrt(eta - 4.*eta2)*(27.927652424857733 - 133.56611389260297*eta + 974.8550901501316*eta2 - 3744.785831952632*eta3 + 5621.897260910284*eta4);
132  double eqSpin = sqrt(eta - 4.*eta2)*S*(7.348313807306079 + eta*(-60.248696675045565 - 37.07212326362276*S) + 5.059236579431119*S + eta2*(159.68630712802727 + 83.33807316873204*S));
133  double uneqSpin = 1412.367880056888*(pWF->dchi)*eta6;
134  total = noSpin + eqSpin + uneqSpin;
135  break;
136  }
137  case 122022:{
138  double eta = pWF->eta;
139  double delta = pWF->delta;
140  double S = pWF->STotR;
141  double chidiff = pWF->dchi_half;
142  double eta1 = eta;
143  double eta2 = eta1 * eta1;
144  double eta3 = eta1 * eta2;
145  double eta4 = eta1 * eta3;
146  double eta5 = eta1 * eta4;
147  double S1 = S;
148  double S2 = S1 * S1;
149  double chidiff1 = chidiff;
150  total = chidiff1*delta*eta1*(-0.3516244197696068*eta1 + 40.425151307421416*eta2 - 148.3162618111991*eta3) + delta*eta1*(26.998512565991778 - 146.29035440932105*eta1 + 914.5350366065115*eta2 - 3047.513201789169*eta3 + 3996.417635728702*eta4) + chidiff1*delta*eta1*(5.575274516197629*eta1 - 44.592719238427094*eta2 + 99.91399033058927*eta3)*S1 + delta*eta1*S1*(-0.5383304368673182*(-7.456619067234563 + 129.36947401891433*eta1 - 843.7897535238325*eta2 + 3507.3655567272644*eta3 - 9675.194644814854*eta4 + 11959.83533107835*eta5) - 0.28042799223829407*(-6.212827413930676 + 266.69059813274475*eta1 - 4241.537539226717*eta2 + 32634.43965039936*eta3 - 119209.70783201039*eta4 + 166056.27237509796*eta5)*S1) + chidiff1*eta5*(199.6863414922219 + 53.36849263931051*S1 + 7.650565415855383*S2);
151  break;
152  }
153  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_33_int1: version %i is not valid.", InterAmpFlag);}
154  }
155  return total;
156 }
157 
158 static double IMRPhenomXHM_Inter_Amp_33_int2(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
159  double total=0;
160  switch (InterAmpFlag){
161  case 122018:{
162  double eta = pWF->eta;
163  double S = pWF->STotR;
164  double eta2,eta6;
165  eta2 = pow(eta,2);
166  eta6 = pow(eta,6);
167  double noSpin = sqrt(eta - 4.*eta2)*(20.162169689041903 - 18.666422946967764*eta + 53.04107631052987*eta2);
168  double eqSpin = sqrt(eta - 4.*eta2)*S*(3.896260108714186 + eta*(-33.707998325000965 - 61.1244771077077*S) + 4.878506403725656*S + eta2*(91.31681057861915 + 196.40535070402336*S));
169  double uneqSpin = 1637.4256048973248*(pWF->dchi)*eta6;
170  total = noSpin + eqSpin + uneqSpin;
171  break;
172  }
173  case 122022:{
174  double eta = pWF->eta;
175  double delta = pWF->delta;
176  double S = pWF->STotR;
177  double chidiff = pWF->dchi_half;
178  double eta1 = eta;
179  double eta2 = eta1 * eta1;
180  double eta3 = eta1 * eta2;
181  double eta4 = eta1 * eta3;
182  double eta5 = eta1 * eta4;
183  double S1 = S;
184  double S2 = S1 * S1;
185  double chidiff1 = chidiff;
186  total = delta*eta1*(17.42562079069636 - 28.970875603981295*eta1 + 50.726220750178435*eta2) + chidiff1*delta*eta1*(-7.861956897615623*eta1 + 93.45476935080045*eta2 - 273.1170921735085*eta3) + chidiff1*delta*eta1*(-0.3265505633310564*eta1 - 9.861644053348053*eta2 + 60.38649425562178*eta3)*S1 + chidiff1*eta5*(234.13476431269862 + 51.2153901931183*S1 - 10.05114600643587*S2) + delta*eta1*S1*(0.3104472390387834*(6.073591341439855 + 169.85423386969634*eta1 - 4964.199967099143*eta2 + 42566.59565666228*eta3 - 154255.3408672655*eta4 + 205525.13910847943*eta5) + 0.2295327944679772*(19.236275867648594 - 354.7914372697625*eta1 + 1876.408148917458*eta2 + 2404.4151687877525*eta3 - 41567.07396803811*eta4 + 79210.33893514868*eta5)*S1 + 0.30983324991828787*(11.302200127272357 - 719.9854052004307*eta1 + 13278.047199998868*eta2 - 104863.50453518033*eta3 + 376409.2335857397*eta4 - 504089.07690692553*eta5)*S2);
187  break;
188  }
189  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_33_int2: version %i is not valid.", InterAmpFlag);}
190  }
191  return total;
192 }
193 
194 static double IMRPhenomXHM_Inter_Amp_32_int1(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
195  double total=0;
196  switch (InterAmpFlag){
197  case 122018:{
198  double eta = pWF->eta;
199  double S = pWF->STotR;
200  double delta=sqrt(1.-4*eta),eta2,eta3,eta4,eta5,eta6;
201  eta2 = pow(eta,2);
202  eta3 = pow(eta,3);
203  eta4 = pow(eta,4);
204  eta5 = pow(eta,5);
205  eta6 = pow(eta,6);
206  double noSpin = sqrt(eta - 3.*eta2)*(6.523612598187996 - 56.93956111746338*eta + 1021.6414686597869*eta2 - 12107.114370361525*eta3 + 76320.90587515048*eta4 - 244144.92645448362*eta5 + 321790.55131499085*eta6);
207  double eqSpin = sqrt(eta - 3.*eta2)*S*(2.9649243713119895 + eta3*(1790.8363334078751 - 5438.911035114849*S) + eta*(-37.87005271181108 - 126.1263286618178*S) + 4.063724538613828*S + eta2*(48.39743086535961 + 1341.2619677741804*S) + eta4*(-5200.659417644607 + 7369.386205324284*S));
208  double uneqSpin = eta2*(-0.4386152975075188*(pow(pWF->chi1L,2) - 2.*pWF->chi1L*pWF->chi2L + pow(pWF->chi2L,2)) + (pWF->chi2L*(3.6527252109313233 - 7.324266404418883*S) + pWF->chi1L*(-3.6527252109313233 + 7.324266404418883*S))*delta);
209  total = noSpin + eqSpin + uneqSpin;
210  break;
211  }
212  case 122022:{
213  double eta = pWF->eta;
214  double delta = pWF->delta;
215  double sqroot = sqrt(eta);
216  double S = pWF->chiPNHat;
217  double chidiff = pWF->dchi_half;
218  double eta1 = eta;
219  double eta2 = eta1 * eta1;
220  double eta3 = eta1 * eta2;
221  double eta4 = eta1 * eta3;
222  double eta5 = eta1 * eta4;
223  double S1 = S;
224  double chidiff1 = chidiff;
225  double chidiff2 = chidiff1 * chidiff1;
226  total = (chidiff2*(-0.2341404256829785*eta1 + 2.606326837996192*eta2 - 8.68296921440857*eta3) + chidiff1*delta*(0.5454562486736877*eta1 - 25.19759222940851*eta2 + 73.40268975811729*eta3))*sqroot + chidiff1*delta*(0.4422257616009941*eta1 - 8.490112284851655*eta2 + 32.22238925527844*eta3)*S1*sqroot + S1*(0.7067243321652764*(0.12885110296881636 + 9.608999847549535*eta1 - 85.46581740280585*eta2 + 325.71940024255775*eta3 + 175.4194342269804*eta4 - 1929.9084724384807*eta5) + 0.1540566313813899*(-0.3261041495083288 + 45.55785402900492*eta1 - 827.591235943271*eta2 + 7184.647314370326*eta3 - 28804.241518798244*eta4 + 43309.69769878964*eta5)*S1)*sqroot + (480.0434256230109*eta1 + 25346.341240810478*eta2 - 99873.4707358776*eta3 + 106683.98302194536*eta4)*sqroot*pow(1 + 1082.6574834474493*eta1 + 10083.297670051445*eta2,-1);
227  break;
228  }
229  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_32_int1: version %i is not valid.", InterAmpFlag);}
230  }
231  return total;
232 }
233 
234 static double IMRPhenomXHM_Inter_Amp_32_int2(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
235  double total=0;
236  switch (InterAmpFlag){
237  case 122018:{
238  double eta = pWF->eta;
239  double S = pWF->STotR;
240  double delta=sqrt(1.-4.*eta),eta2,eta3,eta4,S2;
241  eta2 = pow(eta,2);
242  eta3 = pow(eta,3);
243  eta4 = pow(eta,4);
244  S2 = pow(S,2);
245  double noSpin = sqrt(eta - 3.*eta2)*(5.941845842405418 - 31.905244419036794*eta + 271.105632998832*eta2 - 2113.9652334868965*eta3 + 6214.038393898584*eta4);
246  double eqSpin = sqrt(eta - 3.*eta2)*S*(-2.726472456645038 + 2.9454485454761827*S + eta3*(10581.664858726683 - 8474.190197512324*S - 11680.937129551317*S2) + eta*(98.08119212251981 - 119.88112323140916*S - 145.5079981415436*S2) + 3.5684571473795095*S2 + eta2*(-1595.8027347570667 + 1686.7137359336039*S + 2139.8290160628144*S2) + eta4*(-21488.25117198268 + 13866.428366595079*S + 20863.270079587106*S2));
247  double uneqSpin = 0.0038732029045487884*(pWF->dchi)*eta2*delta;
248  total = noSpin + eqSpin + uneqSpin;
249  break;
250  }
251  case 122022:{
252  double eta = pWF->eta;
253  double delta = pWF->delta;
254  double S = pWF->STotR;
255  double chidiff = pWF->dchi_half;
256  double eta1 = eta;
257  double eta2 = eta1 * eta1;
258  double eta3 = eta1 * eta2;
259  double eta4 = eta1 * eta3;
260  double eta5 = eta1 * eta4;
261  double eta6 = eta1 * eta5;
262  double S1 = S;
263  double chidiff1 = chidiff;
264  double chidiff2 = chidiff1 * chidiff1;
265  total = eta1*(chidiff2*(-4.175680729484314*eta1 + 47.54281549129226*eta2 - 128.88334273588077*eta3) + chidiff1*delta*(-0.18274358639599947*eta1 - 71.01128541687838*eta2 + 208.07105580635888*eta3)) + eta1*(4.760999387359598 - 38.57900689641654*eta1 + 456.2188780552874*eta2 - 4544.076411013166*eta3 + 24956.9592553473*eta4 - 69430.10468748478*eta5 + 77839.74180254337*eta6) + chidiff1*delta*eta1*(1.2198776533959694*eta1 - 26.816651899746475*eta2 + 68.72798751937934*eta3)*S1 + eta1*S1*(1.5098291294292217*(0.4844667556328104 + 9.848766999273414*eta1 - 143.66427232396376*eta2 + 856.9917885742416*eta3 - 1633.3295758142904*eta4) + 0.32413108737204144*(2.835358206961064 - 62.37317183581803*eta1 + 761.6103793011912*eta2 - 3811.5047139343505*eta3 + 6660.304740652403*eta4)*S1);
266  break;
267  }
268  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_32_int2: version %i is not valid.", InterAmpFlag);}
269  }
270  return total;
271 }
272 
273 static double IMRPhenomXHM_Inter_Amp_44_int1(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
274  double total=0;
275  switch (InterAmpFlag){
276  case 122018:{
277  double eta = pWF->eta;
278  double S = pWF->STotR;
279  double delta=sqrt(1.-4.*eta),eta2,eta3,eta4;
280  eta2 = pow(eta,2);
281  eta3 = pow(eta,3);
282  eta4 = pow(eta,4);
283  double noSpin = sqrt(eta - 3.*eta2)*(10.804555518381166 - 72.3834734399584*eta + 540.0541240482852*eta2 - 2612.999845214264*eta3 + 4779.096001663427*eta4);
284  double eqSpin = sqrt(eta - 3.*eta2)*S*(4.26336253142121 + eta*(-47.94914754514519 - 39.31284390368824*S) + 3.0973959822174297*S + eta2*(119.70401520575753 + 106.91295627237112*S));
285  double uneqSpin = 0.7262636326998003*pow(pWF->dchi,2)*(1. - 4.*eta)*eta + 3.001401833124412*(pWF->dchi)*eta2*delta;
286  total = noSpin + eqSpin + uneqSpin;
287  break;
288  }
289  case 122022:{
290  double eta = pWF->eta;
291  double delta = pWF->delta;
292  double S = pWF->chiPNHat;
293  double chidiff = pWF->dchi_half;
294  double eta1 = eta;
295  double eta2 = eta1 * eta1;
296  double eta3 = eta1 * eta2;
297  double eta4 = eta1 * eta3;
298  double S1 = S;
299  double S2 = S1 * S1;
300  double chidiff1 = chidiff;
301  double chidiff2 = chidiff1 * chidiff1;
302  total = eta1*(chidiff1*delta*(1.5378890240544967*eta1 - 3.4499418893734903*eta2 + 16.879953490422782*eta3) + chidiff2*(1.720226708214248*eta1 - 11.87925165364241*eta2 + 23.259283336239545*eta3)) + eta1*(8.790173464969538 - 64.95499142822892*eta1 + 324.1998823562892*eta2 - 1111.9864921907126*eta3 + 1575.602443847111*eta4) + eta1*S1*(-0.062333275821238224*(-21.630297087123807 + 137.4395894877131*eta1 + 64.92115530780129*eta2 - 1013.1110639471394*eta3) - 0.11014697070998722*(4.149721483857751 - 108.6912882442823*eta1 + 831.6073263887092*eta2 - 1828.2527520190122*eta3)*S1 - 0.07704777584463054*(4.581767671445529 - 50.35070009227704*eta1 + 344.9177692251726*eta2 - 858.9168637051405*eta3)*S2);
303  break;
304  }
305  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_44_int1: version %i is not valid.", InterAmpFlag);}
306  }
307  return total;
308 }
309 
310 static double IMRPhenomXHM_Inter_Amp_44_int2(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
311  double total=0;
312  switch (InterAmpFlag){
313  case 122018:{
314  double eta = pWF->eta;
315  double S = pWF->STotR;
316  double delta=sqrt(1.-4.*eta),eta2,eta3,eta4;
317  eta2 = pow(eta,2);
318  eta3 = pow(eta,3);
319  eta4 = pow(eta,4);
320  double noSpin = sqrt(eta - 3.*eta2)*(9.020721305469884 - 53.221883492311235*eta + 508.07176447172264*eta2 - 3194.0620894511508*eta3 + 6769.9274392345915*eta4);
321  double eqSpin = sqrt(eta - 3.*eta2)*S*(3.256591670091969 + eta*(-38.38922554651356 - 25.286684856422735*S) + 2.374434219852751*S + eta2*(96.41777041220982 + 64.74544118094362*S));
322  double uneqSpin = 3.2337593375595417*(pWF->dchi)*eta2*delta;
323  total = noSpin + eqSpin + uneqSpin;
324  break;
325  }
326  case 122022:{
327  double eta = pWF->eta;
328  double delta = pWF->delta;
329  double S = pWF->chiPNHat;
330  double chidiff = pWF->dchi_half;
331  double eta1 = eta;
332  double eta2 = eta1 * eta1;
333  double eta3 = eta1 * eta2;
334  double eta4 = eta1 * eta3;
335  double eta5 = eta1 * eta4;
336  double eta6 = eta1 * eta5;
337  double S1 = S;
338  double chidiff1 = chidiff;
339  double chidiff2 = chidiff1 * chidiff1;
340  total = eta1*(chidiff1*delta*(2.3123974306694057*eta1 - 12.237594841284904*eta2 + 44.78225529547671*eta3) + chidiff2*(2.9282931698944292*eta1 - 25.624210264341933*eta2 + 61.05270871360041*eta3)) + eta1*(6.98072197826729 - 46.81443520117986*eta1 + 236.76146303619544*eta2 - 920.358408667518*eta3 + 1478.050456337336*eta4) + eta1*S1*(-0.07801583359561987*(-28.29972282146242 + 752.1603553640072*eta1 - 10671.072606753183*eta2 + 83447.0461509547*eta3 - 350025.2112501252*eta4 + 760889.6919776166*eta5 - 702172.2934567826*eta6) + 0.013159545629626014*(91.1469833190294 - 3557.5003799977294*eta1 + 52391.684517955284*eta2 - 344254.9973814295*eta3 + 1.0141877915334814e6*eta4 - 1.1505186449682908e6*eta5 + 268756.85659532435*eta6)*S1);
341  break;
342  }
343  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_44_int2: version %i is not valid.", InterAmpFlag);}
344  }
345  return total;
346 }
347 
348 /*
349 Fits for the extra collocation point for EMR cases with 2 intermediate regions
350 */
351 static double IMRPhenomXHM_Inter_Amp_21_int0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
352  double total=0;
353  switch (InterAmpFlag){
354  case 122018:{
355  double eta = pWF->eta;
356  double S = pWF->STotR;
357  double eta2,eta3;
358  eta2 = pow(eta,2);
359  eta3 = pow(eta,3);
360  double noSpin = 0.872895771366973 + 441.76285124642845*eta - 24617.068739152524*eta2 + 518054.9485981792*eta3;
361  double eqSpin = S*(-0.0720494539485585 + eta*(-173.67847091983123 - 113.29725582509889*S) - 0.2687302438646897*S + eta2*(3571.0393588230045 + 2640.919925429635*S));
362  double uneqSpin = 0.;
363  total = noSpin + eqSpin + uneqSpin;
364  break;
365  }
366  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_21_int0: version is not valid.");}
367  }
368  return total;
369 }
370 
371 static double IMRPhenomXHM_Inter_Amp_21_dint0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
372  double total=0;
373  switch (InterAmpFlag){
374  case 122018:{
375  double eta = pWF->eta;
376  double S = pWF->STotR;
377  double eta2,eta3;
378  eta2 = pow(eta,2);
379  eta3 = pow(eta,3);
380  double noSpin = -0.8535048463050732 - 93.1876950411214*eta + 13641.071903017495*eta2 - 337621.44851304166*eta3;
381  double eqSpin = S*(-1.2067842398131878 + eta2*(-1972.284151572111 - 8172.057025783849*S) - 0.26539816223182355*S + eta*(77.26350785961219 + 189.63365484152857*S));
382  double uneqSpin = 0.;
383  total = noSpin + eqSpin + uneqSpin;
384  break;
385  }
386  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_21_dint0: version %i is not valid.", InterAmpFlag);}
387  }
388  return total;
389 }
390 
391 static double IMRPhenomXHM_Inter_Amp_33_int0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
392  double total=0;
393  switch (InterAmpFlag){
394  case 122018:{
395  double eta = pWF->eta;
396  double S = pWF->STotR;
397  double eta2,eta3;
398  eta2 = pow(eta,2);
399  eta3 = pow(eta,3);
400  double noSpin = 1.5852399637975103 + 549.5183711492834*eta - 34257.76380246282*eta2 + 743142.8286902909*eta3;
401  double eqSpin = S*(0.7436306553052219 + eta*(-89.49451655594787 - 174.5730646548662*S) + 0.4253024979725725*S + eta2*(1185.1654325913717 + 6510.983041407191*S));
402  double uneqSpin = 0.;
403  total = noSpin + eqSpin + uneqSpin;
404  break;
405  }
406  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_33_int0: version %i is not valid.", InterAmpFlag);}
407  }
408  return total;
409 }
410 
411 static double IMRPhenomXHM_Inter_Amp_33_dint0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
412  double total=0;
413  switch (InterAmpFlag){
414  case 122018:{
415  double eta = pWF->eta;
416  double S = pWF->STotR;
417  double eta2,eta3;
418  eta2 = pow(eta,2);
419  eta3 = pow(eta,3);
420  double noSpin = -4.691600252198376 + 101.4338937535679*eta + 9262.994550540048*eta2 - 310993.1309846956*eta3;
421  double eqSpin = S*(-4.198232394219111 + eta2*(-28714.904192060643 - 5100.09336069277*S) - 0.40986595512314733*S + eta*(734.7118618746317 + 292.04566260701574*S));
422  double uneqSpin = 0.;
423  total = noSpin + eqSpin + uneqSpin;
424  break;
425  }
426  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_33_dint0: version %i is not valid.", InterAmpFlag);}
427  }
428  return total;
429 }
430 
431 static double IMRPhenomXHM_Inter_Amp_32_int0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
432  double total=0;
433  switch (InterAmpFlag){
434  case 122018:{
435  double eta = pWF->eta;
436  double S = pWF->STotR;
437  double eta2,eta3,S2,S3;
438  eta2 = pow(eta,2);
439  eta3 = pow(eta,3);
440  S2 = pow(S,2);
441  S3 = pow(S,3);
442  double noSpin = 0.24794156582503746 + 115.81823862983131*eta - 6626.167995915723*eta2 + 141004.29332593994*eta3;
443  double eqSpin = (0.21144389781375486 + 35.10041265469983*eta - 1794.2301585086836*eta2)*S + (0.2781735549493081 - 37.038950686633*eta + 1258.628375238807*eta2)*S2 + (0.23428222791962147 - 63.98011009365723*eta + 2118.213562899934*eta2)*S3;
444  double uneqSpin = 0.;
445  total = noSpin + eqSpin + uneqSpin;
446  break;
447  }
448  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_32_int0: version %i is not valid.", InterAmpFlag);}
449  }
450  return total;
451 }
452 
453 static double IMRPhenomXHM_Inter_Amp_32_dint0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
454  double total=0;
455  switch (InterAmpFlag){
456  case 122018:{
457  double eta = pWF->eta;
458  double S = pWF->STotR;
459  double eta2,eta3;
460  eta2 = pow(eta,2);
461  eta3 = pow(eta,3);
462  double noSpin = -0.3391808620221253 - 14.604141885467747*eta + 3694.1706648870427*eta2 - 95482.02951271653*eta3;
463  double eqSpin = S*(-1.2844502090793946 + eta2*(-5018.762853306415 - 6332.389157828062*S) - 1.2356159239385598*S + eta*(149.04865679660233 + 188.2052849646003*S));
464  double uneqSpin = 0.;
465  total = noSpin + eqSpin + uneqSpin;
466  break;
467  }
468  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_32_dint0: version %i is not valid.", InterAmpFlag);}
469  }
470  return total;
471 }
472 
473 static double IMRPhenomXHM_Inter_Amp_44_int0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
474  double total=0;
475  switch (InterAmpFlag){
476  case 122018:{
477  double eta = pWF->eta;
478  double S = pWF->STotR;
479  double eta2,eta3,S2,S3;
480  eta2 = pow(eta,2);
481  eta3 = pow(eta,3);
482  S2 = pow(S,2);
483  S3 = pow(S,3);
484  double noSpin = 0.5664660641971224 + 185.58965113823874*eta - 11458.768824989507*eta2 + 249386.7511724409*eta3;
485  double eqSpin = (0.1741768776210781 - 9.365114803167128*eta + 703.2622732011035*eta2)*S + (0.20169229783048184 - 62.13147149352512*eta + 2833.5738711424974*eta2)*S2 + (0.4423803798742513 - 23.60535149579996*eta - 994.9241585715828*eta2)*S3;
486  double uneqSpin = 0.;
487  total = noSpin + eqSpin + uneqSpin;
488  break;
489  }
490  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_44_int0: version %i is not valid.", InterAmpFlag);}
491  }
492  return total;
493 }
494 
495 static double IMRPhenomXHM_Inter_Amp_44_dint0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag) {
496  double total=0;
497  switch (InterAmpFlag){
498  case 122018:{
499  double eta = pWF->eta;
500  double S = pWF->STotR;
501  double eta2 = pow(eta,2);
502  double noSpin = -1.796444922382065 + 111.51170611049032*eta - 1728.7493675776548*eta2;
503  double eqSpin = S*(-1.842119860613924 + eta2*(-11235.484645624338 - 2927.019210835522*S) - 0.36655273031432567*S + eta*(312.34531117524097 + 128.64488103364167*S));
504  double uneqSpin = 0.;
505  total = noSpin + eqSpin + uneqSpin;
506  break;
507  }
508  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_44_dint0: version %i is not valid.", InterAmpFlag);}
509  }
510  return total;
511 }
512 
513 
514 static double IMRPhenomXHM_Inter_Amp_21_int3(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag){
515  double total=0;
516  switch (InterAmpFlag){
517  case 122022:{
518  double eta = pWF->eta;
519  double delta = pWF->delta;
520  double S = pWF->STotR;
521  double chidiff = pWF->dchi_half;
522  double eta1 = eta;
523  double eta2 = eta1 * eta1;
524  double eta3 = eta1 * eta2;
525  double eta4 = eta1 * eta3;
526  double eta5 = eta1 * eta4;
527  double S1 = S;
528  double S2 = S1 * S1;
529  double chidiff1 = chidiff;
530  total = fabs(delta*eta1*(13.318990196097973 - 21.755549987331054*eta1 + 76.14884211156267*eta2 - 127.62161159798488*eta3) + chidiff1*delta*eta1*(17.704321326939414*eta1 - 434.4390350012534*eta2 + 1366.2408490833282*eta3) + chidiff1*delta*eta1*(11.877985158418596*eta1 - 131.04937626836355*eta2 + 343.79587860999874*eta3)*S1 + chidiff1*eta5*(-1522.8543551416456 - 16.639896279650678*S1 + 3.0053086651515843*S2) + delta*eta1*S1*(-8.665646058245033*(0.7862132291286934 + 8.293609541933655*eta1 - 111.70764910503321*eta2 + 576.7172598056907*eta3 - 1001.2370065269745*eta4) - 0.9459820574514348*(1.309016452198605 + 48.94077040282239*eta1 - 817.7854010574645*eta2 + 4331.56002883546*eta3 - 7518.309520232795*eta4)*S1 - 0.4308267743835775*(9.970654092010587 - 302.9708323417439*eta1 + 3662.099161055873*eta2 - 17712.883990278668*eta3 + 29480.158198408903*eta4)*S2));
531  break;
532  }
533  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_21_int3:version %i is not valid.", InterAmpFlag);}
534  }
535  return total;
536 }
537 
538 static double IMRPhenomXHM_Inter_Amp_21_int4(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag){
539  double total=0;
540  switch (InterAmpFlag){
541  case 122022:{
542  double eta = pWF->eta;
543  double delta = pWF->delta;
544  double S = pWF->STotR;
545  double chidiff = pWF->dchi_half;
546  double eta1 = eta;
547  double eta2 = eta1 * eta1;
548  double eta3 = eta1 * eta2;
549  double eta4 = eta1 * eta3;
550  double eta5 = eta1 * eta4;
551  double S1 = S;
552  double S2 = S1 * S1;
553  double chidiff1 = chidiff;
554  total = fabs(delta*eta1*(13.094382343446163 - 22.831152256559523*eta1 + 83.20619262213437*eta2 - 139.25546924151664*eta3) + chidiff1*delta*eta1*(20.120192352555357*eta1 - 458.2592421214168*eta2 + 1430.3698681181*eta3) + chidiff1*delta*eta1*(12.925363020014743*eta1 - 126.87194512915104*eta2 + 280.6003655502327*eta3)*S1 + chidiff1*eta5*(-1528.956015503355 + 74.44462583487345*S1 - 2.2456928156392197*S2) + delta*eta1*S1*(-9.499741513411829*(0.912120958549489 + 2.400945118514037*eta1 - 33.651192908287236*eta2 + 166.04254881175257*eta3 - 248.5050377498615*eta4) - 0.7850652143322492*(1.534131218043425 + 60.81773903539479*eta1 - 1032.1319480683567*eta2 + 5381.481380750608*eta3 - 9077.037917192794*eta4)*S1 - 0.21540359093306097*(9.42805409480658 - 109.06544597367301*eta1 + 385.8345793110262*eta2 + 1889.9613367802453*eta3 - 9835.416414460055*eta4)*S2));
555  break;
556  }
557  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_21_int4:version %i is not valid.", InterAmpFlag);}
558  }
559  return total;
560 }
561 
562 static double IMRPhenomXHM_Inter_Amp_33_int3(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag){
563  double total=0;
564  switch (InterAmpFlag){
565  case 122022:{
566  double eta = pWF->eta;
567  double delta = pWF->delta;
568  double S = pWF->STotR;
569  double chidiff = pWF->dchi_half;
570  double eta1 = eta;
571  double eta2 = eta1 * eta1;
572  double eta3 = eta1 * eta2;
573  double eta4 = eta1 * eta3;
574  double eta5 = eta1 * eta4;
575  double S1 = S;
576  double S2 = S1 * S1;
577  double chidiff1 = chidiff;
578  total = delta*eta1*(14.555522136327964 - 12.799844096694798*eta1 + 16.79500349318081*eta2) + chidiff1*delta*eta1*(-16.292654447108134*eta1 + 190.3516012682791*eta2 - 562.0936797781519*eta3) + chidiff1*delta*eta1*(-7.048898856045782*eta1 + 49.941617405768135*eta2 - 73.62033985436068*eta3)*S1 + chidiff1*eta5*(263.5151703818307 + 44.408527093031566*S1 + 10.457035444964653*S2) + delta*eta1*S1*(0.4590550434774332*(3.0594364612798635 + 207.74562213604057*eta1 - 5545.0086137386525*eta2 + 50003.94075934942*eta3 - 195187.55422847517*eta4 + 282064.174913521*eta5) + 0.657748992123043*(5.57939137343977 - 124.06189543062042*eta1 + 1276.6209573025596*eta2 - 6999.7659193505915*eta3 + 19714.675715229736*eta4 - 20879.999628681435*eta5)*S1 + 0.3695850566805098*(6.077183107132255 - 498.95526910874986*eta1 + 10426.348944657859*eta2 - 91096.64982858274*eta3 + 360950.6686625352*eta4 - 534437.8832860565*eta5)*S2);
579  break;
580  }
581  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_33_int3:version %i is not valid.", InterAmpFlag);}
582  }
583  return total;
584 }
585 
586 static double IMRPhenomXHM_Inter_Amp_33_int4(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag){
587  double total=0;
588  switch (InterAmpFlag){
589  case 122022:{
590  double eta = pWF->eta;
591  double delta = pWF->delta;
592  double S = pWF->STotR;
593  double chidiff = pWF->dchi_half;
594  double eta1 = eta;
595  double eta2 = eta1 * eta1;
596  double eta3 = eta1 * eta2;
597  double eta4 = eta1 * eta3;
598  double eta5 = eta1 * eta4;
599  double S1 = S;
600  double S2 = S1 * S1;
601  double chidiff1 = chidiff;
602  double chidiff2 = chidiff1 * chidiff1;
603  total = delta*eta1*(13.312095699772305 - 7.449975618083432*eta1 + 17.098576301150125*eta2) + delta*eta1*(chidiff1*(-31.171150896110156*eta1 + 371.1389274783572*eta2 - 1103.1917047361735*eta3) + chidiff2*(32.78644599730888*eta1 - 395.15713118955387*eta2 + 1164.9282236341376*eta3)) + chidiff1*delta*eta1*(-46.85669289852532*eta1 + 522.3965959942979*eta2 - 1485.5134187612182*eta3)*S1 + chidiff1*eta5*(287.90444670305715 - 21.102665129433042*chidiff2 + 7.635582066682054*S1 - 29.471275170013012*S2) + delta*eta1*S1*(0.6893003654021495*(3.1014226377197027 - 44.83989278653052*eta1 + 565.3767256471909*eta2 - 4797.429130246123*eta3 + 19514.812242035154*eta4 - 27679.226582207506*eta5) + 0.7068016563068026*(4.071212304920691 - 118.51094098279343*eta1 + 1788.1730303291356*eta2 - 13485.270489656365*eta3 + 48603.96661003743*eta4 - 65658.74746265226*eta5)*S1 + 0.2181399561677432*(-1.6754158383043574 + 303.9394443302189*eta1 - 6857.936471898544*eta2 + 59288.71069769708*eta3 - 216137.90827404748*eta4 + 277256.38289831823*eta5)*S2);
604  break;
605  }
606  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_33_int4:version %i is not valid.", InterAmpFlag);}
607  }
608  return total;
609 }
610 
611 static double IMRPhenomXHM_Inter_Amp_44_int3(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag){
612  double total=0;
613  switch (InterAmpFlag){
614  case 122022:{
615  double eta = pWF->eta;
616  double delta = pWF->delta;
617  double S = pWF->chiPNHat;
618  double chidiff = pWF->dchi_half;
619  double eta1 = eta;
620  double eta2 = eta1 * eta1;
621  double eta3 = eta1 * eta2;
622  double eta4 = eta1 * eta3;
623  double S1 = S;
624  double S2 = S1 * S1;
625  double chidiff1 = chidiff;
626  double chidiff2 = chidiff1 * chidiff1;
627  total = eta1*(chidiff1*delta*(-0.8765502142143329*eta1 + 22.806632458441996*eta2 - 43.675503209991184*eta3) + chidiff2*(0.48698617426180074*eta1 - 4.302527065360426*eta2 + 16.18571810759235*eta3)) + eta1*(6.379772583015967 - 44.10631039734796*eta1 + 269.44092930942793*eta2 - 1285.7635006711453*eta3 + 2379.538739132234*eta4) + eta1*S1*(-0.23316184683282615*(-1.7279023138971559 - 23.606399143993716*eta1 + 409.3387618483284*eta2 - 1115.4147472977265*eta3) - 0.09653777612560172*(-5.310643306559746 - 2.1852511802701264*eta1 + 541.1248219096527*eta2 - 1815.7529908827103*eta3)*S1 - 0.060477799540741804*(-14.578189130145661 + 175.6116682068523*eta1 - 569.4799973930861*eta2 + 426.0861915646515*eta3)*S2);
628  break;
629  }
630  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_44_int3:version %i is not valid.", InterAmpFlag);}
631  }
632  return total;
633 }
634 
635 static double IMRPhenomXHM_Inter_Amp_44_int4(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag){
636  double total=0;
637  switch (InterAmpFlag){
638  case 122022:{
639  double eta = pWF->eta;
640  double delta = pWF->delta;
641  double S = pWF->chiPNHat;
642  double chidiff = pWF->dchi_half;
643  double eta1 = eta;
644  double eta2 = eta1 * eta1;
645  double eta3 = eta1 * eta2;
646  double eta4 = eta1 * eta3;
647  double eta5 = eta1 * eta4;
648  double S1 = S;
649  double S2 = S1 * S1;
650  double chidiff1 = chidiff;
651  double chidiff2 = chidiff1 * chidiff1;
652  total = eta1*(chidiff1*delta*(-2.461738962276138*eta1 + 45.3240543970684*eta2 - 112.2714974622516*eta3) + chidiff2*(0.9158352037567031*eta1 - 8.724582331021695*eta2 + 28.44633544874233*eta3)) + eta1*(6.098676337298138 - 45.42463610529546*eta1 + 350.97192927929433*eta2 - 2002.2013283876834*eta3 + 4067.1685640401033*eta4) + eta1*S1*(-0.36068516166901304*(-2.120354236840677 - 47.56175350408845*eta1 + 1618.4222330016048*eta2 - 14925.514654896673*eta3 + 60287.45399959349*eta4 - 91269.3745059139*eta5) - 0.09635801207669747*(-11.824692837267394 + 371.7551657959369*eta1 - 4176.398139238679*eta2 + 16655.87939259747*eta3 - 4102.218189945819*eta4 - 67024.98285179552*eta5)*S1 - 0.06565232123453196*(-26.15227471380236 + 1869.0168486099005*eta1 - 33951.35186039629*eta2 + 253694.6032002248*eta3 - 845341.6001856657*eta4 + 1.0442282862506858e6*eta5)*S2);
653  break;
654  }
655  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_44_int4:version %i is not valid.", InterAmpFlag);}
656  }
657  return total;
658 }
659 
660 static double IMRPhenomXHM_Inter_Amp_32_int3(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag){
661  double total=0;
662  switch (InterAmpFlag){
663  case 122022:{
664  double eta = pWF->eta;
665  double delta = pWF->delta;
666  double S = pWF->chiPNHat;
667  double chidiff = pWF->dchi_half;
668  double eta1 = eta;
669  double eta2 = eta1 * eta1;
670  double eta3 = eta1 * eta2;
671  double eta4 = eta1 * eta3;
672  double eta5 = eta1 * eta4;
673  double eta6 = eta1 * eta5;
674  double S1 = S;
675  double S2 = S1 * S1;
676  double chidiff1 = chidiff;
677  double chidiff2 = chidiff1 * chidiff1;
678  total = 3.881450518842405*eta1 - 12.580316392558837*eta2 + 1.7262466525848588*eta3 + chidiff2*(-7.065118823041031*eta2 + 77.97950589523865*eta3 - 203.65975422378446*eta4) - 58.408542930248046*eta4 + chidiff1*delta*(1.924723094787216*eta2 - 90.92716917757797*eta3 + 387.00162600306226*eta4) + 403.5748987560612*eta5 + chidiff1*delta*(-0.2566958540737833*eta2 + 14.488550203412675*eta3 - 26.46699529970884*eta4)*S1 + S1*(0.3650871458400108*(71.57390929624825*eta2 - 994.5272351916166*eta3 + 6734.058809060536*eta4 - 18580.859291282686*eta5 + 16001.318492586077*eta6) + 0.0960146077440495*(451.74917589707513*eta2 - 9719.470997418284*eta3 + 83403.5743434538*eta4 - 318877.43061174755*eta5 + 451546.88775684836*eta6)*S1 - 0.03985156529181297*(-304.92981902871617*eta2 + 3614.518459296278*eta3 - 7859.4784979916085*eta4 - 46454.57664737511*eta5 + 162398.81483375572*eta6)*S2);
679  break;
680  }
681  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_32_int3:version %i is not valid.", InterAmpFlag);}
682  }
683  return total;
684 }
685 
686 static double IMRPhenomXHM_Inter_Amp_32_int4(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag){
687  double total=0;
688  switch (InterAmpFlag){
689  case 122022:{
690  double eta = pWF->eta;
691  double delta = pWF->delta;
692  double S = pWF->STotR;
693  double chidiff = pWF->dchi_half;
694  double eta1 = eta;
695  double eta2 = eta1 * eta1;
696  double eta3 = eta1 * eta2;
697  double eta4 = eta1 * eta3;
698  double S1 = S;
699  double chidiff1 = chidiff;
700  double chidiff2 = chidiff1 * chidiff1;
701  total = eta1*(chidiff2*(-8.572797326909152*eta1 + 92.95723645687826*eta2 - 236.2438921965621*eta3) + chidiff1*delta*(6.674358856924571*eta1 - 171.4826985994883*eta2 + 645.2760206304703*eta3)) + eta1*(3.921660532875504 - 16.57299637423352*eta1 + 25.254017911686333*eta2 - 143.41033155133266*eta3 + 692.926425981414*eta4) + chidiff1*delta*eta1*(-3.582040878719185*eta1 + 57.75888914133383*eta2 - 144.21651114700492*eta3)*S1 + eta1*S1*(1.242750265695504*(-0.522172424518215 + 25.168480118950065*eta1 - 303.5223688400309*eta2 + 1858.1518762309654*eta3 - 3797.3561904195085*eta4) + 0.2927045241764365*(0.5056957789079993 - 15.488754837330958*eta1 + 471.64047356915603*eta2 - 3131.5783196211587*eta3 + 6097.887891566872*eta4)*S1);
702  break;
703  }
704  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Amp_32_int4:version %i is not valid.", InterAmpFlag);}
705  }
706  return total;
707 }
708 
709 /* End of Amp Parameter Space Fits */
710 
711 
712 /* Solves system of equations for 5th order polynomial ansatz */
713 
714 /* We use a 5th order polynomial to connect the inspiral and ringdown regions.
715 This polynomial is built such that it has the same value and derivative at the boundaries of the
716 intermediate region than the inspiral and ringdown part (demanding continuity for the function
717 and its first derivative) and also crosses through the two intermediate collocation points.
718 Sometimes we will not use the 5th order but a lower one to assure a better behaviour.
719 
720 Now we have the functions to get the coefficients of such a polynomial:
721  delta0 + delta1*f + delta2*f^2 + delta3*f^3 + delta4*f^4 + delta5*f^5
722 
723 The different cases inside one function correspond to the different ways of building the polynomial,
724 and it depend on the number of collocation points and derivatives we are using.
725 For example case:105 correspond to the 5th order polynomial we have described before.
726 */
727 
728 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)
729 {
730  double retVal;
731 
732  switch (IntAmpFlag)
733  {
734  case 101: //linear, only v1, v2
735  {
736  double f1mf4 = f1-f4;
737 
738  retVal = (-(f4*v1) + f1*v4)/f1mf4;
739  break;
740  }
741  case 102: //quadratic: v1, v2, d2
742  {
743  double f12 = f1*f1;
744  double f42 = f4*f4;
745  double f1mf4 = f1-f4;
746  double f1mf42 = f1mf4*f1mf4;
747 
748  retVal = (-(d4*f1*f1mf4*f4) + f42*v1 + f12*v4 - 2*f1*f4*v4)/f1mf42;
749  break;
750  }
751  case 1032: // 2 freqs, points and derivatives: v1, v4, d1, d4
752  {
753  double f12 = f1*f1;
754  double f13 = f12*f1;
755  double f42 = f4*f4;
756  double f43 = f42*f4;
757 
758  double f1mf4 = f1-f4;
759  double f1mf42 = f1mf4*f1mf4;
760  double f1mf43 = f1mf42*f1mf4;
761 
762  retVal = (d4*f12*f4*(-f1 + f4) + d1*f1*(-f1 + f4)*f42 + 3*f1*f42*v1 - f43*v1 + f13*v4 - 3*f12*f4*v4)/f1mf43;
763  break;
764  }
765  case 103: // 4 freqs, no boundaries derivatives
766  {
767  double f12 = f1*f1;
768  double f13 = f12*f1;
769 
770  double f22 = f2*f2;
771  double f23 = f22*f2;
772 
773  double f32 = f3*f3;
774  double f33 = f32*f3;
775 
776  double f42 = f4*f4;
777  double f43 = f42*f4;
778 
779  double f1mf2 = f1-f2;
780  double f1mf3 = f1-f3;
781  double f1mf4 = f1-f4;
782  double f2mf3 = f2-f3;
783  double f2mf4 = f2-f4;
784  double f3mf4 = f3-f4;
785 
786  retVal = (f1*f1mf3*f1mf4*f3*f3mf4*f4*v2 + f23*(f1*f1mf4*f4*v3 + f32*(-(f4*v1) + f1*v4) + f3*(f42*v1 - f12*v4)) + f2*(f12*f1mf4*f42*v3 + f33*(-(f42*v1) + f12*v4) + f32*(f43*v1 - f13*v4)) +
787  f22*(f1*f4*(-f12 + f42)*v3 + f33*(f4*v1 - f1*v4) + f3*(-(f43*v1) + f13*v4)))/(f1mf2*f1mf3*f1mf4*f2mf3*f2mf4*f3mf4);
788  break;
789  }
790  case 1043: //no left derivative
791  {
792  double f12 = f1*f1;
793  double f13 = f12*f1;
794  double f14 = f13*f1;
795 
796  double f22 = f2*f2;
797  double f23 = f22*f2;
798  double f24 = f23*f2;
799 
800  double f32 = f3*f3;
801  double f33 = f32*f3;
802  double f34 = f33*f3;
803 
804  double f42 = f4*f4;
805  double f43 = f42*f4;
806  double f44 = f43*f4;
807  double f45 = f44*f4;
808 
809  double f1mf2 = f1-f2;
810  double f1mf3 = f1-f3;
811  double f1mf4 = f1-f4;
812  double f2mf3 = f2-f3;
813  double f2mf4 = f2-f4;
814  double f3mf4 = f3-f4;
815 
816  double f1mf42 = f1mf4*f1mf4;
817  double f3mf42 = f3mf4*f3mf4;
818  double f2mf42 = f2mf4*f2mf4;
819 
820  retVal = (-(d4*f1*f1mf2*f1mf3*f1mf4*f2*f2mf3*f2mf4*f3*f3mf4*f4) - f1*f1mf3*f1mf42*f3*f3mf42*f42*v2 +
821  f24*(-(f1*f1mf42*f42*v3) + f33*(f42*v1 + f12*v4 - 2*f1*f4*v4) + f3*f4*(f43*v1 + 2*f13*v4 - 3*f12*f4*v4) - f32*(2*f43*v1 + f13*v4 - 3*f1*f42*v4)) +
822  f2*f4*(f12*f1mf42*f43*v3 - f34*(f43*v1 + 2*f13*v4 - 3*f12*f4*v4) - f32*f4*(f44*v1 + 3*f14*v4 - 4*f13*f4*v4) + 2*f33*(f44*v1 + f14*v4 - 2*f12*f42*v4)) +
823  f22*(-(f1*f1mf42*(2*f1 + f4)*f43*v3) + f3*f42*(f44*v1 + 3*f14*v4 - 4*f13*f4*v4) + f34*(2*f43*v1 + f13*v4 - 3*f1*f42*v4) - f33*(3*f44*v1 + f14*v4 - 4*f1*f43*v4)) +
824  f23*(f1*f1mf42*(f1 + 2*f4)*f42*v3 - f34*(f42*v1 + f12*v4 - 2*f1*f4*v4) + f32*(3*f44*v1 + f14*v4 - 4*f1*f43*v4) - 2*f3*(f45*v1 + f14*f4*v4 - 2*f12*f43*v4)))/(f1mf2*f1mf3*f1mf42*f2mf3*f2mf42*f3mf42);
825  break;
826  }
827  case 1042: //4th order poly: v1,d1, v4,d4, v3 // used for the first intermediate region
828  {
829 
830  double f12 = f1*f1;
831  double f13 = f12*f1;
832  double f14 = f13*f1;
833  double f15 = f14*f1;
834 
835  double f42 = f4*f4;
836  double f43 = f42*f4;
837  double f44 = f43*f4;
838  double f45 = f44*f4;
839 
840  double f32 = f3*f3;
841  double f33 = f32*f3;
842  double f34 = f33*f3;
843 
844  double f1mf4 = f1-f4;
845  double f1mf3 = f1-f3;
846  double f3mf4 = f3-f4;
847 
848  double f1mf42 = f1mf4*f1mf4;
849  double f1mf32 = f1mf3*f1mf3;
850  double f3mf42 = f3mf4*f3mf4;
851 
852  double f1mf43 = f1mf42*f1mf4;
853 
854  retVal = (-(d4*f12*f1mf32*f1mf4*f3*f3mf4*f4) + d1*f1*f1mf3*f1mf4*f3*f3mf42*f42 - 4*f12*f33*f42*v1 + 3*f1*f34*f42*v1 + 8*f12*f32*f43*v1 - 4*f1*f33*f43*v1 - f34*f43*v1 - 4*f12*f3*f44*v1 - f1*f32*f44*v1 +
855  2*f33*f44*v1 + 2*f1*f3*f45*v1 - f32*f45*v1 + f15*f42*v3 - 3*f14*f43*v3 + 3*f13*f44*v3 - f12*f45*v3 + f15*f32*v4 - 2*f14*f33*v4 + f13*f34*v4 - 2*f15*f3*f4*v4 + f14*f32*f4*v4 + 4*f13*f33*f4*v4 -
856  3*f12*f34*f4*v4 + 4*f14*f3*f42*v4 - 8*f13*f32*f42*v4 + 4*f12*f33*f42*v4)/(f1mf32*f1mf43*f3mf42);
857 
858  break;
859  }
860  case 104: //Geraint's Version, 4th order poly: v1,d1, v4,d4, v2
861  {
862 
863  double f12 = f1*f1;
864 
865  double f42 = f4*f4;
866 
867  double f1mf2 = f1-f2;
868  double f1mf4 = f1-f4;
869  double f2mf4 = f2-f4;
870 
871  double f1mf22 = f1mf2*f1mf2;
872  double f2mf42 = f2mf4*f2mf4;
873  double f1mf43 = f1mf4*f1mf4*f1mf4;
874 
875  retVal = ((-(d4*f12*f1mf22*f1mf4*f2*f2mf4*f4) + d1*f1*f1mf2*f1mf4*f2*f2mf42*f42 + f42*(f2*f2mf42*(-4*f12 + 3*f1*f2 + 2*f1*f4 - f2*f4)*v1 + f12*f1mf43*v2) +
876  f12*f1mf22*f2*(f1*f2 - 2*f1*f4 - 3*f2*f4 + 4*f42)*v4)/(f1mf22*f1mf43*f2mf42));
877  break;
878  }
879  case 105: // Geraint, standard way: v1, v2, v3, v4, d1, d4
880  {
881  double f12 = f1*f1;
882  double f13 = f12*f1;
883  double f14 = f13*f1;
884  double f15 = f14*f1;
885  double f16 = f15*f1;
886  double f17 = f16*f1;
887 
888  double f22 = f2*f2;
889  double f23 = f22*f2;
890  double f24 = f23*f2;
891  double f25 = f24*f2;
892 
893  double f32 = f3*f3;
894  double f33 = f32*f3;
895  double f34 = f33*f3;
896  double f35 = f34*f3;
897 
898  double f42 = f4*f4;
899  double f43 = f42*f4;
900  double f44 = f43*f4;
901  double f45 = f44*f4;
902  double f46 = f45*f4;
903  double f47 = f46*f4;
904 
905  double f1mf2 = f1-f2;
906  double f1mf3 = f1-f3;
907  double f1mf4 = f1-f4;
908  double f2mf3 = f2-f3;
909  double f2mf4 = f2-f4;
910  double f3mf4 = f3-f4;
911 
912  double f1mf22 = f1mf2*f1mf2;
913  double f1mf32 = f1mf3*f1mf3;
914  double f1mf42 = f1mf4*f1mf4;
915  double f2mf42 = f2mf4*f2mf4;
916  double f3mf42 = f3mf4*f3mf4;
917  double f1mf43 = f1mf42*f1mf4;
918 
919  retVal = (
920  (-(d4*f12*f1mf22*f1mf32*f1mf4*f2*f2mf3*f2mf4*f3*f3mf4*f4) - d1*f1*f1mf2*f1mf3*f1mf4*f2*f2mf3*f2mf42*f3*f3mf42*f42 + 5*f13*f24*f33*f42*v1 - 4*f12*f25*f33*f42*v1 - 5*f13*f23*f34*f42*v1 +
921  3*f1*f25*f34*f42*v1 + 4*f12*f23*f35*f42*v1 - 3*f1*f24*f35*f42*v1 - 10*f13*f24*f32*f43*v1 + 8*f12*f25*f32*f43*v1 + 5*f12*f24*f33*f43*v1 - 4*f1*f25*f33*f43*v1 + 10*f13*f22*f34*f43*v1 -
922  5*f12*f23*f34*f43*v1 - f25*f34*f43*v1 - 8*f12*f22*f35*f43*v1 + 4*f1*f23*f35*f43*v1 + f24*f35*f43*v1 + 5*f13*f24*f3*f44*v1 - 4*f12*f25*f3*f44*v1 + 15*f13*f23*f32*f44*v1 -
923  10*f12*f24*f32*f44*v1 - f1*f25*f32*f44*v1 - 15*f13*f22*f33*f44*v1 + 5*f1*f24*f33*f44*v1 + 2*f25*f33*f44*v1 - 5*f13*f2*f34*f44*v1 + 10*f12*f22*f34*f44*v1 - 5*f1*f23*f34*f44*v1 +
924  4*f12*f2*f35*f44*v1 + f1*f22*f35*f44*v1 - 2*f23*f35*f44*v1 - 10*f13*f23*f3*f45*v1 + 5*f12*f24*f3*f45*v1 + 2*f1*f25*f3*f45*v1 - f12*f23*f32*f45*v1 + 2*f1*f24*f32*f45*v1 - f25*f32*f45*v1 +
925  10*f13*f2*f33*f45*v1 + f12*f22*f33*f45*v1 - 3*f24*f33*f45*v1 - 5*f12*f2*f34*f45*v1 - 2*f1*f22*f34*f45*v1 + 3*f23*f34*f45*v1 - 2*f1*f2*f35*f45*v1 + f22*f35*f45*v1 + 5*f13*f22*f3*f46*v1 +
926  2*f12*f23*f3*f46*v1 - 4*f1*f24*f3*f46*v1 - 5*f13*f2*f32*f46*v1 - f1*f23*f32*f46*v1 + 2*f24*f32*f46*v1 - 2*f12*f2*f33*f46*v1 + f1*f22*f33*f46*v1 + 4*f1*f2*f34*f46*v1 - 2*f22*f34*f46*v1 -
927  3*f12*f22*f3*f47*v1 + 2*f1*f23*f3*f47*v1 + 3*f12*f2*f32*f47*v1 - f23*f32*f47*v1 - 2*f1*f2*f33*f47*v1 + f22*f33*f47*v1 - f17*f33*f42*v2 + 2*f16*f34*f42*v2 - f15*f35*f42*v2 + 2*f17*f32*f43*v2 -
928  f16*f33*f43*v2 - 4*f15*f34*f43*v2 + 3*f14*f35*f43*v2 - f17*f3*f44*v2 - 4*f16*f32*f44*v2 + 8*f15*f33*f44*v2 - 3*f13*f35*f44*v2 + 3*f16*f3*f45*v2 - 8*f14*f33*f45*v2 + 4*f13*f34*f45*v2 +
929  f12*f35*f45*v2 - 3*f15*f3*f46*v2 + 4*f14*f32*f46*v2 + f13*f33*f46*v2 - 2*f12*f34*f46*v2 + f14*f3*f47*v2 - 2*f13*f32*f47*v2 + f12*f33*f47*v2 + f17*f23*f42*v3 - 2*f16*f24*f42*v3 +
930  f15*f25*f42*v3 - 2*f17*f22*f43*v3 + f16*f23*f43*v3 + 4*f15*f24*f43*v3 - 3*f14*f25*f43*v3 + f17*f2*f44*v3 + 4*f16*f22*f44*v3 - 8*f15*f23*f44*v3 + 3*f13*f25*f44*v3 - 3*f16*f2*f45*v3 +
931  8*f14*f23*f45*v3 - 4*f13*f24*f45*v3 - f12*f25*f45*v3 + 3*f15*f2*f46*v3 - 4*f14*f22*f46*v3 - f13*f23*f46*v3 + 2*f12*f24*f46*v3 - f14*f2*f47*v3 + 2*f13*f22*f47*v3 - f12*f23*f47*v3 +
932  f12*f1mf22*f1mf32*f2*f2mf3*f3*(f4*(-3*f2*f3 + 4*(f2 + f3)*f4 - 5*f42) + f1*(f2*f3 - 2*(f2 + f3)*f4 + 3*f42))*v4)/(f1mf22*f1mf32*f1mf43*f2mf3*f2mf42*f3mf42)
933  );
934  break;
935  }
936  default:
937  {
938  XLAL_ERROR_REAL8(XLAL_EINVAL, "Error in IMRPhenomXHM_Intermediate_Amp_delta0: IMRPhenomXIntermediateAmpVersion is not valid.\n");
939  }
940  }
941 
942  return retVal;
943 }
944 
945 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)
946 {
947  double retVal;
948  switch ( IntAmpFlag )
949  {
950  case 101: //linear, only v1, v2
951  {
952  double f1mf4 = f1-f4;
953 
954  retVal = (v1 - v4)/f1mf4;
955  break;
956  }
957  case 102: //quadratic: v1, v2, d2
958  {
959  double f12 = f1*f1;
960  double f42 = f4*f4;
961  double f1mf42 = (f1-f4)*(f1-f4);
962 
963  retVal = (d4*(f12 - f42) + 2*f4*(-v1 + v4))/f1mf42;
964  break;
965  }
966  case 1032: // 2 freqs, points and derivatives: v1, v4, d1, d4
967  {
968  double f12 = f1*f1;
969  double f42 = f4*f4;
970 
971  double f1mf4 = f1-f4;
972  double f1mf42 = f1mf4*f1mf4;
973  double f1mf43 = f1mf42*f1mf4;
974 
975  retVal = (d4*f1*f1mf4*(f1 + 2*f4) - f4*(d1*(-2*f12 + f1*f4 + f42) + 6*f1*(v1 - v4)))/f1mf43;
976  break;
977  }
978  case 103: // 4 freqs, no boundaries derivatives
979  {
980  double f12 = f1*f1;
981  double f13 = f12*f1;
982 
983  double f22 = f2*f2;
984  double f23 = f22*f2;
985 
986  double f32 = f3*f3;
987  double f33 = f32*f3;
988 
989  double f42 = f4*f4;
990  double f43 = f42*f4;
991 
992  double f1mf2 = f1-f2;
993  double f1mf3 = f1-f3;
994  double f1mf4 = f1-f4;
995  double f2mf3 = f2-f3;
996  double f2mf4 = f2-f4;
997  double f3mf4 = f3-f4;
998 
999  retVal = (f12*f1mf4*f42*(v2 - v3) + f33*(f42*(v1 - v2) + f12*(v2 - v4)) + f22*(f43*(v1 - v3) + f13*(v3 - v4) + f33*(-v1 + v4)) + f32*(f43*(-v1 + v2) + f13*(-v2 + v4)) +
1000  f23*(f42*(-v1 + v3) + f32*(v1 - v4) + f12*(-v3 + v4)))/(f1mf2*f1mf3*f1mf4*f2mf3*f2mf4*f3mf4);
1001  break;
1002  }
1003  case 1043: //no left derivative
1004  {
1005  double f12 = f1*f1;
1006  double f13 = f12*f1;
1007  double f14 = f13*f1;
1008 
1009  double f22 = f2*f2;
1010  double f23 = f22*f2;
1011  double f24 = f23*f2;
1012 
1013  double f32 = f3*f3;
1014  double f33 = f32*f3;
1015  double f34 = f33*f3;
1016 
1017  double f42 = f4*f4;
1018  double f43 = f42*f4;
1019  double f44 = f43*f4;
1020 
1021  double f1mf2 = f1-f2;
1022  double f1mf3 = f1-f3;
1023  double f1mf4 = f1-f4;
1024  double f2mf3 = f2-f3;
1025  double f2mf4 = f2-f4;
1026  double f3mf4 = f3-f4;
1027 
1028  double f1mf42 = f1mf4*f1mf4;
1029  double f2mf42 = f2mf4*f2mf4;
1030  double f3mf42 = f3mf4*f3mf4;
1031 
1032  double v1mv2 = v1-v2;
1033  double v2mv3 = v2-v3;
1034  double v2mv4 = v2-v4;
1035  double v1mv3 = v1-v3;
1036  double v1mv4 = v1-v4;
1037  double v3mv4 = v3-v4;
1038 
1039  retVal =(d4*f1mf4*f2mf4*f3mf4*(f1*f2*f3 + f2*f3*f4 + f1*(f2 + f3)*f4) + (f4*(f12*f1mf42*f43*v2mv3 + f34*(f43*v1mv2 +
1040  3*f12*f4*v2mv4 + 2*f13*(-v2 + v4)) + f32*f4*(f44*v1mv2 + 4*f13*f4*v2mv4 + 3*f14*(-v2 + v4)) + 2*f33*(f44*(-v1 + v2)
1041  + f14*v2mv4 + 2*f12*f42*(-v2 + v4)) + 2*f23*(f44*v1mv3 + f34*v1mv4 + 2*f12*f42*v3mv4 + 2*f32*f42*(-v1 + v4) + f14*(-v3 + v4))
1042  + f24*(3*f32*f4*v1mv4 + f43*(-v1 + v3) + 2*f13*v3mv4 + 2*f33*(-v1 + v4) + 3*f12*f4*(-v3 + v4)) + f22*f4*(4*f33*f4*v1mv4 + f44*(-v1 + v3)
1043  + 3*f14*v3mv4 + 3*f34*(-v1 + v4) + 4*f13*f4*(-v3 + v4))))/(f1mf2*f1mf3*f2mf3))/(f1mf42*f2mf42*f3mf42);
1044 
1045  break;
1046  }
1047  case 1042: //4th order poly: v1,d1, v4,d4, v3 // used for the first intermediate region
1048  {
1049  double f12 = f1*f1;
1050  double f13 = f12*f1;
1051  double f14 = f13*f1;
1052  double f15 = f14*f1;
1053 
1054  double f42 = f4*f4;
1055  double f43 = f42*f4;
1056  double f44 = f43*f4;
1057  double f45 = f44*f4;
1058 
1059  double f32 = f3*f3;
1060  double f33 = f32*f3;
1061  double f34 = f33*f3;
1062 
1063  double f1mf4 = f1-f4;
1064  double f1mf3 = f1-f3;
1065  double f3mf4 = f3-f4;
1066 
1067  double f1mf42 = f1mf4*f1mf4;
1068  double f1mf32 = f1mf3*f1mf3;
1069  double f3mf42 = f3mf4*f3mf4;
1070 
1071  double f1mf43 = f1mf42*f1mf4;
1072 
1073  retVal = (d4*f15*f32 - 2*d4*f14*f33 + d4*f13*f34 + d4*f14*f32*f4 - 2*d1*f13*f33*f4 - 2*d4*f13*f33*f4 + 2*d1*f12*f34*f4 + d4*f12*f34*f4 - d4*f15*f42 + 3*d1*f13*f32*f42 + d4*f13*f32*f42 - 2*d1*f12*f33*f42 +
1074  2*d4*f12*f33*f42 - d1*f1*f34*f42 - 2*d4*f1*f34*f42 + d4*f14*f43 - d1*f12*f32*f43 - 3*d4*f12*f32*f43 + 2*d1*f1*f33*f43 + 2*d4*f1*f33*f43 - d1*f34*f43 - d1*f13*f44 - d1*f1*f32*f44 + 2*d1*f33*f44 +
1075  d1*f12*f45 - d1*f32*f45 + 8*f12*f33*f4*v1 - 6*f1*f34*f4*v1 - 12*f12*f32*f42*v1 + 8*f1*f33*f42*v1 + 4*f12*f44*v1 - 2*f1*f45*v1 - 2*f15*f4*v3 + 4*f14*f42*v3 - 4*f12*f44*v3 + 2*f1*f45*v3 +
1076  2*f15*f4*v4 - 8*f12*f33*f4*v4 + 6*f1*f34*f4*v4 - 4*f14*f42*v4 + 12*f12*f32*f42*v4 - 8*f1*f33*f42*v4)/(f1mf32*f1mf43*f3mf42);
1077 
1078  #if DEBUG == 1
1079  printf("\ndelta1 = %.16f", retVal);
1080  printf("\nf1 = %.16f", f1);
1081  printf("\nf2 = %.16f", f2);
1082  printf("\nf3 = %.16f", f3);
1083  printf("\nf4 = %.16f", f4);
1084  printf("\nv1 = %.16f", v1);
1085  printf("\nv2 = %.16f", v2);
1086  printf("\nv3 = %.16f", v3);
1087  printf("\nv4 = %.16f", v4);
1088  printf("\nd1 = %.16f", d1);
1089  printf("\nd4 = %.16f", d4);
1090  #endif
1091 
1092  break;
1093  }
1094  case 104: //Geraint's Version, 4th order poly: v1,d1, v2,d2, v3
1095  {
1096 
1097  double f12 = f1*f1;
1098  double f13 = f12*f1;
1099  double f14 = f13*f1;
1100 
1101  double f22 = f2*f2;
1102  double f23 = f22*f2;
1103  double f24 = f23*f2;
1104 
1105  double f42 = f4*f4;
1106  double f43 = f42*f4;
1107  double f44 = f43*f4;
1108 
1109  double f1mf2 = f1-f2;
1110  double f1mf4 = f1-f4;
1111  double f2mf4 = f2-f4;
1112 
1113  double f1mf22 = f1mf2*f1mf2;
1114  double f2mf42 = f2mf4*f2mf4;
1115  double f1mf43 = f1mf4*f1mf4*f1mf4;
1116 
1117  retVal = ((d4*f1*f1mf22*f1mf4*f2mf4*(2*f2*f4 + f1*(f2 + f4)) + f4*(-(d1*f1mf2*f1mf4*f2mf42*(2*f1*f2 + (f1 + f2)*f4)) -
1118  2*f1*(f44*(v1 - v2) + 3*f24*(v1 - v4) + f14*(v2 - v4) + 4*f23*f4*(-v1 + v4)
1119  + 2*f13*f4*(-v2 + v4) + f1*(2*f43*(-v1 + v2) + 6*f22*f4*(v1 - v4) + 4*f23*(-v1 + v4)))))/(f1mf22*f1mf43*f2mf42));
1120  break;
1121  }
1122  case 105: // Geraint, standard way: v1, v2, v3, v4, d1, d4
1123  {
1124 
1125  double f12 = f1*f1;
1126  double f13 = f12*f1;
1127  double f14 = f13*f1;
1128  double f15 = f14*f1;
1129  double f16 = f15*f1;
1130 
1131  double f22 = f2*f2;
1132  double f23 = f22*f2;
1133  double f24 = f23*f2;
1134  double f25 = f24*f2;
1135 
1136  double f32 = f3*f3;
1137  double f33 = f32*f3;
1138  double f34 = f33*f3;
1139  double f35 = f34*f3;
1140 
1141  double f42 = f4*f4;
1142  double f43 = f42*f4;
1143  double f44 = f43*f4;
1144  double f45 = f44*f4;
1145 
1146  double f1mf2 = f1-f2;
1147  double f1mf3 = f1-f3;
1148  double f1mf4 = f1-f4;
1149  double f2mf3 = f2-f3;
1150  double f2mf4 = f2-f4;
1151  double f3mf4 = f3-f4;
1152 
1153  double f1mf22 = f1mf2*f1mf2;
1154  double f1mf32 = f1mf3*f1mf3;
1155  double f2mf42 = f2mf4*f2mf4;
1156  double f3mf42 = f3mf4*f3mf4;
1157  double f1mf43 = f1mf4*f1mf4*f1mf4;
1158 
1159  retVal = (
1160  (d4*f1*f1mf22*f1mf32*f1mf4*f2mf3*f2mf4*f3mf4*(f1*f2*f3 + 2*f2*f3*f4 + f1*(f2 + f3)*f4) +
1161  f4*(d1*f1mf2*f1mf3*f1mf4*f2mf3*f2mf42*f3mf42*(2*f1*f2*f3 + f2*f3*f4 + f1*(f2 + f3)*f4) +
1162  f1*(f16*(f43*(v2 - v3) + 2*f33*(v2 - v4) + 3*f22*f4*(v3 - v4) + 3*f32*f4*(-v2 + v4) + 2*f23*(-v3 + v4)) +
1163  f13*f4*(f45*(-v2 + v3) + 5*f34*f4*(v2 - v4) + 4*f25*(v3 - v4) + 4*f35*(-v2 + v4) + 5*f24*f4*(-v3 + v4)) +
1164  f14*(3*f45*(v2 - v3) + 2*f35*(v2 - v4) + 5*f34*f4*(v2 - v4) + 10*f23*f42*(v3 - v4) + 10*f33*f42*(-v2 + v4) + 2*f25*(-v3 + v4) + 5*f24*f4*(-v3 + v4)) +
1165  f15*(3*f44*(-v2 + v3) + 2*f33*f4*(v2 - v4) + 5*f32*f42*(v2 - v4) + 4*f24*(v3 - v4) + 4*f34*(-v2 + v4) + 2*f23*f4*(-v3 + v4) + 5*f22*f42*(-v3 + v4)) -
1166  5*f12*(-(f32*f3mf42*f43*(v1 - v2)) + 2*f23*(f44*(-v1 + v3) + 2*f32*f42*(v1 - v4) + f34*(-v1 + v4)) + f24*(f43*(v1 - v3) + 2*f33*(v1 - v4) + 3*f32*f4*(-v1 + v4)) +
1167  f22*f4*(f44*(v1 - v3) + 3*f34*(v1 - v4) + 4*f33*f4*(-v1 + v4))) +
1168  f1*(-(f32*f3mf42*(4*f3 + 3*f4)*f43*(v1 - v2)) + 2*f23*(f45*(-v1 + v3) + 5*f34*f4*(v1 - v4) + 4*f35*(-v1 + v4)) + 4*f25*(f43*(v1 - v3) + 2*f33*(v1 - v4) + 3*f32*f4*(-v1 + v4)) -
1169  5*f24*f4*(f43*(v1 - v3) + 2*f33*(v1 - v4) + 3*f32*f4*(-v1 + v4)) + 3*f22*f4*(f45*(v1 - v3) + 4*f35*(v1 - v4) + 5*f34*f4*(-v1 + v4))) -
1170  2*(-(f33*f3mf42*f44*(v1 - v2)) + f24*(2*f45*(-v1 + v3) + 5*f33*f42*(v1 - v4) + 3*f35*(-v1 + v4)) + f25*(f44*(v1 - v3) + 3*f34*(v1 - v4) + 4*f33*f4*(-v1 + v4)) +
1171  f23*f4*(f45*(v1 - v3) + 4*f35*(v1 - v4) + 5*f34*f4*(-v1 + v4))))))/(f1mf22*f1mf32*f1mf43*f2mf3*f2mf42*f3mf42)
1172  );
1173  break;
1174  }
1175  default:
1176  {
1177  XLAL_ERROR_REAL8(XLAL_EINVAL, "Error in IMRPhenomXHM_Intermediate_Amp_delta1: IMRPhenomXIntermediateAmpVersion is not valid.\n");
1178  }
1179  }
1180  return retVal;
1181  }
1182 
1183 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)
1184 {
1185  double retVal;
1186 
1187  switch ( IntAmpFlag )
1188  {
1189  case 101: //linear, only v1, v2
1190  {
1191  retVal = 0.;
1192  break;
1193  }
1194  case 102: //quadratic: v1, v2, d2
1195  {
1196  double f1mf4 = f1-f4;
1197  double f1mf42 = f1mf4*f1mf4;
1198 
1199  retVal = (-(d4*f1mf4) + v1 - v4)/f1mf42;
1200  break;
1201  }
1202  case 1032: // 2 freqs, points and derivatives: v1, v4, d1, d4
1203  {
1204  double f12 = f1*f1;
1205  double f42 = f4*f4;
1206 
1207  double f1mf4 = f1-f4;
1208  double f1mf42 = f1mf4*f1mf4;
1209  double f1mf43 = f1mf42*f1mf4;
1210 
1211  retVal = (-(d1*(f12 + f1*f4 - 2*f42)) + d4*(-2*f12 + f1*f4 + f42) + 3*(f1 + f4)*(v1 - v4))/f1mf43;
1212  break;
1213  }
1214  case 103: // 4 freqs, no boundaries derivatives
1215  {
1216  double f12 = f1*f1;
1217  double f13 = f12*f1;
1218 
1219  double f22 = f2*f2;
1220  double f23 = f22*f2;
1221 
1222  double f32 = f3*f3;
1223  double f33 = f32*f3;
1224 
1225  double f42 = f4*f4;
1226  double f43 = f42*f4;
1227 
1228  double f1mf2 = f1-f2;
1229  double f1mf3 = f1-f3;
1230  double f1mf4 = f1-f4;
1231  double f2mf3 = f2-f3;
1232  double f2mf4 = f2-f4;
1233  double f3mf4 = f3-f4;
1234 
1235  retVal = (-(f1*f4*(f12 - f42)*(v2 - v3)) + f3*(f43*(v1 - v2) + f13*(v2 - v4)) + f23*(f4*(v1 - v3) + f1*(v3 - v4) + f3*(-v1 + v4)) + f33*(f4*(-v1 + v2) + f1*(-v2 + v4)) +
1236  f2*(f43*(-v1 + v3) + f33*(v1 - v4) + f13*(-v3 + v4)))/(f1mf2*f1mf3*f1mf4*f2mf3*f2mf4*f3mf4);
1237  break;
1238  }
1239  case 1043: //no left derivative: v1, v2, v3, v4, d4
1240  {
1241  double f12 = f1*f1;
1242  double f13 = f12*f1;
1243  double f14 = f13*f1;
1244 
1245  double f22 = f2*f2;
1246  double f23 = f22*f2;
1247  double f24 = f23*f2;
1248 
1249  double f32 = f3*f3;
1250  double f33 = f32*f3;
1251  double f34 = f33*f3;
1252 
1253  double f42 = f4*f4;
1254  double f43 = f42*f4;
1255  double f44 = f43*f4;
1256  double f46 = f44*f42;
1257 
1258  double f1mf2 = f1-f2;
1259  double f1mf3 = f1-f3;
1260  double f1mf4 = f1-f4;
1261  double f2mf3 = f2-f3;
1262  double f2mf4 = f2-f4;
1263  double f3mf4 = f3-f4;
1264 
1265  double f1mf42 = f1mf4*f1mf4;
1266  double f2mf42 = f2mf4*f2mf4;
1267  double f3mf42 = f3mf4*f3mf4;
1268 
1269  double v1mv3 = v1-v3;
1270  double v1mv4 = v1-v4;
1271  double v3mv4 = v3-v4;
1272 
1273  retVal = (-(d4*f1mf2*f1mf3*f1mf4*f2mf3*f2mf4*f3mf4*(f3*f4 + f2*(f3 + f4) + f1*(f2 + f3 + f4))) - 2*f34*f43*v1 + 3*f33*f44*v1 - f3*f46*v1 - f14*f33*v2 + f13*f34*v2 + 3*f14*f3*f42*v2 - 3*f1*f34*f42*v2 -
1274  2*f14*f43*v2 - 4*f13*f3*f43*v2 + 4*f1*f33*f43*v2 + 2*f34*f43*v2 + 3*f13*f44*v2 - 3*f33*f44*v2 - f1*f46*v2 + f3*f46*v2 + 2*f14*f43*v3 - 3*f13*f44*v3 + f1*f46*v3 +
1275  f2*f42*(f44*v1mv3 + 3*f34*v1mv4 - 4*f33*f4*v1mv4 - 3*f14*v3mv4 + 4*f13*f4*v3mv4) + f24*(2*f43*v1mv3 + f33*v1mv4 - 3*f3*f42*v1mv4 - f13*v3mv4 + 3*f1*f42*v3mv4) +
1276  f23*(-3*f44*v1mv3 - f34*v1mv4 + 4*f3*f43*v1mv4 + f14*v3mv4 - 4*f1*f43*v3mv4) + f14*f33*v4 - f13*f34*v4 - 3*f14*f3*f42*v4 + 3*f1*f34*f42*v4 + 4*f13*f3*f43*v4 - 4*f1*f33*f43*v4)/
1277  (f1mf2*f1mf3*f1mf42*f2mf3*f2mf42*f3mf42);
1278  break;
1279  }
1280  case 1042: //4th order poly: v1,d1, v2,d2, v3 // used for the first intermediate region
1281  {
1282  double f12 = f1*f1;
1283  double f13 = f12*f1;
1284  double f14 = f13*f1;
1285  double f15 = f14*f1;
1286 
1287  double f42 = f4*f4;
1288  double f43 = f42*f4;
1289  double f44 = f43*f4;
1290  double f45 = f44*f4;
1291 
1292  double f32 = f3*f3;
1293  double f33 = f32*f3;
1294  double f34 = f33*f3;
1295 
1296  double f1mf4 = f1-f4;
1297  double f1mf3 = f1-f3;
1298  double f3mf4 = f3-f4;
1299 
1300  double f1mf42 = f1mf4*f1mf4;
1301  double f1mf32 = f1mf3*f1mf3;
1302  double f3mf42 = f3mf4*f3mf4;
1303 
1304  double f1mf43 = f1mf42*f1mf4;
1305 
1306  retVal = (-(d4*f1mf32*f1mf4*f3mf4*(f12 + f3*f4 + 2*f1*(f3 + f4))) + d1*f1mf3*f1mf4*f3mf42*(f1*f3 + 2*(f1 + f3)*f4 + f42) - 4*f12*f33*v1 + 3*f1*f34*v1 - 4*f1*f33*f4*v1 + 3*f34*f4*v1 + 12*f12*f3*f42*v1 -
1307  4*f33*f42*v1 - 8*f12*f43*v1 + f1*f44*v1 + f45*v1 + f15*v3 + f14*f4*v3 - 8*f13*f42*v3 + 8*f12*f43*v3 - f1*f44*v3 - f45*v3 -
1308  f1mf32*(f13 + f3*(3*f3 - 4*f4)*f4 + f12*(2*f3 + f4) + f1*(3*f3 - 4*f4)*(f3 + 2*f4))*v4)/(f1mf32*f1mf43*f3mf42);
1309  break;
1310  }
1311  case 104: //Geraint's Version, 4th order poly: v1,d1, v2,d2, v3
1312  {
1313  double f12 = f1*f1;
1314  double f13 = f12*f1;
1315  double f14 = f13*f1;
1316  double f15 = f14*f1;
1317 
1318  double f22 = f2*f2;
1319  double f23 = f22*f2;
1320  double f24 = f23*f2;
1321 
1322  double f42 = f4*f4;
1323  double f43 = f42*f4;
1324  double f44 = f43*f4;
1325  double f45 = f44*f4;
1326 
1327  double f1mf2 = f1-f2;
1328  double f1mf4 = f1-f4;
1329  double f2mf4 = f2-f4;
1330 
1331  double f1mf22 = f1mf2*f1mf2;
1332  double f2mf42 = f2mf4*f2mf4;
1333  double f1mf43 = f1mf4*f1mf4*f1mf4;
1334 
1335  retVal = ((-(d4*f1mf22*f1mf4*f2mf4*(f12 + f2*f4 + 2*f1*(f2 + f4))) + d1*f1mf2*f1mf4*f2mf42*(f1*f2 + 2*(f1 + f2)*f4 + f42)
1336  - 4*f12*f23*v1 + 3*f1*f24*v1 - 4*f1*f23*f4*v1 + 3*f24*f4*v1 + 12*f12*f2*f42*v1 -
1337  4*f23*f42*v1 - 8*f12*f43*v1 + f1*f44*v1 + f45*v1 + f15*v2 + f14*f4*v2 - 8*f13*f42*v2 + 8*f12*f43*v2 - f1*f44*v2 - f45*v2 -
1338  f1mf22*(f13 + f2*(3*f2 - 4*f4)*f4 + f12*(2*f2 + f4) + f1*(3*f2 - 4*f4)*(f2 + 2*f4))*v4)/(f1mf22*f1mf43*f2mf42));
1339 
1340  break;
1341  }
1342  case 105: // Geraint, standard way: v1, v2, v3, v4, d1, d4
1343  {
1344  double f12 = f1*f1;
1345  double f13 = f12*f1;
1346  double f14 = f13*f1;
1347  double f15 = f14*f1;
1348  double f16 = f15*f1;
1349  double f17 = f16*f1;
1350 
1351  double f22 = f2*f2;
1352  double f23 = f22*f2;
1353  double f24 = f23*f2;
1354  double f25 = f24*f2;
1355 
1356  double f32 = f3*f3;
1357  double f33 = f32*f3;
1358  double f34 = f33*f3;
1359  double f35 = f34*f3;
1360 
1361  double f42 = f4*f4;
1362  double f43 = f42*f4;
1363  double f44 = f43*f4;
1364  double f45 = f44*f4;
1365  double f46 = f45*f4;
1366  double f47 = f46*f4;
1367 
1368  double f1mf2 = f1-f2;
1369  double f1mf3 = f1-f3;
1370  double f1mf4 = f1-f4;
1371  double f2mf3 = f2-f3;
1372  double f2mf4 = f2-f4;
1373  double f3mf4 = f3-f4;
1374 
1375  double f1mf22 = f1mf2*f1mf2;
1376  double f1mf32 = f1mf3*f1mf3;
1377  double f2mf42 = f2mf4*f2mf4;
1378  double f3mf42 = f3mf4*f3mf4;
1379  double f1mf43 = f1mf4*f1mf4*f1mf4;
1380 
1381  retVal = (
1382  (-(d4*f1mf22*f1mf32*f1mf4*f2mf3*f2mf4*f3mf4*(f2*f3*f4 + f12*(f2 + f3 + f4) + 2*f1*(f2*f3 + (f2 + f3)*f4))) -
1383  d1*f1mf2*f1mf3*f1mf4*f2mf3*f2mf42*f3mf42*(f1*f2*f3 + 2*(f2*f3 + f1*(f2 + f3))*f4 + (f1 + f2 + f3)*f42) + 5*f13*f24*f33*v1 - 4*f12*f25*f33*v1 - 5*f13*f23*f34*v1 + 3*f1*f25*f34*v1 +
1384  4*f12*f23*f35*v1 - 3*f1*f24*f35*v1 + 5*f12*f24*f33*f4*v1 - 4*f1*f25*f33*f4*v1 - 5*f12*f23*f34*f4*v1 + 3*f25*f34*f4*v1 + 4*f1*f23*f35*f4*v1 - 3*f24*f35*f4*v1 - 15*f13*f24*f3*f42*v1 +
1385  12*f12*f25*f3*f42*v1 + 5*f1*f24*f33*f42*v1 - 4*f25*f33*f42*v1 + 15*f13*f2*f34*f42*v1 - 5*f1*f23*f34*f42*v1 - 12*f12*f2*f35*f42*v1 + 4*f23*f35*f42*v1 + 10*f13*f24*f43*v1 - 8*f12*f25*f43*v1 +
1386  20*f13*f23*f3*f43*v1 - 15*f12*f24*f3*f43*v1 - 20*f13*f2*f33*f43*v1 + 5*f24*f33*f43*v1 - 10*f13*f34*f43*v1 + 15*f12*f2*f34*f43*v1 - 5*f23*f34*f43*v1 + 8*f12*f35*f43*v1 - 15*f13*f23*f44*v1 +
1387  10*f12*f24*f44*v1 + f1*f25*f44*v1 + 15*f13*f33*f44*v1 - 10*f12*f34*f44*v1 - f1*f35*f44*v1 + f12*f23*f45*v1 - 2*f1*f24*f45*v1 + f25*f45*v1 - f12*f33*f45*v1 + 2*f1*f34*f45*v1 - f35*f45*v1 +
1388  5*f13*f2*f46*v1 + f1*f23*f46*v1 - 2*f24*f46*v1 - 5*f13*f3*f46*v1 - f1*f33*f46*v1 + 2*f34*f46*v1 - 3*f12*f2*f47*v1 + f23*f47*v1 + 3*f12*f3*f47*v1 - f33*f47*v1 - f17*f33*v2 + 2*f16*f34*v2 -
1389  f15*f35*v2 - f16*f33*f4*v2 + 2*f15*f34*f4*v2 - f14*f35*f4*v2 + 3*f17*f3*f42*v2 - f15*f33*f42*v2 - 10*f14*f34*f42*v2 + 8*f13*f35*f42*v2 - 2*f17*f43*v2 - 5*f16*f3*f43*v2 + 15*f14*f33*f43*v2 -
1390  8*f12*f35*f43*v2 + 4*f16*f44*v2 - 15*f13*f33*f44*v2 + 10*f12*f34*f44*v2 + f1*f35*f44*v2 + f12*f33*f45*v2 - 2*f1*f34*f45*v2 + f35*f45*v2 - 4*f14*f46*v2 + 5*f13*f3*f46*v2 + f1*f33*f46*v2 -
1391  2*f34*f46*v2 + 2*f13*f47*v2 - 3*f12*f3*f47*v2 + f33*f47*v2 + f17*f23*v3 - 2*f16*f24*v3 + f15*f25*v3 + f16*f23*f4*v3 - 2*f15*f24*f4*v3 + f14*f25*f4*v3 - 3*f17*f2*f42*v3 + f15*f23*f42*v3 +
1392  10*f14*f24*f42*v3 - 8*f13*f25*f42*v3 + 2*f17*f43*v3 + 5*f16*f2*f43*v3 - 15*f14*f23*f43*v3 + 8*f12*f25*f43*v3 - 4*f16*f44*v3 + 15*f13*f23*f44*v3 - 10*f12*f24*f44*v3 - f1*f25*f44*v3 -
1393  f12*f23*f45*v3 + 2*f1*f24*f45*v3 - f25*f45*v3 + 4*f14*f46*v3 - 5*f13*f2*f46*v3 - f1*f23*f46*v3 + 2*f24*f46*v3 - 2*f13*f47*v3 + 3*f12*f2*f47*v3 - f23*f47*v3 -
1394  f1mf22*f1mf32*f2mf3*(f13*(f22 + f2*f3 + f32 - 3*f42) + f2*f3*f4*(3*f2*f3 - 4*(f2 + f3)*f4 + 5*f42) + f1*(f2*f3 + 2*(f2 + f3)*f4)*(3*f2*f3 - 4*(f2 + f3)*f4 + 5*f42) +
1395  f12*(2*f2*f3*(f2 + f3) + (f22 + f2*f3 + f32)*f4 - 6*(f2 + f3)*f42 + 5*f43))*v4)/(f1mf22*f1mf32*f1mf43*f2mf3*f2mf42*f3mf42)
1396  );
1397 
1398  break;
1399  }
1400  default:
1401  {
1402  XLAL_ERROR_REAL8(XLAL_EINVAL, "Error in IMRPhenomXHM_Intermediate_Amp_delta2: IMRPhenomXIntermediateAmpVersion is not valid.\n");
1403  }
1404  }
1405 
1406  return retVal;
1407 }
1408 
1409 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)
1410 {
1411  double retVal;
1412 
1413  switch ( IntAmpFlag )
1414  {
1415  case 101: //linear, only v1, v2
1416  {
1417  retVal = 0.;
1418  break;
1419  }
1420  case 102: //quadratic: v1, v2, d2
1421  {
1422  retVal = 0.;
1423  break;
1424  }
1425  case 1032: // 2 freqs, points and derivatives: v1, v4, d1, d4
1426  {
1427  double f1mf4 = f1-f4;
1428  double f1mf42 = f1mf4*f1mf4;
1429  double f1mf43 = f1mf42*f1mf4;
1430 
1431  retVal = (d1*f1mf4 + d4*f1mf4 - 2*v1 + 2*v4)/f1mf43;
1432  break;
1433  }
1434  case 103: // 4 freqs, no boundaries derivatives
1435  {
1436  double f12 = f1*f1;
1437 
1438  double f22 = f2*f2;
1439 
1440  double f32 = f3*f3;
1441 
1442  double f42 = f4*f4;
1443 
1444  double f1mf2 = f1-f2;
1445  double f1mf3 = f1-f3;
1446  double f1mf4 = f1-f4;
1447  double f2mf3 = f2-f3;
1448  double f2mf4 = f2-f4;
1449  double f3mf4 = f3-f4;
1450 
1451  retVal = (f1*f1mf4*f4*(v2 - v3) + f32*(f4*(v1 - v2) + f1*(v2 - v4)) + f2*(f42*(v1 - v3) + f12*(v3 - v4) + f32*(-v1 + v4)) + f3*(f42*(-v1 + v2) + f12*(-v2 + v4)) +
1452  f22*(f4*(-v1 + v3) + f3*(v1 - v4) + f1*(-v3 + v4)))/(f1mf2*f1mf3*f1mf4*f2mf3*f2mf4*f3mf4);
1453  break;
1454  }
1455  case 1043: //no left derivative: v1, v2, v3, v4, d4
1456  {
1457  double f12 = f1*f1;
1458  double f13 = f12*f1;
1459  double f14 = f13*f1;
1460 
1461  double f22 = f2*f2;
1462  double f23 = f22*f2;
1463  double f24 = f23*f2;
1464 
1465  double f32 = f3*f3;
1466  double f34 = f32*f32;
1467 
1468  double f42 = f4*f4;
1469  double f43 = f42*f4;
1470  double f44 = f43*f4;
1471  double f45 = f44*f4;
1472 
1473  double f1mf2 = f1-f2;
1474  double f1mf3 = f1-f3;
1475  double f1mf4 = f1-f4;
1476  double f2mf3 = f2-f3;
1477  double f2mf4 = f2-f4;
1478  double f3mf4 = f3-f4;
1479 
1480  double f1mf42 = f1mf4*f1mf4;
1481  double f2mf42 = f2mf4*f2mf4;
1482  double f3mf42 = f3mf4*f3mf4;
1483 
1484  double v1mv3 = v1-v3;
1485  double v1mv4 = v1-v4;
1486  double v3mv4 = v3-v4;
1487 
1488  retVal = (d4*f1mf2*f1mf3*f1mf4*f2mf3*f2mf4*f3mf4*(f1 + f2 + f3 + f4) + f34*f42*v1 - 3*f32*f44*v1 + 2*f3*f45*v1 + f14*f32*v2 - f12*f34*v2 - 2*f14*f3*f4*v2 + 2*f1*f34*f4*v2 + f14*f42*v2 - f34*f42*v2 +
1489  4*f12*f3*f43*v2 - 4*f1*f32*f43*v2 - 3*f12*f44*v2 + 3*f32*f44*v2 + 2*f1*f45*v2 - 2*f3*f45*v2 - f14*f42*v3 + 3*f12*f44*v3 - 2*f1*f45*v3 +
1490  f24*(-(f42*v1mv3) - f32*v1mv4 + 2*f3*f4*v1mv4 + f12*v3mv4 - 2*f1*f4*v3mv4) - 2*f2*f4*(f44*v1mv3 + f34*v1mv4 - 2*f32*f42*v1mv4 - f14*v3mv4 + 2*f12*f42*v3mv4) +
1491  f22*(3*f44*v1mv3 + f34*v1mv4 - 4*f3*f43*v1mv4 - f14*v3mv4 + 4*f1*f43*v3mv4) - f14*f32*v4 + f12*f34*v4 + 2*f14*f3*f4*v4 - 2*f1*f34*f4*v4 - 4*f12*f3*f43*v4 + 4*f1*f32*f43*v4)/
1492  (f1mf2*f1mf3*f1mf42*f2mf3*f2mf42*f3mf42);
1493  break;
1494  }
1495  case 1042: //4th order poly: v1,d1, v2,d2, v3 // used for the first intermediate region
1496  {
1497 
1498  double f12 = f1*f1;
1499  double f13 = f12*f1;
1500  double f14 = f13*f1;
1501 
1502  double f42 = f4*f4;
1503  double f43 = f42*f4;
1504  double f44 = f43*f4;
1505 
1506  double f32 = f3*f3;
1507  double f33 = f32*f3;
1508  double f34 = f33*f3;
1509 
1510  double f1mf4 = f1-f4;
1511  double f1mf3 = f1-f3;
1512  double f3mf4 = f3-f4;
1513 
1514  double f1mf42 = f1mf4*f1mf4;
1515  double f1mf32 = f1mf3*f1mf3;
1516  double f3mf42 = f3mf4*f3mf4;
1517 
1518  double f1mf43 = f1mf42*f1mf4;
1519 
1520  retVal = (2*d4*f14*f3 - d1*f13*f32 - 3*d4*f13*f32 + d1*f1*f34 + d4*f1*f34 - 2*d4*f14*f4 + 2*d1*f13*f3*f4 + 2*d4*f13*f3*f4 - d1*f12*f32*f4 + d4*f12*f32*f4 - d1*f34*f4 - d4*f34*f4 - d1*f13*f42 + d4*f13*f42 +
1521  2*d1*f12*f3*f42 - 2*d4*f12*f3*f42 - d1*f1*f32*f42 + d4*f1*f32*f42 - d1*f12*f43 + d4*f12*f43 - 2*d1*f1*f3*f43 - 2*d4*f1*f3*f43 + 3*d1*f32*f43 + d4*f32*f43 + 2*d1*f1*f44 - 2*d1*f3*f44 +
1522  4*f12*f32*v1 - 2*f34*v1 - 8*f12*f3*f4*v1 + 4*f1*f32*f4*v1 + 4*f12*f42*v1 - 8*f1*f3*f42*v1 + 4*f32*f42*v1 + 4*f1*f43*v1 - 2*f44*v1 - 2*f14*v3 + 4*f13*f4*v3 - 4*f1*f43*v3 + 2*f44*v3 + 2*f14*v4 -
1523  4*f12*f32*v4 + 2*f34*v4 - 4*f13*f4*v4 + 8*f12*f3*f4*v4 - 4*f1*f32*f4*v4 - 4*f12*f42*v4 + 8*f1*f3*f42*v4 - 4*f32*f42*v4)/(f1mf32*f1mf43*f3mf42);
1524  break;
1525  }
1526  case 104: //Geraint's Version, 4th order poly: v1,d1, v2,d2, v3
1527  {
1528 
1529  double f12 = f1*f1;
1530  double f13 = f12*f1;
1531  double f14 = f13*f1;
1532 
1533  double f22 = f2*f2;
1534  double f23 = f22*f2;
1535  double f24 = f23*f2;
1536 
1537  double f42 = f4*f4;
1538  double f43 = f42*f4;
1539  double f44 = f43*f4;
1540 
1541  double f1mf2 = f1-f2;
1542  double f1mf4 = f1-f4;
1543  double f2mf4 = f2-f4;
1544 
1545  double f1mf22 = f1mf2*f1mf2;
1546  double f2mf42 = f2mf4*f2mf4;
1547  double f1mf43 = f1mf4*f1mf4*f1mf4;
1548 
1549  retVal = ((d4*f1mf22*f1mf4*f2mf4*(2*f1 + f2 + f4) - d1*f1mf2*f1mf4*f2mf42*(f1 + f2 + 2*f4)
1550  + 2*(f44*(-v1 + v2) + 2*f12*f2mf42*(v1 - v4) + 2*f22*f42*(v1 - v4)
1551  + 2*f13*f4*(v2 - v4) + f24*(-v1 + v4) + f14*(-v2 + v4) + 2*f1*f4*(f42*(v1 - v2) + f22*(v1 - v4) + 2*f2*f4*(-v1 + v4)))) / (f1mf22*f1mf43*f2mf42));
1552 
1553 
1554  break;
1555  }
1556  case 105: // Geraint, standard way: v1, v2, v3, v4, d1, d4
1557  {
1558  double f12 = f1*f1;
1559  double f13 = f12*f1;
1560  double f14 = f13*f1;
1561  double f15 = f14*f1;
1562  double f16 = f15*f1;
1563  double f17 = f16*f1;
1564 
1565  double f22 = f2*f2;
1566  double f23 = f22*f2;
1567  double f24 = f23*f2;
1568  double f25 = f24*f2;
1569 
1570  double f32 = f3*f3;
1571  double f33 = f32*f3;
1572  double f34 = f33*f3;
1573  double f35 = f34*f3;
1574 
1575  double f42 = f4*f4;
1576  double f43 = f42*f4;
1577  double f44 = f43*f4;
1578  double f45 = f44*f4;
1579  double f46 = f45*f4;
1580  double f47 = f46*f4;
1581 
1582  double f1mf2 = f1-f2;
1583  double f1mf3 = f1-f3;
1584  double f1mf4 = f1-f4;
1585  double f2mf3 = f2-f3;
1586  double f2mf4 = f2-f4;
1587  double f3mf4 = f3-f4;
1588 
1589  double f1mf22 = f1mf2*f1mf2;
1590  double f1mf32 = f1mf3*f1mf3;
1591  double f2mf42 = f2mf4*f2mf4;
1592  double f3mf42 = f3mf4*f3mf4;
1593  double f1mf43 = f1mf4*f1mf4*f1mf4;
1594 
1595  retVal = (
1596  (d4*f1mf22*f1mf32*f1mf4*f2mf3*f2mf4*f3mf4*(f12 + f2*f3 + (f2 + f3)*f4 + 2*f1*(f2 + f3 + f4)) + d1*f1mf2*f1mf3*f1mf4*f2mf3*f2mf42*f3mf42*(f1*f2 + f1*f3 + f2*f3 + 2*(f1 + f2 + f3)*f4 + f42) -
1597  5*f13*f24*f32*v1 + 4*f12*f25*f32*v1 + 5*f13*f22*f34*v1 - 2*f25*f34*v1 - 4*f12*f22*f35*v1 + 2*f24*f35*v1 + 10*f13*f24*f3*f4*v1 - 8*f12*f25*f3*f4*v1 - 5*f12*f24*f32*f4*v1 + 4*f1*f25*f32*f4*v1 -
1598  10*f13*f2*f34*f4*v1 + 5*f12*f22*f34*f4*v1 + 8*f12*f2*f35*f4*v1 - 4*f1*f22*f35*f4*v1 - 5*f13*f24*f42*v1 + 4*f12*f25*f42*v1 + 10*f12*f24*f3*f42*v1 - 8*f1*f25*f3*f42*v1 - 5*f1*f24*f32*f42*v1 +
1599  4*f25*f32*f42*v1 + 5*f13*f34*f42*v1 - 10*f12*f2*f34*f42*v1 + 5*f1*f22*f34*f42*v1 - 4*f12*f35*f42*v1 + 8*f1*f2*f35*f42*v1 - 4*f22*f35*f42*v1 - 5*f12*f24*f43*v1 + 4*f1*f25*f43*v1 -
1600  20*f13*f22*f3*f43*v1 + 10*f1*f24*f3*f43*v1 + 20*f13*f2*f32*f43*v1 - 5*f24*f32*f43*v1 + 5*f12*f34*f43*v1 - 10*f1*f2*f34*f43*v1 + 5*f22*f34*f43*v1 - 4*f1*f35*f43*v1 + 15*f13*f22*f44*v1 -
1601  5*f1*f24*f44*v1 - 2*f25*f44*v1 - 15*f13*f32*f44*v1 + 5*f1*f34*f44*v1 + 2*f35*f44*v1 - 10*f13*f2*f45*v1 - f12*f22*f45*v1 + 3*f24*f45*v1 + 10*f13*f3*f45*v1 + f12*f32*f45*v1 - 3*f34*f45*v1 +
1602  2*f12*f2*f46*v1 - f1*f22*f46*v1 - 2*f12*f3*f46*v1 + f1*f32*f46*v1 + 2*f1*f2*f47*v1 - f22*f47*v1 - 2*f1*f3*f47*v1 + f32*f47*v1 + f17*f32*v2 - 3*f15*f34*v2 + 2*f14*f35*v2 - 2*f17*f3*f4*v2 +
1603  f16*f32*f4*v2 + 5*f14*f34*f4*v2 - 4*f13*f35*f4*v2 + f17*f42*v2 - 2*f16*f3*f42*v2 + f15*f32*f42*v2 + f16*f43*v2 + 10*f15*f3*f43*v2 - 15*f14*f32*f43*v2 + 4*f1*f35*f43*v2 - 8*f15*f44*v2 +
1604  15*f13*f32*f44*v2 - 5*f1*f34*f44*v2 - 2*f35*f44*v2 + 8*f14*f45*v2 - 10*f13*f3*f45*v2 - f12*f32*f45*v2 + 3*f34*f45*v2 - f13*f46*v2 + 2*f12*f3*f46*v2 - f1*f32*f46*v2 - f12*f47*v2 +
1605  2*f1*f3*f47*v2 - f32*f47*v2 - f17*f22*v3 + 3*f15*f24*v3 - 2*f14*f25*v3 + 2*f17*f2*f4*v3 - f16*f22*f4*v3 - 5*f14*f24*f4*v3 + 4*f13*f25*f4*v3 - f17*f42*v3 + 2*f16*f2*f42*v3 - f15*f22*f42*v3 -
1606  f16*f43*v3 - 10*f15*f2*f43*v3 + 15*f14*f22*f43*v3 - 4*f1*f25*f43*v3 + 8*f15*f44*v3 - 15*f13*f22*f44*v3 + 5*f1*f24*f44*v3 + 2*f25*f44*v3 - 8*f14*f45*v3 + 10*f13*f2*f45*v3 + f12*f22*f45*v3 -
1607  3*f24*f45*v3 + f13*f46*v3 - 2*f12*f2*f46*v3 + f1*f22*f46*v3 + f12*f47*v3 - 2*f1*f2*f47*v3 + f22*f47*v3 +
1608  f1mf22*f1mf32*f2mf3*(2*f22*f32 + f13*(f2 + f3 - 2*f4) + f12*(f2 + f3 - 2*f4)*(2*(f2 + f3) + f4) - 4*(f22 + f2*f3 + f32)*f42 + 5*(f2 + f3)*f43 +
1609  f1*(4*f2*f3*(f2 + f3) - 4*(f22 + f2*f3 + f32)*f4 - 3*(f2 + f3)*f42 + 10*f43))*v4)/(f1mf22*f1mf32*f1mf43*f2mf3*f2mf42*f3mf42)
1610  );
1611 
1612  break;
1613  }
1614  default:
1615  {
1616  XLAL_ERROR_REAL8(XLAL_EINVAL, "Error in IMRPhenomXHM_Intermediate_Amp_delta3: IMRPhenomXIntermediateAmpVersion is not valid.\n");
1617  }
1618  }
1619 
1620  return retVal;
1621 }
1622 
1623 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)
1624 {
1625  double retVal;
1626 
1627  switch ( IntAmpFlag )
1628  {
1629  case 101: //linear, only v1, v2
1630  {
1631  retVal = 0.;
1632  break;
1633  }
1634  case 102: //quadratic: v1, v2, d2
1635  {
1636  retVal = 0.;
1637  break;
1638  }
1639  case 1032: // 2 freqs, points and derivatives: v1, v4, d1, d4
1640  {
1641  retVal = 0.;
1642  break;
1643  }
1644  case 103: // 4 freqs, no boundaries derivatives
1645  {
1646  retVal = 0.;
1647  break;
1648  }
1649  case 1043: //no left derivative: v1, v2, v3, v4, d4
1650  {
1651  double f12 = f1*f1;
1652  double f13 = f12*f1;
1653 
1654  double f22 = f2*f2;
1655  double f23 = f22*f2;
1656 
1657  double f32 = f3*f3;
1658  double f33 = f32*f3;
1659 
1660  double f42 = f4*f4;
1661  double f43 = f42*f4;
1662  double f44 = f43*f4;
1663 
1664  double f1mf2 = f1-f2;
1665  double f1mf3 = f1-f3;
1666  double f1mf4 = f1-f4;
1667  double f2mf3 = f2-f3;
1668  double f2mf4 = f2-f4;
1669  double f3mf4 = f3-f4;
1670 
1671  double f1mf42 = f1mf4*f1mf4;
1672  double f2mf42 = f2mf4*f2mf4;
1673  double f3mf42 = f3mf4*f3mf4;
1674 
1675  double v1mv3 = v1-v3;
1676  double v1mv4 = v1-v4;
1677  double v3mv4 = v3-v4;
1678 
1679  retVal = (-(d4*f1mf2*f1mf3*f1mf4*f2mf3*f2mf4*f3mf4) - f33*f42*v1 + 2*f32*f43*v1 - f3*f44*v1 - f13*f32*v2 + f12*f33*v2 + 2*f13*f3*f4*v2 - 2*f1*f33*f4*v2 - f13*f42*v2 - 3*f12*f3*f42*v2 + 3*f1*f32*f42*v2 +
1680  f33*f42*v2 + 2*f12*f43*v2 - 2*f32*f43*v2 - f1*f44*v2 + f3*f44*v2 + f13*f42*v3 - 2*f12*f43*v3 + f1*f44*v3 + f23*(f42*v1mv3 + f32*v1mv4 - 2*f3*f4*v1mv4 - f12*v3mv4 + 2*f1*f4*v3mv4) +
1681  f2*f4*(f43*v1mv3 + 2*f33*v1mv4 - 3*f32*f4*v1mv4 - 2*f13*v3mv4 + 3*f12*f4*v3mv4) + f22*(-2*f43*v1mv3 - f33*v1mv4 + 3*f3*f42*v1mv4 + f13*v3mv4 - 3*f1*f42*v3mv4) + f13*f32*v4 - f12*f33*v4 -
1682  2*f13*f3*f4*v4 + 2*f1*f33*f4*v4 + 3*f12*f3*f42*v4 - 3*f1*f32*f42*v4)/(f1mf2*f1mf3*f1mf42*f2mf3*f2mf42*f3mf42);
1683  break;
1684  }
1685  case 1042: //4th order poly: v1,d1, v2,d2, v3 // used for the first intermediate region
1686  {
1687  double f12 = f1*f1;
1688  double f13 = f12*f1;
1689 
1690  double f42 = f4*f4;
1691  double f43 = f42*f4;
1692 
1693  double f32 = f3*f3;
1694  double f33 = f32*f3;
1695 
1696  double f1mf4 = f1-f4;
1697  double f1mf3 = f1-f3;
1698  double f3mf4 = f3-f4;
1699 
1700  double f1mf42 = f1mf4*f1mf4;
1701  double f1mf32 = f1mf3*f1mf3;
1702  double f3mf42 = f3mf4*f3mf4;
1703 
1704  double f1mf43 = f1mf42*f1mf4;
1705 
1706  retVal = (-(d4*f1mf32*f1mf4*f3mf4) + d1*f1mf3*f1mf4*f3mf42 - 3*f1*f32*v1 + 2*f33*v1 + 6*f1*f3*f4*v1 - 3*f32*f4*v1 - 3*f1*f42*v1 + f43*v1 + f13*v3 - 3*f12*f4*v3 + 3*f1*f42*v3 - f43*v3 -
1707  f1mf32*(f1 + 2*f3 - 3*f4)*v4)/(f1mf32*f1mf43*f3mf42);
1708  break;
1709  }
1710  case 104: //Geraint's Version, 4th order poly: v1,d1, v2, v4,d4
1711  {
1712  double f12 = f1*f1;
1713  double f13 = f12*f1;
1714 
1715  double f22 = f2*f2;
1716  double f23 = f22*f2;
1717 
1718  double f42 = f4*f4;
1719  double f43 = f42*f4;
1720 
1721  double f1mf2 = f1-f2;
1722  double f1mf4 = f1-f4;
1723  double f2mf4 = f2-f4;
1724 
1725  double f1mf22 = f1mf2*f1mf2;
1726  double f2mf42 = f2mf4*f2mf4;
1727  double f1mf43 = f1mf4*f1mf4*f1mf4;
1728 
1729  retVal = ((-(d4*f1mf22*f1mf4*f2mf4) + d1*f1mf2*f1mf4*f2mf42 - 3*f1*f22*v1 + 2*f23*v1 + 6*f1*f2*f4*v1 - 3*f22*f4*v1
1730  - 3*f1*f42*v1 + f43*v1 + f13*v2 - 3*f12*f4*v2 + 3*f1*f42*v2 - f43*v2 - f1mf22*(f1 + 2*f2 - 3*f4)*v4)/(f1mf22*f1mf43*f2mf42));
1731 
1732  break;
1733  }
1734  case 105: // Geraint, standard way: v1, v2, v3, v4, d1, d4
1735  {
1736  double f12 = f1*f1;
1737  double f13 = f12*f1;
1738  double f14 = f13*f1;
1739  double f15 = f14*f1;
1740  double f16 = f15*f1;
1741 
1742  double f22 = f2*f2;
1743  double f23 = f22*f2;
1744  double f24 = f23*f2;
1745  double f25 = f24*f2;
1746 
1747  double f32 = f3*f3;
1748  double f33 = f32*f3;
1749  double f34 = f33*f3;
1750  double f35 = f34*f3;
1751 
1752  double f42 = f4*f4;
1753  double f43 = f42*f4;
1754  double f44 = f43*f4;
1755  double f45 = f44*f4;
1756  double f46 = f45*f4;
1757 
1758  double f1mf2 = f1-f2;
1759  double f1mf3 = f1-f3;
1760  double f1mf4 = f1-f4;
1761  double f2mf3 = f2-f3;
1762  double f2mf4 = f2-f4;
1763  double f3mf4 = f3-f4;
1764 
1765  double f1mf22 = f1mf2*f1mf2;
1766  double f1mf32 = f1mf3*f1mf3;
1767  double f2mf42 = f2mf4*f2mf4;
1768  double f3mf42 = f3mf4*f3mf4;
1769  double f1mf43 = f1mf4*f1mf4*f1mf4;
1770 
1771  retVal = (
1772  (-(d4*f1mf22*f1mf32*f1mf4*f2mf3*f2mf4*f3mf4*(2*f1 + f2 + f3 + f4)) - d1*f1mf2*f1mf3*f1mf4*f2mf3*f2mf42*f3mf42*(f1 + f2 + f3 + 2*f4) + 5*f13*f23*f32*v1 - 3*f1*f25*f32*v1 - 5*f13*f22*f33*v1 +
1773  2*f25*f33*v1 + 3*f1*f22*f35*v1 - 2*f23*f35*v1 - 10*f13*f23*f3*f4*v1 + 6*f1*f25*f3*f4*v1 + 5*f12*f23*f32*f4*v1 - 3*f25*f32*f4*v1 + 10*f13*f2*f33*f4*v1 - 5*f12*f22*f33*f4*v1 -
1774  6*f1*f2*f35*f4*v1 + 3*f22*f35*f4*v1 + 5*f13*f23*f42*v1 - 3*f1*f25*f42*v1 + 15*f13*f22*f3*f42*v1 - 10*f12*f23*f3*f42*v1 - 15*f13*f2*f32*f42*v1 + 5*f1*f23*f32*f42*v1 - 5*f13*f33*f42*v1 +
1775  10*f12*f2*f33*f42*v1 - 5*f1*f22*f33*f42*v1 + 3*f1*f35*f42*v1 - 10*f13*f22*f43*v1 + 5*f12*f23*f43*v1 + f25*f43*v1 + 15*f12*f22*f3*f43*v1 - 10*f1*f23*f3*f43*v1 + 10*f13*f32*f43*v1 -
1776  15*f12*f2*f32*f43*v1 + 5*f23*f32*f43*v1 - 5*f12*f33*f43*v1 + 10*f1*f2*f33*f43*v1 - 5*f22*f33*f43*v1 - f35*f43*v1 + 5*f13*f2*f44*v1 - 10*f12*f22*f44*v1 + 5*f1*f23*f44*v1 - 5*f13*f3*f44*v1 +
1777  10*f12*f32*f44*v1 - 5*f1*f33*f44*v1 + 5*f12*f2*f45*v1 + 2*f1*f22*f45*v1 - 3*f23*f45*v1 - 5*f12*f3*f45*v1 - 2*f1*f32*f45*v1 + 3*f33*f45*v1 - 4*f1*f2*f46*v1 + 2*f22*f46*v1 + 4*f1*f3*f46*v1 -
1778  2*f32*f46*v1 - 2*f16*f32*v2 + 3*f15*f33*v2 - f13*f35*v2 + 4*f16*f3*f4*v2 - 2*f15*f32*f4*v2 - 5*f14*f33*f4*v2 + 3*f12*f35*f4*v2 - 2*f16*f42*v2 - 5*f15*f3*f42*v2 + 10*f14*f32*f42*v2 -
1779  3*f1*f35*f42*v2 + 4*f15*f43*v2 - 5*f14*f3*f43*v2 + f35*f43*v2 + 5*f13*f3*f44*v2 - 10*f12*f32*f44*v2 + 5*f1*f33*f44*v2 - 4*f13*f45*v2 + 5*f12*f3*f45*v2 + 2*f1*f32*f45*v2 - 3*f33*f45*v2 +
1780  2*f12*f46*v2 - 4*f1*f3*f46*v2 + 2*f32*f46*v2 + 2*f16*f22*v3 - 3*f15*f23*v3 + f13*f25*v3 - 4*f16*f2*f4*v3 + 2*f15*f22*f4*v3 + 5*f14*f23*f4*v3 - 3*f12*f25*f4*v3 + 2*f16*f42*v3 +
1781  5*f15*f2*f42*v3 - 10*f14*f22*f42*v3 + 3*f1*f25*f42*v3 - 4*f15*f43*v3 + 5*f14*f2*f43*v3 - f25*f43*v3 - 5*f13*f2*f44*v3 + 10*f12*f22*f44*v3 - 5*f1*f23*f44*v3 + 4*f13*f45*v3 - 5*f12*f2*f45*v3 -
1782  2*f1*f22*f45*v3 + 3*f23*f45*v3 - 2*f12*f46*v3 + 4*f1*f2*f46*v3 - 2*f22*f46*v3 -
1783  f1mf22*f1mf32*f2mf3*(2*f2*f3*(f2 + f3) + 2*f12*(f2 + f3 - 2*f4) - 3*(f22 + f2*f3 + f32)*f4 + f1*(f22 + 5*f2*f3 + f32 - 6*(f2 + f3)*f4 + 5*f42) + 5*f43)*v4)/
1784  (f1mf22*f1mf32*f1mf43*f2mf3*f2mf42*f3mf42)
1785  );
1786 
1787  break;
1788  }
1789  default:
1790  {
1791  XLAL_ERROR_REAL8(XLAL_EINVAL, "Error in IMRPhenomXHM_Intermediate_Amp_delta4: IMRPhenomXIntermediateAmpVersion is not valid.\n");
1792  }
1793  }
1794 
1795  return retVal;
1796 }
1797 
1798 static double IMRPhenomXHM_Intermediate_Amp_delta5(double d1, double d4, double v1, double v2, double v3, double v4, double f1, double f2, double f3, double f4, int IntAmpFlag)
1799 {
1800  double retVal;
1801 
1802  switch ( IntAmpFlag )
1803  {
1804  case 101:
1805  {
1806  retVal = 0.;
1807  break;
1808  }
1809  case 102: //quadratic: v1, v2, d2
1810  {
1811  retVal = 0.;
1812  break;
1813  }
1814  case 1032: // 2 freqs, points and derivatives: v1, v4, d1, d4
1815  {
1816  //printf("\nIMRPhenomXHM_Intermediate_Amp_delta5 = 0 \r\n");
1817  retVal = 0.;
1818  break;
1819  }
1820  case 103: // 4 freqs, no boundaries derivatives
1821  {
1822  retVal = 0.;
1823  break;
1824  }
1825  case 1043: //no left derivative: v1, v2, v3, v4, d4
1826  {
1827  retVal = 0.;
1828  break;
1829  }
1830  case 1042: //4th order poly: v1,d1, v2,d2, v3 // used for the first intermediate region
1831  {
1832  retVal = 0.0;
1833  break;
1834  }
1835  case 104: //Geraint's Version, 4th order poly: v1,d1, v2,d2, v3
1836  {
1837  retVal = 0.0;
1838  break;
1839  }
1840  case 105: // Geraint, standard way: v1, v2, v3, v4, d1, d4
1841  {
1842  double f12 = f1*f1;
1843  double f13 = f12*f1;
1844  double f14 = f13*f1;
1845  double f15 = f14*f1;
1846 
1847  double f22 = f2*f2;
1848  double f23 = f22*f2;
1849  double f24 = f23*f2;
1850 
1851  double f32 = f3*f3;
1852  double f33 = f32*f3;
1853  double f34 = f33*f3;
1854 
1855  double f42 = f4*f4;
1856  double f43 = f42*f4;
1857  double f44 = f43*f4;
1858  double f45 = f44*f4;
1859 
1860  double f1mf2 = f1-f2;
1861  double f1mf3 = f1-f3;
1862  double f1mf4 = f1-f4;
1863  double f2mf3 = f2-f3;
1864  double f2mf4 = f2-f4;
1865  double f3mf4 = f3-f4;
1866 
1867  double f1mf22 = f1mf2*f1mf2;
1868  double f1mf32 = f1mf3*f1mf3;
1869  double f2mf42 = f2mf4*f2mf4;
1870  double f3mf42 = f3mf4*f3mf4;
1871  double f1mf43 = f1mf4*f1mf4*f1mf4;
1872 
1873  retVal = (
1874  (d4*f1mf22*f1mf32*f1mf4*f2mf3*f2mf4*f3mf4 + d1*f1mf2*f1mf3*f1mf4*f2mf3*f2mf42*f3mf42 - 4*f12*f23*f32*v1 + 3*f1*f24*f32*v1 + 4*f12*f22*f33*v1 - 2*f24*f33*v1 - 3*f1*f22*f34*v1 + 2*f23*f34*v1 +
1875  8*f12*f23*f3*f4*v1 - 6*f1*f24*f3*f4*v1 - 4*f1*f23*f32*f4*v1 + 3*f24*f32*f4*v1 - 8*f12*f2*f33*f4*v1 + 4*f1*f22*f33*f4*v1 + 6*f1*f2*f34*f4*v1 - 3*f22*f34*f4*v1 - 4*f12*f23*f42*v1 +
1876  3*f1*f24*f42*v1 - 12*f12*f22*f3*f42*v1 + 8*f1*f23*f3*f42*v1 + 12*f12*f2*f32*f42*v1 - 4*f23*f32*f42*v1 + 4*f12*f33*f42*v1 - 8*f1*f2*f33*f42*v1 + 4*f22*f33*f42*v1 - 3*f1*f34*f42*v1 +
1877  8*f12*f22*f43*v1 - 4*f1*f23*f43*v1 - f24*f43*v1 - 8*f12*f32*f43*v1 + 4*f1*f33*f43*v1 + f34*f43*v1 - 4*f12*f2*f44*v1 - f1*f22*f44*v1 + 2*f23*f44*v1 + 4*f12*f3*f44*v1 + f1*f32*f44*v1 -
1878  2*f33*f44*v1 + 2*f1*f2*f45*v1 - f22*f45*v1 - 2*f1*f3*f45*v1 + f32*f45*v1 + f15*f32*v2 - 2*f14*f33*v2 + f13*f34*v2 - 2*f15*f3*f4*v2 + f14*f32*f4*v2 + 4*f13*f33*f4*v2 - 3*f12*f34*f4*v2 +
1879  f15*f42*v2 + 4*f14*f3*f42*v2 - 8*f13*f32*f42*v2 + 3*f1*f34*f42*v2 - 3*f14*f43*v2 + 8*f12*f32*f43*v2 - 4*f1*f33*f43*v2 - f34*f43*v2 + 3*f13*f44*v2 - 4*f12*f3*f44*v2 - f1*f32*f44*v2 +
1880  2*f33*f44*v2 - f12*f45*v2 + 2*f1*f3*f45*v2 - f32*f45*v2 - f15*f22*v3 + 2*f14*f23*v3 - f13*f24*v3 + 2*f15*f2*f4*v3 - f14*f22*f4*v3 - 4*f13*f23*f4*v3 + 3*f12*f24*f4*v3 - f15*f42*v3 -
1881  4*f14*f2*f42*v3 + 8*f13*f22*f42*v3 - 3*f1*f24*f42*v3 + 3*f14*f43*v3 - 8*f12*f22*f43*v3 + 4*f1*f23*f43*v3 + f24*f43*v3 - 3*f13*f44*v3 + 4*f12*f2*f44*v3 + f1*f22*f44*v3 - 2*f23*f44*v3 +
1882  f12*f45*v3 - 2*f1*f2*f45*v3 + f22*f45*v3 + f1mf22*f1mf32*f2mf3*(2*f2*f3 + f1*(f2 + f3 - 2*f4) - 3*(f2 + f3)*f4 + 4*f42)*v4)/(f1mf22*f1mf32*f1mf43*f2mf3*f2mf42*f3mf42)
1883  );
1884  break ;
1885  }
1886  default:
1887  {
1888  XLAL_ERROR_REAL8(XLAL_EINVAL, "Error in IMRPhenomXHM_Intermediate_Amp_delta5: IMRPhenomXIntermediateAmpVersion is not valid.\n");
1889  }
1890  }
1891  return retVal;
1892 }
1893 
1894 
1895 /*********************************************/
1896 /* INTERMEDIATE AMPLITUDE ANSATZ */
1897 /*********************************************/
1898 
1899 // Build the polynomial with the coefficients given and return the inverse of the polynomial (this is the ansatz)
1900 static double IMRPhenomXHM_Intermediate_Amp_Ansatz(IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
1901 {
1902  if(pWFHM->IMRPhenomXHMReleaseVersion != 122019){
1903  double result = 0., fpower = 1.;
1904  /* Ansatz = f^(-7/6) * polynomial */
1905  for (UINT2 i = 0; i < pAmp->nCoefficientsInter; i++){
1906  //printf("pAmp->nCoefficientsInter = %i, pAmp->InterCoefficient[%i] = %.16e\n", pAmp->nCoefficientsInter, i, pAmp->InterCoefficient[i]);
1907  result += (pAmp->InterCoefficient[i] * fpower);
1908  fpower *= powers_of_f->itself;
1909  }
1910  result *= powers_of_f->m_seven_sixths;
1911  return result;
1912  }
1913  else{
1914  double a0 = pAmp->delta0;
1915  double a1 = pAmp->delta1;
1916  double a2 = pAmp->delta2;
1917  double a3 = pAmp->delta3;
1918  double a4 = pAmp->delta4;
1919  double a5 = pAmp->delta5;
1920  double polynomial;
1921  int InterAmpPolOrder = pAmp->InterAmpPolOrder;
1922 
1923  switch ( InterAmpPolOrder )
1924  {
1925  case 101: // linear order
1926  {
1927  double ff = powers_of_f->itself;
1928  polynomial = a0 + a1*ff ;
1929  break;
1930  }
1931  case 102: // quadratic order
1932  {
1933  double ff = powers_of_f->itself;
1934  double ff2 = powers_of_f->two;
1935  polynomial = a0 + a1*ff + a2*ff2;
1936  break;
1937  }
1938  case 103: //cubic order
1939  {
1940  double ff = powers_of_f->itself;
1941  double ff2 = powers_of_f->two;
1942  double ff3 = powers_of_f->three;
1943  polynomial = a0 + a1*ff + a2*ff2 + a3*ff3;
1944  break;
1945  }
1946  case 1042: // 4th order, used for the first intermediate region
1947  {
1948  double ff = powers_of_f->itself;
1949  double ff2 = powers_of_f->two;
1950  double ff3 = powers_of_f->three;
1951  double ff4 = powers_of_f->four;
1952  a0 = pAmp->alpha0;
1953  a1 = pAmp->alpha1;
1954  a2 = pAmp->alpha2;
1955  a3 = pAmp->alpha3;
1956  a4 = pAmp->alpha4;
1957  polynomial = a0 + a1*ff + a2*ff2 + a3*ff3 + a4*ff4;
1958  break;
1959  }
1960  case 104: // 4th order
1961  {
1962  double ff = powers_of_f->itself;
1963  double ff2 = powers_of_f->two;
1964  double ff3 = powers_of_f->three;
1965  double ff4 = powers_of_f->four;
1966  polynomial = a0 + a1*ff + a2*ff2 + a3*ff3 + a4*ff4;
1967  break;
1968  }
1969  case 105: // 5th order
1970  {
1971  double ff = powers_of_f->itself;
1972  double ff2 = powers_of_f->two;
1973  double ff3 = powers_of_f->three;
1974  double ff4 = powers_of_f->four;
1975  double ff5 = powers_of_f->five;
1976  polynomial = a0 + a1*ff + a2*ff2 + a3*ff3 + a4*ff4 + a5*ff5;
1977  break;
1978  }
1979  default:
1980  {
1981  XLAL_ERROR_REAL8(XLAL_EINVAL, "Error in IMRPhenomXHM_Intermediate_Amp_Ansatz: InterAmpPolOrder is not valid.\n");
1982  }
1983  }
1984  return pAmp->ampNorm / polynomial;
1985  }
1986 }
1987 
1988 
1989 /**** VETO functions ****/
1990 
1991 // Utility functions used to decide how many collocation points and which kind of reconstruction we do in the intermediate
1992 
1993 
1994 // Remove too low collocation points (heuristic)
1995 void IMRPhenomXHM_Intermediate_Amplitude_Veto(double *int1, double *int2, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22){
1996  double threshold = 0.2/(pWF22->ampNorm);
1997  if( 1./(*int1) < threshold )
1998  {
1999  *int1 = 1.;
2000  pWFHM->IMRPhenomXHMIntermediateAmpVersion = 1042;
2001  if( 1./(*int2) < threshold ){
2002  *int2 = 1.;
2003  pWFHM->IMRPhenomXHMIntermediateAmpVersion = 1032;
2004  }
2005  }else if( 1./(*int2) < threshold )
2006  {
2007  *int2 = 1.;
2008  pWFHM->IMRPhenomXHMIntermediateAmpVersion = 1042;
2009  }
2010 }
2011 
2012 // Check if a particular frequency belong to an interval
2013 int InsideInterval(double ftest, double fmin, double fmax){
2014 
2015  if(ftest >= fmin && ftest <= fmax){
2016  return 1;
2017  }else{
2018  return 0;
2019  }
2020 }
2021 
2022 // Check if the 2nd order polynomial crosses zero
2023 int CrossZeroP2(double a0, double a1, double a2, double fstart, double fend){
2024 
2025  double complex f1, f2;
2026  double discriminant = -4*a0*a2 + a1*a1;
2027 
2028  f1 = creall((-a1 - sqrt(discriminant))/(2.*a2));
2029  f2 = creall((-a1 + sqrt(discriminant))/(2.*a2));
2030 
2031  if( discriminant >= 0 && (InsideInterval(f1, fstart, fend) || InsideInterval(f2, fstart, fend) )){
2032  return 1;
2033  }else{
2034  return 0;
2035  }
2036 }
2037 
2038 // Check if the 3th order polynomial crosses zero
2039 int CrossZeroP3(double a0, double a1, double a2, double a3, double fstart, double fend)
2040 {
2041 
2042  double q1 = -a2/(3.*a3);
2043  double q2 = -a2*a2 + 3*a1*a3;
2044  double q3 = -2*a2*a2*a2 + 9*a1*a2*a3-27*a0*a3*a3;
2045  double discri = 4*q2*q2*q2 + q3*q3;
2046  double complex onethirdpower = cpow( q3 + csqrt(discri) , 1/3.);
2047  double complex f1, f2, f3;
2048  double twothird = pow(2,1/3.);
2049  double complex z1 = (1. + I*sqrt(3.)), z2 = (1. - I*sqrt(3.));
2050  int i1=0, i2=0, i3=0;
2051 
2052  // These are the points where the polynomial is zero
2053  f1 = q1 - q2*twothird/(3*a3*onethirdpower) + onethirdpower/(3.*a3*twothird);
2054  f2 = q1 + q2*z1/(3*twothird*twothird*a3*onethirdpower) - onethirdpower*z2/(6*twothird*a3);
2055  f3 = q1 + q2*z2/(3*twothird*twothird*a3*onethirdpower) - onethirdpower*z1/(6*twothird*a3);
2056 
2057  // If the solution is real and lay inside the interval then we return true.
2058  // Sometimes due to finite precission the imaginary part can be not exactly zero and we set the threshold 10^-15 as limiting value.
2059  double threshold = pow(10.,-15);
2060  if (fabs(cimag(f1)) < threshold ){
2061  i1 = InsideInterval(creall(f1), fstart, fend);
2062  }
2063  if (fabs(cimag(f2)) < threshold ){
2064  i2 = InsideInterval(creall(f2), fstart, fend);
2065  }
2066  if (fabs(cimag(f3)) < threshold ){
2067  i3 = InsideInterval(creall(f3), fstart, fend);
2068  }
2069  #if DEBUG == 1
2070  printf("\nCrossZeroP3: Coefficients = %.16f %.16f %.16f %.16f",a0, a1, a2, a3);
2071  printf("\nCrossZeroP3: discri = %.16f ", discri);
2072  printf("\nCrossZeroP3: Imag f1, f2, f3 = %.16e %.16e %.16e ", fabs(cimag(f1)), fabs(cimag(f2)), fabs(cimag(f3)));
2073  printf("\nCrossZeroP3: a3 = %.16f %.16f", creal(a3), cimag(a3));
2074  printf("\nCrossZeroP3: onethirdpower = %.16f %.16f", creal(onethirdpower), cimag(onethirdpower));
2075  printf("\nCrossZeroP3: f1 = %.16e %.16e", creal(f1), cimag(f1));
2076  printf("\nCrossZeroP3: f2 = %.16e %.16e", creal(f2), cimag(f2));
2077  printf("\nCrossZeroP3: f3 = %.16e %.16e", creal(f3), cimag(f3));
2078  printf("\nCrossZeroP3: i1, i2, i3 = %i %i %i", i1, i2, i3);
2079  #endif
2080  if (i1 == 0 && i2 == 0 && i3 == 0 ){
2081  return 0;
2082  }else{
2083  return 1;
2084  }
2085 }
2086 
2087 // Check if the 4th order polynomial crosses zero
2088 int CrossZeroP4(double a0, double a1, double a2, double a3, double a4, double fstart, double fend){
2089 
2090  double q0 = -a3/(4*a4);
2091  double q1 = a3*a3/(4*a4*a4) -2*a2/(3*a4);
2092  double q1b = 2*q1;
2093  double q2 = a2*a2 - 3*a1*a3 + 12*a0*a4;
2094  double q3 = 2*a2*a2*a2 - 9*a1*a2*a3 + 27*a0*a3*a3 + 27*a1*a1*a4 - 72*a0*a2*a4;
2095  double complex squareroot = csqrt(-4*q2*q2*q2 + q3*q3);
2096  double complex onethird = cpow(q3 + squareroot , 1/3.);
2097  double twothird = pow(2,1/3.);
2098  double complex frac1 = twothird*q2/(3*a4*onethird);
2099  double complex frac2 = onethird/(3*twothird*a4);
2100  double complex bigdenom = 4*sqrt(q1 + frac1 + frac2);
2101  double bignum = -a3*a3*a3/(a4*a4*a4) + 4*a2*a3/(a4*a4) - 8*a1/a4;
2102  double threshold = pow(10.,-15);
2103 
2104  complex double f1, f2, f3, f4;
2105 
2106  // These are the solutions
2107  f1 = q0 - 0.5*csqrt(q1 + frac1 + frac2) - 0.5*csqrt(q1b - frac1 - frac2 - bignum/bigdenom);
2108  f2 = q0 - 0.5*csqrt(q1 + frac1 + frac2) + 0.5*csqrt(q1b - frac1 - frac2 - bignum/bigdenom);
2109  f3 = q0 + 0.5*csqrt(q1 + frac1 + frac2) - 0.5*csqrt(q1b - frac1 - frac2 + bignum/bigdenom);
2110  f4 = q0 + 0.5*csqrt(q1 + frac1 + frac2) + 0.5*csqrt(q1b - frac1 - frac2 + bignum/bigdenom);
2111 
2112  #if DEBUG == 1
2113  printf("\n***** CrossZeroP4 *********\n");
2114  printf("q0, q1, q1b, q2, q3 %.16f %.16f %.16f %.16f %.16f\n", q0, q1, q1b, q2, q3);
2115  printf("bigdenom bignum %.16f %.16f %.16f %.16f\n", creal(bigdenom), cimag(bigdenom), creal(bignum), cimag(bignum));
2116  printf("squareroot %.16f %.16f \n", creal(squareroot), cimag(squareroot));
2117  printf("frac1 frac2 %.16f %.16f %.16f %.16f \n", creal(frac1), cimag(frac1), creal(frac2), cimag(frac2));
2118  printf("twothird onethird %.16f %.16f %.16f %.16f", creal(twothird), cimag(twothird), creal(onethird), cimag(onethird));
2119  printf("\nfstart, fend = %.16f %.16f\n", fstart, fend);
2120  printf("\nf1 = %.16f %.16f", creal(f1), cimag(f1));
2121  printf("\nf2 = %.16f %.16f", creal(f2), cimag(f2));
2122  printf("\nf3 = %.16f %.16f", creal(f3), cimag(f3));
2123  printf("\nf4 = %.16f %.16f", creal(f4), cimag(f4));
2124  #endif
2125 
2126  int i1=0, i2=0, i3=0, i4=0;
2127 
2128  // Check if the soultions are real and lay in the interval
2129  if (fabs(cimag(f1)) < threshold ){
2130  i1 = InsideInterval(creall(f1), fstart, fend);
2131  }
2132  if (fabs(cimag(f2)) < threshold){
2133  i2 = InsideInterval(creall(f2), fstart, fend);
2134  }
2135  if (fabs(cimag(f3)) < threshold ){
2136  i3 = InsideInterval(creall(f3), fstart, fend);
2137  }
2138  if (fabs(cimag(f4)) < threshold ){
2139  i4 = InsideInterval(creall(f4), fstart, fend);
2140  }
2141 
2142  if (i1 == 0 && i2 == 0 && i3 == 0 && i4 == 0 ){
2143  return 0;
2144  }else{
2145  return 1;
2146  }
2147 }
2148 
2149 
2150 // Check if the 5th order polynomial crosses zero.
2151 // In this case there is no analytical solution and we have to check numerically
2152 int CrossZeroP5(double a0, double a1, double a2, double a3, double a4, double a5, double fstart, double fend)
2153 {//https://www.gnu.org/software/gsl/doc/html/poly.html
2154  double threshold = pow(10.,-15);
2155  /* coefficients of P(x) = -1 + x^5 */
2156  double a[6] = { a0, a1, a2, a3, a4, a5 };
2157  double z[10];
2158 
2159  gsl_poly_complex_workspace * w
2160  = gsl_poly_complex_workspace_alloc (6);
2161 
2162  gsl_poly_complex_solve (a, 6, w, z);
2163 
2164  gsl_poly_complex_workspace_free (w);
2165 
2166  #if DEBUG == 1
2167  for (int i = 0; i < 5; i++)
2168  {
2169 
2170  printf ("z%d = %+.18f %+.18f\n",
2171  i, z[2*i], z[2*i+1]);
2172  }
2173  #endif
2174 
2175  int i1=0, i2=0, i3=0, i4=0, i5=0;
2176 
2177  // Check if the soultions are real and lay in the interval
2178  if (fabs(z[1]) < threshold ){
2179  i1 = InsideInterval(z[0], fstart, fend);
2180  }
2181  if (fabs(z[3]) < threshold){
2182  i2 = InsideInterval(z[2], fstart, fend);
2183  }
2184  if (fabs(z[4]) < threshold ){
2185  i3 = InsideInterval(z[4], fstart, fend);
2186  }
2187  if (fabs(z[5]) < threshold ){
2188  i4 = InsideInterval(z[6], fstart, fend);
2189  }
2190  if (fabs(z[7]) < threshold ){
2191  i5 = InsideInterval(z[8], fstart, fend);
2192  }
2193 
2194  if (i1 == 0 && i2 == 0 && i3 == 0 && i4 == 0 && i5 == 0 ){
2195  return 0;
2196  }else{
2197  return 1;
2198  }
2199 
2200  return 0;
2201 }
2202 
2203 
2204 
2205 
2206 // Get the coefficients of the polynomial for a particular reconstruction indicated with IntAmpFlag
2207 void Update_Intermediate_Amplitude_Coefficients(IMRPhenomXHMAmpCoefficients *pAmp, int IntAmpFlag)
2208 {
2209  double d1 = pAmp->d1;
2210  double d4 = pAmp->d4;
2211  double v1 = pAmp->v1;
2212  double v2 = pAmp->v2;
2213  double v3 = pAmp->v3;
2214  double v4 = pAmp->v4;
2215  double f1 = pAmp->f1;
2216  double f2 = pAmp->f2;
2217  double f3 = pAmp->f3;
2218  double f4 = pAmp->f4;
2219  #if DEBUG == 1
2220  printf("\n UpdateCoeff IMRPhenomXHMIntermediateAmpVersion = %i", IntAmpFlag);
2221  #endif
2222  pAmp->delta0 = IMRPhenomXHM_Intermediate_Amp_delta0(d1,d4,v1,v2,v3,v4,f1,f2,f3,f4,IntAmpFlag);
2223  pAmp->delta1 = IMRPhenomXHM_Intermediate_Amp_delta1(d1,d4,v1,v2,v3,v4,f1,f2,f3,f4,IntAmpFlag);
2224  pAmp->delta2 = IMRPhenomXHM_Intermediate_Amp_delta2(d1,d4,v1,v2,v3,v4,f1,f2,f3,f4,IntAmpFlag);
2225  pAmp->delta3 = IMRPhenomXHM_Intermediate_Amp_delta3(d1,d4,v1,v2,v3,v4,f1,f2,f3,f4,IntAmpFlag);
2226  pAmp->delta4 = IMRPhenomXHM_Intermediate_Amp_delta4(d1,d4,v1,v2,v3,v4,f1,f2,f3,f4,IntAmpFlag);
2227  pAmp->delta5 = IMRPhenomXHM_Intermediate_Amp_delta5(d1,d4,v1,v2,v3,v4,f1,f2,f3,f4,IntAmpFlag);
2228 }
2229 
2230 // Check if the polynomial crosses zero wrapper.
2231 // If it crosses zero, then remove one collocation point and lower the order of the polynomial.
2232 // The order can be as low as linear order, when it is certain that it does not cross zero.
2233 void ChoosePolOrder(IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp){
2234  switch(pWFHM->IMRPhenomXHMIntermediateAmpVersion)
2235  {
2236  case 105: // v1, v2, v3, v4, d1, d4
2237  {
2238  #if DEBUG == 1
2239  printf("\nChoosePolOrder 105\n");
2240  #endif
2241  // struct pol5_params params;
2242  // params.a0 = pAmp->delta0;
2243  // params.a1 = pAmp->delta1;
2244  // params.a2 = pAmp->delta2;
2245  // params.a3 = pAmp->delta3;
2246  // params.a4 = pAmp->delta4;
2247  // params.a5 = pAmp->delta5;
2248 
2249  double fstart, fend;
2250  fstart = pAmp->fAmpMatchIN;
2251  fend = pAmp->fAmpMatchIM;
2252  #if DEBUG == 1
2253  printf("\n In ChoosePolOrder\n");
2254  #endif
2255  if(CrossZeroP5(pAmp->delta0, pAmp->delta1, pAmp->delta2, pAmp->delta3, pAmp->delta4, pAmp->delta5, fstart, fend)==1){
2256  //if (RootPol5_finder_gsl(params, fstart, fend) == 1){
2257  #if DEBUG == 1
2258  printf("\n Pol5 crosses zero \n");
2259  #endif
2260  //int dummy = RootPol5_finder_gsl(params, fstart, fend);
2261  //update Coefficients and intermediate version
2262  Update_Intermediate_Amplitude_Coefficients(pAmp, 1042);
2263  pWFHM->IMRPhenomXHMIntermediateAmpVersion=104;
2264  //check if the lower order cross zero
2265  ChoosePolOrder(pWFHM, pAmp);
2266  }else{
2267  pAmp->InterAmpPolOrder = 105;
2268  }
2269  break;
2270  }
2271  case 1043: //no left derivative: v1, v2, v3, v4, d4
2272  {
2273  #if DEBUG == 1
2274  printf("\nChoosePolOrder 1043\n");
2275  #endif
2276  double a0, a1, a2, a3, a4, fstart, fend;
2277  a0 = pAmp->delta0;
2278  a1 = pAmp->delta1;
2279  a2 = pAmp->delta2;
2280  a3 = pAmp->delta3;
2281  a4 = pAmp->delta4;
2282  fstart = pAmp->fAmpMatchIN;
2283  fend = pAmp->fAmpMatchIM;
2284 
2285  if(CrossZeroP4(a0, a1, a2, a3, a4, fstart, fend) == 1){
2286  //update Coefficients and intermediate version
2287  Update_Intermediate_Amplitude_Coefficients(pAmp, 1032);
2288  pWFHM->IMRPhenomXHMIntermediateAmpVersion=1032;
2289  //check if the lower order cross zero
2290  ChoosePolOrder(pWFHM, pAmp);
2291  }else{
2292  pAmp->InterAmpPolOrder = 104;
2293  }
2294  break;
2295  }
2296  case 1042: //Remove one inter collocation point: v1, d1, v3, v4, d4.
2297  {
2298  #if DEBUG == 1
2299  printf("\nChoosePolOrder 1042\n");
2300  #endif
2301  double a0, a1, a2, a3, a4, fstart, fend;
2302  a0 = pAmp->delta0;
2303  a1 = pAmp->delta1;
2304  a2 = pAmp->delta2;
2305  a3 = pAmp->delta3;
2306  a4 = pAmp->delta4;
2307  fstart = pAmp->fAmpMatchIN;
2308  fend = pAmp->fAmpMatchIM;
2309 
2310  if(CrossZeroP4(a0, a1, a2, a3, a4, fstart, fend) == 1){
2311  //update Coefficients and intermediate version
2312  Update_Intermediate_Amplitude_Coefficients(pAmp, 1032);
2313  pWFHM->IMRPhenomXHMIntermediateAmpVersion=1032;
2314  //check if the lower order cross zero
2315  ChoosePolOrder(pWFHM, pAmp);
2316  }else{
2317  pAmp->InterAmpPolOrder = 104;
2318  }
2319  break;
2320  }
2321  case 104: //Remove one inter collocation point: v1, d1, v2, v4,d4,
2322  {
2323  #if DEBUG == 1
2324  printf("\nChoosePolOrder 104\n");
2325  #endif
2326  double a0, a1, a2, a3, a4, fstart, fend;
2327  a0 = pAmp->delta0;
2328  a1 = pAmp->delta1;
2329  a2 = pAmp->delta2;
2330  a3 = pAmp->delta3;
2331  a4 = pAmp->delta4;
2332  fstart = pAmp->fAmpMatchIN;
2333  fend = pAmp->fAmpMatchIM;
2334 
2335  if(CrossZeroP4(a0, a1, a2, a3, a4, fstart, fend) == 1){
2336  //update Coefficients and intermediate version
2337  Update_Intermediate_Amplitude_Coefficients(pAmp, 1032);
2338  pWFHM->IMRPhenomXHMIntermediateAmpVersion=1032;
2339  //check if the lower order cross zero
2340  ChoosePolOrder(pWFHM, pAmp);
2341  }else{
2342  pAmp->InterAmpPolOrder = 104;
2343  }
2344  break;
2345  }
2346  case 1032: // 2 freqs, points and derivatives: v1, v4, d1, d4
2347  {
2348  #if DEBUG == 1
2349  printf("\nChoosePolOrder 1032\n");
2350  #endif
2351  double a0, a1, a2, a3, fstart, fend;
2352  a0 = pAmp->delta0;
2353  a1 = pAmp->delta1;
2354  a2 = pAmp->delta2;
2355  a3 = pAmp->delta3;
2356  fstart = pAmp->fAmpMatchIN;
2357  fend = pAmp->fAmpMatchIM;
2358 
2359  if(CrossZeroP3(a0, a1, a2, a3, fstart, fend) == 1){
2360  //update Coefficients and intermediate version
2361  #if DEBUG == 1
2362  printf("\nCrossZeroP3 True\n");
2363  #endif
2364  Update_Intermediate_Amplitude_Coefficients(pAmp, 102);
2365  pWFHM->IMRPhenomXHMIntermediateAmpVersion = 102;
2366  //check if the lower order cross zero
2367  ChoosePolOrder(pWFHM, pAmp);
2368  }else{
2369  #if DEBUG == 1
2370  printf("\nCrossZeroP3 False\n");
2371  #endif
2372  pAmp->InterAmpPolOrder = 103;
2373  }
2374  break;
2375  }
2376  case 103: // 4 freqs, no boundaries derivatives
2377  {
2378  #if DEBUG == 1
2379  printf("\nChoosePolOrder 103\n");
2380  #endif
2381  double a0, a1, a2, a3, fstart, fend;
2382  a0 = pAmp->delta0;
2383  a1 = pAmp->delta1;
2384  a2 = pAmp->delta2;
2385  a3 = pAmp->delta3;
2386  fstart = pAmp->fAmpMatchIN;
2387  fend = pAmp->fAmpMatchIM;
2388 
2389  if(
2390  CrossZeroP3(a0, a1, a2, a3, fstart, fend) == 1){
2391  //update Coefficients and intermediate version
2392  Update_Intermediate_Amplitude_Coefficients(pAmp, 102);
2393  pWFHM->IMRPhenomXHMIntermediateAmpVersion=102;
2394  //check if the lower order cross zero
2395  ChoosePolOrder(pWFHM, pAmp);
2396  }else{
2397  pAmp->InterAmpPolOrder = 103;
2398  }
2399  break;
2400  }
2401  case 102: //quadratic: v1, v2, d2
2402  {
2403  #if DEBUG == 1
2404  printf("\nChoosePolOrder 102\n");
2405  #endif
2406  double a0, a1, a2, fstart, fend;
2407  a0 = pAmp->delta0;
2408  a1 = pAmp->delta1;
2409  a2 = pAmp->delta2;
2410  fstart = pAmp->fAmpMatchIN;
2411  fend = pAmp->fAmpMatchIM;
2412 
2413  if(CrossZeroP2(a0, a1, a2, fstart, fend) == 1){
2414  //update Coefficients and intermediate version
2415  Update_Intermediate_Amplitude_Coefficients(pAmp, 101);
2416  pWFHM->IMRPhenomXHMIntermediateAmpVersion=101;
2417  //check if the lower order cross zero
2418  ChoosePolOrder(pWFHM, pAmp);
2419  }else{
2420  pAmp->InterAmpPolOrder = 102;
2421  }
2422  break;
2423  }
2424  case 101: //linear, only v1, v2
2425  {
2426  #if DEBUG == 1
2427  printf("\nChoosePolOrder 101\n");
2428  #endif
2429  //The linear reconstruction is not going to cross zero, since both boundaries to connect are positive
2430  pAmp->InterAmpPolOrder = 101;
2431  break;
2432  }
2433  }
2434 }
2435 
2436 static void IMRPhenomXHM_Intermediate_Amp_CollocationPoints(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22){
2437 
2438 
2439  /* Define collocation points frequencies */
2440  switch(pWFHM->IMRPhenomXHMIntermediateAmpFreqsVersion){
2441  case 0:{ // Equispaced. Get boundaries too
2442  REAL8 deltaf = (pAmp->fAmpMatchIM - pAmp->fAmpMatchIN) / (pWFHM->nCollocPtsInterAmp - 1);
2443  UINT2 idx = 0;
2444  for (UINT2 i = 0; i < pWFHM->nCollocPtsInterAmp; i++){
2445  if(pAmp->VersionCollocPtsInter[i] == 1){
2446  // Add point
2447  pAmp->CollocationPointsFreqsAmplitudeInter[idx] = pAmp->fAmpMatchIN + deltaf * i;
2448  }
2449  else if (pAmp->VersionCollocPtsInter[i] == 2){
2450  // Add point + derivative
2451  pAmp->CollocationPointsFreqsAmplitudeInter[idx] = pAmp->fAmpMatchIN + deltaf * i;
2452  pAmp->CollocationPointsFreqsAmplitudeInter[idx + 1] = pAmp->CollocationPointsFreqsAmplitudeInter[idx];
2453  }
2454  idx += pAmp->VersionCollocPtsInter[i];
2455  }
2456  break;
2457  }
2458  case 1:{ // Chebyshev. Get boundaries too
2459  REAL8 semisum = 0.5 * (pAmp->fAmpMatchIN + pAmp->fAmpMatchIM);
2460  REAL8 semidif = 0.5 * (pAmp->fAmpMatchIM - pAmp->fAmpMatchIN);
2461  for (INT4 i = pWFHM->nCollocPtsInterAmp + 1; i >=0; i--){
2462  pAmp->CollocationPointsFreqsAmplitudeInter[i] = semisum + semidif * cos( i * LAL_PI / pWFHM->nCollocPtsInterAmp );
2463  }
2464  break;
2465  }
2466  default: {XLAL_ERROR_VOID(XLAL_EDOM, "Error in IMRPhenomXHM_Intermediate_Amp_CollocationPoints: IMRPhenomXHMIntermediateAmpFreqsVersion = %i is not valid. \n", pWFHM->IMRPhenomXHMInspiralAmpFreqsVersion);}
2467  }
2468  /* Define values */
2469 
2470  IMRPhenomX_UsefulPowers powers_of_finsp;
2471  IMRPhenomX_Initialize_Powers(&powers_of_finsp, pAmp->fAmpMatchIN);
2472  // Be careful with TGR thing. Inspiral affecting intermediate region.
2473  UINT2 tmp_factor = pAmp->InspRescaleFactor;
2474  pAmp->InspRescaleFactor = pAmp->InterRescaleFactor;
2475  UINT2 tmpnCollocPts = 0;
2476  if (pAmp->VersionCollocPtsInter[0] == 1){
2477  pAmp->CollocationPointsValuesAmplitudeInter[0] = IMRPhenomXHM_Inspiral_Amp_Ansatz(&powers_of_finsp, pWFHM, pAmp);//FIXME: Rescale with InterFactor
2478  tmpnCollocPts += 1;
2479  }
2480  else if (pAmp->VersionCollocPtsInter[0] == 2){
2481  pAmp->CollocationPointsValuesAmplitudeInter[0] = IMRPhenomXHM_Inspiral_Amp_Ansatz(&powers_of_finsp, pWFHM, pAmp);
2482  pAmp->CollocationPointsValuesAmplitudeInter[1] = IMRPhenomXHM_Inspiral_Amp_NDAnsatz(&powers_of_finsp, pWFHM, pAmp);//pAmp->CollocationPointsValuesAmplitudeInsp[2];
2483  tmpnCollocPts += 2;
2484  }
2485  pAmp->InspRescaleFactor = tmp_factor;
2486 
2487  /* Call parameter space fits */
2488  /* pWFHM->nCollocPtsInterAmp also includes the boundaries, for the parameter space fits we just need the in between points */
2489  UINT2 idx = 0;
2490  for(UINT2 i = 1; i < pWFHM->nCollocPtsInterAmp - 1; i++){
2491  if(i <= 2)
2492  idx = pWFHM->modeInt * 2 + i - 1;
2493  else
2494  idx = pWFHM->modeInt * 2 + i - 3 + 16; //FIXME
2495  if(pAmp->VersionCollocPtsInter[i] == 1){
2496  pAmp->CollocationPointsValuesAmplitudeInter[tmpnCollocPts] = fabs(pAmp->IntermediateAmpFits[idx](pWF22, pWFHM->IMRPhenomXHMIntermediateAmpFitsVersion)); //FIXME: negative fits?
2497  tmpnCollocPts++;
2498  }
2499  // FIXME: throw error here if pAmp->VersionCollocPtsInter[i] == 2, cannot use derivatives if you don't have the parameter space fit
2500  }
2501 
2502  tmp_factor = pAmp->RDRescaleFactor;
2503  pAmp->RDRescaleFactor = pAmp->InterRescaleFactor;
2504  IMRPhenomX_UsefulPowers powers_of_fRD;
2505  IMRPhenomX_Initialize_Powers(&powers_of_fRD, pAmp->fAmpMatchIM);
2506  switch(pAmp->VersionCollocPtsInter[pWFHM->nCollocPtsInterAmp - 1]){
2507  case 1:{ // Add point
2508  if (pWFHM->MixingOn == 0){
2509  pAmp->CollocationPointsValuesAmplitudeInter[tmpnCollocPts] = IMRPhenomXHM_RD_Amp_Ansatz(&powers_of_fRD, pWFHM, pAmp);//pAmp->CollocationPointsValuesAmplitudeRD[0];
2510  }
2511  else{
2512  pAmp->CollocationPointsValuesAmplitudeInter[tmpnCollocPts] = cabs(SpheroidalToSpherical(&powers_of_fRD, pAmp22, pPhase22, pAmp, pPhase, pWFHM, pWF22));
2513  }
2514  tmpnCollocPts++;
2515  break;
2516  }
2517  case 2:{ // Add point + derivative
2518  if (pWFHM->MixingOn == 0){
2519  pAmp->CollocationPointsValuesAmplitudeInter[tmpnCollocPts] = IMRPhenomXHM_RD_Amp_Ansatz(&powers_of_fRD, pWFHM, pAmp);//pAmp->CollocationPointsValuesAmplitudeRD[0];
2520  pAmp->CollocationPointsValuesAmplitudeInter[tmpnCollocPts + 1] = IMRPhenomXHM_RD_Amp_DAnsatz(&powers_of_fRD, pWFHM, pAmp);//pAmp->CollocationPointsValuesAmplitudeRD[0];
2521  }
2522  else{
2523  pAmp->CollocationPointsValuesAmplitudeInter[tmpnCollocPts] = cabs(SpheroidalToSpherical(&powers_of_fRD, pAmp22, pPhase22, pAmp, pPhase, pWFHM, pWF22));
2524  pAmp->CollocationPointsValuesAmplitudeInter[tmpnCollocPts + 1] = IMRPhenomXHM_RD_Amp_NDAnsatz(&powers_of_fRD, pAmp, pPhase, pWFHM, pAmp22, pPhase22, pWF22);
2525  }
2526  tmpnCollocPts += 2;
2527  break;
2528  }
2529  default: {XLALPrintError("Error in IMRPhenomXHM_Intermediate_Amp_Coefficients: version %i is not valid.", pAmp->VersionCollocPtsInter[pAmp->nCoefficientsInter - 1]);}
2530  }
2531  pAmp->RDRescaleFactor = tmp_factor;
2532 
2533  /* tmpnCollocPts must be the same tahn pWFHM->nCollocPtsInterAmp + 2 (=number of free coefficients in intermediate ansatz) */
2534  if(tmpnCollocPts != pAmp->nCoefficientsInter)
2535  XLAL_ERROR_VOID(XLAL_EFUNC, "IMRPhenomXHM_Intermediate_Amp_CollocationPoints failed. Inconsistent number of free parameters %i, %i.", tmpnCollocPts, pAmp->nCoefficientsInter);
2536 }
2537 
2538 void IMRPhenomXHM_Intermediate_Amp_Coefficients(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22){
2539 
2540  /* Previously we checked that nCollocPtsInterAmp read from the IMRPhenomXHMIntermediateAmpVersion is equal to the number of free coefficients in the ansatz. */
2541  UINT2 nCollocPtsInterAmp = pAmp->nCoefficientsInter;
2542 
2543  /* Define set of collocation points */
2544  IMRPhenomXHM_Intermediate_Amp_CollocationPoints(pAmp, pWFHM, pWF22, pPhase, pAmp22, pPhase22);
2545 
2546  /* GSL objects for solving system of equations via LU decomposition */
2547  gsl_vector *b, *x;
2548  gsl_matrix *A;
2549  gsl_permutation *p;
2550  int signum; // No need to set, used internally by gsl_linalg_LU_decomp
2551 
2552  /* Initialize gsl objects */
2553  p = gsl_permutation_alloc(nCollocPtsInterAmp);
2554  b = gsl_vector_alloc(nCollocPtsInterAmp);
2555  x = gsl_vector_alloc(nCollocPtsInterAmp);
2556  A = gsl_matrix_alloc(nCollocPtsInterAmp, nCollocPtsInterAmp);
2557 
2558 
2559  /* Define linear system of equations: A x = b */
2560  /* x is the solution vector: the coefficients of the intermediate ansatz */
2561  /* b is the vector of collocation points for a set of frequencies */
2562  /* A is the matrix of multiplicative factors to each coefficient of the ansatz.
2563  Each row gives the ansatz evaluated at a collocation point frequency. */
2564  UINT2 tmpnCollocPts = 0;
2565  for(UINT2 i = 0; i < pWFHM->nCollocPtsInterAmp; i++){
2566  /* Skip the 0 cases, means that collocation point is not used */
2567  if(pAmp->VersionCollocPtsInter[i] > 0){
2568  // Set b vector
2569  gsl_vector_set(b, tmpnCollocPts, pAmp->CollocationPointsValuesAmplitudeInter[tmpnCollocPts]);
2570  //FIXME: distinguish InterAmp ansatzaes versions
2571  // Set system matrix: Polynomial/f^7/6 at the collocation points frequencies.
2572  /* A = (1, f1, f1^2, f1^3, f1^4, ...) * f1^(-7/6)
2573  (1, f2, f2^2, f2^3, f2^4, ...) * f2^(-7/6)
2574  ....
2575  Until number of collocation points
2576  If one of the conditions is a derivative we substitute by
2577  (0, 1, 2*fi, 3*fi^2, 4*fi^3, ...)
2578  */
2579  REAL8 fcollpoint = pAmp->CollocationPointsFreqsAmplitudeInter[tmpnCollocPts];
2580  REAL8 seven_sixths = 7/6.;
2581  REAL8 fcollpoint_m_seven_sixths = pow(fcollpoint, -seven_sixths);
2582  REAL8 fpower = 1.; // 1, f, f^2, f^3, f^4, ...
2583 
2584  // Add equation for Point. Assuming I will always use point, not the derivative alone
2585  for(INT4 j = 0; j < nCollocPtsInterAmp; j++){
2586  gsl_matrix_set(A, tmpnCollocPts, j, fpower * fcollpoint_m_seven_sixths);
2587  fpower *= fcollpoint;
2588  }
2589  tmpnCollocPts++;
2590  // Add equation for Derivative
2591  if (pAmp->VersionCollocPtsInter[i] == 2 ){
2592  gsl_vector_set(b, tmpnCollocPts, pAmp->CollocationPointsValuesAmplitudeInter[tmpnCollocPts]);
2593  fpower = 1./fcollpoint;
2594  for(INT4 j = 0; j < nCollocPtsInterAmp; j++){
2595  REAL8 derivative = (j - seven_sixths) * fpower * fcollpoint_m_seven_sixths;
2596  gsl_matrix_set(A, tmpnCollocPts, j, derivative);
2597  fpower *= fcollpoint;
2598  }
2599  tmpnCollocPts++;
2600  }
2601  } /* End non-zero if statement */
2602  } /* End of loop over number of free coefficients. System of equations setup. */
2603 
2604  /* tmpnCollocPts must be the same than the number of free coefficients and to pAmp->nCoefficientsInter */
2605  if (tmpnCollocPts != pAmp->nCoefficientsInter ){
2606  /* Free gsl variables */
2607  gsl_vector_free(b);
2608  gsl_vector_free(x);
2609  gsl_matrix_free(A);
2610  gsl_permutation_free(p);
2611  XLAL_ERROR_VOID(XLAL_EFUNC, "IMRPhenomXHM_Intermediate_Amp_Coefficients failed. Inconsistent number of collocation points (%i) and free parameters (%i).", tmpnCollocPts, pAmp->nCoefficientsInter);
2612  }
2613 
2614  /* We now solve the system A x = b via an LU decomposition. x is the solution vector */
2615  gsl_linalg_LU_decomp(A, p, &signum);
2616  gsl_linalg_LU_solve(A, p, b, x);
2617 
2618  /* The solution corresponds to the coefficients of the ansatz */
2619  for (UINT2 i = 0; i < pAmp->nCoefficientsInter; i++){
2620  pAmp->InterCoefficient[i] = gsl_vector_get(x, i);
2621  }
2622 
2623  /* Free gsl variables */
2624  gsl_vector_free(b);
2625  gsl_vector_free(x);
2626  gsl_matrix_free(A);
2627  gsl_permutation_free(p);
2628 
2629 }
2630 
2631 
2632 /***********************************************/
2633 /* */
2634 /* PHASE */
2635 /* */
2636 /***********************************************/
2637 
2638 // Fits of phase derivatives at collocation points for each mode
2639 
2640 /* Start of Phase Parameter Space Fits */
2641 
2642 static double IMRPhenomXHM_Inter_Phase_21_p1(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2643  double total=0;
2644  switch (InterPhaseFlag){
2645  case 122019:{
2646  double eta = pWF->eta;
2647  double S = pWF->STotR;
2648  double eta2,eta3,eta4,eta5,eta6,S2,S3,S4;
2649  eta2 = pow(eta,2);
2650  eta3 = pow(eta,3);
2651  eta4 = pow(eta,4);
2652  eta5 = pow(eta,5);
2653  eta6 = pow(eta,6);
2654  S2 = pow(S,2);
2655  S3 = pow(S,3);
2656  S4 = pow(S,4);
2657  double noSpin = 4045.84 + 7.63226/eta - 1956.93*eta - 23428.1*eta2 + 369153.*eta3 - 2.28832e6*eta4 + 6.82533e6*eta5 - 7.86254e6*eta6;
2658  double eqSpin = - 347.273*S + 83.5428*S2 - 355.67*S3 + (4.44457*S + 16.5548*S2 + 13.6971*S3)/eta + eta*( - 79.761*S - 355.299*S2 + 1114.51*S3 - 1077.75*S4) + 92.6654*S4 + eta2*(- 619.837*S - 722.787*S2 + 2392.73*S3 + 2689.18*S4);
2659  double uneqSpin = ( 918.976*pWF->chi1L*sqrt(1.-4.*eta) - 918.976*pWF->chi2L*sqrt(1.-4.*eta))*eta + ( 91.7679*pWF->chi1L*sqrt(1.-4.*eta) - 91.7679*pWF->chi2L*sqrt(1.-4.*eta))*eta2;
2660  total = noSpin + eqSpin + uneqSpin;
2661  break;
2662  }
2663  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_21_p1: version is not valid. Recommended version is 122019.");}
2664  }
2665  return total;
2666 }
2667 
2668 static double IMRPhenomXHM_Inter_Phase_33_p1(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2669  double total=0;
2670  switch (InterPhaseFlag){
2671  case 122019:{
2672  double eta = pWF->eta;
2673  double S = pWF->STotR;
2674  double eta2,eta3,eta4,eta5,eta6,S2;
2675  eta2 = pow(eta,2);
2676  eta3 = pow(eta,3);
2677  eta4 = pow(eta,4);
2678  eta5 = pow(eta,5);
2679  eta6 = pow(eta,6);
2680  S2 = pow(S,2);
2681  double noSpin = 4360.19 + 4.27128/eta - 8727.4*eta + 18485.9*eta2 + 371303.00000000006*eta3 - 3.22792e6*eta4 + 1.01799e7*eta5 - 1.15659e7*eta6;
2682  double eqSpin = ((11.6635 - 251.579*eta - 3255.6400000000003*eta2 + 19614.6*eta3 - 34860.2*eta4)*S + (14.8017 + 204.025*eta - 5421.92*eta2 + 36587.3*eta3 - 74299.5*eta4)*S2)/(eta);
2683  double uneqSpin = eta*(223.65100000000004*pWF->chi1L*sqrt(1.-4.*eta)*(3.9201300240106223 + 1.*eta) - 223.65100000000004*pWF->chi2L*sqrt(1.-4.*eta)*(3.9201300240106223 + 1.*eta));
2684  total = noSpin + eqSpin + uneqSpin;
2685  break;
2686  }
2687  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_33_p1: version is not valid. Recommended version is 122019.");}
2688  }
2689  return total;
2690 }
2691 
2692 static double IMRPhenomXHM_Inter_Phase_32_p1(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2693  double total=0;
2694  switch (InterPhaseFlag){
2695  case 122019:{
2696  double eta = pWF->eta;
2697  double S = pWF->STotR;
2698  double eta2,eta3,eta4,eta5,eta6,S2,S3,S4;
2699  eta2 = pow(eta,2);
2700  eta3 = pow(eta,3);
2701  eta4 = pow(eta,4);
2702  eta5 = pow(eta,5);
2703  eta6 = pow(eta,6);
2704  S2 = pow(S,2);
2705  S3 = pow(S,3);
2706  S4 = pow(S,4);
2707  double noSpin = 4414.11 + 4.21564/eta - 10687.8*eta + 58234.6*eta2 - 64068.40000000001*eta3 - 704442.*eta4 + 2.86393e6*eta5 - 3.26362e6*eta6;
2708  double eqSpin = ((6.39833 - 610.267*eta + 2095.72*eta2 - 3970.89*eta3)*S + (22.956700000000005 - 99.1551*eta + 331.593*eta2 - 794.79*eta3)*S2 + (10.4333 + 43.8812*eta - 541.261*eta2 + 294.289*eta3)*S3 + eta*(106.047 - 1569.0299999999997*eta + 4810.61*eta2)*S4)/(eta);
2709  double uneqSpin = 132.244*sqrt(1.-4.*eta)*eta*(pWF->chi1L*(6.227738120444028 - 1.*eta) + pWF->chi2L*(-6.227738120444028 + 1.*eta));
2710  total = noSpin + eqSpin + uneqSpin;
2711  break;
2712  }
2713  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_32_p1: version is not valid. Recommended version is 122019.");}
2714  }
2715  return total;
2716 }
2717 
2718 static double IMRPhenomXHM_Inter_Phase_44_p1(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2719  double total=0;
2720  switch (InterPhaseFlag){
2721  case 122019:{
2722  double eta = pWF->eta;
2723  double S = pWF->STotR;
2724  double eta2,eta3,eta4,eta5,eta6,S2,S3;
2725  eta2 = pow(eta,2);
2726  eta3 = pow(eta,3);
2727  eta4 = pow(eta,4);
2728  eta5 = pow(eta,5);
2729  eta6 = pow(eta,6);
2730  S2 = pow(S,2);
2731  S3 = pow(S,3);
2732  double noSpin = 4349.66 + 4.34125/eta - 8202.33*eta + 5534.1*eta2 + 536500.*eta3 - 4.33197e6*eta4 + 1.37792e7*eta5 - 1.60802e7*eta6;
2733  double eqSpin = ((12.0704 - 528.098*eta + 1822.9100000000003*eta2 - 9349.73*eta3 + 17900.9*eta4)*S + (10.4092 + 253.334*eta - 5452.04*eta2 + 35416.6*eta3 - 71523.*eta4)*S2 + eta*(492.60300000000007 - 9508.5*eta + 57303.4*eta2 - 109418.*eta3)*S3)/(eta);
2734  double uneqSpin = -262.143*sqrt(1.-4.*eta)*eta*(pWF->chi1L*(-3.0782778864970646 - 1.*eta) + pWF->chi2L*(3.0782778864970646 + 1.*eta));
2735  total = noSpin + eqSpin + uneqSpin;
2736  break;
2737  }
2738  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_44_p1: version is not valid. Recommended version is 122019.");}
2739  }
2740  return total;
2741 }
2742 
2743 
2744 static double IMRPhenomXHM_Inter_Phase_21_p2(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2745  double total=0;
2746  switch (InterPhaseFlag){
2747  case 122019:{
2748  double eta = pWF->eta;
2749  double S = pWF->STotR;
2750  double eta2,eta3,eta4,eta5,eta6,S2,S3,S4;
2751  eta2 = pow(eta,2);
2752  eta3 = eta2*eta;
2753  eta4 = eta3*eta;
2754  eta5 = eta3*eta2;
2755  eta6 = eta4*eta2;
2756  S2 = pow(S,2);
2757  S3 = pow(S,3);
2758  S4= S3*S;
2759  double noSpin = 3509.09 + 0.91868/eta + 194.72*eta - 27556.2*eta2 + 369153.*eta3 - 2.28832e6*eta4 + 6.82533e6*eta5 - 7.86254e6*eta6;
2760  double eqSpin = ((0.7083999999999999 - 60.1611*eta + 131.815*eta2 - 619.837*eta3)*S + (6.104720000000001 - 59.2068*eta + 278.588*eta2 - 722.787*eta3)*S2 + (5.7791 + 117.913*eta - 1180.4*eta2 + 2392.73*eta3)*S3 + eta*(92.6654 - 1077.75*eta + 2689.18*eta2)*S4)/(eta);
2761  double uneqSpin = -91.7679*sqrt(1.-4.*eta)*eta*(pWF->chi1L*(-1.6012352903357276 - 1.*eta) + pWF->chi2L*(1.6012352903357276 + 1.*eta));
2762  total = noSpin + eqSpin + uneqSpin;
2763  break;
2764  }
2765  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_21_p2: version is not valid.Recommended version is 122019.");}
2766  }
2767  return total;
2768 }
2769 
2770 static double IMRPhenomXHM_Inter_Phase_33_p2(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2771  double total=0;
2772  switch (InterPhaseFlag){
2773  case 122019:{
2774  double eta = pWF->eta;
2775  double S = pWF->STotR;
2776  double eta2,eta3,eta4,eta5,eta6,S2;
2777  eta2 = pow(eta,2);
2778  eta3 = pow(eta,3);
2779  eta4 = pow(eta,4);
2780  eta5 = pow(eta,5);
2781  eta6 = pow(eta,6);
2782  S2 = pow(S,2);
2783  double noSpin = 3797.06 + 0.786684/eta - 2397.09*eta - 25514.*eta2 + 518314.99999999994*eta3 - 3.41708e6*eta4 + 1.01799e7*eta5 - 1.15659e7*eta6;
2784  double eqSpin = ((6.7812399999999995 + 39.4668*eta - 3520.37*eta2 + 19614.6*eta3 - 34860.2*eta4)*S + (4.80384 + 293.215*eta - 5914.61*eta2 + 36587.3*eta3 - 74299.5*eta4)*S2)/(eta);
2785  double uneqSpin = -223.65100000000004*sqrt(1.-4.*eta)*eta*(pWF->chi1L*(-1.3095134830606614 - 1.*eta) + pWF->chi2L*(1.3095134830606614 + 1.*eta));
2786  total = noSpin + eqSpin + uneqSpin;
2787  break;
2788  }
2789  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_33_p2: version is not valid.Recommended version is 122019.");}
2790  }
2791  return total;
2792 }
2793 
2794 static double IMRPhenomXHM_Inter_Phase_32_p2(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2795  double total=0;
2796  switch (InterPhaseFlag){
2797  case 122019:{
2798  double eta = pWF->eta;
2799  double S = pWF->STotR;
2800  double eta2,eta3,eta4,eta5,eta6,S2,S3,S4;
2801  eta2 = pow(eta,2);
2802  eta3 = pow(eta,3);
2803  eta4 = pow(eta,4);
2804  eta5 = pow(eta,5);
2805  eta6 = pow(eta,6);
2806  S2 = pow(S,2);
2807  S3 = pow(S,3);
2808  S4 = pow(S,4);
2809  double noSpin = 3980.7 + 0.956703/eta - 6202.38*eta + 29218.1*eta2 + 24484.2*eta3 - 807629.*eta4 + 2.86393e6*eta5 - 3.26362e6*eta6;
2810  double eqSpin = ((1.92692 - 226.825*eta + 75.246*eta2 + 1291.56*eta3)*S + (15.328700000000001 - 99.1551*eta + 608.328*eta2 - 2402.94*eta3)*S2 + (10.4333 + 43.8812*eta - 541.261*eta2 + 294.289*eta3)*S3 + eta*(106.047 - 1569.0299999999997*eta + 4810.61*eta2)*S4)/(eta);
2811  double uneqSpin = 132.244*sqrt(1.-4.*eta)*eta*(pWF->chi1L*(2.5769789177580837 - 1.*eta) + pWF->chi2L*(-2.5769789177580837 + 1.*eta));
2812  total = noSpin + eqSpin + uneqSpin;
2813  break;
2814  }
2815  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_32_p2: version is not valid. Recommended version is 122019.");}
2816  }
2817  return total;
2818 }
2819 
2820 static double IMRPhenomXHM_Inter_Phase_44_p2(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2821  double total=0;
2822  switch (InterPhaseFlag){
2823  case 122019:{
2824  double eta = pWF->eta;
2825  double S = pWF->STotR;
2826  double eta2,eta3,eta4,eta5,eta6,S2,S3;
2827  eta2 = pow(eta,2);
2828  eta3 = pow(eta,3);
2829  eta4 = pow(eta,4);
2830  eta5 = pow(eta,5);
2831  eta6 = pow(eta,6);
2832  S2 = pow(S,2);
2833  S3 = pow(S,3);
2834  double noSpin = 3804.19 + 0.66144/eta - 2421.77*eta - 33475.8*eta2 + 665951.*eta3 - 4.50145e6*eta4 + 1.37792e7*eta5 - 1.60802e7*eta6;
2835  double eqSpin = ((5.83038 - 172.047*eta + 926.576*eta2 - 7676.87*eta3 + 17900.9*eta4)*S + (6.17601 + 253.334*eta - 5672.02*eta2 + 35722.1*eta3 - 71523.*eta4)*S2 + eta*(492.60300000000007 - 9508.5*eta + 57303.4*eta2 - 109418.*eta3)*S3)/(eta);
2836  double uneqSpin = -262.143*sqrt(1.-4.*eta)*eta*(pWF->chi1L*(-1.0543062374352932 - 1.*eta) + pWF->chi2L*(1.0543062374352932 + 1.*eta));
2837  total = noSpin + eqSpin + uneqSpin;
2838  break;
2839  }
2840  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_44_p2: version is not valid. Recommended version is 122019.");}
2841  }
2842  return total;
2843 }
2844 
2845 
2846 static double IMRPhenomXHM_Inter_Phase_21_p3(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2847  double total=0;
2848  switch (InterPhaseFlag){
2849  case 122019:{
2850  double eta = pWF->eta;
2851  double S = pWF->STotR;
2852  double delta=sqrt(1.-4.*eta);
2853  double eta2,eta3,eta4,eta5,eta6,S2,S3,S4;
2854  eta2 = pow(eta,2);
2855  eta3 = eta2*eta;
2856  eta4 = eta2*eta2;
2857  eta5 = eta2*eta3;
2858  eta6 = eta3*eta3;
2859  S2 = pow(S,2);
2860  S3 = pow(S,3);
2861  S4 = pow(S,4);
2862  double noSpin = 3241.68 + 890.016*eta - 28651.9*eta2 + 369153.*eta3 - 2.28832e6*eta4 + 6.82533e6*eta5 - 7.86254e6*eta6;
2863  double eqSpin = (-2.2484 + 187.641*eta - 619.837*eta2)*S + (3.22603 + 166.323*eta - 722.787*eta2)*S2 + (117.913 - 1094.59*eta + 2392.73*eta2)*S3 + (92.6654 - 1077.75*eta + 2689.18*eta2)*S4;
2864  double uneqSpin = 91.7679*(pWF->dchi)*delta*eta2;
2865  total = noSpin + eqSpin + uneqSpin;
2866  break;
2867  }
2868  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_21_p3: version is not valid. Recommended version is 122019.");}
2869  }
2870  return total;
2871 }
2872 
2873 static double IMRPhenomXHM_Inter_Phase_33_p3(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2874  double total=0;
2875  switch (InterPhaseFlag){
2876  case 122019:{
2877  double eta = pWF->eta;
2878  double S = pWF->STotR;
2879  double eta2,eta3,eta4,eta5,eta6,S2;
2880  eta2 = pow(eta,2);
2881  eta3 = pow(eta,3);
2882  eta4 = pow(eta,4);
2883  eta5 = pow(eta,5);
2884  eta6 = pow(eta,6);
2885  S2 = pow(S,2);
2886  double noSpin = 3321.83 + 1796.03*eta - 52406.1*eta2 + 605028.*eta3 - 3.52532e6*eta4 + 1.01799e7*eta5 - 1.15659e7*eta6;
2887  double eqSpin = (223.601 - 3714.77*eta + 19614.6*eta2 - 34860.2*eta3)*S + (314.317 - 5906.46*eta + 36587.3*eta2 - 74299.5*eta3)*S2;
2888  double uneqSpin = 223.651*(pWF->dchi)*sqrt(1.-4.*eta)*eta2;
2889  total = noSpin + eqSpin + uneqSpin;
2890  break;
2891  }
2892  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_33_p3: version is not valid. Recommended version is 122019.");}
2893  }
2894  return total;
2895 }
2896 
2897 static double IMRPhenomXHM_Inter_Phase_32_p3(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2898  double total=0;
2899  switch (InterPhaseFlag){
2900  case 122019:{
2901  double eta = pWF->eta;
2902  double S = pWF->STotR;
2903  double eta2,eta3,eta4,eta5,eta6,S2,S3,S4;
2904  eta2 = pow(eta,2);
2905  eta3 = pow(eta,3);
2906  eta4 = pow(eta,4);
2907  eta5 = pow(eta,5);
2908  eta6 = pow(eta,6);
2909  S2 = pow(S,2);
2910  S3 = pow(S,3);
2911  S4 = pow(S,4);
2912  double noSpin = 3416.57 + 2308.63*eta - 84042.9*eta2 + 1.01936e6*eta3 - 6.0644e6*eta4 + 1.76399e7*eta5 - 2.0065e7*eta6;
2913  double eqSpin = (24.6295 - 282.354*eta - 2582.55*eta2 + 12750.*eta3)*S + (433.675 - 8775.86*eta + 56407.8*eta2 - 114798.*eta3)*S2 + (559.705 - 10627.4*eta + 61581.*eta2 - 114029.*eta3)*S3 + (106.047 - 1569.03*eta + 4810.61*eta2)*S4;
2914  double uneqSpin = 63.9466*pWF->dchi*sqrt(1.-4.*eta)*eta2;
2915  total = noSpin + eqSpin + uneqSpin;
2916  break;
2917  }
2918  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_32_p3: version is not valid. Recommended version is 122019.");}
2919  }
2920  return total;
2921 }
2922 
2923 static double IMRPhenomXHM_Inter_Phase_44_p3(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2924  double total=0;
2925  switch (InterPhaseFlag){
2926  case 122019:{
2927  double eta = pWF->eta;
2928  double S = pWF->STotR;
2929  double eta2,eta3,eta4, eta5, eta6, S2, S3;
2930  eta2 = pow(eta,2);
2931  eta3 = pow(eta,3);
2932  eta4 = pow(eta,4);
2933  eta5 = pow(eta,5);
2934  eta6 = pow(eta,6);
2935  S2 = pow(S,2);
2936  S3 = pow(S,3);
2937  double noSpin = 3308.97 + 2353.58*eta - 66340.1*eta2 + 777272.*eta3 - 4.64438e6*eta4 + 1.37792e7*eta5 - 1.60802e7*eta6;
2938  double eqSpin = (-21.5697 + 926.576*eta - 7989.26*eta2 + 17900.9*eta3)*S + (353.539 - 6403.24*eta + 37599.5*eta2 - 71523.*eta3)*S2 + (492.603 - 9508.5*eta + 57303.4*eta2 - 109418.*eta3)*S3;
2939  double uneqSpin = 262.143*(pWF->dchi)*sqrt(1.-4.*eta)*eta2;
2940  total = noSpin + eqSpin + uneqSpin;
2941  break;
2942  }
2943  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_44_p3: version is not valid.Recommended version is 122019.");}
2944  }
2945  return total;
2946 }
2947 
2948 
2949 static double IMRPhenomXHM_Inter_Phase_21_p4(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2950  double total=0;
2951  switch (InterPhaseFlag){
2952  case 122019:{
2953  double eta = pWF->eta;
2954  double S = pWF->STotR;
2955  double eta2,eta3,eta4,eta5,eta6,S2,S3,S4;
2956  eta2 = pow(eta,2);
2957  eta3 = pow(eta,3);
2958  eta4 = pow(eta,4);
2959  eta5 = pow(eta,5);
2960  eta6 = pow(eta,6);
2961  S2 = pow(S,2);
2962  S3 = pow(S,3);
2963  S4 = pow(S,4);
2964  double noSpin = 3160.88 + 974.355*eta - 28932.5*eta2 + 369780.*eta3 - 2.28832e6*eta4 + 6.82533e6*eta5 - 7.86254e6*eta6;
2965  double eqSpin = (26.3355 - 196.851*eta + 438.401*eta2)*S + (45.9957 - 256.248*eta + 117.563*eta2)*S2 + (-20.0261 + 467.057*eta - 1613.*eta2)*S3 + (-61.7446 + 577.057*eta - 1096.81*eta2)*S4;
2966  double uneqSpin = 65.3326*pWF->dchi*sqrt(1.-4.*eta)*eta2;
2967  total = noSpin + eqSpin + uneqSpin;
2968  break;
2969  }
2970  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_21_p4: version is not valid.Recommended version is 122019.");}
2971  }
2972  return total;
2973 }
2974 
2975 static double IMRPhenomXHM_Inter_Phase_33_p4(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
2976  double total=0;
2977  switch (InterPhaseFlag){
2978  case 122019:{
2979  double eta = pWF->eta;
2980  double S = pWF->STotR;
2981  double eta2,eta3,eta4,eta5, eta6,S2,S3;
2982  eta2 = pow(eta,2);
2983  eta3 = pow(eta,3);
2984  eta4 = pow(eta,4);
2985  eta5 = pow(eta,5);
2986  eta6 = pow(eta,6);
2987  S2 = pow(S,2);
2988  S3 = pow(S,3);
2989  double noSpin = 3239.44 - 661.15*eta + 5139.79*eta2 + 3456.2*eta3 - 248477.*eta4 + 1.17255e6*eta5 - 1.70363e6*eta6;
2990  double eqSpin = (225.859 - 4150.09*eta + 24364.*eta2 - 46537.3*eta3)*S + (35.2439 - 994.971*eta + 8953.98*eta2 - 23603.5*eta3)*S2 + (-310.489 + 5946.15*eta - 35337.1*eta2 + 67102.4*eta3)*S3;
2991  double uneqSpin = 30.484*pWF->dchi*sqrt(1.-4.*eta)*eta2;
2992  total = noSpin + eqSpin + uneqSpin;
2993  break;
2994  }
2995  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_33_p4: version is not valid. Recommended version is 122019.");}
2996  }
2997  return total;
2998 }
2999 
3000 static double IMRPhenomXHM_Inter_Phase_32_p4(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
3001  double total=0;
3002  switch (InterPhaseFlag){
3003  case 122019:{
3004  double eta = pWF->eta;
3005  double S = pWF->STotR;
3006  double eta2,eta3,eta4,eta5,eta6,S2,S3,S4;
3007  eta2 = pow(eta,2);
3008  eta3 = pow(eta,3);
3009  eta4 = pow(eta,4);
3010  eta5 = pow(eta,5);
3011  eta6 = pow(eta,6);
3012  S2 = pow(S,2);
3013  S3 = pow(S,3);
3014  S4 = pow(S,4);
3015  double noSpin = 3307.49 - 476.909*eta - 5980.37*eta2 + 127610.*eta3 - 919108.*eta4 + 2.86393e6*eta5 - 3.26362e6*eta6;
3016  double eqSpin = (-5.02553 - 282.354*eta + 1291.56*eta2)*S + (-43.8823 + 740.123*eta - 2402.94*eta2)*S2 + (43.8812 - 370.362*eta + 294.289*eta2)*S3 + (106.047 - 1569.03*eta + 4810.61*eta2)*S4;
3017  double uneqSpin = -132.244*(pWF->dchi)*sqrt(1.-4.*eta)*eta2;
3018  total = noSpin + eqSpin + uneqSpin;
3019  break;
3020  }
3021  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_32_p4: version is not valid. Recommended version is 122019.");}
3022  }
3023  return total;
3024 }
3025 
3026 static double IMRPhenomXHM_Inter_Phase_44_p4(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
3027  double total=0;
3028  switch (InterPhaseFlag){
3029  case 122019:{
3030  double eta = pWF->eta;
3031  double S = pWF->STotR;
3032  double eta2,eta3,eta4,eta5,eta6,S2,S3;
3033  eta2 = pow(eta,2);
3034  eta3 = pow(eta,3);
3035  eta4 = pow(eta,4);
3036  eta5 = pow(eta,5);
3037  eta6 = pow(eta,6);
3038  S2 = pow(S,2);
3039  S3 = pow(S,3);
3040  double noSpin = 3245.63 - 928.56*eta + 8463.89*eta2 - 17422.6*eta3 - 165169.*eta4 + 908279.*eta5 - 1.31138e6*eta6;
3041  double eqSpin = (32.506 - 590.293*eta + 3536.61*eta2 - 6758.52*eta3)*S + (-25.7716 + 738.141*eta - 4867.87*eta2 + 9129.45*eta3)*S2 + (-15.7439 + 620.695*eta - 4679.24*eta2 + 9582.58*eta3)*S3;
3042  double uneqSpin = 87.0832*pWF->dchi*sqrt(1.-4.*eta)*eta2;
3043  total = noSpin + eqSpin + uneqSpin;
3044  break;
3045  }
3046  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_44_p4: version is not valid. Recommended version is 122019.");}
3047  }
3048  return total;
3049 }
3050 
3051 
3052 static double IMRPhenomXHM_Inter_Phase_21_p5(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
3053  double total=0;
3054  switch (InterPhaseFlag){
3055  case 122019:{
3056  double eta = pWF->eta;
3057  double S = pWF->STotR;
3058  double eta2,eta3,S2,S3,S4;
3059  eta2 = pow(eta,2);
3060  eta3 = pow(eta,3);
3061  S2 = pow(S,2);
3062  S3 = pow(S,3);
3063  S4 = pow(S,4);
3064  double noSpin = 3102.36 + 315.911*eta - 1688.26*eta2 + 3635.76*eta3;
3065  double eqSpin = (-23.0959 + 320.93*eta - 1029.76*eta2)*S + (-49.5435 + 826.816*eta - 3079.39*eta2)*S2 + (40.7054 - 365.842*eta + 1094.11*eta2)*S3 + (81.8379 - 1243.26*eta + 4689.22*eta2)*S4;
3066  double uneqSpin = 119.014*(pWF->dchi)*sqrt(1.-4.*eta)*eta2;
3067  total = noSpin + eqSpin + uneqSpin;
3068  break;
3069  }
3070  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_21_p5: version is not valid. Recommended version is 122019.");}
3071  }
3072  return total;
3073 }
3074 
3075 static double IMRPhenomXHM_Inter_Phase_33_p5(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
3076  double total=0;
3077  switch (InterPhaseFlag){
3078  case 122019:{
3079  double eta = pWF->eta;
3080  double S = pWF->STotR;
3081  double eta2,eta3,eta4,eta5,eta6,S2,S3;
3082  eta2 = pow(eta,2);
3083  eta3 = pow(eta,3);
3084  eta4 = pow(eta,4);
3085  eta5 = pow(eta,5);
3086  eta6 = pow(eta,6);
3087  S2 = pow(S,2);
3088  S3 = pow(S,3);
3089  double noSpin = 3114.3 + 2143.06*eta - 49428.3*eta2 + 563997.*eta3 - 3.35991e6*eta4 + 9.99745e6*eta5 - 1.17123e7*eta6;
3090  double eqSpin = (190.051 - 3705.08*eta + 23046.2*eta2 - 46537.3*eta3)*S + (63.6615 - 1414.2*eta + 10166.1*eta2 - 23603.5*eta3)*S2 + (-257.524 + 5179.97*eta - 33001.4*eta2 + 67102.4*eta3)*S3;
3091  double uneqSpin = 54.9833*pWF->dchi*sqrt(1.-4.*eta)*eta2;
3092  total = noSpin + eqSpin + uneqSpin;
3093  break;
3094  }
3095  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_33_p5: version is not valid. Recommended version is 122019.");}
3096  }
3097  return total;
3098 }
3099 
3100 static double IMRPhenomXHM_Inter_Phase_32_p5(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
3101  double total=0;
3102  switch (InterPhaseFlag){
3103  case 122019:{
3104  double eta = pWF->eta;
3105  double S = pWF->STotR;
3106  double eta2,eta3,eta4,eta5,eta6,eta7,S2,S3,S4;
3107  eta2 = pow(eta,2);
3108  eta3 = pow(eta,3);
3109  eta4 = pow(eta,4);
3110  eta5 = pow(eta,5);
3111  eta6 = pow(eta,6);
3112  eta7 = pow(eta,7);
3113  S2 = pow(S,2);
3114  S3 = pow(S,3);
3115  S4 = pow(S,4);
3116  double noSpin = 3259.03 - 3967.58*eta + 111203.*eta2 - 1.81883e6*eta3 + 1.73811e7*eta4 - 9.56988e7*eta5 + 2.75056e8*eta6 - 3.15866e8*eta7;
3117  double eqSpin = (19.7509 - 1104.53*eta + 3810.18*eta2)*S + (-230.07 + 2314.51*eta - 5944.49*eta2)*S2 + (-201.633 + 2183.43*eta - 6233.99*eta2)*S3 + (106.047 - 1569.03*eta + 4810.61*eta2)*S4;
3118  double uneqSpin = 112.714*pWF->dchi*sqrt(1.-4.*eta)*eta2;
3119  total = noSpin + eqSpin + uneqSpin;
3120  break;
3121  }
3122  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_32_p5: version is not valid. Recommended version is 122019.");}
3123  }
3124  return total;
3125 }
3126 
3127 static double IMRPhenomXHM_Inter_Phase_44_p5(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
3128  double total=0;
3129  switch (InterPhaseFlag){
3130  case 122019:{
3131  double eta = pWF->eta;
3132  double S = pWF->STotR;
3133  double eta2,eta3,eta4,eta5,eta6,eta7,S2,S3;
3134  eta2 = pow(eta,2);
3135  eta3 = pow(eta,3);
3136  eta4 = pow(eta,4);
3137  eta5 = pow(eta,5);
3138  eta6 = pow(eta,6);
3139  eta7 = pow(eta,7);
3140  S2 = pow(S,2);
3141  S3 = pow(S,3);
3142  double noSpin = 3108.38 + 3722.46*eta - 119588.*eta2 + 1.92148e6*eta3 - 1.69796e7*eta4 + 8.39194e7*eta5 - 2.17143e8*eta6 + 2.2829700000000003e8*eta7;
3143  double eqSpin = (118.319 - 529.854*eta)*eta*S + (21.0314 - 240.648*eta + 516.333*eta2)*S2 + (20.3384 - 356.241*eta + 999.417*eta2)*S3;
3144  double uneqSpin = 97.1364*(pWF->dchi)*sqrt(1.-4.*eta)*eta2;
3145  total = noSpin + eqSpin + uneqSpin;
3146  break;
3147  }
3148  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_44_p5: version is not valid. Recommended version is 122019.");}
3149  }
3150  return total;
3151 }
3152 
3153 
3154 static double IMRPhenomXHM_Inter_Phase_21_p6(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
3155  double total=0;
3156  switch (InterPhaseFlag){
3157  case 122019:{
3158  double eta = pWF->eta;
3159  double S = pWF->STotR;
3160  double eta2,eta3,S2,S3,S4;
3161  eta2 = pow(eta,2);
3162  eta3 = pow(eta,3);
3163  S2 = pow(S,2);
3164  S3 = pow(S,3);
3165  S4 = pow(S,4);
3166  double noSpin = 3089.18 + 4.89194*eta + 190.008*eta2 - 255.245*eta3;
3167  double eqSpin = (2.96997 + 57.1612*eta - 432.223*eta2)*S + (-18.8929 + 630.516*eta - 2804.66*eta2)*S2 + (-24.6193 + 549.085*eta2)*S3 + (-12.8798 - 722.674*eta + 3967.43*eta2)*S4;
3168  double uneqSpin = 74.0984*(pWF->dchi)*sqrt(1.-4.*eta)*eta2;
3169  total = noSpin + eqSpin + uneqSpin;
3170  break;
3171  }
3172  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_21_p6: version is not valid. Recommended version is 122019.");}
3173  }
3174  return total;
3175 }
3176 
3177 static double IMRPhenomXHM_Inter_Phase_33_p6(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
3178  double total=0;
3179  switch (InterPhaseFlag){
3180  case 122019:{
3181  double eta = pWF->eta;
3182  double S = pWF->STotR;
3183  double eta2,eta3, eta4, eta5, eta6,S2,S3;
3184  eta2 = pow(eta,2);
3185  eta3 = pow(eta,3);
3186  eta4 = pow(eta,4);
3187  eta5 = pow(eta,5);
3188  eta6 = pow(eta,6);
3189  S2 = pow(S,2);
3190  S3 = pow(S,3);
3191  double noSpin = 3111.46 + 384.121*eta - 13003.6*eta2 + 179537.*eta3 - 1.19313e6*eta4 + 3.79886e6*eta5 - 4.64858e6*eta6;
3192  double eqSpin = (182.864 - 3834.22*eta + 24532.9*eta2 - 50165.9*eta3)*S + (21.0158 - 746.957*eta + 6701.33*eta2 - 17842.3*eta3)*S2 + (-292.855 + 5886.62*eta - 37382.4*eta2 + 75501.8*eta3)*S3;
3193  double uneqSpin = 75.5162*pWF->dchi*sqrt(1.-4.*eta)*eta2;
3194  total = noSpin + eqSpin + uneqSpin;
3195  break;
3196  }
3197  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_33_p6: version is not valid. Recommended version is 122019.");}
3198  }
3199  return total;
3200 }
3201 
3202 static double IMRPhenomXHM_Inter_Phase_32_p6(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
3203  double total=0;
3204  switch (InterPhaseFlag){
3205  case 122019:{
3206  double eta = pWF->eta;
3207  double S = pWF->STotR;
3208  double eta2,eta3,eta4,eta5,eta6,eta7,S2,S3,S4;
3209  eta2 = pow(eta,2);
3210  eta3 = pow(eta,3);
3211  eta4 = pow(eta,4);
3212  eta5 = pow(eta,5);
3213  eta6 = pow(eta,6);
3214  eta7 = pow(eta,7);
3215  S2 = pow(S,2);
3216  S3 = pow(S,3);
3217  S4 = pow(S,4);
3218  double noSpin = 3259.03 - 3967.58*eta + 111203.*eta2 - 1.81883e6*eta3 + 1.73811e7*eta4 - 9.56988e7*eta5 + 2.75056e8*eta6 - 3.15866e8*eta7;
3219  double eqSpin = (19.7509 - 1104.53*eta + 3810.18*eta2)*S + (-230.07 + 2314.51*eta - 5944.49*eta2)*S2 + (-201.633 + 2183.43*eta - 6233.99*eta2)*S3 + (106.047 - 1569.03*eta + 4810.61*eta2)*S4;
3220  double uneqSpin = 112.714*pWF->dchi*sqrt(1.-4.*eta)*eta2;
3221  total = noSpin + eqSpin + uneqSpin;
3222  break;
3223  }
3224  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_32_p6: version is not valid. Recommended version is 122019.");}
3225  }
3226  return total;
3227 }
3228 
3229 static double IMRPhenomXHM_Inter_Phase_44_p6(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag) {
3230  double total=0.;
3231  switch (InterPhaseFlag){
3232  case 122019:{
3233  double eta = pWF->eta;
3234  double S = pWF->STotR;
3235  double eta2,eta3,eta4,eta5,eta6,S2,S3;
3236  eta2 = pow(eta,2);
3237  eta3 = pow(eta,3);
3238  eta4 = pow(eta,4);
3239  eta5 = pow(eta,5);
3240  eta6 = pow(eta,6);
3241  S2 = pow(S,2);
3242  S3 = pow(S,3);
3243  double noSpin = 3096.03 + 986.752*eta - 20371.1*eta2 + 220332.*eta3 - 1.31523e6*eta4 + 4.29193e6*eta5 - 6.01179e6*eta6;
3244  double eqSpin = (-9.96292 - 118.526*eta + 2255.76*eta2 - 6758.52*eta3)*S + (-14.4869 + 370.039*eta - 3605.8*eta2 + 9129.45*eta3)*S2 + (17.0209 + 70.1931*eta - 3070.08*eta2 + 9582.58*eta3)*S3;
3245  double uneqSpin = 23.0759*pWF->dchi*sqrt(1.-4.*eta)*eta2;
3246  total = noSpin + eqSpin + uneqSpin;
3247  break;
3248  }
3249  default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Inter_Phase_44_p6: version is not valid. Recommended version is 122019.");}
3250  }
3251  return total;
3252 }
3253 
3254 /* End of Phase Parameter Space Fits */
3255 
3256 /************* PHASE ANSATZ ****************/
3257 
3258 static double IMRPhenomXHM_Inter_Phase_AnsatzInt(double ff, IMRPhenomX_UsefulPowers *powers_of_f,IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMPhaseCoefficients *pPhase){
3259 
3260  int modeTag = pWFHM->modeTag;
3261  double invf = powers_of_f->m_one;
3262  double invf2 = powers_of_f->m_two;
3263  double invf3 = powers_of_f->m_three;
3264  double logfv=powers_of_f->log;
3265  double phaseIR;
3266  double fda=pWFHM->fDAMP;
3267  double frd=pWFHM->fRING;
3268 
3269  if(modeTag!=32)
3270  /* 5 coefficients */
3271 
3272  phaseIR = pPhase->c0 *ff + (pPhase->c1)*logfv - (pPhase->c2)*invf -1./3.* (pPhase->c4)*invf3 + (pPhase->cL)* atan((ff - frd)/fda);
3273  else { /* 6 coefficients */
3274 
3275  invf3 = powers_of_f->m_three;
3276  /*(c0 + c1 /f + c2 /(f)^2 + c4 /(f)^4 + c3 /f^3 +
3277  cL fdamp/((fdamp)^2 + (f - fRD )^2))*/
3278  phaseIR = pPhase->c0 *ff + (pPhase->c1)*logfv - (pPhase->c2)*invf -1./3.* (pPhase->c4)*invf3 -0.5*pPhase->c3*invf2+ (pPhase->cL)* atan((ff - frd)/fda);
3279  }
3280 
3281  return phaseIR;
3282 
3283 }
3284 
3285 static double IMRPhenomXHM_Inter_Phase_Ansatz(double ff, IMRPhenomX_UsefulPowers *powers_of_f,IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMPhaseCoefficients *pPhase){
3286 
3287  int modeTag = pWFHM->modeTag;
3288  double invf = powers_of_f->m_one;
3289  double invf2 = powers_of_f->m_two;
3290  double invf4 = powers_of_f->m_four;
3291  double dphaseIR;
3292  double fda=pWFHM->fDAMP;
3293  double frd=pWFHM->fRING;
3294 
3295  if(modeTag!=32)
3296  /* 5 coefficients */
3297 
3298  /*(c0 + c1 /f + c2 /(f)^2 + c4 /(f)^4 +
3299  cL fdamp/((fdamp)^2 + (f - fRD )^2))*/
3300  dphaseIR = ( pPhase->c0 + (pPhase->c1)*invf + (pPhase->c2)*invf2 + (pPhase->c4)*invf4 + ( (pPhase->cL)* fda/(fda*fda +(ff - frd)*(ff - frd)) ) );
3301  else { /* 6 coefficients */
3302 
3303  double invf3 = powers_of_f->m_three;
3304  /*(c0 + c1 /f + c2 /(f)^2 + c4 /(f)^4 + c3 /f^3 +
3305  cL fdamp/((fdamp)^2 + (f - fRD )^2))*/
3306  dphaseIR = ( pPhase->c0 + (pPhase->c1)*invf + (pPhase->c2)*invf2 + (pPhase->c4)*invf4 +
3307  (pPhase->c3)*invf3 + ( (pPhase->cL)* fda/(fda*fda +(ff - frd)*(ff - frd)) ) );
3308  }
3309 
3310  return dphaseIR;
3311 
3312 }
3313 
3314 static double IMRPhenomXHM_Inter_Phase_dAnsatz(double ff, IMRPhenomX_UsefulPowers *powers_of_f,IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMPhaseCoefficients *pPhase){
3315 
3316  int modeTag = pWFHM->modeTag;
3317  double invf2 = powers_of_f->m_two;
3318  double invf3 = powers_of_f->m_three;
3319  double invf5 = powers_of_f->m_five;
3320  double d2phaseIR;
3321  double fda=pWFHM->fDAMP;
3322  double frd=pWFHM->fRING;
3323 
3324  if(modeTag!=32)
3325  /* 5 coefficients */
3326 
3327  /*(-((4 c4)/f^5) - (2 c2)/f^3 - c1/f^2 - (
3328  2 cL fdamp (f - fRD))/(fdamp^2 + (f - fRD)^2)^2)*/
3329  d2phaseIR = -(pPhase->c1)*invf2 -2. *(pPhase->c2)*invf3 -4. *(pPhase->c4)*invf5 -2.* (pPhase->cL)* (ff - frd)*fda/pow((fda*fda +(ff - frd)*(ff - frd)),2);
3330  else { /* 6 coefficients */
3331 
3332  double invf4 = powers_of_f->m_four;
3333  /*(-((4 c4)/f^5) - (2 c2)/f^3 - c1/f^2 -3 c3/ f^4- (
3334  2 cL fdamp (f - fRD))/(fdamp^2 + (f - fRD)^2)^2)*/
3335  d2phaseIR =-(pPhase->c1)*invf2 -3.*pPhase->c3*invf4-2. *(pPhase->c2)*invf3 -4. *(pPhase->c4)*invf5 -2.* (pPhase->cL)* (ff - frd)*fda/pow((fda*fda +(ff - frd)*(ff - frd)),2);
3336  }
3337 
3338  return d2phaseIR;
3339 
3340 }
delta
static double double delta
Definition: LALSimIMRPhenomTHM_fits.h:105
IMRPhenomX_Initialize_Powers
int IMRPhenomX_Initialize_Powers(IMRPhenomX_UsefulPowers *p, REAL8 number)
Definition: LALSimIMRPhenomX_internals.c:54
IMRPhenomXHM_Inspiral_Amp_NDAnsatz
static double IMRPhenomXHM_Inspiral_Amp_NDAnsatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
Definition: LALSimIMRPhenomXHM_inspiral.c:777
IMRPhenomXHM_Inspiral_Amp_Ansatz
static double IMRPhenomXHM_Inspiral_Amp_Ansatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
Definition: LALSimIMRPhenomXHM_inspiral.c:742
IMRPhenomXHM_Intermediate_Amp_delta0
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)
Definition: LALSimIMRPhenomXHM_intermediate.c:728
IMRPhenomXHM_Inter_Phase_21_p3
static double IMRPhenomXHM_Inter_Phase_21_p3(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2846
IMRPhenomXHM_Inter_Amp_21_int4
static double IMRPhenomXHM_Inter_Amp_21_int4(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:538
InsideInterval
int InsideInterval(double ftest, double fmin, double fmax)
Definition: LALSimIMRPhenomXHM_intermediate.c:2013
IMRPhenomXHM_Intermediate_Amp_delta5
static double IMRPhenomXHM_Intermediate_Amp_delta5(double d1, double d4, double v1, double v2, double v3, double v4, double f1, double f2, double f3, double f4, int IntAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:1798
CrossZeroP2
int CrossZeroP2(double a0, double a1, double a2, double fstart, double fend)
Definition: LALSimIMRPhenomXHM_intermediate.c:2023
ChoosePolOrder
void ChoosePolOrder(IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
Definition: LALSimIMRPhenomXHM_intermediate.c:2233
IMRPhenomXHM_Inter_Amp_44_int2
static double IMRPhenomXHM_Inter_Amp_44_int2(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:310
IMRPhenomXHM_Inter_Amp_33_int1
static double IMRPhenomXHM_Inter_Amp_33_int1(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:120
IMRPhenomXHM_Inter_Amp_33_int2
static double IMRPhenomXHM_Inter_Amp_33_int2(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:158
IMRPhenomXHM_Inter_Phase_dAnsatz
static double IMRPhenomXHM_Inter_Phase_dAnsatz(double ff, IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMPhaseCoefficients *pPhase)
Definition: LALSimIMRPhenomXHM_intermediate.c:3314
IMRPhenomXHM_Inter_Amp_32_int2
static double IMRPhenomXHM_Inter_Amp_32_int2(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:234
IMRPhenomXHM_Inter_Amp_21_dint0
static double IMRPhenomXHM_Inter_Amp_21_dint0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:371
IMRPhenomXHM_Inter_Phase_44_p5
static double IMRPhenomXHM_Inter_Phase_44_p5(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:3127
IMRPhenomXHM_Intermediate_Amp_CollocationPoints
static void IMRPhenomXHM_Intermediate_Amp_CollocationPoints(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22)
Definition: LALSimIMRPhenomXHM_intermediate.c:2436
IMRPhenomXHM_Inter_Amp_32_int0
static double IMRPhenomXHM_Inter_Amp_32_int0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:431
IMRPhenomXHM_Inter_Phase_44_p6
static double IMRPhenomXHM_Inter_Phase_44_p6(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:3229
IMRPhenomXHM_Inter_Phase_32_p6
static double IMRPhenomXHM_Inter_Phase_32_p6(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:3202
IMRPhenomXHM_Inter_Phase_21_p6
static double IMRPhenomXHM_Inter_Phase_21_p6(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:3154
IMRPhenomXHM_Inter_Amp_33_int3
static double IMRPhenomXHM_Inter_Amp_33_int3(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:562
IMRPhenomXHM_Inter_Amp_44_int4
static double IMRPhenomXHM_Inter_Amp_44_int4(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:635
IMRPhenomXHM_Inter_Phase_Ansatz
static double IMRPhenomXHM_Inter_Phase_Ansatz(double ff, IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMPhaseCoefficients *pPhase)
Definition: LALSimIMRPhenomXHM_intermediate.c:3285
IMRPhenomXHM_Intermediate_Amplitude_Veto
void IMRPhenomXHM_Intermediate_Amplitude_Veto(double *int1, double *int2, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
Definition: LALSimIMRPhenomXHM_intermediate.c:1995
IMRPhenomXHM_Inter_Amp_33_int4
static double IMRPhenomXHM_Inter_Amp_33_int4(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:586
IMRPhenomXHM_Inter_Phase_33_p3
static double IMRPhenomXHM_Inter_Phase_33_p3(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2873
IMRPhenomXHM_Inter_Amp_21_int1
static double IMRPhenomXHM_Inter_Amp_21_int1(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:49
IMRPhenomXHM_Inter_Amp_32_dint0
static double IMRPhenomXHM_Inter_Amp_32_dint0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:453
Update_Intermediate_Amplitude_Coefficients
void Update_Intermediate_Amplitude_Coefficients(IMRPhenomXHMAmpCoefficients *pAmp, int IntAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2207
IMRPhenomXHM_Inter_Phase_44_p1
static double IMRPhenomXHM_Inter_Phase_44_p1(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2718
IMRPhenomXHM_Inter_Phase_33_p4
static double IMRPhenomXHM_Inter_Phase_33_p4(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2975
IMRPhenomXHM_Inter_Phase_32_p4
static double IMRPhenomXHM_Inter_Phase_32_p4(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:3000
IMRPhenomXHM_Inter_Amp_32_int3
static double IMRPhenomXHM_Inter_Amp_32_int3(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:660
IMRPhenomXHM_Inter_Phase_44_p2
static double IMRPhenomXHM_Inter_Phase_44_p2(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2820
IMRPhenomXHM_Inter_Phase_44_p4
static double IMRPhenomXHM_Inter_Phase_44_p4(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:3026
IMRPhenomXHM_Inter_Phase_21_p5
static double IMRPhenomXHM_Inter_Phase_21_p5(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:3052
IMRPhenomXHM_Inter_Amp_21_int0
static double IMRPhenomXHM_Inter_Amp_21_int0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:351
IMRPhenomXHM_Inter_Phase_32_p1
static double IMRPhenomXHM_Inter_Phase_32_p1(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2692
IMRPhenomXHM_Intermediate_Amp_delta3
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)
Definition: LALSimIMRPhenomXHM_intermediate.c:1409
IMRPhenomXHM_Inter_Amp_21_int2
static double IMRPhenomXHM_Inter_Amp_21_int2(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:86
IMRPhenomXHM_Inter_Amp_44_int0
static double IMRPhenomXHM_Inter_Amp_44_int0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:473
IMRPhenomXHM_Intermediate_Amp_delta2
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)
Definition: LALSimIMRPhenomXHM_intermediate.c:1183
IMRPhenomXHM_Inter_Amp_44_int3
static double IMRPhenomXHM_Inter_Amp_44_int3(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:611
IMRPhenomXHM_Inter_Phase_AnsatzInt
static double IMRPhenomXHM_Inter_Phase_AnsatzInt(double ff, IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMPhaseCoefficients *pPhase)
Definition: LALSimIMRPhenomXHM_intermediate.c:3258
IMRPhenomXHM_Inter_Amp_32_int1
static double IMRPhenomXHM_Inter_Amp_32_int1(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:194
IMRPhenomXHM_Inter_Amp_32_int4
static double IMRPhenomXHM_Inter_Amp_32_int4(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:686
IMRPhenomXHM_Intermediate_Amp_delta4
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)
Definition: LALSimIMRPhenomXHM_intermediate.c:1623
IMRPhenomXHM_Inter_Phase_44_p3
static double IMRPhenomXHM_Inter_Phase_44_p3(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2923
IMRPhenomXHM_Inter_Phase_33_p1
static double IMRPhenomXHM_Inter_Phase_33_p1(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2668
IMRPhenomXHM_Inter_Amp_44_dint0
static double IMRPhenomXHM_Inter_Amp_44_dint0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:495
IMRPhenomXHM_Intermediate_Amp_Ansatz
static double IMRPhenomXHM_Intermediate_Amp_Ansatz(IMRPhenomX_UsefulPowers *powers_of_f, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
Definition: LALSimIMRPhenomXHM_intermediate.c:1900
CrossZeroP3
int CrossZeroP3(double a0, double a1, double a2, double a3, double fstart, double fend)
Definition: LALSimIMRPhenomXHM_intermediate.c:2039
IMRPhenomXHM_Inter_Phase_32_p5
static double IMRPhenomXHM_Inter_Phase_32_p5(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:3100
IMRPhenomXHM_Inter_Amp_21_int3
static double IMRPhenomXHM_Inter_Amp_21_int3(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:514
IMRPhenomXHM_Inter_Phase_33_p5
static double IMRPhenomXHM_Inter_Phase_33_p5(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:3075
CrossZeroP4
int CrossZeroP4(double a0, double a1, double a2, double a3, double a4, double fstart, double fend)
Definition: LALSimIMRPhenomXHM_intermediate.c:2088
IMRPhenomXHM_Inter_Amp_33_dint0
static double IMRPhenomXHM_Inter_Amp_33_dint0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:411
IMRPhenomXHM_Inter_Phase_33_p2
static double IMRPhenomXHM_Inter_Phase_33_p2(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2770
IMRPhenomXHM_Inter_Phase_21_p2
static double IMRPhenomXHM_Inter_Phase_21_p2(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2744
IMRPhenomXHM_Inter_Phase_33_p6
static double IMRPhenomXHM_Inter_Phase_33_p6(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:3177
IMRPhenomXHM_Inter_Amp_33_int0
static double IMRPhenomXHM_Inter_Amp_33_int0(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:391
IMRPhenomXHM_Inter_Amp_44_int1
static double IMRPhenomXHM_Inter_Amp_44_int1(IMRPhenomXWaveformStruct *pWF, int InterAmpFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:273
IMRPhenomXHM_Inter_Phase_32_p3
static double IMRPhenomXHM_Inter_Phase_32_p3(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2897
IMRPhenomXHM_Inter_Phase_21_p4
static double IMRPhenomXHM_Inter_Phase_21_p4(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2949
IMRPhenomXHM_Inter_Phase_32_p2
static double IMRPhenomXHM_Inter_Phase_32_p2(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2794
IMRPhenomXHM_Intermediate_Amp_delta1
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)
Definition: LALSimIMRPhenomXHM_intermediate.c:945
IMRPhenomXHM_Intermediate_Amp_Coefficients
void IMRPhenomXHM_Intermediate_Amp_Coefficients(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22)
Definition: LALSimIMRPhenomXHM_intermediate.c:2538
CrossZeroP5
int CrossZeroP5(double a0, double a1, double a2, double a3, double a4, double a5, double fstart, double fend)
Definition: LALSimIMRPhenomXHM_intermediate.c:2152
IMRPhenomXHM_Inter_Phase_21_p1
static double IMRPhenomXHM_Inter_Phase_21_p1(IMRPhenomXWaveformStruct *pWF, int InterPhaseFlag)
Definition: LALSimIMRPhenomXHM_intermediate.c:2642
LALSimIMRPhenomXHM_intermediate.h
SpheroidalToSpherical
COMPLEX16 SpheroidalToSpherical(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22, IMRPhenomXHMAmpCoefficients *pAmplm, IMRPhenomXHMPhaseCoefficients *pPhaselm, IMRPhenomXHMWaveformStruct *pWFlm, IMRPhenomXWaveformStruct *pWF22)
Definition: LALSimIMRPhenomXHM_internals.c:2653
IMRPhenomXHM_RD_Amp_DAnsatz
static double IMRPhenomXHM_RD_Amp_DAnsatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWF, IMRPhenomXHMAmpCoefficients *pAmp)
Definition: LALSimIMRPhenomXHM_ringdown.c:955
IMRPhenomXHM_RD_Amp_Ansatz
static double IMRPhenomXHM_RD_Amp_Ansatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
Definition: LALSimIMRPhenomXHM_ringdown.c:897
IMRPhenomXHM_RD_Amp_NDAnsatz
static double IMRPhenomXHM_RD_Amp_NDAnsatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMPhaseCoefficients *pPhase, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXAmpCoefficients *pAmp22, IMRPhenomXPhaseCoefficients *pPhase22, IMRPhenomXWaveformStruct *pWF22)
Definition: LALSimIMRPhenomXHM_ringdown.c:1002
f14
static const REAL8 f14[]
Definition: LALSimIMRTEOBResumS_Internals.c:1920
q3
static REAL8 UNUSED q3(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
Definition: LALSimInspiralOptimizedCoefficientsEccentricityFD.c:1195
q2
static REAL8 UNUSED q2(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
Definition: LALSimInspiralOptimizedCoefficientsEccentricityFD.c:1175
q1
static REAL8 UNUSED q1(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
Definition: LALSimInspiralOptimizedCoefficientsEccentricityFD.c:1155
z1
static REAL8 UNUSED z1(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
Definition: LALSimInspiralOptimizedCoefficientsEccentricityFD.c:505
z2
static REAL8 UNUSED z2(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
Definition: LALSimInspiralOptimizedCoefficientsEccentricityFD.c:527
i
double i
Definition: bh_ringdown.c:118
a4
const double a4
Definition: exact_derivatives-HrealHeader.c:34
a2
const double a2
Definition: exact_derivatives-HrealHeader.c:33
w
const double w
Definition: exact_derivatives-HrealHeader.c:106
LAL_PI
#define LAL_PI
REAL8
double REAL8
UINT2
uint16_t UINT2
INT4
int32_t INT4
a
static const INT4 a
XLAL_ERROR_VOID
#define XLAL_ERROR_VOID(...)
XLAL_ERROR_REAL8
#define XLAL_ERROR_REAL8(...)
XLALPrintError
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_EFUNC
XLAL_EFUNC
XLAL_EDOM
XLAL_EDOM
XLAL_EINVAL
XLAL_EINVAL
x
list x
p
list p
IMRPhenomX_UsefulPowers
Definition: LALSimIMRPhenomX_internals.h:249
IMRPhenomX_UsefulPowers::m_five
REAL8 m_five
Definition: LALSimIMRPhenomX_internals.h:269
IMRPhenomX_UsefulPowers::itself
REAL8 itself
Definition: LALSimIMRPhenomX_internals.h:263
IMRPhenomX_UsefulPowers::three
REAL8 three
Definition: LALSimIMRPhenomX_internals.h:260
IMRPhenomX_UsefulPowers::m_seven_sixths
REAL8 m_seven_sixths
Definition: LALSimIMRPhenomX_internals.h:278
IMRPhenomX_UsefulPowers::two
REAL8 two
Definition: LALSimIMRPhenomX_internals.h:261
IMRPhenomX_UsefulPowers::m_four
REAL8 m_four
Definition: LALSimIMRPhenomX_internals.h:268
IMRPhenomX_UsefulPowers::log
REAL8 log
Definition: LALSimIMRPhenomX_internals.h:280
IMRPhenomX_UsefulPowers::four
REAL8 four
Definition: LALSimIMRPhenomX_internals.h:259
IMRPhenomX_UsefulPowers::five
REAL8 five
Definition: LALSimIMRPhenomX_internals.h:258
IMRPhenomX_UsefulPowers::m_two
REAL8 m_two
Definition: LALSimIMRPhenomX_internals.h:266
IMRPhenomX_UsefulPowers::m_one
REAL8 m_one
Definition: LALSimIMRPhenomX_internals.h:265
IMRPhenomX_UsefulPowers::m_three
REAL8 m_three
Definition: LALSimIMRPhenomX_internals.h:267
IMRPhenomXHMAmpCoefficients::CollocationPointsFreqsAmplitudeInter
REAL8 CollocationPointsFreqsAmplitudeInter[N_MAX_COEFFICIENTS_AMPLITUDE_INTER]
Definition: LALSimIMRPhenomXHM_structs.h:214
IMRPhenomXHMAmpCoefficients::VersionCollocPtsInter
UINT2 VersionCollocPtsInter[N_MAX_COEFFICIENTS_AMPLITUDE_INTER]
Definition: LALSimIMRPhenomXHM_structs.h:218
IMRPhenomXHMAmpCoefficients::IntermediateAmpFits
ParameterSpaceFit IntermediateAmpFits[N_HIGHERMODES_IMPLEMENTED *N_MAX_COEFFICIENTS_AMPLITUDE_INTER]
Definition: LALSimIMRPhenomXHM_structs.h:205
IMRPhenomXHMAmpCoefficients::InterCoefficient
REAL8 InterCoefficient[N_MAX_COEFFICIENTS_AMPLITUDE_INTER]
Definition: LALSimIMRPhenomXHM_structs.h:216
IMRPhenomXHMAmpCoefficients::CollocationPointsValuesAmplitudeInter
REAL8 CollocationPointsValuesAmplitudeInter[N_MAX_COEFFICIENTS_AMPLITUDE_INTER]
Definition: LALSimIMRPhenomXHM_structs.h:215