LALSimulation  5.4.0.1-fe68b98
LALSimIMRSpinEOBFactorizedWaveformCoefficientsPrec.c
Go to the documentation of this file.
1 /*
2  * \author Stas Babak, Prayush Kumar, Andrea Taracchini
3  */
4 #ifndef _LALSIMIMRSPINPRECEOBFACTORIZEDWAVEFORMCOEFFICIENTS_C
5 #define _LALSIMIMRSPINPRECEOBFACTORIZEDWAVEFORMCOEFFICIENTS_C
6 
7 #define v4Pwave 451
8 
9 
10 #include <complex.h>
11 #include <lal/LALSimInspiral.h>
12 //#include <lal/LALSimIMR.h>
13 
14 #include "LALSimIMREOBNRv2.h"
15 #include "LALSimIMRSpinEOB.h"
16 
17 /**
18  *
19  * Waveform expressions are given by
20  * Taracchini et al. ( PRD 86, 024011 (2012), arXiv 1202.0790 ).
21  * All equation numbers in this file refer to equations of this paper,
22  * unless otherwise specified.
23  * Coefficients of the so-called "deltalm" terms are given by
24  * Damour et al. PRD 79, 064004 (2009) [arXiv:0811.2069] and Pan et al. PRD 83, 064003 (2011) [arXiv:1006.0431],
25  * henceforth DIN and PBFRT.
26  *
27  * Functions to compute the factorized waveform for the purpose of computing the
28  * flux, i.e. returning only the absolute value of the multipoles. The tail term
29  * Tlm is used in its resummed form, given by Eq. (42) of Damour, Nagar and
30  * Bernuzzi, PRD 87 084035 (2013) [arxiv:1212.4357], called DNB here.
31  *
32  */
33 
34 
35 /*--------------------------------------------------------------*/
36 /**
37  * Spin Factors
38  */
39 
40 static int
41 XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients(
42  FacWaveformCoeffs * const coeffs, /**< OUTPUT, pre-computed waveform coefficients */
43  const REAL8 m1, /**< mass 1 */
44  const REAL8 m2, /**< mass 2 */
45  const REAL8 eta, /**< symmetric mass ratio */
46  const REAL8 a, /**< Kerr spin parameter for test-particle terms */
47  const REAL8 chiS, /**< (chi1+chi2)/2 */
48  const REAL8 chiA, /**< (chi1-chi2)/2 */
49  const UINT4 SpinAlignedEOBversion /**< 1 for SEOBNRv1; 2 for SEOBNRv2; 3 for a small special treatment for SEOBNRv3*/
50 );
51 
52 
53 /****************************************************
54  * Definition of Functions
55  * **********************************************/
56 
57 
58 /**
59  * This function calculates coefficients for hlm mode factorized-resummed waveform.
60  * The coefficients are pre-computed and stored in the SpinEOBParams structure.
61  * Eq. 17 and the entire Appendix of PRD 86, 024011 (2012) + changes
62  * described in the section "Factorized waveforms" of https://dcc.ligo.org/T1400476
63  */
64 
65 static int
66 XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients(
67  FacWaveformCoeffs * const coeffs, /**< OUTPUT, pre-computed waveform coefficients */
68  const REAL8 m1, /**< mass 1 */
69  const REAL8 m2, /**< mass 2 */
70  const REAL8 eta, /**< symmetric mass ratio */
71  const REAL8 tmpa, /**< Kerr spin parameter for test-particle terms */
72  const REAL8 chiS, /**< (chi1+chi2)/2 */
73  const REAL8 chiA, /**< (chi1-chi2)/2 */
74  UINT4 SpinAlignedEOBversion /**< 1 for SEOBNRv1; 2 for SEOBNRv2; 4 for the coefficients
75  in the flux of v4P and v4Pwave for the coefficients in the waveform of v4P */
76  )
77 {
78  int debugPK = 0;
79  if (debugPK) {
80  XLAL_PRINT_INFO("In XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients: Renewing hlm coefficients.\n");
81  XLAL_PRINT_INFO("PK:: chiS = %.12e, chiA = %.12e, a = %.12e (UNUSED), EOBVERSION = %d\n",
82  chiS, chiA, tmpa, SpinAlignedEOBversion);
83  }
84  REAL8 a = tmpa;
85  INT4 waveform = 0;
86  REAL8 eta2 = eta * eta;
87  REAL8 eta3 = eta2 * eta;
88 
89  if (SpinAlignedEOBversion == v4Pwave) {
90  SpinAlignedEOBversion = 4;
91  waveform = 1;
92  }
93  //RC: Here I am taking care of the fact that the flux used to calibrate SEOBNRv4 is missing some terms
94  //in the 21 mode that I have added later in SEOBNRv4HM. For this reason when computing the Coefficients
95  //for the 21 mode one has to specify if they are needed for the flux or for the waveform.
96 
97  REAL8 dM , dM2, chiA2, chiS2;
98  //dM3;
99  REAL8 aDelta , a2, a3;
100 
101  /* Combination which appears a lot */
102  REAL8 m1Plus3eta, m1Plus3eta2;
103  REAL8 UNUSED m1Plus3eta3;
104 
105  dM2 = 1. - 4. * eta;
106  chiA2 = chiA*chiA;
107  chiS2 = chiS*chiS;
108  REAL8 chiS3 = chiS2*chiS;
109  REAL8 chiA3 = chiA2*chiA;
110 
111 
112  //XLAL_PRINT_INFO("****************************** a = %e *********************************\n", a);
113 
114  /* Check that deltaM has a reasonable value */
115  if (dM2 < 0) {
116  XLALPrintError("XLAL Error - %s: eta seems to be > 0.25 - this isn't allowed!\n", __func__);
117  XLAL_ERROR(XLAL_EINVAL);
118  }
119  dM = sqrt(dM2);
120  if (m1 < m2) {
121  dM = -dM;
122  }
123  //dM3 = dM2 * dM;
124 
125  aDelta = 0.;
126  //a value in delta_lm is 0 in both SEOBNRv1 and SEOBNRv2
127  a2 = a * a;
128  a3 = a2 * a;
129 
130  m1Plus3eta = -1. + 3. * eta;
131  m1Plus3eta2 = m1Plus3eta * m1Plus3eta;
132  m1Plus3eta3 = m1Plus3eta * m1Plus3eta2;
133 
134  /* Initialize all coefficients to zero */
135  /* This is important, as we will not set some if dM is zero */
136  memset(coeffs, 0, sizeof(FacWaveformCoeffs));
137 
138 
139  /*
140  * l = 2, Eqs. A8a and A8b for rho, Eq. A15a for f, Eqs. 20 and 21 of
141  * DIN and Eqs. 27a and 27b of PBFRT for delta as well as eqns 28-29 of PBFRT
142  */
143 
144  coeffs->delta22vh3 = 7. / 3.;
145  coeffs->delta22vh6 = (-4. * aDelta) / 3. + (428. * LAL_PI) / 105.;
146  if (SpinAlignedEOBversion == 4)
147  {
148  coeffs->delta22vh6 =
149  -4. / 3. * (dM * chiA + chiS * (1 - 2 * eta)) +
150  (428. * LAL_PI) / 105.;
151  }
152  coeffs->delta22v8 = (20. * aDelta) / 63.;
153  coeffs->delta22vh9 = -2203. / 81. + (1712. * LAL_PI * LAL_PI) / 315.;
154  coeffs->delta22v5 = -24. * eta;
155  coeffs->delta22v6 = 0.0;
156  if (SpinAlignedEOBversion == 2 && chiS + chiA * dM / (1. - 2. * eta) > 0.) {
157  coeffs->delta22v6 = -540. * eta * (chiS + chiA * dM / (1. - 2. * eta));
158 // double chi = (chiS + chiA * dM / (1. - 2. * eta));
159 // coeffs->delta22v6 = eta*(1./4.*(1. - 1080.*chi - sqrt((1. - 1080.*chi)*(1. - 1080*chi) + 8.*(13.5 +270.*chi +13.5*chi*chi))));
160  }
161  if (SpinAlignedEOBversion == 3 /*&& chiS + chiA * dM / (1. - 2. * eta) > 0.*/) {
162 // coeffs->delta22v6 = -540. * eta * (chiS + chiA * dM / (1. - 2. * eta));
163  double chi = (chiS + chiA * dM / (1. - 2. * eta));
164  coeffs->delta22v6 = eta*(1./4.*(1. - 1080.*chi - sqrt((1. - 1080.*chi)*(1. - 1080*chi) + 8.*(13.5 +270.*chi +13.5*chi*chi))));
165  }
166  coeffs->rho22v2 = -43. / 42. + (55. * eta) / 84.;
167  coeffs->rho22v3 = (-2. * (chiS + chiA * dM - chiS * eta)) / 3.;
168  switch (SpinAlignedEOBversion) {
169  case 1:
170  coeffs->rho22v4 = -20555. / 10584. + 0.5 * a2
171  - (33025. * eta) / 21168. + (19583. * eta2) / 42336.;
172  break;
173  case 2:
174  coeffs->rho22v4 = -20555. / 10584. + 0.5 * (chiS + chiA * dM) * (chiS + chiA * dM)
175  - (33025. * eta) / 21168. + (19583. * eta2) / 42336.;
176  break;
177  case 3:
178  coeffs->rho22v4 = -20555. / 10584. + 0.5 * (chiS + chiA * dM) * (chiS + chiA * dM)
179  - (33025. * eta) / 21168. + (19583. * eta2) / 42336.;
180  break;
181  case 4:
182  coeffs->rho22v4 =
183  -20555. / 10584. + 0.5 * (chiS + chiA * dM) * (chiS + chiA * dM) -
184  (33025. * eta) / 21168. + (19583. * eta2) / 42336.;
185  break;
186  default:
187  XLALPrintError("XLAL Error - %s: wrong SpinAlignedEOBversion value, must be 1, 2 or 3!\n", __func__);
188  XLAL_ERROR(XLAL_EINVAL);
189  break;
190  }
191  coeffs->rho22v5 = (-34. * a) / 21.;
192 if (SpinAlignedEOBversion == 4)
193  {
194  coeffs->rho22v5 =
195  (-34. / 21. + 49. * eta / 18. + 209. * eta2 / 126.) * chiS +
196  (-34. / 21. - 19. * eta / 42.) * dM * chiA;
197  }
198  coeffs->rho22v6 = 1556919113. / 122245200. + (89. * a2) / 252. - (48993925. * eta) / 9779616.
199  - (6292061. * eta2) / 3259872. + (10620745. * eta3) / 39118464.
200  + (41. * eta * LAL_PI * LAL_PI) / 192.;
201  coeffs->rho22v6l = -428. / 105.;
202  coeffs->rho22v7 = (18733. * a) / 15876. + a * a2 / 3.;
203  if (SpinAlignedEOBversion == 4)
204  {
205  coeffs->rho22v7 =
206  a3 / 3. + chiA * dM * (18733. / 15876. + (50140. * eta) / 3969. +
207  (97865. * eta2) / 63504.) +
208  chiS * (18733. / 15876. + (74749. * eta) / 5292. -
209  (245717. * eta2) / 63504. + (50803. * eta3) / 63504.);
210  }
211  switch (SpinAlignedEOBversion) {
212  case 1:
213  coeffs->rho22v8 = eta * (-5.6 - 117.6 * eta) - 387216563023. / 160190110080. + (18353. * a2) / 21168. - a2 * a2 / 8.;
214  break;
215  case 2:
216  coeffs->rho22v8 = -387216563023. / 160190110080. + (18353. * a2) / 21168. - a2 * a2 / 8.;
217  break;
218  case 3:
219  coeffs->rho22v8 = -387216563023. / 160190110080. + (18353. * a2) / 21168. - a2 * a2 / 8.;
220  break;
221  case 4:
222  coeffs->rho22v8 =
223  -387216563023. / 160190110080. + (18353. * a2) / 21168. -
224  a2 * a2 / 8.;
225  break;
226  default:
227  XLALPrintError("XLAL Error - %s: wrong SpinAlignedEOBversion value, must be 1, 2 or 3!\n", __func__);
228  XLAL_ERROR(XLAL_EINVAL);
229  break;
230  }
231  coeffs->rho22v8l = 9202. / 2205.;
232  coeffs->rho22v10 = -16094530514677. / 533967033600.;
233  coeffs->rho22v10l = 439877. / 55566.;
234 
235  if (debugPK) {
236  XLAL_PRINT_INFO("\nPK:: dM, eta, chiS, chiA while renewing hlm coeffs: %e, %e, %e, %e\n",
237  dM, eta, chiS, chiA);
238  XLAL_PRINT_INFO("PK:: Renewed rho-lm coeffs: v2 = %.16e, v3 = %.16e, v4 = %.16e, v5 = %.16e\nv6 = %.16e, v6l = %.16e v7 = %.16e v8 = %.16e, v8l = %.16e v10 = %.16e v10l = %.16e\n",
239  coeffs->rho22v2, coeffs->rho22v3, coeffs->rho22v4, coeffs->rho22v5,
240  coeffs->rho22v6, coeffs->rho22v6l, coeffs->rho22v7, coeffs->rho22v8,
241  coeffs->rho22v8l, coeffs->rho22v10, coeffs->rho22v10l);
242  }
243  if(waveform == 1){
244  a = 0;
245  a2 = 0;
246  a3 = 0;
247  }
248  coeffs->delta21vh3 = 2. / 3.;
249  coeffs->delta21vh6 = (-17. * aDelta) / 35. + (107. * LAL_PI) / 105.;
250  coeffs->delta21vh7 = (3. * aDelta * aDelta) / 140.;
251  coeffs->delta21vh9 = -272. / 81. + (214. * LAL_PI * LAL_PI) / 315.;
252  coeffs->delta21v5 = -493. * eta / 42.;
253  coeffs->delta21v7 = 0.0;
254  if (dM2) {
255 
256  //coeffs->rho21v1 = (-3. * (chiS + chiA / dM)) / (4.);
257  coeffs->rho21v1 = 0.0;
258  //coeffs->rho21v2 = -59. / 56 - (9. * chiAPlusChiSdM * chiAPlusChiSdM) / (32. * dM2) + (23. * eta) / 84.;
259  switch (SpinAlignedEOBversion) {
260  case 1:
261  coeffs->rho21v2 = -59. / 56 + (23. * eta) / 84. - 9. / 32. * a2;
262  coeffs->rho21v3 = 1177. / 672. * a - 27. / 128. * a3;
263  break;
264  case 2:
265  coeffs->rho21v2 = -59. / 56 + (23. * eta) / 84.;
266  coeffs->rho21v3 = 0.0;
267  break;
268  case 3:
269  coeffs->rho21v2 = -59. / 56 + (23. * eta) / 84.;
270  coeffs->rho21v3 = 0.0;
271  break;
272  case 4:
273  coeffs->rho21v2 = -59. / 56 + (23. * eta) / 84.;
274  coeffs->rho21v3 = 0.0;
275  break;
276  default:
277  XLALPrintError("XLAL Error - %s: wrong SpinAlignedEOBversion value, must be 1 or 2!\n", __func__);
278  XLAL_ERROR(XLAL_EINVAL);
279  break;
280  }
281  /*
282  * coeffs->rho21v3 = (-567.*chiA*chiA*chiA -
283  * 1701.*chiA*chiA*chiS*dM + chiA*(-4708. + 1701.*chiS*chiS -
284  * 2648.*eta)*(-1. + 4.*eta) + chiS* dM3 *(4708. -
285  * 567.*chiS*chiS + 1816.*eta))/(2688.*dM3);
286  */
287  coeffs->rho21v4 = -47009. / 56448. - (865. * a2) / 1792. - (405. * a2 * a2) / 2048. - (10993. * eta) / 14112.
288  + (617. * eta2) / 4704.;
289  coeffs->rho21v5 = (-98635. * a) / 75264. + (2031. * a * a2) / 7168. - (1701. * a2 * a3) / 8192.;
290  coeffs->rho21v6 = 7613184941. / 2607897600. + (9032393. * a2) / 1806336. + (3897. * a2 * a2) / 16384.
291  - (15309. * a3 * a3) / 65536.;
292  coeffs->rho21v6l = -107. / 105.;
293  coeffs->rho21v7 = (-3859374457. * a) / 1159065600. - (55169. * a3) / 16384.
294  + (18603. * a2 * a3) / 65536. - (72171. * a2 * a2 * a3) / 262144.;
295  coeffs->rho21v7l = 107. * a / 140.;
296  coeffs->rho21v8 = -1168617463883. / 911303737344.;
297  coeffs->rho21v8l = 6313. / 5880.;
298  coeffs->rho21v10 = -63735873771463. / 16569158860800.;
299  coeffs->rho21v10l = 5029963. / 5927040.;
300 
301  coeffs->f21v1 = (-3. * (chiS + chiA / dM)) / (2.);
302  switch (SpinAlignedEOBversion) {
303  case 1:
304  coeffs->f21v3 = 0.0;
305  break;
306  case 2:
307  coeffs->f21v3 = (chiS * dM * (427. + 79. * eta) + chiA * (147. + 280. * dM * dM + 1251. * eta)) / 84. / dM;
308  break;
309  case 3:
310  coeffs->f21v3 = (chiS * dM * (427. + 79. * eta) + chiA * (147. + 280. * dM * dM + 1251. * eta)) / 84. / dM;
311  break;
312  case 4:
313  coeffs->f21v3 =
314  (chiS * dM * (427. + 79. * eta) +
315  chiA * (147. + 280. * dM * dM + 1251. * eta)) / 84. / dM;
316  /*RC: New terms for SEOBNRv4HM, they are put to zero if use_hm == 0 */
317  if (waveform == 0)
318  {
319  coeffs->f21v4 = 0.0;
320  coeffs->f21v5 = 0.0;
321  coeffs->f21v6 = 0.0;
322  coeffs->f21v7c = 0;
323  }
324  else{
325  coeffs->f21v4 = (-3.-2.*eta)*chiA2 + (-3.+ eta/2.)*chiS2 + (-6.+21.*eta/2.)*chiS*chiA/dM;
326  coeffs->f21v5 = (3./4.-3.*eta)*chiA3/dM + (-81./16. +1709.*eta/1008. + 613.*eta2/1008.+(9./4.-3*eta)*chiA2)*chiS + 3./4.*chiS3
327  + (-81./16. - 703.*eta2/112. + 8797.*eta/1008.+(9./4. - 6.*eta)*chiS2)*chiA/dM;
328  coeffs->f21v6 = (4163./252.-9287.*eta/1008. - 85.*eta2/112.)*chiA2 + (4163./252. - 2633.*eta/1008. + 461.*eta2/1008.)*chiS2 + (4163./126.-1636.*eta/21. + 1088.*eta2/63.)*chiS*chiA/dM;
329  }
330  /* End new terms for SEOBNRv4HM */
331  break;
332  default:
333  XLALPrintError("XLAL Error - %s: wrong SpinAlignedEOBversion value, must be 1, 2 or 3!\n", __func__);
334  XLAL_ERROR(XLAL_EINVAL);
335  break;
336  }
337  } else {
338  coeffs->f21v1 = -3. * chiA / 2.; // Odd modes in dm->0 limit, note that there is a typo in Taracchini et.al. discussion after A7 (even->odd)
339  switch (SpinAlignedEOBversion) {
340  case 1:
341  coeffs->f21v3 = 0.0;
342  break;
343  case 2:
344  coeffs->f21v3 = (chiS * dM * (427. + 79. * eta) + chiA * (147. + 280. * dM * dM + 1251. * eta)) / 84.;
345  break;
346  case 3:
347  coeffs->f21v3 = (chiS * dM * (427. + 79. * eta) + chiA * (147. + 280. * dM * dM + 1251. * eta)) / 84.;
348  break;
349  case 4:
350  coeffs->f21v3 =
351  (chiS * dM * (427. + 79. * eta) +
352  chiA * (147. + 280. * dM * dM + 1251. * eta)) / 84.;
353  /* New terms for SEOBNRv4HM, they are put to zero if use_hm == 0 */
354  if (waveform == 0)
355  {
356  coeffs->f21v4 = 0.0;
357  coeffs->f21v5 = 0.0;
358  coeffs->f21v6 = 0.0;
359  }
360  else{
361  coeffs->f21v4 = (-6+21*eta/2.)*chiS*chiA;
362  coeffs->f21v5 = (3./4.-3.*eta)*chiA3 + (-81./16. - 703.*eta2/112. + 8797.*eta/1008. + (9./4. - 6.*eta)*chiS2)*chiA;
363  coeffs->f21v6 = (4163./126.-1636.*eta/21. + 1088.*eta2/63.)*chiS*chiA;
364  // printf("f21v1 = %.16f\n f21v3 = %.16f\n f21v4 = %.16f\n f21v5 = %.16f\n f21v6 = %.16f\n", coeffs->f21v1, coeffs->f21v3, coeffs->f21v4, coeffs->f21v5, coeffs->f21v6);
365  }
366  /* End new terms for SEOBNRv4HM */
367  break;
368  default:
369  XLALPrintError("XLAL Error - %s: wrong SpinAlignedEOBversion value, must be 1 or 2!\n", __func__);
370  XLAL_ERROR(XLAL_EINVAL);
371  break;
372  }
373  }
374 
375  /*
376  * l = 3, Eqs. A9a - A9c for rho, Eqs. A15b and A15c for f, Eqs. 22 -
377  * 24 of DIN and Eqs. 27c - 27e of PBFRT for delta
378  */
379  coeffs->delta33vh3 = 13. / 10.;
380  coeffs->delta33vh6 = (-81. * aDelta) / 20. + (39. * LAL_PI) / 7.;
381  coeffs->delta33vh9 = -227827. / 3000. + (78. * LAL_PI * LAL_PI) / 7.;
382  coeffs->delta33v5 = -80897. * eta / 2430.;
383  if (dM2) {
384  coeffs->rho33v2 = -7. / 6. + (2. * eta) / 3.;
385  //coeffs->rho33v3 = (chiS * dM * (-4. + 5. * eta) + chiA * (-4. + 19. * eta)) / (6. * dM);
386  coeffs->rho33v3 = 0.0;
387  coeffs->rho33v4 = -6719. / 3960. + a2 / 2. - (1861. * eta) / 990. + (149. * eta2) / 330.;
388  coeffs->rho33v5 = (-4. * a) / 3.;
389  coeffs->rho33v6 = 3203101567. / 227026800. + (5. * a2) / 36.;
390  coeffs->rho33v6l = -26. / 7.;
391  coeffs->rho33v7 = (5297. * a) / 2970. + a * a2 / 3.;
392  coeffs->rho33v8 = -57566572157. / 8562153600.;
393  coeffs->rho33v8l = 13. / 3.;
394  if(waveform == 1){
395  //RC: This terms are in Eq.A6 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
396  coeffs->rho33v6 = 3203101567. / 227026800. + (5. * a2) / 36. + (-129509./25740. + 41./192. * LAL_PI*LAL_PI)*eta - 274621./154440.*eta2+ 12011./46332.*eta3;
397  coeffs->rho33v10 = -903823148417327./30566888352000.;
398  coeffs->rho33v10l = 87347./13860.;
399  }
400 
401  coeffs->f33v3 = (chiS * dM * (-4. + 5. * eta) + chiA * (-4. + 19. * eta)) / (2. * dM);
402  coeffs->f33v4 = 0;
403  coeffs->f33v5 = 0;
404  coeffs->f33v6 = 0;
405  coeffs->f33vh6 = 0;
406  if(waveform == 1){
407  //RC: This terms are in Eq.A10 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
408  coeffs->f33v4 = (3./2. * chiS2 * dM + (3. - 12 * eta) * chiA * chiS + dM * (3./2. -6. * eta) * chiA2)/(dM);
409  coeffs->f33v5 = (dM * (241./30. * eta2 + 11./20. * eta + 2./3.) * chiS + (407./30. * eta2 - 593./60. * eta + 2./3.)* chiA)/(dM);
410  coeffs->f33v6 = (dM * (6. * eta2 -27. / 2. * eta - 7./ 4.) * chiS2 + (44. * eta2 - 1. * eta - 7./2.) * chiA * chiS + dM * (-12 * eta2 + 11./2. * eta - 7./4.) * chiA2)/dM;
411  coeffs->f33vh6 = (dM * (593. / 108. * eta - 81./20.) * chiS + (7339./540. * eta - 81./20.) * chiA)/(dM);
412  }
413  }
414  else {
415  coeffs->f33v3 = chiA * 3. / 8.;
416  if(waveform == 1){
417  //RC: This terms are in Eq.A10 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
418  coeffs->f33v4 = ((3. - 12 * eta) * chiA * chiS);
419  coeffs->f33v5 = ((407./30. * eta2 - 593./60. * eta + 2./3.)* chiA);
420  coeffs->f33v6 = ((44. * eta2 - 1. * eta - 7./2.) * chiA * chiS);
421  coeffs->f33vh6 = ((7339./540. * eta - 81./20.) * chiA);
422  }
423  }
424 
425  coeffs->delta32vh3 = (10. + 33. * eta) / (-15. * m1Plus3eta);
426  coeffs->delta32vh4 = 4. * aDelta;
427  coeffs->delta32vh6 = (-136. * aDelta) / 45. + (52. * LAL_PI) / 21.;
428  coeffs->delta32vh9 = -9112. / 405. + (208. * LAL_PI * LAL_PI) / 63.;
429 
430  coeffs->rho32v = (4. * chiS * eta) / (-3. * m1Plus3eta);
431  /** TODO The term proportional to eta^2 a^2 in coeffs->rho32v2 is wrong, but it was used in the calibration of SEOBNRv2 */
432  coeffs->rho32v2 = (-4. * a2 * eta2) / (9. * m1Plus3eta2) + (328. - 1115. * eta
433  + 320. * eta2) / (270. * m1Plus3eta);
434  if (SpinAlignedEOBversion == 4) {
435  coeffs->rho32v2 = (328. - 1115. * eta +
436  320. * eta2) / (270. * m1Plus3eta);
437  }
438  //coeffs->rho32v3 = (2. * (45. * a * m1Plus3eta3
439  // - a * eta * (328. - 2099. * eta + 5. * (733. + 20. * a2) * eta2
440  // - 960. * eta3))) / (405. * m1Plus3eta3);
441  coeffs->rho32v3 = 2. / 9. * a;
442  coeffs->rho32v4 = a2 / 3. + (-1444528.
443  + 8050045. * eta - 4725605. * eta2 - 20338960. * eta3
444  + 3085640. * eta2 * eta2) / (1603800. * m1Plus3eta2);
445  coeffs->rho32v5 = (-2788. * a) / 1215.;
446  coeffs->rho32v6 = 5849948554. / 940355325. + (488. * a2) / 405.;
447  coeffs->rho32v6l = -104. / 63.;
448  coeffs->rho32v8 = -10607269449358. / 3072140846775.;
449  coeffs->rho32v8l = 17056. / 8505.;
450 
451  if (dM2) {
452  coeffs->delta31vh3 = 13. / 30.;
453  coeffs->delta31vh6 = (61. * aDelta) / 20. + (13. * LAL_PI) / 21.;
454  coeffs->delta31vh7 = (-24. * aDelta * aDelta) / 5.;
455  coeffs->delta31vh9 = -227827. / 81000. + (26. * LAL_PI * LAL_PI) / 63.;
456  coeffs->delta31v5 = -17. * eta / 10.;
457 
458  coeffs->rho31v2 = -13. / 18. - (2. * eta) / 9.;
459  //coeffs->rho31v3 = (chiA * (-4. + 11. * eta) + chiS * dM * (-4. + 13. * eta)) / (6. * dM);
460  coeffs->rho31v3 = 0.0;
461  coeffs->rho31v4 = 101. / 7128.
462  - (5. * a2) / 6. - (1685. * eta) / 1782. - (829. * eta2) / 1782.;
463  coeffs->rho31v5 = (4. * a) / 9.;
464  coeffs->rho31v6 = 11706720301. / 6129723600. - (49. * a2) / 108.;
465  coeffs->rho31v6l = -26. / 63.;
466  coeffs->rho31v7 = (-2579. * a) / 5346. + a * a2 / 9.;
467  coeffs->rho31v8 = 2606097992581. / 4854741091200.;
468  coeffs->rho31v8l = 169. / 567.;
469 
470  coeffs->f31v3 = (chiA * (-4. + 11. * eta) + chiS * dM * (-4. + 13. * eta)) / (2. * dM);
471  } else {
472  coeffs->f31v3 = -chiA * 5. / 8.;
473  }
474 
475  /*
476  * l = 4, Eqs. A10a - A10d for rho, Eq. A15d for f Eqs. 25 - 28 of
477  * DIN and Eqs. 27f - 27i of PBFRT for delta
478  */
479 
480  coeffs->delta44vh3 = (112. + 219. * eta) / (-120. * m1Plus3eta);
481  coeffs->delta44vh6 = (-464. * aDelta) / 75. + (25136. * LAL_PI) / 3465.;
482  coeffs->delta44vh9 = 0.;
483  if(waveform == 1){
484  //RC: This terms are in Eq.A15 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
485  coeffs->delta44vh9 = -55144./375. + 201088.*LAL_PI*LAL_PI/10395.;
486  }
487 
488  coeffs->rho44v2 = (1614. - 5870. * eta + 2625. * eta2) / (1320. * m1Plus3eta);
489  coeffs->rho44v3 = (chiA * (10. - 39. * eta) * dM + chiS * (10. - 41. * eta
490  + 42. * eta2)) / (15. * m1Plus3eta);
491  coeffs->rho44v4 = a2 / 2. + (-511573572.
492  + 2338945704. * eta - 313857376. * eta2 - 6733146000. * eta3
493  + 1252563795. * eta2 * eta2) / (317116800. * m1Plus3eta2);
494  coeffs->rho44v5 = (-69. * a) / 55.;
495  coeffs->rho44v8 = 0.;
496  coeffs->rho44v8l = 0.;
497  coeffs->rho44v10 = 0.;
498  coeffs->rho44v10l = 0;
499  if(waveform == 1){
500  //RC: This terms are in Eq.A8 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
501  coeffs->rho44v4 =
502  (-511573572. + 2338945704. * eta - 313857376. * eta2 -
503  6733146000. * eta3 +
504  1252563795. * eta2 * eta2) / (317116800. * m1Plus3eta2)
505  + chiS2/2. + dM*chiS*chiA + dM2*chiA2/2.;
506  coeffs->rho44v5 = chiA*dM*(-8280. + 42716.*eta - 57990.*eta2 + 8955*eta3)/(6600.*m1Plus3eta2)
507  + chiS*(-8280. + 66284.*eta-176418.*eta2+128085.*eta3 + 88650*eta2*eta2)/(6600.*m1Plus3eta2);
508  coeffs->rho44v8 = -172066910136202271./19426955708160000.;
509  coeffs->rho44v8l = 845198./190575.;
510  coeffs->rho44v10 = - 17154485653213713419357./568432724020761600000.;
511  coeffs->rho44v10l = 22324502267./3815311500;
512  }
513  coeffs->rho44v6 = 16600939332793. / 1098809712000. + (217. * a2) / 3960.;
514  coeffs->rho44v6l = -12568. / 3465.;
515 
516  if (dM2) {
517  coeffs->delta43vh3 = (486. + 4961. * eta) / (810. * (1. - 2. * eta));
518  coeffs->delta43vh4 = (11. * aDelta) / 4.;
519  coeffs->delta43vh6 = 1571. * LAL_PI / 385.;
520 
521  //coeffs->rho43v = (5. * (chiA - chiS * dM) * eta) / (8. * dM * (-1. + 2. * eta));
522  coeffs->rho43v = 0.0;
523  coeffs->rho43v2 = (222. - 547. * eta + 160. * eta2) / (176. * (-1. + 2. * eta));
524  coeffs->rho43v4 = -6894273. / 7047040. + (3. * a2) / 8.;
525  coeffs->rho43v5 = (-12113. * a) / 6160.;
526  coeffs->rho43v6 = 1664224207351. / 195343948800.;
527  coeffs->rho43v6l = -1571. / 770.;
528 
529  coeffs->f43v = (5. * (chiA - chiS * dM) * eta) / (2. * dM * (-1. + 2. * eta));
530  } else {
531  coeffs->f43v = -5. * chiA / 4.;
532  }
533 
534  coeffs->delta42vh3 = (7. * (1. + 6. * eta)) / (-15. * m1Plus3eta);
535  coeffs->delta42vh6 = (212. * aDelta) / 75. + (6284. * LAL_PI) / 3465.;
536 
537  coeffs->rho42v2 = (1146. - 3530. * eta + 285. * eta2) / (1320. * m1Plus3eta);
538  coeffs->rho42v3 = (chiA * (10. - 21. * eta) * dM + chiS * (10. - 59. * eta
539  + 78. * eta2)) / (15. * m1Plus3eta);
540  coeffs->rho42v4 = a2 / 2. + (-114859044. + 295834536. * eta + 1204388696. * eta2 - 3047981160. * eta3
541  - 379526805. * eta2 * eta2) / (317116800. * m1Plus3eta2);
542  coeffs->rho42v5 = (-7. * a) / 110.;
543  coeffs->rho42v6 = 848238724511. / 219761942400. + (2323. * a2) / 3960.;
544  coeffs->rho42v6l = -3142. / 3465.;
545 
546  if (dM2) {
547  coeffs->delta41vh3 = (2. + 507. * eta) / (10. * (1. - 2. * eta));
548  coeffs->delta41vh4 = (11. * aDelta) / 12.;
549  coeffs->delta41vh6 = 1571. * LAL_PI / 3465.;
550 
551  //coeffs->rho41v = (5. * (chiA - chiS * dM) * eta) / (8. * dM * (-1. + 2. * eta));
552  coeffs->rho41v = 0.0;
553  coeffs->rho41v2 = (602. - 1385. * eta + 288. * eta2) / (528. * (-1. + 2. * eta));
554  coeffs->rho41v4 = -7775491. / 21141120. + (3. * a2) / 8.;
555  coeffs->rho41v5 = (-20033. * a) / 55440. - (5 * a * a2) / 6.;
556  coeffs->rho41v6 = 1227423222031. / 1758095539200.;
557  coeffs->rho41v6l = -1571. / 6930.;
558 
559  coeffs->f41v = (5. * (chiA - chiS * dM) * eta) / (2. * dM * (-1. + 2. * eta));
560  } else {
561  coeffs->f41v = -5. * chiA / 4.;
562  }
563 
564  /*
565  * l = 5, Eqs. A11a - A11e for rho, Eq. 29 of DIN and Eqs. E1a and
566  * E1b of PBFRT for delta
567  */
568  coeffs->delta55vh3 =
569  (96875. + 857528. * eta) / (131250. * (1 - 2 * eta));
570  coeffs->delta55vh6 = 0;
571  coeffs->delta55vh9 = 0;
572  if(waveform == 1){
573  //RC: This terms are in Eq.A16 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
574  coeffs->delta55vh6 = 3865./429.*LAL_PI;
575  coeffs->delta55vh9 = (-7686949127. + 954500400.*LAL_PI*LAL_PI)/31783752.;
576  }
577  if (dM2) {
578 
579  coeffs->rho55v2 = (487. - 1298. * eta + 512. * eta2) / (390. * (-1. + 2. * eta));
580  coeffs->rho55v3 = (-2. * a) / 3.;
581  coeffs->rho55v4 = -3353747. / 2129400. + a2 / 2.;
582  coeffs->rho55v5 = -241. * a / 195.;
583  coeffs->rho55v6 = 0.;
584  coeffs->rho55v6l = 0.;
585  coeffs->rho55v8 = 0.;
586  coeffs->rho55v8l = 0.;
587  coeffs->rho55v10 = 0.;
588  coeffs->rho55v10l = 0.;
589  coeffs->f55v3 = 0.;
590  coeffs->f55v4 = 0.;
591  coeffs->f55v5c = 0;
592  if(waveform == 1){
593  //RC: This terms are in Eq.A9 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
594  coeffs->rho55v6 = 190606537999247./11957879934000.;
595  coeffs->rho55v6l = - 1546./429.;
596  coeffs->rho55v8 = - 1213641959949291437./118143853747920000.;
597  coeffs->rho55v8l = 376451./83655.;
598  coeffs->rho55v10 = -150082616449726042201261./4837990810977324000000.;
599  coeffs->rho55v10l = 2592446431./456756300.;
600 
601  coeffs->f55v3 = chiA/dM *(10./(3.*(-1.+2.*eta)) - 70.*eta/(3.*(-1.+2.*eta)) + 110.*eta2/(3.*(-1.+2.*eta)) ) +
602  chiS*(10./(3.*(-1.+2.*eta)) -10.*eta/(-1.+2.*eta) + 10*eta2/(-1.+2.*eta));
603  coeffs->f55v4 = chiS2*(-5./(2.*(-1.+2.*eta)) + 5.*eta/(-1.+2.*eta)) +
604  chiA*chiS/dM *(-5./(-1.+2.*eta) + 30.*eta/(-1.+2.*eta) - 40.*eta2/(-1.+2.*eta)) +
605  chiA2*(-5./(2.*(-1.+2.*eta)) + 15.*eta/(-1.+2.*eta) - 20.*eta2/(-1.+2.*eta));
606  coeffs->f55v5c = 0; //RC: this is the calibration parameter which is initially set to 0.
607  }
608  }
609  else{
610  coeffs->f55v3 = 0;
611  coeffs->f55v4 = 0;
612  coeffs->f55v5c = 0;
613  if(waveform == 1){
614  //RC: This terms are in Eq.A12 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
615  coeffs->f55v3 = chiA *(10./(3.*(-1.+2.*eta)) - 70.*eta/(3.*(-1.+2.*eta)) + 110.*eta2/(3.*(-1.+2.*eta)) );
616  coeffs->f55v4 = chiA*chiS *(-5./(-1.+2.*eta) + 30.*eta/(-1.+2.*eta) - 40.*eta2/(-1.+2.*eta));
617  coeffs->f55v5c = 0;
618  }
619  }
620  coeffs->delta54vh3 = 8. / 15.;
621  coeffs->delta54vh4 = 12. * aDelta / 5.;
622 
623  coeffs->rho54v2 = (-17448. + 96019. * eta - 127610. * eta2
624  + 33320. * eta3) / (13650. * (1. - 5. * eta + 5. * eta2));
625  coeffs->rho54v3 = (-2. * a) / 15.;
626  coeffs->rho54v4 = -16213384. / 15526875. + (2. * a2) / 5.;
627 
628  if (dM2) {
629  coeffs->delta53vh3 = 31. / 70.;
630 
631  coeffs->rho53v2 = (375. - 850. * eta + 176. * eta2) / (390. * (-1. + 2. * eta));
632  coeffs->rho53v3 = (-2. * a) / 3.;
633  coeffs->rho53v4 = -410833. / 709800. + a2 / 2.;
634  coeffs->rho53v5 = -103. * a / 325.;
635  }
636  coeffs->delta52vh3 = 4. / 15.;
637  coeffs->delta52vh4 = 6. * aDelta / 5.;
638 
639  coeffs->rho52v2 = (-15828. + 84679. * eta - 104930. * eta2
640  + 21980. * eta3) / (13650. * (1. - 5. * eta + 5. * eta2));
641  coeffs->rho52v3 = (-2. * a) / 15.;
642  coeffs->rho52v4 = -7187914. / 15526875. + (2. * a2) / 5.;
643 
644  if (dM2) {
645  coeffs->delta51vh3 = 31. / 210.;
646 
647  coeffs->rho51v2 = (319. - 626. * eta + 8. * eta2) / (390. * (-1. + 2. * eta));
648  coeffs->rho51v3 = (-2. * a) / 3.;
649  coeffs->rho51v4 = -31877. / 304200. + a2 / 2.;
650  coeffs->rho51v5 = 139. * a / 975.;
651  }
652  /*
653  * l = 6, Eqs. A12a - A12f for rho, Eqs. E1c and E1d of PBFRT for
654  * delta
655  */
656 
657  coeffs->delta66vh3 = 43. / 70.;
658 
659  coeffs->rho66v2 = (-106. + 602. * eta - 861. * eta2
660  + 273. * eta3) / (84. * (1. - 5. * eta + 5. * eta2));
661  coeffs->rho66v3 = (-2. * a) / 3.;
662  coeffs->rho66v4 = -1025435. / 659736. + a2 / 2.;
663 
664  if (dM2) {
665  coeffs->delta65vh3 = 10. / 21.;
666 
667  coeffs->rho65v2 = (-185. + 838. * eta - 910. * eta2
668  + 220. * eta3) / (144. * (dM2 + 3. * eta2));
669  coeffs->rho65v3 = -2. * a / 9.;
670  }
671  coeffs->delta64vh3 = 43. / 105.;
672 
673  coeffs->rho64v2 = (-86. + 462. * eta - 581. * eta2
674  + 133. * eta3) / (84. * (1. - 5. * eta + 5. * eta2));
675  coeffs->rho64v3 = (-2. * a) / 3.;
676  coeffs->rho64v4 = -476887. / 659736. + a2 / 2.;
677 
678  if (dM2) {
679  coeffs->delta63vh3 = 2. / 7.;
680 
681  coeffs->rho63v2 = (-169. + 742. * eta - 750. * eta2
682  + 156. * eta3) / (144. * (dM2 + 3. * eta2));
683  coeffs->rho63v3 = -2. * a / 9.;
684  }
685  coeffs->delta62vh3 = 43. / 210.;
686 
687  coeffs->rho62v2 = (-74. + 378. * eta - 413. * eta2
688  + 49. * eta3) / (84. * (1. - 5. * eta + 5. * eta2));
689  coeffs->rho62v3 = (-2. * a) / 3.;
690  coeffs->rho62v4 = -817991. / 3298680. + a2 / 2.;
691 
692  if (dM2) {
693  coeffs->delta61vh3 = 2. / 21.;
694 
695  coeffs->rho61v2 = (-161. + 694. * eta - 670. * eta2
696  + 124. * eta3) / (144. * (dM2 + 3. * eta2));
697  coeffs->rho61v3 = -2. * a / 9.;
698  }
699  /*
700  * l = 7, Eqs. A13a - A13g for rho, Eqs. E1e and E1f of PBFRT for
701  * delta
702  */
703  if (dM2) {
704  coeffs->delta77vh3 = 19. / 36.;
705 
706  coeffs->rho77v2 = (-906. + 4246. * eta - 4963. * eta2
707  + 1380. * eta3) / (714. * (dM2 + 3. * eta2));
708  coeffs->rho77v3 = -2. * a / 3.;
709  }
710  coeffs->rho76v2 = (2144. - 16185. * eta + 37828. * eta2 - 29351. * eta3
711  + 6104. * eta2 * eta2) / (1666. * (-1 + 7 * eta - 14 * eta2
712  + 7 * eta3));
713 
714  if (dM2) {
715  coeffs->delta75vh3 = 95. / 252.;
716 
717  coeffs->rho75v2 = (-762. + 3382. * eta - 3523. * eta2
718  + 804. * eta3) / (714. * (dM2 + 3. * eta2));
719  coeffs->rho75v3 = -2. * a / 3.;
720  }
721  coeffs->rho74v2 = (17756. - 131805. * eta + 298872. * eta2 - 217959. * eta3
722  + 41076. * eta2 * eta2) / (14994. * (-1. + 7. * eta - 14. * eta2
723  + 7. * eta3));
724 
725  if (dM2) {
726  coeffs->delta73vh3 = 19. / 84.;
727 
728  coeffs->rho73v2 = (-666. + 2806. * eta - 2563. * eta2
729  + 420. * eta3) / (714. * (dM2 + 3. * eta2));
730  coeffs->rho73v3 = -2. * a / 3.;
731  }
732  coeffs->rho72v2 = (16832. - 123489. * eta + 273924. * eta2 - 190239. * eta3
733  + 32760. * eta2 * eta2) / (14994. * (-1. + 7. * eta - 14. * eta2
734  + 7. * eta3));
735 
736  if (dM2) {
737  coeffs->delta71vh3 = 19. / 252.;
738 
739  coeffs->rho71v2 = (-618. + 2518. * eta - 2083. * eta2
740  + 228. * eta3) / (714. * (dM2 + 3. * eta2));
741  coeffs->rho71v3 = -2. * a / 3.;
742  }
743  /* l = 8, Eqs. A14a - A14h */
744 
745  coeffs->rho88v2 = (3482. - 26778. * eta + 64659. * eta2 - 53445. * eta3
746  + 12243. * eta2 * eta2) / (2736. * (-1. + 7. * eta - 14. * eta2
747  + 7. * eta3));
748 
749  if (dM2) {
750  coeffs->rho87v2 = (23478. - 154099. * eta + 309498. * eta2 - 207550. * eta3
751  + 38920 * eta2 * eta2) / (18240. * (-1 + 6 * eta - 10 * eta2
752  + 4 * eta3));
753  }
754  coeffs->rho86v2 = (1002. - 7498. * eta + 17269. * eta2 - 13055. * eta3
755  + 2653. * eta2 * eta2) / (912. * (-1. + 7. * eta - 14. * eta2
756  + 7. * eta3));
757 
758  if (dM2) {
759  coeffs->rho85v2 = (4350. - 28055. * eta + 54642. * eta2 - 34598. * eta3
760  + 6056. * eta2 * eta2) / (3648. * (-1. + 6. * eta - 10. * eta2
761  + 4. * eta3));
762  }
763  coeffs->rho84v2 = (2666. - 19434. * eta + 42627. * eta2 - 28965. * eta3
764  + 4899. * eta2 * eta2) / (2736. * (-1. + 7. * eta - 14. * eta2
765  + 7. * eta3));
766 
767  if (dM2) {
768  coeffs->rho83v2 = (20598. - 131059. * eta + 249018. * eta2 - 149950. * eta3
769  + 24520. * eta2 * eta2) / (18240. * (-1. + 6. * eta - 10. * eta2
770  + 4. * eta3));
771  }
772  coeffs->rho82v2 = (2462. - 17598. * eta + 37119. * eta2 - 22845. * eta3
773  + 3063. * eta2 * eta2) / (2736. * (-1. + 7. * eta - 14. * eta2
774  + 7. * eta3));
775 
776  if (dM2) {
777  coeffs->rho81v2 = (20022. - 126451. * eta + 236922. * eta2 - 138430. * eta3
778  + 21640. * eta2 * eta2) / (18240. * (-1. + 6. * eta - 10. * eta2
779  + 4. * eta3));
780  }
781  /* All relevant coefficients should be set, so we return */
782 
783  return XLAL_SUCCESS;
784 }
785 
786 #endif
XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients
static int XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients(FacWaveformCoeffs *const coeffs, const REAL8 m1, const REAL8 m2, const REAL8 eta, const REAL8 a, const REAL8 chiS, const REAL8 chiA, const UINT4 SpinAlignedEOBversion)
Waveform expressions are given by Taracchini et al.
Definition: LALSimIMRSpinEOBFactorizedWaveformCoefficientsPrec.c:66
a2
const double a2
Definition: exact_derivatives-HrealHeader.c:33
LAL_PI
#define LAL_PI
REAL8
double REAL8
UINT4
uint32_t UINT4
INT4
int32_t INT4
a
static const INT4 a
XLAL_PRINT_INFO
#define XLAL_PRINT_INFO(...)
XLAL_ERROR
#define XLAL_ERROR(...)
XLALPrintError
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EINVAL
XLAL_EINVAL
FacWaveformCoeffs
Structure containing the coefficients for calculating the factorized waveform.
Definition: LALSimIMREOBNRv2.h:97