Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALSimIMRPhenomXHM_inspiral.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2019 Marta Colleoni, Cecilio Garcia Quiros
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with with program; see the file COPYING. If not, write to the
16 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17 * MA 02110-1301 USA
18 */
19 //
20// LALSimIMRPhenomXHM_inspiral.c
21//
22//
23// Created by Marta on 06/02/2019.
24//
25
27
28/*************************************/
29/* */
30/* AMPLITUDE */
31/* */
32/*************************************/
33
34/* Fits of the collocation points in the inspiral across parameter space. 3 for each mode: iv1, iv2, iv3
35 [iv1,iv2,iv3]=[0.5*f^Ins_lm,0.75*f^Ins_lm,f^Ins_lm]
36 For more details about the reconstruction procedure, see Sec. IV.A
37 */
38
39 /* These parameter-space fits are documented in the supplementary material of https://dcc.ligo.org/P2000011-v2.
40 There are 2 Mathematica notebooks (one for amplitude and one for phase) that read the fits data and automatically generate the C-code below.
41 For more information read https://git.ligo.org/waveforms/reviews/imrphenomx/blob/master/documentation/ParspaceFits/README and the documentation in the notebooks. */
42
43// Spin parameter S = chi_PN = chi_eff - 38/113 * eta * (chi1 + chi2)
44// chi_eff = (m1*chi1 + m2*chi2)/(m1 + m2)
45
46/* Start of Amp Parameter Space Fits */
47
48static double IMRPhenomXHM_Insp_Amp_21_iv1(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag) {
49 double total=0;
50 switch (InspAmpFlag){
51 case 122018:{
52 double eta = pWF->eta;
53 double S = pWF->chiPNHat;
54 double eta2,eta3,eta4,eta5,S2;
55 eta2 = pow(eta,2);
56 eta3 = pow(eta,3);
57 eta4 = pow(eta,4);
58 eta5 = pow(eta,5);
59 S2 = pow(S,2);
60 double noSpin = sqrt(1.-4.*eta)*(0.037868557189995156 + 0.10740090317702103*eta + 1.963812986867654*eta2 - 16.706455229589558*eta3 + 69.75910808095745*eta4 - 98.3062466823662*eta5);
61 double eqSpin = sqrt(1.-4.*eta)*S*(-0.007963757232702219 + 0.10627108779259965*eta - 0.008044970210401218*S + eta2*(-0.4735861262934258 - 0.5985436493302649*S - 0.08217216660522082*S2));
62 double uneqSpin = -0.257787704938017*pWF->dchi*eta2*(1. + 8.75928187268504*eta2) - 0.2597503605427412*pWF->dchi*eta2*S;
63 total = noSpin + eqSpin + uneqSpin;
64 break;
65 }
66 case 122022:{
67 double eta = pWF->eta;
68 double delta = pWF->delta;
69 double sqroot = sqrt(eta);
70 double S = pWF->chiPNHat;
71 double chidiff = pWF->dchi_half;
72 double eta1 = eta;
73 double eta2 = eta1 * eta1;
74 double eta3 = eta1 * eta2;
75 double eta4 = eta1 * eta3;
76 double eta5 = eta1 * eta4;
77 double S1 = S;
78 double S2 = S1 * S1;
79 double chidiff1 = chidiff;
80 total = fabs(chidiff1*eta5*(-3962.5020052272976 + 987.635855365408*S1 - 134.98527058315528*S2) + delta*(19.30531354642419 + 16.6640319856064*eta1 - 120.58166037019478*eta2 + 220.77233521626252*eta3)*sqroot + chidiff1*delta*(31.364509907424765*eta1 - 843.6414532232126*eta2 + 2638.3077554662905*eta3)*sqroot + chidiff1*delta*(32.374226994179054*eta1 - 202.86279451816662*eta2 + 347.1621871204769*eta3)*S1*sqroot + delta*S1*(-16.75726972301224*(1.1787350890261943 - 7.812073811917883*eta1 + 99.47071002831267*eta2 - 500.4821414428368*eta3 + 876.4704270866478*eta4) + 2.3439955698372663*(0.9373952326655807 + 7.176140122833879*eta1 - 279.6409723479635*eta2 + 2178.375177755584*eta3 - 4768.212511142035*eta4)*S1)*sqroot);
81 break;
82 }
83 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Insp_Amp_21_iv1: version %i is not valid.", InspAmpFlag);}
84 }
85 return total;
86}
87
88static double IMRPhenomXHM_Insp_Amp_21_iv2(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag) {
89 double total=0;
90 switch (InspAmpFlag){
91 case 122018:{
92 double eta = pWF->eta;
93 double S = pWF->chiPNHat;
94 double delta=sqrt(1.-4.*eta),eta2,eta3,eta4,S2;
95 eta2 = pow(eta,2);
96 eta3 = pow(eta,3);
97 eta4 = pow(eta,4);
98 S2 = pow(S,2);
99 double noSpin = sqrt(1.-4.*eta)*(0.05511628628738656 - 0.12579599745414977*eta + 2.831411618302815*eta2 - 14.27268643447161*eta3 + 28.3307320191161*eta4);
100 double eqSpin = sqrt(1.-4.*eta)*S*(-0.008692738851491525 + eta*(0.09512553997347649 + 0.116470975986383*S) - 0.009520793625590234*S + eta2*(-0.3409769288480959 - 0.8321002363767336*S - 0.13099477081654226*S2) - 0.006383232900211555*S2);
101 double uneqSpin = -0.2962753588645467*pWF->dchi*eta2*(1. + 1.3993978458830476*eta2) - 0.17100612756133535*pWF->dchi*eta2*S*(1. + 18.974303741922743*eta2*delta);
102 total = noSpin + eqSpin + uneqSpin;
103 break;
104 }
105 case 122022:{
106 double eta = pWF->eta;
107 double delta = pWF->delta;
108 double sqroot = sqrt(eta);
109 double S = pWF->chiPNHat;
110 double chidiff = pWF->dchi_half;
111 double eta1 = eta;
112 double eta2 = eta1 * eta1;
113 double eta3 = eta1 * eta2;
114 double eta4 = eta1 * eta3;
115 double eta5 = eta1 * eta4;
116 double S1 = S;
117 double S2 = S1 * S1;
118 double chidiff1 = chidiff;
119 double chidiff2 = chidiff1 * chidiff1;
120 total = fabs(chidiff1*eta5*(-2898.9172078672705 + 580.9465034962822*S1 + 22.251142639924076*S2) + delta*(chidiff2*(-18.541685007214625*eta1 + 166.7427445020744*eta2 - 417.5186332459383*eta3) + chidiff1*(41.61457952037761*eta1 - 779.9151607638761*eta2 + 2308.6520892707795*eta3))*sqroot + delta*(11.414934585404561 + 30.883118528233638*eta1 - 260.9979123967537*eta2 + 1046.3187137392433*eta3 - 1556.9475493549746*eta4)*sqroot + delta*S1*(-10.809007068469844*(1.1408749895922659 - 18.140470190766937*eta1 + 368.25127088896744*eta2 - 3064.7291458207815*eta3 + 11501.848278358668*eta4 - 16075.676528787526*eta5) + 1.0088254664333147*(1.2322739396680107 - 192.2461213084741*eta1 + 4257.760834055382*eta2 - 35561.24587952242*eta3 + 130764.22485304279*eta4 - 177907.92440833704*eta5)*S1)*sqroot + delta*(chidiff1*(36.88578491943111*eta1 - 321.2569602623214*eta2 + 748.6659668096737*eta3)*S1 + chidiff1*(-95.42418611585117*eta1 + 1217.338674959742*eta2 - 3656.192371615541*eta3)*S2)*sqroot);
121 break;
122 }
123 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Insp_Amp_21_iv2: version %i is not valid.", InspAmpFlag);}
124 }
125 return total;
126}
127
128static double IMRPhenomXHM_Insp_Amp_21_iv3(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag) {
129 double total=0;
130 switch (InspAmpFlag){
131 case 122018:{
132 double eta = pWF->eta;
133 double S = pWF->chiPNHat;
134 double delta=sqrt(1.-4.*eta),eta2,S2;
135 eta2 = pow(eta,2);
136 S2 = pow(S,2);
137 double noSpin = sqrt(1.-4.*eta)*(0.059110044024271766 - 0.0024538774422098405*eta + 0.2428578654261086*eta2);
138 double eqSpin = sqrt(1.-4.*eta)*S*(-0.007044339356171243 - 0.006952154764487417*S + eta2*(-0.016643018304732624 - 0.12702579620537421*S + 0.004623467175906347*S2) - 0.007685497720848461*S2);
139 double uneqSpin = -0.3172310538516028*pWF->dchi*(1. - 2.9155919835488024*eta2)*eta2 - 0.11975485688200693*pWF->dchi*eta2*S*(1. + 17.27626751837825*eta2*delta);
140 total = noSpin + eqSpin + uneqSpin;
141 break;
142 }
143 case 122022:{
144 double eta = pWF->eta;
145 double delta = pWF->delta;
146 double sqroot = sqrt(eta);
147 double S = pWF->chiPNHat;
148 double chidiff = pWF->dchi_half;
149 double eta1 = eta;
150 double eta2 = eta1 * eta1;
151 double eta3 = eta1 * eta2;
152 double eta4 = eta1 * eta3;
153 double eta5 = eta1 * eta4;
154 double S1 = S;
155 double S2 = S1 * S1;
156 double chidiff1 = chidiff;
157 total = fabs(chidiff1*eta5*(-2282.9983216879655 + 157.94791186394787*S1 + 16.379731479465033*S2) + chidiff1*delta*(21.935833431534224*eta1 - 460.7130131927895*eta2 + 1350.476411541137*eta3)*sqroot + delta*(5.390240326328237 + 69.01761987509603*eta1 - 568.0027716789259*eta2 + 2435.4098320959706*eta3 - 3914.3390484239667*eta4)*sqroot + chidiff1*delta*(29.731007410186827*eta1 - 372.09609843131386*eta2 + 1034.4897198648962*eta3)*S1*sqroot + delta*S1*(-7.1976397556450715*(0.7603360145475428 - 6.587249958654174*eta1 + 120.87934060776237*eta2 - 635.1835857158857*eta3 + 1109.0598539312573*eta4) - 0.0811847192323969*(7.951454648295709 + 517.4039644814231*eta1 - 9548.970156895082*eta2 + 52586.63520999897*eta3 - 93272.17990295641*eta4)*S1 - 0.28384547935698246*(-0.8870770459576875 + 180.0378964169756*eta1 - 2707.9572896559484*eta2 + 14158.178124971111*eta3 - 24507.800226675925*eta4)*S2)*sqroot);
158 break;
159 }
160 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Insp_Amp_21_iv3: version %i is not valid.", InspAmpFlag);}
161 }
162 return total;
163}
164
165static double IMRPhenomXHM_Insp_Amp_33_iv1(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag) {
166 double total=0;
167 switch (InspAmpFlag){
168 case 122018:{
169 double eta = pWF->eta;
170 double S = pWF->chiPNHat;
171 double eta2,S2;
172 eta2 = pow(eta,2);
173 S2 = pow(S,2);
174 double noSpin = (sqrt(1.-4.*eta)*(-0.056586690934283326 - 0.14374841547279146*eta + 0.5584776628959615*eta2))/(-0.3996185676368123 + eta);
175 double eqSpin = sqrt(1.-4.*eta)*S*((0.056042044149691175 + 0.12482426029674777*S)*S + eta*(2.1108074577110343 - 1.7827773156978863*S2) + eta2*(-7.657635515668849 - 0.07646730296478217*S + 5.343277927456605*S2));
176 double uneqSpin = 0.45866449225302536*pWF->dchi*(1. - 9.603750707244906*eta2)*eta2;
177 total = noSpin + eqSpin + uneqSpin;
178 break;
179 }
180 case 122022:{
181 double eta = pWF->eta;
182 double delta = pWF->delta;
183 double sqroot = sqrt(eta);
184 double S = pWF->chiPNHat;
185 double chidiff = pWF->dchi_half;
186 double eta1 = eta;
187 double eta2 = eta1 * eta1;
188 double eta3 = eta1 * eta2;
189 double eta4 = eta1 * eta3;
190 double eta5 = eta1 * eta4;
191 double S1 = S;
192 double S2 = S1 * S1;
193 double chidiff1 = chidiff;
194 total = chidiff1*eta5*(155.1434307076563 + 26.852777193715088*S1 + 1.4157230717300835*S2) + chidiff1*delta*(6.296698171560171*eta1 + 15.81328761563562*eta2 - 141.85538063933927*eta3)*sqroot + delta*(20.94372147101354 + 68.14577638017842*eta1 - 898.470298591732*eta2 + 4598.64854748635*eta3 - 8113.199260593833*eta4)*sqroot + chidiff1*delta*(29.221863857271703*eta1 - 348.1658322276406*eta2 + 965.4670353331536*eta3)*S1*sqroot + delta*S1*(-9.753610761811967*(1.7819678168496158 - 44.07982999150369*eta1 + 750.8933447725581*eta2 - 5652.44754829634*eta3 + 19794.855873435758*eta4 - 26407.40988450443*eta5) + 0.014210376114848208*(-196.97328616330392 + 7264.159472864562*eta1 - 125763.47850622259*eta2 + 1.1458022059130718e6*eta3 - 4.948175330328345e6*eta4 + 7.911048294733888e6*eta5)*S1 - 0.26859293613553986*(-8.029069605349488 + 888.7768796633982*eta1 - 16664.276483466252*eta2 + 128973.72291098491*eta3 - 462437.2690007375*eta4 + 639989.1197424605*eta5)*S2)*sqroot;
195 break;
196 }
197 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Insp_Amp_33_iv1: version %i is not valid.", InspAmpFlag);}
198 }
199 return total;
200}
201
202static double IMRPhenomXHM_Insp_Amp_33_iv2(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag) {
203 double total=0;
204 switch (InspAmpFlag){
205 case 122018:{
206 double eta = pWF->eta;
207 double S = pWF->chiPNHat;
208 double eta2,eta3,eta4,eta5,eta6,S2;
209 eta2 = pow(eta,2);
210 eta3 = pow(eta,3);
211 eta4 = pow(eta,4);
212 eta5 = pow(eta,5);
213 eta6 = pow(eta,6);
214 S2 = pow(S,2);
215 double noSpin = sqrt(1.-4.*eta)*(0.2137734510411439 - 0.7692194209223682*eta + 26.10570221351058*eta2 - 316.0643979123107*eta3 + 2090.9063511488234*eta4 - 6897.3285171507105*eta5 + 8968.893362362503*eta6);
216 double eqSpin = sqrt(1.-4.*eta)*S*(0.018546836505210842 + 0.05924304311104228*S + eta*(1.6484440612224325 - 0.4683932646001618*S - 2.110311135456494*S2) + 0.10701786057882816*S2 + eta2*(-6.51575737684721 + 1.6692205620001157*S + 8.351789152096782*S2));
217 double uneqSpin = 0.3929315188124088*pWF->dchi*(1. - 11.289452844364227*eta2)*eta2;
218 total = noSpin + eqSpin + uneqSpin;
219 break;
220 }
221 case 122022:{
222 double eta = pWF->eta;
223 double delta = pWF->delta;
224 double sqroot = sqrt(eta);
225 double S = pWF->chiPNHat;
226 double chidiff = pWF->dchi_half;
227 double eta1 = eta;
228 double eta2 = eta1 * eta1;
229 double eta3 = eta1 * eta2;
230 double eta4 = eta1 * eta3;
231 double eta5 = eta1 * eta4;
232 double S1 = S;
233 double S2 = S1 * S1;
234 double chidiff1 = chidiff;
235 total = chidiff1*eta5*(161.62678370819597 + 37.141092711336846*S1 - 0.16889712161410445*S2) + chidiff1*delta*(3.4895829486899825*eta1 + 51.07954458810889*eta2 - 249.71072528701757*eta3)*sqroot + delta*(12.501397517602173 + 35.75290806646574*eta1 - 357.6437296928763*eta2 + 1773.8883882162215*eta3 - 3100.2396041211605*eta4)*sqroot + chidiff1*delta*(13.854211287141906*eta1 - 135.54916401086845*eta2 + 327.2467193417936*eta3)*S1*sqroot + delta*S1*(-5.2580116732827085*(1.7794900975289085 - 48.20753331991333*eta1 + 861.1650630146937*eta2 - 6879.681319382729*eta3 + 25678.53964955809*eta4 - 36383.824902258915*eta5) + 0.028627002336747746*(-50.57295946557892 + 734.7581857539398*eta1 - 2287.0465658878725*eta2 + 15062.821881048358*eta3 - 168311.2370167227*eta4 + 454655.37836367317*eta5)*S1 - 0.15528289788512326*(-12.738184090548508 + 1129.44485109116*eta1 - 25091.14888164863*eta2 + 231384.03447562453*eta3 - 953010.5908118751*eta4 + 1.4516597366230418e6*eta5)*S2)*sqroot;
236 break;
237 }
238 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Insp_Amp_33_iv2: version %i is not valid.", InspAmpFlag);}
239 }
240 return total;
241}
242
243static double IMRPhenomXHM_Insp_Amp_33_iv3(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag) {
244 double total=0;
245 switch (InspAmpFlag){
246 case 122018:{
247 double eta = pWF->eta;
248 double S = pWF->chiPNHat;
249 double eta2,eta3,eta4,eta5,eta6,S2;
250 eta2 = pow(eta,2);
251 eta3 = pow(eta,3);
252 eta4 = pow(eta,4);
253 eta5 = pow(eta,5);
254 eta6 = pow(eta,6);
255 S2 = pow(S,2);
256 double noSpin = sqrt(1.-4.*eta)*(0.2363760327127446 + 0.2855410252403732*eta - 10.159877125359897*eta2 + 162.65372389693505*eta3 - 1154.7315106095564*eta4 + 3952.61320206691*eta5 - 5207.67472857814*eta6);
257 double eqSpin = sqrt(1.-4.*eta)*S*(0.04573095188775319 + 0.048249943132325494*S + eta*(0.15922377052827502 - 0.1837289613228469*S - 0.2834348500565196*S2) + 0.052963737236081304*S2);
258 double uneqSpin = 0.25187274502769835*pWF->dchi*(1. - 12.172961866410864*eta2)*eta2;
259 total = noSpin + eqSpin + uneqSpin;
260 break;
261 }
262 case 122022:{
263 double eta = pWF->eta;
264 double delta = pWF->delta;
265 double sqroot = sqrt(eta);
266 double S = pWF->chiPNHat;
267 double chidiff = pWF->dchi_half;
268 double eta1 = eta;
269 double eta2 = eta1 * eta1;
270 double eta3 = eta1 * eta2;
271 double eta4 = eta1 * eta3;
272 double eta5 = eta1 * eta4;
273 double eta6 = eta1 * eta5;
274 double eta7 = eta1 * eta6;
275 double S1 = S;
276 double S2 = S1 * S1;
277 double chidiff1 = chidiff;
278 double chidiff2 = chidiff1 * chidiff1;
279 total = chidiff1*delta*(-0.5869777957488564*eta1 + 32.65536124256588*eta2 - 110.10276573567405*eta3) + chidiff1*delta*(3.524800489907584*eta1 - 40.26479860265549*eta2 + 113.77466499598913*eta3)*S1 + delta*S1*(-1.2846335585108297*(0.09991079016763821 + 1.37856806162599*eta1 + 23.26434219690476*eta2 - 34.842921754693386*eta3 - 70.83896459998664*eta4) - 0.03496714763391888*(-0.230558571912664 + 188.38585449575902*eta1 - 3736.1574640444287*eta2 + 22714.70643022915*eta3 - 43221.0453556626*eta4)*S1) + chidiff1*eta7*(2667.3441342894776 + 47.94869769580204*chidiff2 + 793.5988192446642*S1 + 293.89657731755483*S2) + delta*(5.148353856800232 + 148.98231189649468*eta1 - 2774.5868652930294*eta2 + 29052.156454239772*eta3 - 162498.31493332976*eta4 + 460912.76402476896*eta5 - 521279.50781871413*eta6)*sqroot;
280 break;
281 }
282 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Insp_Amp_33_iv3: version %i is not valid.", InspAmpFlag);}
283 }
284 return total;
285}
286
287static double IMRPhenomXHM_Insp_Amp_32_iv1(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag) {
288 double total=0;
289 switch (InspAmpFlag){
290 case 122018:{
291 double eta = pWF->eta;
292 double S = pWF->chiPNHat;
293 double delta=sqrt(1.-4.*eta),eta2,eta3,eta4,eta5,eta6,eta8,S2,S3;
294 eta2 = pow(eta,2);
295 eta3 = pow(eta,3);
296 eta4 = pow(eta,4);
297 eta5 = pow(eta,5);
298 eta6 = pow(eta,6);
299 eta8 = pow(eta,8);
300 S2 = pow(S,2);
301 S3 = pow(S,3);
302 double noSpin = sqrt(1. - 3.*eta)*(0.019069933430190773 - 0.19396651989685837*eta + 11.95224600241255*eta2 - 158.90113442757382*eta3 + 1046.65239329071*eta4 - 3476.940285294999*eta5 + 4707.249209858949*eta6);
303 double eqSpin = sqrt(1. - 3.*eta)*S*(0.0046910348789512895 + 0.40231360805609434*eta - 0.0038263656140933152*S + 0.018963579407636953*S2 + eta2*(-1.955352354930108 + 2.3753413452420133*S - 0.9085620866763245*S3) + 0.02738043801805805*S3 + eta3*(7.977057990568723 - 7.9259853291789515*S + 0.49784942656123987*S2 + 5.2255665027119145*S3));
304 double uneqSpin = 0.058560321425018165*pow(pWF->dchi,2)*(1. - 19.936477485971217*eta2)*eta2 + 1635.4240644598524*pWF->dchi*eta8*delta + 0.2735219358839411*pWF->dchi*eta2*S*delta;
305 total = noSpin + eqSpin + uneqSpin;
306 break;
307 }
308 case 122022:{
309 double eta = pWF->eta;
310 double delta = pWF->delta;
311 double sqroot = sqrt(eta);
312 double S = pWF->chiPNHat;
313 double chidiff = pWF->dchi_half;
314 double eta1 = eta;
315 double eta2 = eta1 * eta1;
316 double eta3 = eta1 * eta2;
317 double eta4 = eta1 * eta3;
318 double eta5 = eta1 * eta4;
319 double eta6 = eta1 * eta5;
320 double S1 = S;
321 double chidiff1 = chidiff;
322 double chidiff2 = chidiff1 * chidiff1;
323 total = (chidiff1*delta*(-0.739317114582042*eta1 - 47.473246070362634*eta2 + 278.9717709112207*eta3 - 566.6420939162068*eta4) + chidiff2*(-0.5873680378268906*eta1 + 6.692187014925888*eta2 - 24.37776782232888*eta3 + 23.783684827838247*eta4))*sqroot + (3.2940434453819694 + 4.94285331708559*eta1 - 343.3143244815765*eta2 + 3585.9269057886418*eta3 - 19279.186145681153*eta4 + 51904.91007211022*eta5 - 55436.68857586653*eta6)*sqroot + chidiff1*delta*(12.488240781993923*eta1 - 209.32038774208385*eta2 + 1160.9833883184604*eta3 - 2069.5349737049073*eta4)*S1*sqroot + S1*(0.6343034651912586*(-2.5844888818001737 + 78.98200041834092*eta1 - 1087.6241783616488*eta2 + 7616.234910399297*eta3 - 24776.529123239357*eta4 + 30602.210950069973*eta5) - 0.062088720220899465*(6.5586380356588565 + 36.01386705325694*eta1 - 3124.4712274775407*eta2 + 33822.437731298516*eta3 - 138572.93700180828*eta4 + 198366.10615196894*eta5)*S1)*sqroot;
324 break;
325 }
326 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Insp_Amp_32_iv1: version %i is not valid.", InspAmpFlag);}
327 }
328 return total;
329}
330
331static double IMRPhenomXHM_Insp_Amp_32_iv2(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag) {
332 double total=0;
333 switch (InspAmpFlag){
334 case 122018:{
335 double eta = pWF->eta;
336 double S = pWF->chiPNHat;
337 double delta=sqrt(1.-4.*eta),eta2,eta3,eta4,eta8,S2,S3;
338 eta2 = pow(eta,2);
339 eta3 = pow(eta,3);
340 eta4 = pow(eta,4);
341 eta8 = pow(eta,8);
342 S2 = pow(S,2);
343 S3 = pow(S,3);
344 double noSpin = sqrt(1. - 3.*eta)*(0.024621376891809633 - 0.09692699636236377*eta + 2.7200998230836158*eta2 - 16.160563094841066*eta3 + 32.930430889650836*eta4);
345 double eqSpin = sqrt(1. - 3.*eta)*S*(0.008522695567479373 - 1.1104639098529456*eta2 - 0.00362963820787208*S + 0.016978054142418417*S2 + eta*(0.24280554040831698 + 0.15878436411950506*S - 0.1470288177047577*S3) + 0.029465887557447824*S3 + eta3*(4.649438233164449 - 0.7550771176087877*S + 0.3381436950547799*S2 + 2.5663386135613093*S3));
346 double uneqSpin = -0.007061187955941243*pow(pWF->dchi,2)*(1. - 2.024701925508361*eta2)*eta2 + 215.06940561269835*pWF->dchi*eta8*delta + 0.1465612311350642*pWF->dchi*eta2*S*delta;
347 total = noSpin + eqSpin + uneqSpin;
348 break;
349 }
350 case 122022:{
351 double eta = pWF->eta;
352 double delta = pWF->delta;
353 double sqroot = sqrt(eta);
354 double S = pWF->chiPNHat;
355 double chidiff = pWF->dchi_half;
356 double eta1 = eta;
357 double eta2 = eta1 * eta1;
358 double eta3 = eta1 * eta2;
359 double eta4 = eta1 * eta3;
360 double eta5 = eta1 * eta4;
361 double eta6 = eta1 * eta5;
362 double S1 = S;
363 double chidiff1 = chidiff;
364 double chidiff2 = chidiff1 * chidiff1;
365 total = (chidiff2*(-0.03940151060321499*eta1 + 1.9034209537174116*eta2 - 8.78587250202154*eta3) + chidiff1*delta*(-1.704299788495861*eta1 - 4.923510922214181*eta2 + 0.36790005839460627*eta3))*sqroot + (2.2911849711339123 - 5.1846950040514335*eta1 + 60.10368251688146*eta2 - 1139.110227749627*eta3 + 7970.929280907627*eta4 - 25472.73682092519*eta5 + 30950.67053883646*eta6)*sqroot + S1*(0.7718201508695763*(-1.3012906461000349 + 26.432880113146012*eta1 - 186.5001124789369*eta2 + 712.9101229418721*eta3 - 970.2126139442341*eta4) + 0.04832734931068797*(-5.9999628512498315 + 78.98681284391004*eta1 + 1.8360177574514709*eta2 - 2537.636347529708*eta3 + 6858.003573909322*eta4)*S1)*sqroot;
366 break;
367 }
368 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Insp_Amp_32_iv2: version %i is not valid.", InspAmpFlag);}
369 }
370 return total;
371}
372
373static double IMRPhenomXHM_Insp_Amp_32_iv3(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag) {
374 double total=0;
375 switch (InspAmpFlag){
376 case 122018:{
377 double eta = pWF->eta;
378 double S = pWF->chiPNHat;
379 double delta=sqrt(1.-4.*eta),eta2,eta3,eta8,S2,S3;
380 eta2 = pow(eta,2);
381 eta3 = pow(eta,3);
382 eta8 = pow(eta,8);
383 S2 = pow(S,2);
384 S3 = pow(S,3);
385 double noSpin = (sqrt(1. - 3.*eta)*(-0.006150151041614737 + 0.017454430190035*eta + 0.02620962593739105*eta2 - 0.019043090896351363*eta3))/(-0.2655505633361449 + eta);
386 double eqSpin = sqrt(1. - 3.*eta)*S*(0.011073381681404716 + 0.00347699923233349*S + eta*S*(0.05592992411391443 - 0.15666140197050316*S2) + 0.012079324401547036*S2 + eta2*(0.5440307361144313 - 0.008730335213434078*S + 0.04615964369925028*S2 + 0.6703688097531089*S3) + 0.016323101357296865*S3);
387 double uneqSpin = -0.020140175824954427*pow(pWF->dchi,2)*(1. - 12.675522774051249*eta2)*eta2 - 417.3604094454253*pWF->dchi*eta8*delta + 0.10464021067936538*pWF->dchi*eta2*S*delta;
388 total = noSpin + eqSpin + uneqSpin;
389 break;
390 }
391 case 122022:{
392 double eta = pWF->eta;
393 double delta = pWF->delta;
394 double sqroot = sqrt(eta);
395 double S = pWF->chiPNHat;
396 double chidiff = pWF->dchi_half;
397 double eta1 = eta;
398 double eta2 = eta1 * eta1;
399 double eta3 = eta1 * eta2;
400 double eta4 = eta1 * eta3;
401 double eta5 = eta1 * eta4;
402 double eta6 = eta1 * eta5;
403 double S1 = S;
404 double chidiff1 = chidiff;
405 double chidiff2 = chidiff1 * chidiff1;
406 total = (chidiff2*(-0.6358511175987503*eta1 + 5.555088747533164*eta2 - 14.078156877577733*eta3) + chidiff1*delta*(0.23205448591711159*eta1 - 19.46049432345157*eta2 + 36.20685853857613*eta3))*sqroot + (1.1525594672495008 + 7.380126197972549*eta1 - 17.51265776660515*eta2 - 976.9940395257111*eta3 + 8880.536804741967*eta4 - 30849.228936891763*eta5 + 38785.53683146884*eta6)*sqroot + chidiff1*delta*(1.904350804857431*eta1 - 25.565242391371093*eta2 + 80.67120303906654*eta3)*S1*sqroot + S1*(0.785171689871352*(-0.4634745514643032 + 18.70856733065619*eta1 - 167.9231114864569*eta2 + 744.7699462372949*eta3 - 1115.008825153004*eta4) + 0.13469300326662165*(-2.7311391326835133 + 72.17373498208947*eta1 - 483.7040402103785*eta2 + 1136.8367114738041*eta3 - 472.02962341590774*eta4)*S1)*sqroot;
407 break;
408 }
409 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Insp_Amp_32_iv3: version %i is not valid.", InspAmpFlag);}
410 }
411 return total;
412}
413
414static double IMRPhenomXHM_Insp_Amp_44_iv1(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag) {
415 double total=0;
416 switch (InspAmpFlag){
417 case 122018:{
418 double eta = pWF->eta;
419 double S = pWF->chiPNHat;
420 double eta2,eta3,eta4,S2;
421 eta2 = pow(eta,2);
422 eta3 = pow(eta,3);
423 eta4 = pow(eta,4);
424 S2 = pow(S,2);
425 double noSpin = sqrt(1. - 3.*eta)*(0.06190013067931406 + 0.1928897813606222*eta + 1.9024723168424225*eta2 - 15.988716302668415*eta3 + 35.21461767354364*eta4);
426 double eqSpin = sqrt(1. - 3.*eta)*S*(0.011454874900772544 + 0.044702230915643903*S + eta*(0.6600413908621988 + 0.12149520289658673*S - 0.4482406547006759*S2) + 0.07327810908370004*S2 + eta2*(-2.1705970511116486 - 0.6512813450832168*S + 1.1237234702682313*S2));
427 double uneqSpin = 0.4766851579723911*pWF->dchi*(1. - 15.950025762198988*eta2)*eta2 + 0.127900699645338*pow(pWF->dchi,2)*(1. - 15.79329306044842*eta2)*eta2;
428 total = noSpin + eqSpin + uneqSpin;
429 break;
430 }
431 case 122022:{
432 double eta = pWF->eta;
433 double delta = pWF->delta;
434 double sqroot = sqrt(eta);
435 double S = pWF->chiPNHat;
436 double chidiff = pWF->dchi_half;
437 double eta1 = eta;
438 double eta2 = eta1 * eta1;
439 double eta3 = eta1 * eta2;
440 double eta4 = eta1 * eta3;
441 double eta5 = eta1 * eta4;
442 double eta6 = eta1 * eta5;
443 double S1 = S;
444 double S2 = S1 * S1;
445 double chidiff1 = chidiff;
446 double chidiff2 = chidiff1 * chidiff1;
447 total = (chidiff1*delta*(0.5697308729057493*eta1 + 8.895576813118867*eta2 - 34.98399465240273*eta3) + chidiff2*(1.6370346538130884*eta1 - 14.597095790380884*eta2 + 33.182723737396294*eta3))*sqroot + (5.2601381002242595 - 3.557926105832778*eta1 - 138.9749850448088*eta2 + 603.7453704122706*eta3 - 923.5495700703648*eta4)*sqroot + S1*(-0.41839636169678796*(5.143510231379954 + 104.62892421207803*eta1 - 4232.508174045782*eta2 + 50694.024801783446*eta3 - 283097.33358214336*eta4 + 758333.2655404843*eta5 - 788783.0559069642*eta6) - 0.05653522061311774*(5.605483124564013 + 694.00652410087*eta1 - 17551.398321516353*eta2 + 165236.6480734229*eta3 - 761661.9645651339*eta4 + 1.7440315410044065e6*eta5 - 1.6010489769238676e6*eta6)*S1 - 0.023693246676754775*(16.437107575918503 - 2911.2154288136217*eta1 + 89338.32554683842*eta2 - 1.0803340811860575e6*eta3 + 6.255666490084672e6*eta4 - 1.7434160932177313e7*eta5 + 1.883460394974573e7*eta6)*S2)*sqroot;
448 break;
449 }
450 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Insp_Amp_44_iv1: version %i is not valid.", InspAmpFlag);}
451 }
452 return total;
453}
454
455static double IMRPhenomXHM_Insp_Amp_44_iv2(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag) {
456 double total=0;
457 switch (InspAmpFlag){
458 case 122018:{
459 double eta = pWF->eta;
460 double S = pWF->chiPNHat;
461 double delta=sqrt(1.-4.*eta),eta2,eta3,S2,S3;
462 eta2 = pow(eta,2);
463 eta3 = pow(eta,3);
464 S2 = pow(S,2);
465 S3 = pow(S,3);
466 double noSpin = 0.08406011695496626 - 0.1469952725049322*eta + 0.2997223283799925*eta2 - 1.2910560244510723*eta3;
467 double eqSpin = (0.023924074703897662 + 0.26110236039648027*eta - 1.1536009170220438*eta2)*S + (0.04479727299752669 - 0.1439868858871802*eta + 0.05736387085230215*eta2)*S2 + (0.06028104440131858 - 0.4759412992529712*eta + 1.1090751649419717*eta2)*S3;
468 double uneqSpin = 0.10346324686812074*pow(pWF->dchi,2)*(1. - 16.135903382018213*eta2)*eta2 + 0.2648241309154185*pWF->dchi*eta2*delta;
469 total = noSpin + eqSpin + uneqSpin;
470 break;
471 }
472 case 122022:{
473 double eta = pWF->eta;
474 double delta = pWF->delta;
475 double sqroot = sqrt(eta);
476 double S = pWF->chiPNHat;
477 double chidiff = pWF->dchi_half;
478 double eta1 = eta;
479 double eta2 = eta1 * eta1;
480 double eta3 = eta1 * eta2;
481 double eta4 = eta1 * eta3;
482 double eta5 = eta1 * eta4;
483 double eta6 = eta1 * eta5;
484 double S1 = S;
485 double S2 = S1 * S1;
486 double S3 = S1 * S2;
487 double chidiff1 = chidiff;
488 double chidiff2 = chidiff1 * chidiff1;
489 total = (chidiff2*(-0.8318312659717388*eta1 + 7.6541168007977864*eta2 - 16.648660653220123*eta3) + chidiff1*delta*(2.214478316304753*eta1 - 7.028104574328955*eta2 + 5.56587823143958*eta3))*sqroot + (3.173191054680422 + 6.707695566702527*eta1 - 155.22519772642607*eta2 + 604.0067075996933*eta3 - 876.5048298377644*eta4)*sqroot + chidiff1*delta*(4.749663394334708*eta1 - 42.62996105525792*eta2 + 97.01712147349483*eta3)*S1*sqroot + S1*(-0.2627203100303006*(6.460396349297595 - 52.82425783851536*eta1 - 552.1725902144143*eta2 + 12546.255587592654*eta3 - 81525.50289542897*eta4 + 227254.37897941095*eta5 - 234487.3875219032*eta6) - 0.008424003742397579*(-109.26773035716548 + 15514.571912666677*eta1 - 408022.6805482195*eta2 + 4.620165968920881e6*eta3 - 2.6446950627957724e7*eta4 + 7.539643948937692e7*eta5 - 8.510662871580401e7*eta6)*S1 - 0.008830881730801855*(-37.49992494976597 + 1359.7883958101172*eta1 - 23328.560285901796*eta2 + 260027.4121353132*eta3 - 1.723865744472182e6*eta4 + 5.858455766230802e6*eta5 - 7.756341721552802e6*eta6)*S2 - 0.027167813927224657*(34.281932237450256 - 3312.7658728016568*eta1 + 84126.14531363266*eta2 - 956052.0170024392*eta3 + 5.570748509263883e6*eta4 - 1.6270212243584689e7*eta5 + 1.8855858173287075e7*eta6)*S3)*sqroot;
490 break;
491 }
492 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Insp_Amp_44_iv2: version %i is not valid.", InspAmpFlag);}
493 }
494 return total;
495}
496
497static double IMRPhenomXHM_Insp_Amp_44_iv3(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag) {
498 double total=0;
499 switch (InspAmpFlag){
500 case 122018:{
501 double eta = pWF->eta;
502 double S = pWF->chiPNHat;
503 double delta=sqrt(1.-4.*eta),eta2,eta3,eta4,eta5,S2;
504 eta2 = pow(eta,2);
505 eta3 = pow(eta,3);
506 eta4 = pow(eta,4);
507 eta5 = pow(eta,5);
508 S2 = pow(S,2);
509 double noSpin = 0.08212436946985402 - 0.025332770704783136*eta - 3.2466088293309885*eta2 + 28.404235115663706*eta3 - 111.36325359782991*eta4 + 157.05954559045156*eta5;
510 double eqSpin = S*(0.03488890057062679 + 0.039491331923244756*S + eta*(-0.08968833480313292 - 0.12754920943544915*S - 0.11199012099701576*S2) + 0.034468577523793176*S2);
511 double uneqSpin = 0.2062291124580944*pWF->dchi*eta2*delta;
512 total = noSpin + eqSpin + uneqSpin;
513 break;
514 }
515 case 122022:{
516 double eta = pWF->eta;
517 double delta = pWF->delta;
518 double sqroot = sqrt(eta);
519 double S = pWF->chiPNHat;
520 double chidiff = pWF->dchi_half;
521 double eta1 = eta;
522 double eta2 = eta1 * eta1;
523 double eta3 = eta1 * eta2;
524 double eta4 = eta1 * eta3;
525 double eta5 = eta1 * eta4;
526 double eta6 = eta1 * eta5;
527 double S1 = S;
528 double S2 = S1 * S1;
529 double chidiff1 = chidiff;
530 double chidiff2 = chidiff1 * chidiff1;
531 total = (chidiff1*delta*(1.4739380748149558*eta1 + 0.06541707987699942*eta2 - 9.473290540936633*eta3) + chidiff2*(-0.3640838331639651*eta1 + 3.7369795937033756*eta2 - 8.709159662885131*eta3))*sqroot + (1.7335503724888923 + 12.656614578053683*eta1 - 139.6610487470118*eta2 + 456.78649322753824*eta3 - 599.2709938848282*eta4)*sqroot + chidiff1*delta*(2.3532739003216254*eta1 - 21.37216554136868*eta2 + 53.35003268489743*eta3)*S1*sqroot + S1*(-0.15782329022461472*(6.0309399412954345 - 229.16361598098678*eta1 + 3777.477006415653*eta2 - 31109.307191210424*eta3 + 139319.8239886073*eta4 - 324891.4001578353*eta5 + 307714.3954026392*eta6) - 0.03050157254864058*(4.232861441291087 + 1609.4251694451375*eta1 - 51213.27604422822*eta2 + 612317.1751155312*eta3 - 3.5589766538499263e6*eta4 + 1.0147654212772278e7*eta5 - 1.138861230369246e7*eta6)*S1 - 0.026407497690308382*(-17.184685557542196 + 744.4743953122965*eta1 - 10494.512487701073*eta2 + 66150.52694069289*eta3 - 184787.79377504133*eta4 + 148102.4257785174*eta5 + 128167.89151782403*eta6)*S2)*sqroot;
532 break;
533 }
534 default:{XLAL_ERROR_REAL8(XLAL_EINVAL,"Error in IMRPhenomXHM_Insp_Amp_44_iv3: version %i is not valid.", InspAmpFlag);}
535 }
536 return total;
537}
538
539/* End of Amp Parameter Space Fits */
540
541
542/***********************************************/
543/* Pseudo PN Amplitude Coefficients */
544/***********************************************/
545
546/* The whole inspiral ansatz is given by
547
548 PNAnsatz(f) + pseudo-PN(f), with pseudo-PN(f) = rho1 *(f/fcutInsp)^(7/3) + rho2(f/fcutInsp)^(8/3) + rho3(f/fcutInsp)^(9/3).
549
550 The coefficients are computed demanding that the pseudo-PN part at the three collocation points frequencies returns the actual collocation points: iv1, iv2, iv3.
551
552*/
553
554/* The input arguments for these rho coefficients are the value of the collocation points: v1, v2, v3, useful powers of the frequencies of these points f1, f2, f3 and the cutting frequency for the
555inspiral fcutInsp. The values v1, v2, v3 are given by the fit of the collocation points minus the value of the PNAnsatz at that frequency. */
556
557/* There are 3 functions, one for each pseudoPN coefficient rho1,2,3.
558Every function contains three possible results, depending on if you reconstruct with 3, 2 or 1 collocation points.
559With 2 colloc points rho3 = 0, and with 1 colloc point rho2=rho3=0 */
560
561/* Get Pseudo PN Amp coefficient rho1 at f^(7/3) */
562static double IMRPhenomXHM_Inspiral_Amp_rho1(double v1, double v2, double v3, IMRPhenomX_UsefulPowers *powers_of_fcutInsp, IMRPhenomX_UsefulPowers *powers_of_f1, IMRPhenomX_UsefulPowers *powers_of_f2, IMRPhenomX_UsefulPowers *powers_of_f3, IMRPhenomXHMWaveformStruct *pWFHM){
563 double retVal = 0.;
564 switch(pWFHM->IMRPhenomXHMInspiralAmpVersion){
565 case 3: // 3 inspiral collocation points
566 {
567 // default version
568 retVal = (powers_of_fcutInsp->seven_thirds*(-(powers_of_f1->three*powers_of_f3->eight_thirds*v2) + powers_of_f1->eight_thirds*powers_of_f3->three*v2 + powers_of_f2->three*(powers_of_f3->eight_thirds*v1 - powers_of_f1->eight_thirds*v3) + powers_of_f2->eight_thirds*(-(powers_of_f3->three*v1) + powers_of_f1->three*v3)))/(powers_of_f1->seven_thirds*(powers_of_f1->one_third - powers_of_f2->one_third)*powers_of_f2->seven_thirds*(powers_of_f1->one_third - powers_of_f3->one_third)*(powers_of_f2->one_third - powers_of_f3->one_third)*powers_of_f3->seven_thirds);
569 break;
570 }
571 case 2: // 2 inspiral collocation points
572 {
573 retVal=(powers_of_fcutInsp->seven_thirds*(-(powers_of_f2->eight_thirds*v1) + powers_of_f1->eight_thirds*v2))/(powers_of_f1->seven_thirds*(powers_of_f1->one_third - powers_of_f2->one_third)*powers_of_f2->seven_thirds);
574 break;
575 }
576 case 1: // 1 inspiral collocation points
577 {
578 retVal=(powers_of_fcutInsp->seven_thirds*v1)/(powers_of_f1->seven_thirds);
579 break;
580 }
581 case 0: // No collocation points, reconstruct with PN only
582 {
583 retVal=0;
584 break;
585 }
586 default: {
587 XLALPrintError("Error in IMRPhenomXHM_Inspiral_Amp_rho1: version is not valid. Versions available are 0,1,2,3.\n");
588 XLAL_ERROR(XLAL_EDOM, "IMRPhenomXHM_Inspiral_Amp_rho1 version is not valid. Aborting. Versions available are 0,1,2,3.\n");
589 }
590 }
591 return retVal;
592
593}
594
595/* Get Pseudo PN Amp coefficient rho2 at f^(8/3) */
596static double IMRPhenomXHM_Inspiral_Amp_rho2(double v1, double v2, double v3, IMRPhenomX_UsefulPowers *powers_of_fcutInsp, IMRPhenomX_UsefulPowers *powers_of_f1, IMRPhenomX_UsefulPowers *powers_of_f2, IMRPhenomX_UsefulPowers *powers_of_f3, IMRPhenomXHMWaveformStruct *pWFHM){
597
598 double retVal = 0;
599 switch(pWFHM->IMRPhenomXHMInspiralAmpVersion){
600 case 3: // 3 inspiral collocation points
601 {
602 // default version
603 retVal=(-(powers_of_fcutInsp->eight_thirds*(-(powers_of_f1->three*powers_of_f3->seven_thirds*v2) + powers_of_f1->seven_thirds*powers_of_f3->three*v2 + powers_of_f2->three*(powers_of_f3->seven_thirds*v1 - powers_of_f1->seven_thirds*v3) + powers_of_f2->seven_thirds*(-(powers_of_f3->three*v1) + powers_of_f1->three*v3))))/(powers_of_f1->seven_thirds*(powers_of_f1->one_third - powers_of_f2->one_third)*powers_of_f2->seven_thirds*(powers_of_f1->one_third - powers_of_f3->one_third)*(powers_of_f2->one_third - powers_of_f3->one_third)*powers_of_f3->seven_thirds);
604 break;
605 }
606 case 2: // 2 inspiral collocation points
607 {
608
609 retVal=(-(powers_of_fcutInsp->eight_thirds*(-(powers_of_f2->seven_thirds*v1) + powers_of_f1->seven_thirds*v2)))/(powers_of_f1->seven_thirds*(powers_of_f1->one_third - powers_of_f2->one_third)*powers_of_f2->seven_thirds);
610 break;
611 }
612 case 1: // 1 inspiral collocation points
613 {
614 retVal = 0;
615 break;
616 }
617 case 0: // No collocation points, reconstruct with PN only
618 {
619 retVal=0;
620 break;
621 }
622 default: {XLALPrintError("Error in IMRPhenomXHM_Inspiral_Amp_rho2: version is not valid. Versions avilable are 0,1,2,3.\n");XLAL_ERROR(XLAL_EDOM, "IMRPhenomXHM_Inspiral_Amp_rho2 version is not valid. Aborting. Versions available are 0,1,2,3.\n");}
623 }
624 return retVal;
625
626}
627
628/* Get Pseudo PN Amp coefficient rho3 at f^(9/3) */
629static double IMRPhenomXHM_Inspiral_Amp_rho3(double v1, double v2, double v3, IMRPhenomX_UsefulPowers *powers_of_fcutInsp, IMRPhenomX_UsefulPowers *powers_of_f1, IMRPhenomX_UsefulPowers *powers_of_f2, IMRPhenomX_UsefulPowers *powers_of_f3, IMRPhenomXHMWaveformStruct *pWFHM){
630
631 double retVal = 0;
632 switch(pWFHM->IMRPhenomXHMInspiralAmpVersion){
633 case 3: // 3 inspiral collocation points
634 {
635 // default version
636 retVal=(powers_of_fcutInsp->three*(powers_of_f1->seven_thirds*(-powers_of_f1->one_third + powers_of_f3->one_third)*powers_of_f3->seven_thirds*v2 + powers_of_f2->seven_thirds*(-(powers_of_f3->eight_thirds*v1) + powers_of_f1->eight_thirds*v3) + powers_of_f2->eight_thirds*(powers_of_f3->seven_thirds*v1 - powers_of_f1->seven_thirds*v3)))/(powers_of_f1->seven_thirds*(powers_of_f1->one_third - powers_of_f2->one_third)*powers_of_f2->seven_thirds*(powers_of_f1->one_third - powers_of_f3->one_third)*(powers_of_f2->one_third - powers_of_f3->one_third)*powers_of_f3->seven_thirds);
637 break;
638 }
639 case 2: // 2 inspiral collocation points
640 {
641 retVal = 0;
642 break;
643 }
644 case 1: // 1 inspiral collocation point
645 {
646 retVal = 0;
647 break;
648 }
649 case 0: // No collocation points, reconstruct with PN only
650 {
651 retVal=0;
652 break;
653 }
654 default: {XLALPrintError("Error in IMRPhenomXHM_Inspiral_Amp_rho3: version is not valid.\n");XLAL_ERROR(XLAL_EDOM, "IMRPhenomXHM_Inspiral_Amp_rho3 version is not valid. Aborting. Versions available: 0,1,2,3.\n");}
655 }
656 return retVal;
657}
658
659
660/************* INSPIRAL AMPLITUDE ANSATZ ******************/
661
662/*
663 The ansatz is built by performing the Stationary Phase Approximation of the Time-Domain Amplitude up to 3PN.
664 The result is rexpanded in power series up to 3PN, except for the 21 mode, which is better behaved without the reexpansion.
665*/
666
667//Return the Fourier Domain Post-Newtonian ansatz up to 3PN without the pseudoPN terms for a particular frequency
669
670 // The 21 mode is special, is not a power series
671 if(pWFHM->useFAmpPN==1){
672 return IMRPhenomXHM_Inspiral_PNAmp_21Ansatz(powers_of_Mf, pWFHM, pAmp);
673 }
674
675 //This returns the amplitude strain rescaled with the prefactor of the 22 mode: divided by sqrt(2*eta/3.)/pi^(1/6)*Mf^(-7/6)
676 double pnAmp;
677 pnAmp = cabs(pAmp->pnInitial
678 + powers_of_Mf->one_third * pAmp->pnOneThird
679 + powers_of_Mf->two_thirds * pAmp->pnTwoThirds
680 + powers_of_Mf->itself * pAmp->pnThreeThirds
681 + powers_of_Mf->four_thirds * pAmp->pnFourThirds
682 + powers_of_Mf->five_thirds * pAmp->pnFiveThirds
683 + powers_of_Mf->two * pAmp->pnSixThirds
684 );// * (pAmp->PNglobalfactor); // Added
685 if (pAmp->InspRescaleFactor == 0){
686 pnAmp *= pAmp->PNglobalfactor * powers_of_Mf->m_seven_sixths * pWFHM->ampNorm;
687 }
688 else if (pAmp->InspRescaleFactor == 1){
689 pnAmp *= pAmp->PNglobalfactor;
690 }
691
692 return pnAmp;
693}
694
695//Return the 21 mode Fourier Domain Post-Newtonian ansatz up to 3PN without the pseudoPN terms for a particular frequency
697
698 double complex CpnAmpTD; //This is not the real strain of the lm mode. It is the strain rescaled with the prefactor of the 22 mode: divided by sqrt(2*eta/3.)/pi^(1/6)
699 double pnAmp, XdotT4, x_to_m_one_four, two_to_m_one_sixths = 0.8908987181403393, three_to_m_one_second = 0.5773502691896257;
700 x_to_m_one_four = two_to_m_one_sixths * powers_of_lalpiHM.m_one_sixth * powers_of_Mf->m_one_sixth;
701 int InsAmpFlag = pWFHM->IMRPhenomXHMInspiralAmpFreqsVersion;
702 switch(InsAmpFlag)
703 {
704 case 122018:
705 {
706 // Complex time-domain Post-Newtonina amplitude, power series
707 CpnAmpTD = (
708 powers_of_Mf->one_third * pAmp->x05
709 + powers_of_Mf->two_thirds * pAmp->x1
710 + powers_of_Mf->itself * pAmp->x15
711 + powers_of_Mf->four_thirds * pAmp->x2
712 + powers_of_Mf->five_thirds * pAmp->x25
713 + powers_of_Mf->two * pAmp->x3
714 );
715 CpnAmpTD = CpnAmpTD*powers_of_Mf->two_thirds*pAmp->PNTDfactor;
716
717 // Value of the the derivative of the PN expansion parameter X given by TaylorT4
718 XdotT4 = powers_of_Mf->five_thirds * powers_of_Mf->five_thirds * pAmp->xdot5
719 + powers_of_Mf->four * pAmp->xdot6
720 + powers_of_Mf->eight_thirds*powers_of_Mf->five_thirds * pAmp->xdot65
721 + powers_of_Mf->seven_thirds*powers_of_Mf->seven_thirds * pAmp->xdot7
722 + powers_of_Mf->eight_thirds*powers_of_Mf->seven_thirds * pAmp->xdot75
723 + powers_of_Mf->eight_thirds * powers_of_Mf->eight_thirds * pAmp->xdot8
724 + (powers_of_Mf->log*2./3. + pAmp->log2pi_two_thirds )*powers_of_Mf->eight_thirds * powers_of_Mf->eight_thirds * pAmp->xdot8Log
725 + powers_of_Mf->eight_thirds * powers_of_Mf->eight_thirds * powers_of_Mf->one_third * pAmp->xdot85;
726
727 // Perform the SPA, multiply time-domain by the phasing factor
728 pnAmp = 2. * powers_of_lalpiHM.sqrt * three_to_m_one_second * cabs(CpnAmpTD) * x_to_m_one_four / sqrt(XdotT4);///powers_of_Mf->m_seven_sixths/pWFHM->ampNorm;//Added last two
729 if (pAmp->InspRescaleFactor != 0){
730 pnAmp /= RescaleFactor(powers_of_Mf, pAmp, pAmp->InspRescaleFactor);
731 }
732
733 break;
734 }
735 default:{XLAL_ERROR_REAL8(XLAL_EINVAL, "Error in IMRPhenomXHM_Inspiral_PNAmp_Ansatz: IMRPhenomXInspiralAmpVersion is not valid. Recommended version is 122018.\n");}
736 }
737 return pnAmp;
738}
739
740
741//This is the complete Inspiral Amplitude Ansatz: PN ansatz + pseudoPN terms
743{
744 double InspAmp; //This is the amplitude strain rescaled with the prefactor of the 22 mode: divided by [sqrt(2*eta/3.)/pi^(1/6) * f^(-7/6)]
745 double pseudoterms = 0.;
746 if(pWFHM->IMRPhenomXHMInspiralAmpFreqsVersion == 122018){
747 INT2 tmp = pAmp->InspRescaleFactor;
748 pAmp->InspRescaleFactor = 1;
749 InspAmp = IMRPhenomXHM_Inspiral_PNAmp_Ansatz(powers_of_Mf, pWFHM, pAmp)
750 + powers_of_Mf->seven_thirds / pAmp->fcutInsp_seven_thirds * pAmp->rho1
751 + powers_of_Mf->eight_thirds / pAmp->fcutInsp_eight_thirds * pAmp->rho2
752 + powers_of_Mf->three / pAmp->fcutInsp_three * pAmp->rho3
753 ;
754 pAmp->InspRescaleFactor = tmp;
755 if (pAmp->InspRescaleFactor == 0){
756 InspAmp *= powers_of_Mf->m_seven_sixths * pWFHM->ampNorm;
757 }
758 }
759 else{
760 /* New release. FIXME: improve how the ansatz is built */
761 InspAmp = IMRPhenomXHM_Inspiral_PNAmp_Ansatz(powers_of_Mf, pWFHM, pAmp);
762 pseudoterms = powers_of_Mf->seven_thirds / pAmp->fcutInsp_seven_thirds * pAmp->InspiralCoefficient[0]
763 + powers_of_Mf->eight_thirds / pAmp->fcutInsp_eight_thirds * pAmp->InspiralCoefficient[1]
764 + powers_of_Mf->three / pAmp->fcutInsp_three * pAmp->InspiralCoefficient[2];
765 pseudoterms *= powers_of_Mf->m_seven_sixths * pAmp->PNdominant;
766 InspAmp += pseudoterms;
767 if (pAmp->InspRescaleFactor != 0){
768 InspAmp /= RescaleFactor(powers_of_Mf, pAmp, pAmp->InspRescaleFactor);
769 }
770 }
771
772 return InspAmp ;
773}
774
775/* Numerical derivative of the inspiral (4th order finite differences)
776 It is used for reconstructing the intermediate region */
778
779 double df = 10e-10;
780 double Nder;
781 double fun2R, funR, funL, fun2L;
782 double centralfreq = powers_of_Mf->itself;
783
784 IMRPhenomX_UsefulPowers powers_of_Mf2R, powers_of_MfR, powers_of_MfL, powers_of_Mf2L;
785
786 IMRPhenomX_Initialize_Powers(&powers_of_Mf2R, centralfreq + 2*df);
787 IMRPhenomX_Initialize_Powers(&powers_of_MfR, centralfreq + df);
788 IMRPhenomX_Initialize_Powers(&powers_of_MfL, centralfreq - df);
789 IMRPhenomX_Initialize_Powers(&powers_of_Mf2L, centralfreq - 2*df);
790
791 fun2R = IMRPhenomXHM_Inspiral_Amp_Ansatz(&powers_of_Mf2R, pWFHM, pAmp);
792 funR = IMRPhenomXHM_Inspiral_Amp_Ansatz(&powers_of_MfR, pWFHM, pAmp);
793 funL = IMRPhenomXHM_Inspiral_Amp_Ansatz(&powers_of_MfL, pWFHM, pAmp);
794 fun2L = IMRPhenomXHM_Inspiral_Amp_Ansatz(&powers_of_Mf2L, pWFHM, pAmp);
795 Nder = (-fun2R + 8*funR - 8*funL + fun2L )/(12*df);
796
797 return Nder;
798}
799
800/* VETO function: when extrapolating the model outside the calibration region the collocation points can be bad behaved.
801 With this function we select those that will be used in the reconstruction. */
803 double *iv1, double *iv2, double *iv3,
804 IMRPhenomX_UsefulPowers *powers_of_f1,
805 IMRPhenomX_UsefulPowers *powers_of_f2,
806 IMRPhenomX_UsefulPowers *powers_of_f3,
809)
810{
811 double threshold = 0.2/(pWFHM->ampNorm);
812 #if DEBUG == 1
813 printf("\n\nf1, f2, f3 = %.16f %.16f %.16f\n\n",powers_of_f1->itself,powers_of_f2->itself,powers_of_f3->itself);
814 printf("\n\nf1, f2, f3 = %.16f %.16f %.16f\n\n",threshold*powers_of_f1->seven_sixths,threshold*powers_of_f2->seven_sixths,threshold*powers_of_f3->seven_sixths);
815 printf("\nInspiral Veto: AmpVersion = %i",pWFHM->IMRPhenomXHMInspiralAmpVersion);
816 #endif
817 // Remove too low collocation points (heuristic).
818 if(pAmp->CollocationPointsValuesAmplitudeInsp[0] < threshold*powers_of_f1->seven_sixths){
819 *iv1 = 0;
821 }
822 if(pAmp->CollocationPointsValuesAmplitudeInsp[1] < threshold*powers_of_f2->seven_sixths){
823 *iv2 = 0;
825 }
826 if(pAmp->CollocationPointsValuesAmplitudeInsp[2] < threshold*powers_of_f3->seven_sixths){
827 *iv3 = 0;
829 }
830}
831
832// Check if the three collocation points are wavy
833int WavyPoints(double p1, double p2, double p3){
834 if((p1>p2 && p2<p3) || (p1<p2 && p2>p3)){
835 return 1;
836 }else{
837 return 0;
838 }
839}
840
842
843 // Get CollocationPointsFreqsAmplitudeInsp and CollocationPointsValuesAmplitudeInsp
844 IMRPhenomX_UsefulPowers *powers_of_Mf_inspcollpoints = (IMRPhenomX_UsefulPowers *)XLALMalloc((pWFHM->nCollocPtsInspAmp + 1) * sizeof(IMRPhenomX_UsefulPowers));
846
847 for(UINT2 i = 0; i < pWFHM->nCollocPtsInspAmp; i++){
848 int status = IMRPhenomX_Initialize_Powers(&(powers_of_Mf_inspcollpoints[i]), pAmp->CollocationPointsFreqsAmplitudeInsp[i]);
849 if(status != XLAL_SUCCESS)
850 XLALPrintError("IMRPhenomXHM_Get_Inspiral_Amp_Coefficients failed for Mf, initial_status=%d",status);
851 }
852 int status = IMRPhenomX_Initialize_Powers(&(powers_of_Mf_inspcollpoints[pWFHM->nCollocPtsInspAmp]), pAmp->fAmpMatchIN);
853 if(status != XLAL_SUCCESS)
854 XLALPrintError("IMRPhenomXHM_Get_Inspiral_Amp_Coefficients failed for Mf, initial_status=%d",status);
855 pAmp->fcutInsp_seven_thirds = powers_of_Mf_inspcollpoints[pWFHM->nCollocPtsInspAmp].seven_thirds;
856 pAmp->fcutInsp_eight_thirds = powers_of_Mf_inspcollpoints[pWFHM->nCollocPtsInspAmp].eight_thirds;
857 pAmp->fcutInsp_three = powers_of_Mf_inspcollpoints[pWFHM->nCollocPtsInspAmp].three;
858
859
860 // Do Vetos? set final pWFHM->IMRPhenomXHMInspiralAmpVersion
861
862 // Get PN values at collocation points frequencies
863 for (INT4 i=0; i < pWFHM->nCollocPtsInspAmp; i++){
864 pAmp->PNAmplitudeInsp[i] = IMRPhenomXHM_Inspiral_PNAmp_Ansatz(&(powers_of_Mf_inspcollpoints[i]), pWFHM, pAmp);
865 }
866
867 IMRPhenomXHM_Inspiral_Amp_Coefficients(pAmp, powers_of_Mf_inspcollpoints, pWFHM);
868
869 LALFree(powers_of_Mf_inspcollpoints);
870}
871
874 case 122022:
875 case 102021:{
877 pAmp->CollocationPointsFreqsAmplitudeInsp[1] = 0.75 * pAmp->fAmpMatchIN;
879 break;
880 }
881 // FIXME: Add cases for equispaced, Chebyshev
882 default: {XLAL_ERROR_VOID(XLAL_EDOM, "Error in IMRPhenomXHM_Inspiral_CollocationPoints: IMRPhenomXHMInspiralAmpFreqsVersion = %i is not valid. Recommneded version is 102021.\n", pWFHM->IMRPhenomXHMInspiralAmpFreqsVersion);}
883 }
884 for(UINT2 i = 0; i < pWFHM->nCollocPtsInspAmp; i++){
886 }
887
888 if (pWFHM->InspiralAmpVeto == 1){
890 }
891}
892
894
896 pAmp->InspiralCoefficient[i] = 0;
897 }
898
899 IMRPhenomX_UsefulPowers *powers_of_f1 = &(powers_of_Mf_inspcollpoints[0]);
900 IMRPhenomX_UsefulPowers *powers_of_f2 = &(powers_of_Mf_inspcollpoints[1]);
901 IMRPhenomX_UsefulPowers *powers_of_f3 = &(powers_of_Mf_inspcollpoints[2]);
902 IMRPhenomX_UsefulPowers *powers_of_finsp = &(powers_of_Mf_inspcollpoints[3]);
903
904 REAL8 v1 = (pAmp->CollocationPointsValuesAmplitudeInsp[0] - pAmp->PNAmplitudeInsp[0]) / pAmp->PNdominant / powers_of_f1->m_seven_sixths;
905 REAL8 v2 = (pAmp->CollocationPointsValuesAmplitudeInsp[1] - pAmp->PNAmplitudeInsp[1]) / pAmp->PNdominant / powers_of_f2->m_seven_sixths;
906 REAL8 v3 = (pAmp->CollocationPointsValuesAmplitudeInsp[2] - pAmp->PNAmplitudeInsp[2]) / pAmp->PNdominant / powers_of_f3->m_seven_sixths;
907
908 switch(pWFHM->IMRPhenomXHMInspiralAmpVersion){
909 case 1:{
910 pAmp->InspiralCoefficient[0] = (powers_of_finsp->seven_thirds*v1)/powers_of_f1->seven_thirds;
911 break;
912 }
913 case 2:{
914 pAmp->InspiralCoefficient[0] = (powers_of_finsp->seven_thirds*v2)/powers_of_f2->seven_thirds;
915 break;
916 }
917 case 3:{
918 pAmp->InspiralCoefficient[0] = (powers_of_finsp->seven_thirds*v3)/powers_of_f3->seven_thirds;
919 break;
920 }
921 case 12:{
922 pAmp->InspiralCoefficient[0] = (powers_of_finsp->seven_thirds*(-(powers_of_f2->eight_thirds*v1) + powers_of_f1->eight_thirds*v2))/(powers_of_f1->seven_thirds*(powers_of_f1->one_third - powers_of_f2->one_third)*powers_of_f2->seven_thirds);
923 pAmp->InspiralCoefficient[1] = (powers_of_finsp->eight_thirds*(powers_of_f1->m_seven_thirds*v1 - powers_of_f2->m_seven_thirds*v2))/(powers_of_f1->one_third - powers_of_f2->one_third);
924 break;
925 }
926 case 13:{
927 pAmp->InspiralCoefficient[0] = (powers_of_finsp->seven_thirds*(-(powers_of_f3->eight_thirds*v1) + powers_of_f1->eight_thirds*v3))/(powers_of_f1->seven_thirds*(powers_of_f1->one_third - powers_of_f3->one_third)*powers_of_f3->seven_thirds);
928 pAmp->InspiralCoefficient[1] = (powers_of_finsp->eight_thirds*(powers_of_f1->m_seven_thirds*v1 - powers_of_f3->m_seven_thirds*v3))/(powers_of_f1->one_third - powers_of_f3->one_third);
929 break;
930 }
931 case 23:{
932 pAmp->InspiralCoefficient[0] = (powers_of_finsp->seven_thirds*(-(powers_of_f3->eight_thirds*v2) + powers_of_f2->eight_thirds*v3))/(powers_of_f2->seven_thirds*(powers_of_f2->one_third - powers_of_f3->one_third)*powers_of_f3->seven_thirds);
933 pAmp->InspiralCoefficient[1] = (powers_of_finsp->eight_thirds*(powers_of_f1->m_seven_thirds*v1 - powers_of_f3->m_seven_thirds*v3))/(powers_of_f1->one_third - powers_of_f3->one_third);
934 break;
935 }
936 case 123:{
937 pAmp->InspiralCoefficient[0] = (powers_of_finsp->seven_thirds*(-(powers_of_f1->three*powers_of_f3->eight_thirds*v2) + powers_of_f1->eight_thirds*powers_of_f3->three*v2 + powers_of_f2->three*(powers_of_f3->eight_thirds*v1 - powers_of_f1->eight_thirds*v3) + powers_of_f2->eight_thirds*(-(powers_of_f3->three*v1) + powers_of_f1->three*v3)))/(powers_of_f1->seven_thirds*(powers_of_f1->one_third - powers_of_f2->one_third)*powers_of_f2->seven_thirds*(powers_of_f1->one_third - powers_of_f3->one_third)*(powers_of_f2->one_third - powers_of_f3->one_third)*powers_of_f3->seven_thirds);
938 pAmp->InspiralCoefficient[1] = (powers_of_finsp->eight_thirds*(powers_of_f1->three*powers_of_f3->seven_thirds*v2 - powers_of_f1->seven_thirds*powers_of_f3->three*v2 + powers_of_f2->three*(-(powers_of_f3->seven_thirds*v1) + powers_of_f1->seven_thirds*v3) + powers_of_f2->seven_thirds*(powers_of_f3->three*v1 - powers_of_f1->three*v3)))/(powers_of_f1->seven_thirds*(powers_of_f1->one_third - powers_of_f2->one_third)*powers_of_f2->seven_thirds*(powers_of_f1->one_third - powers_of_f3->one_third)*(powers_of_f2->one_third - powers_of_f3->one_third)*powers_of_f3->seven_thirds);
939 pAmp->InspiralCoefficient[2] = (powers_of_finsp->three*(powers_of_f1->seven_thirds*(-powers_of_f1->one_third + powers_of_f3->one_third)*powers_of_f3->seven_thirds*v2 + powers_of_f2->seven_thirds*(-(powers_of_f3->eight_thirds*v1) + powers_of_f1->eight_thirds*v3) + powers_of_f2->eight_thirds*(powers_of_f3->seven_thirds*v1 - powers_of_f1->seven_thirds*v3)))/(powers_of_f1->seven_thirds*(powers_of_f1->one_third - powers_of_f2->one_third)*powers_of_f2->seven_thirds*(powers_of_f1->one_third - powers_of_f3->one_third)*(powers_of_f2->one_third - powers_of_f3->one_third)*powers_of_f3->seven_thirds);
940 break;
941 }
942 }
943}
944/*************************************/
945/* */
946/* PHASE */
947/* */
948/*************************************/
949
950// Spin parameter S = (m1^2*chi1 + m2^2*chi2)/(m1^2 + m2^2)
951
952// Below we give paramater-space fits for the weighted difference between each mode's phase and the 22-phase: phi_lm-m/2 phi_22(2/m f), see Eqs. (4.10-4.12)
953
954/* Start of Phase Parameter Space Fits */
955
956static double IMRPhenomXHM_Insp_Phase_21_lambda(IMRPhenomXWaveformStruct *pWF, int InspPhaseFlag) {
957 double total;
958 switch (InspPhaseFlag){
959 case 122019:{
960 double eta = pWF->eta;
961 double S = pWF->STotR;
962 double eta2,eta3,eta4,eta5,S2;
963 eta2 = pow(eta,2);
964 eta3 = pow(eta,3);
965 eta4 = pow(eta,4);
966 eta5 = pow(eta,5);
967 S2 = pow(S,2);
968 double noSpin = 13.664473636545068 - 170.08866400251395*eta + 3535.657736681598*eta2 - 26847.690494515424*eta3 + 96463.68163125668*eta4 - 133820.89317471132*eta5;
969 double eqSpin = (S*(18.52571430563905 - 41.55066592130464*S + eta3*(83493.24265292779 + 16501.749243703132*S - 149700.4915210766*S2) + eta*(3642.5891077598003 + 1198.4163078715173*S - 6961.484805326852*S2) + 33.8697137964237*S2 + eta2*(-35031.361998480075 - 7233.191207000735*S + 62149.00902591944*S2)))/(6.880288191574696 + 1.*S);
970 double uneqSpin = -134.27742343186577*pWF->dchi*sqrt(1.-4.*eta)*eta2;
971 total = noSpin + eqSpin + uneqSpin;
972 break;
973 }
974 default:{XLAL_ERROR(XLAL_EDOM, "Error in IMRPhenomXHM_Insp_Phase_21_lambda: version is not valid. Recommended version is 122019.");}
975 }
976 return total;
977}
978
979static double IMRPhenomXHM_Insp_Phase_33_lambda(IMRPhenomXWaveformStruct *pWF, int InspPhaseFlag) {
980 double total=0;
981 switch (InspPhaseFlag){
982 case 122019:{
983 double eta = pWF->eta;
984 double S = pWF->STotR;
985 double delta=sqrt(1.-4.*eta);
986 double eta2,eta3;
987 eta2 = pow(eta,2);
988 eta3 = eta2*eta;
989 double noSpin = 4.1138398568400705 + 9.772510519809892*eta - 103.92956504520747*eta2 + 242.3428625556764*eta3;
990 double eqSpin = ((-0.13253553909611435 + 26.644159828590055*eta - 105.09339163109497*eta2)*S)/(1. + 0.11322426762297967*S);
991 double uneqSpin = -19.705359163581168*pWF->dchi*eta2*delta;
992 total = noSpin + eqSpin + uneqSpin;
993 break;
994 }
995 default:{XLAL_ERROR(XLAL_EDOM, "Error in IMRPhenomXHM_Insp_Phase_32_lambda: version is not valid. Recommended version is 122019.");}
996 }
997 return total;
998}
999
1000static double IMRPhenomXHM_Insp_Phase_32_lambda(IMRPhenomXWaveformStruct *pWF, int InspPhaseFlag) {
1001 double total;
1002 switch (InspPhaseFlag){
1003 case 122019:{
1004 double eta = pWF->eta;
1005 double S = pWF->STotR;
1006 double eta2,eta3,eta4,S2,S3,S4;
1007 eta2 = pow(eta,2);
1008 eta3 = pow(eta,3);
1009 eta4 = pow(eta,4);
1010 S2 = pow(S,2);
1011 S3 = pow(S,3);
1012 S4 = pow(S,4);
1013 double noSpin = (9.913819875501506 + 18.424900617803107*eta - 574.8672384388947*eta2 + 2671.7813055097877*eta3 - 6244.001932443913*eta4)/(1. - 0.9103118343073325*eta);
1014 double eqSpin = (-4.367632806613781 + 245.06757304950986*eta - 2233.9319708029775*eta2 + 5894.355429022858*eta3)*S + (-1.375112297530783 - 1876.760129419146*eta + 17608.172965575013*eta2 - 40928.07304790013*eta3)*S2 + (-1.28324755577382 - 138.36970336658558*eta + 708.1455154504333*eta2 - 273.23750933544176*eta3)*S3 + (1.8403161863444328 + 2009.7361967331492*eta - 18636.271414571278*eta2 + 42379.205045791656*eta3)*S4;
1015 double uneqSpin = pWF->dchi*sqrt(1.-4.*eta)*eta2*(-105.34550407768225 - 1566.1242344157668*pWF->chi1L*eta + 1566.1242344157668*pWF->chi2L*eta + 2155.472229664981*eta*S);
1016 total = noSpin + eqSpin + uneqSpin;
1017 break;
1018 }
1019 default:{XLAL_ERROR(XLAL_EDOM, "Error in IMRPhenomXHM_Insp_Phase_32_lambda: version is not valid. Recommended version is 122019.");}
1020 }
1021 return total;
1022}
1023
1024
1025static double IMRPhenomXHM_Insp_Phase_44_lambda(IMRPhenomXWaveformStruct *pWF, int InspPhaseFlag) {
1026 double total;
1027 switch (InspPhaseFlag){
1028 case 122019:{
1029 double eta = pWF->eta;
1030 double S = pWF->STotR;
1031 double eta2,eta3,eta4,eta5,S2;
1032 eta2 = pow(eta,2);
1033 eta3 = pow(eta,3);
1034 eta4 = pow(eta,4);
1035 eta5 = pow(eta,5);
1036 S2 = pow(S,2);
1037 double noSpin = 5.254484747463392 - 21.277760168559862*eta + 160.43721442910618*eta2 - 1162.954360723399*eta3 + 1685.5912722190276*eta4 - 1538.6661348106031*eta5;
1038 double eqSpin = (0.007067861615983771 - 10.945895160727437*eta + 246.8787141453734*eta2 - 810.7773268493444*eta3)*S + (0.17447830920234977 + 4.530539154777984*eta - 176.4987316167203*eta2 + 621.6920322846844*eta3)*S2;
1039 double uneqSpin = -8.384066369867833*pWF->dchi*sqrt(1.-4.*eta)*eta2;
1040 total = noSpin + eqSpin + uneqSpin;
1041 break;
1042 }
1043 default:{XLAL_ERROR(XLAL_EDOM, "Error in IMRPhenomXHM_Insp_Phase_44_lambda: version is not valid. Recommended version is 122019.");}
1044 }
1045 return total;
1046}
1047
1048/* End of Phase Parameter Space Fits */
1049
1050// Returns linear-in-f term coming from the complex part of each mode's amplitude that will be added to the orbital phase, this is an expansion of Eq. (4.9) truncated at linear order
1051static double IMRPhenomXHM_Insp_Phase_LambdaPN(double eta, int modeInt){
1052
1053 double output;
1054 const double PI = powers_of_lalpiHM.itself;
1055
1056 switch(modeInt){
1057 case 21:
1058 {//2 f \[Pi] (-(1/2) - Log[16]/2)
1059 output=(2.*PI*(-0.5-2.*log(2)));
1060 break;
1061 }
1062 case 33:
1063 {//2/3 f \[Pi] (-(21/5) + 6 Log[3/2])
1064 output=2./3.*PI*(-21./5.+6.*log(1.5));
1065 break;
1066 }
1067 case 32:
1068 { //-((2376 f \[Pi] (-5 + 22 \[Eta]))/(-3960 + 11880 \[Eta]))
1069 output=-((2376.*PI*(-5.+22.*eta))/(-3960.+11880*eta));
1070 break;
1071 }
1072 case 44:
1073 { //(45045 f \[powers_of_lalpiHM.itself] (336 - 1193 \[Eta] +320 (-1 + 3 \[Eta]) Log[2]))/(2 (-1801800 + 5405400 \[Eta]))
1074 output=45045.*PI*(336.-1193.*eta+320.*(-1.+3.*eta)*log(2))/(2.*(-1801800. + 5405400.*eta));
1075 break;
1076 }
1077 default: output=0.;
1078 }
1079
1080 return -1.*output;
1081}
1082
1083
1084/************* INSPIRAL PHASE ANSATZ *************/
1085
1086/*
1087 The function loads the rescaled phenX coefficients for the phase and correct the result with a linear-in-f term coming from the complex part of the PN amplitude
1088
1089 The rescaling of the 22-coefficients is explained in App. D
1090*/
1092{
1093 //compute the orbital phase by laoding the rescaled phenX coefficients
1094 double philm=0.;
1095 double freqs[N_MAX_COEFFICIENTS_PHASE_INS]={powers_of_Mf->m_five_thirds,powers_of_Mf->m_four_thirds,powers_of_Mf->m_one,powers_of_Mf->m_two_thirds,powers_of_Mf->m_one_third,1.,powers_of_Mf->one_third,powers_of_Mf->two_thirds, Mf, powers_of_Mf->four_thirds, powers_of_Mf->five_thirds, powers_of_Mf->two,powers_of_Mf->seven_thirds};
1096 double logMf=powers_of_Mf->log;
1097 for(int i=0; i<N_MAX_COEFFICIENTS_PHASE_INS; i++)
1098 philm+=(pPhase->phi[i]+pPhase->phiL[i]*(logMf))*freqs[i];
1099 return philm;
1100}
1101
1103{
1104 //compute the orbital phase by laoding the rescaled phenX coefficients
1105 double dphilm=0.;
1106 double coeffs[N_MAX_COEFFICIENTS_PHASE_INS]={-5./3,-4./3,-1.,-2./3,-1./3,0.,1./3, 2./3, 1., 4./3, 5./3, 2., 7./3};
1107 double freqs[N_MAX_COEFFICIENTS_PHASE_INS]={powers_of_Mf->m_eight_thirds,powers_of_Mf->m_seven_thirds,powers_of_Mf->m_two,powers_of_Mf->m_five_thirds,powers_of_Mf->m_four_thirds,powers_of_Mf->m_one,powers_of_Mf->m_two_thirds,powers_of_Mf->m_one_third,1.,powers_of_Mf->one_third, powers_of_Mf->two_thirds, Mf,powers_of_Mf->four_thirds};
1108 double logMf=powers_of_Mf->log;
1109
1110 for(int i=0; i<N_MAX_COEFFICIENTS_PHASE_INS; i++)
1111 dphilm+=((pPhase->phi[i]+pPhase->phiL[i]*(logMf))*coeffs[i]+pPhase->phiL[i])*freqs[i];
1112 return dphilm;
1113}
#define LALFree(p)
static double double delta
int IMRPhenomX_Initialize_Powers(IMRPhenomX_UsefulPowers *p, REAL8 number)
IMRPhenomX_UsefulPowers powers_of_lalpiHM
void IMRPhenomXHM_Get_Inspiral_Amp_Coefficients(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
static double IMRPhenomXHM_Inspiral_Phase_AnsatzInt(double Mf, IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMPhaseCoefficients *pPhase)
static double IMRPhenomXHM_Insp_Amp_21_iv2(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Inspiral_Amp_rho3(double v1, double v2, double v3, IMRPhenomX_UsefulPowers *powers_of_fcutInsp, IMRPhenomX_UsefulPowers *powers_of_f1, IMRPhenomX_UsefulPowers *powers_of_f2, IMRPhenomX_UsefulPowers *powers_of_f3, IMRPhenomXHMWaveformStruct *pWFHM)
static double IMRPhenomXHM_Inspiral_Amp_NDAnsatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
static double IMRPhenomXHM_Inspiral_Amp_Ansatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
static double IMRPhenomXHM_Insp_Amp_44_iv2(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Inspiral_Amp_rho2(double v1, double v2, double v3, IMRPhenomX_UsefulPowers *powers_of_fcutInsp, IMRPhenomX_UsefulPowers *powers_of_f1, IMRPhenomX_UsefulPowers *powers_of_f2, IMRPhenomX_UsefulPowers *powers_of_f3, IMRPhenomXHMWaveformStruct *pWFHM)
static double IMRPhenomXHM_Insp_Amp_33_iv1(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Insp_Phase_33_lambda(IMRPhenomXWaveformStruct *pWF, int InspPhaseFlag)
static double IMRPhenomXHM_Insp_Amp_33_iv2(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Insp_Amp_32_iv3(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static void IMRPhenomXHM_Inspiral_Amp_Coefficients(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomX_UsefulPowers *powers_of_Mf_inspcollpoints, IMRPhenomXHMWaveformStruct *pWFHM)
static double IMRPhenomXHM_Insp_Phase_21_lambda(IMRPhenomXWaveformStruct *pWF, int InspPhaseFlag)
static double IMRPhenomXHM_Insp_Phase_44_lambda(IMRPhenomXWaveformStruct *pWF, int InspPhaseFlag)
static double IMRPhenomXHM_Inspiral_PNAmp_21Ansatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
static double IMRPhenomXHM_Insp_Amp_32_iv1(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Insp_Phase_32_lambda(IMRPhenomXWaveformStruct *pWF, int InspPhaseFlag)
static void IMRPhenomXHM_Inspiral_Amp_CollocationPoints(IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXWaveformStruct *pWF22)
static double IMRPhenomXHM_Insp_Amp_33_iv3(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Insp_Amp_44_iv1(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
void IMRPhenomXHM_Inspiral_Amplitude_Veto(double *iv1, double *iv2, double *iv3, IMRPhenomX_UsefulPowers *powers_of_f1, IMRPhenomX_UsefulPowers *powers_of_f2, IMRPhenomX_UsefulPowers *powers_of_f3, IMRPhenomXHMAmpCoefficients *pAmp, IMRPhenomXHMWaveformStruct *pWFHM)
static double IMRPhenomXHM_Insp_Amp_21_iv3(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
int WavyPoints(double p1, double p2, double p3)
static double IMRPhenomXHM_Insp_Amp_44_iv3(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Inspiral_PNAmp_Ansatz(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMWaveformStruct *pWFHM, IMRPhenomXHMAmpCoefficients *pAmp)
static double IMRPhenomXHM_Insp_Amp_21_iv1(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
static double IMRPhenomXHM_Inspiral_Phase_Ansatz(double Mf, IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMPhaseCoefficients *pPhase)
static double IMRPhenomXHM_Inspiral_Amp_rho1(double v1, double v2, double v3, IMRPhenomX_UsefulPowers *powers_of_fcutInsp, IMRPhenomX_UsefulPowers *powers_of_f1, IMRPhenomX_UsefulPowers *powers_of_f2, IMRPhenomX_UsefulPowers *powers_of_f3, IMRPhenomXHMWaveformStruct *pWFHM)
static double IMRPhenomXHM_Insp_Phase_LambdaPN(double eta, int modeInt)
static double IMRPhenomXHM_Insp_Amp_32_iv2(IMRPhenomXWaveformStruct *pWF, int InspAmpFlag)
double RescaleFactor(IMRPhenomX_UsefulPowers *powers_of_Mf, IMRPhenomXHMAmpCoefficients *pAmp, UINT2 rescalefactor)
#define N_MAX_COEFFICIENTS_PHASE_INS
#define N_MAX_COEFFICIENTS_AMPLITUDE_INS
double i
Definition: bh_ringdown.c:118
double e
Definition: bh_ringdown.c:117
double REAL8
int16_t INT2
uint16_t UINT2
int32_t INT4
void * XLALMalloc(size_t n)
#define XLAL_ERROR_VOID(...)
#define XLAL_ERROR_REAL8(...)
#define XLAL_ERROR(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_SUCCESS
XLAL_EDOM
XLAL_EINVAL
string status
REAL8 CollocationPointsFreqsAmplitudeInsp[N_MAX_COEFFICIENTS_AMPLITUDE_INS]
REAL8 InspiralCoefficient[N_MAX_COEFFICIENTS_AMPLITUDE_INS]
REAL8 PNAmplitudeInsp[N_MAX_COEFFICIENTS_AMPLITUDE_INS]
REAL8 CollocationPointsValuesAmplitudeInsp[N_MAX_COEFFICIENTS_AMPLITUDE_INS]
ParameterSpaceFit InspiralAmpFits[N_HIGHERMODES_IMPLEMENTED *N_MAX_COEFFICIENTS_AMPLITUDE_INS]
REAL8 phiL[N_MAX_COEFFICIENTS_PHASE_INS]
REAL8 phi[N_MAX_COEFFICIENTS_PHASE_INS]
char output[FILENAME_MAX]
Definition: unicorn.c:26