LALSimulation  5.4.0.1-fe68b98
LALSimIMRSpinEOBFactorizedWaveformPrec.c
Go to the documentation of this file.
1 /**
2  * \author Craig Robinson, Yi Pan, Prayush Kumar, Stas Babak, Andrea Taracchini
3  *
4  * \brief Function to compute the factorized waveform as used in the SEOBNRv1 model.
5  * Waveform expressions are given by
6  * Taracchini et al. ( PRD 86, 024011 (2012), arXiv 1202.0790 ).
7  * All equation numbers in this file refer to equations of this paper,
8  * unless otherwise specified.
9  * Coefficients of the so-called "deltalm" terms are given by
10  * Damour et al. PRD 79, 064004 (2009) and Pan et al. PRD 83, 064003 (2011),
11  * henceforth DIN and PBFRT.
12  *
13  * Functions to compute the factorized waveform for the purpose of computing the
14  * flux, i.e. returning only the absolute value of the multipoles. The tail term
15  * Tlm is used in its resummed form, given by Eq. (42) of Damour, Nagar and
16  * Bernuzzi, PRD 87 084035 (2013), called DNB here.
17  *
18  */
19 
20 #ifndef _LALSIMIMRSPINPRECEOBFACTORIZEDWAVEFORM_C
21 #define _LALSIMIMRSPINPRECEOBFACTORIZEDWAVEFORM_C
22 #include <gsl/gsl_deriv.h>
23 #include <complex.h>
24 #include <lal/LALSimInspiral.h>
25 #include <lal/LALSimIMR.h>
26 #include <gsl/gsl_spline.h>
27 #include <gsl/gsl_linalg.h>
28 
29 #include "LALSimIMREOBNRv2.h"
30 #include "LALSimIMRSpinEOB.h"
31 
32 #include "LALSimIMREOBNewtonianMultipole.c"
33 #include "LALSimIMRSpinEOBAuxFuncs.c"
34 #include "LALSimIMREOBNQCCorrection.c"
35 #include "LALSimIMRSpinEOBHamiltonian.c"
36 #include "LALSimIMRSpinEOBHamiltonianPrec.c"
37 #include "LALSimIMRSpinEOBFactorizedWaveform.c"
38 #define v4Pwave 451
39 /* #include "LALSimIMRSpinEOBFactorizedWaveform.c" */
40 
41 
42 /*------------------------------------------------------------------------------------------
43  *
44  * Prototypes of functions defined in this code.
45  *
46  *------------------------------------------------------------------------------------------
47  */
48 static INT4
49 XLALSimIMRSpinEOBGetPrecSpinFactorizedWaveform(
50  COMPLEX16 * restrict hlm, /**< OUTPUT, hlm waveforms */
51  REAL8Vector * restrict values, /**< dyanmical variables: \f$(r,\phi,p_r,p_\phi)\f$ */
52  REAL8Vector * restrict cartvalues, /**< dyanmical variables */
53  const REAL8 v, /**< velocity */
54  const REAL8 Hreal, /**< real Hamiltonian */
55  const INT4 l, /**< l mode index */
56  const INT4 m, /**< m mode index */
57  SpinEOBParams * restrict params /**< Spin EOB parameters */
58 );
59 
60 static INT4
61 XLALSimIMRSpinEOBFluxGetPrecSpinFactorizedWaveform(
62  COMPLEX16 * restrict hlmTab, /**< OUTPUT, hlm waveforms */
63  REAL8Vector * restrict values, /**< dyanmical variables: \f$(r,\phi,p_r,p_\phi)\f$ */
64  REAL8Vector * restrict cartvalues, /**< dyanmical variables */
65  const REAL8 v, /**< velocity */
66  const REAL8 Hreal, /**< real Hamiltonian */
67  const INT4 lMax, /**< maximum l mode to compute, compute 0 < m <= lMax */
68  SpinEOBParams * restrict params /**< Spin EOB parameters */
69 );
70 
71 UNUSED static INT4 XLALSimIMRSpinEOBGetAmplitudeResidualPrec (
72  COMPLEX16 * restrict rholmpwrl,
73  const REAL8 v,
74  const REAL8 Hreal,
75  const INT4 modeL,
76  const INT4 modeM,
77  SpinEOBParams * restrict params
78 );
79 
80 UNUSED static INT4 XLALSimIMREOBCalcCalibCoefficientHigherModesPrec (
81  SpinEOBParams * restrict UNUSED params, /**<< Output */
82  const UINT4 modeL, /**<< Mode index L */
83  const UINT4 modeM, /**<< Mode index M */
84  SEOBdynamics *seobdynamics /** Dynamics array*/,
85  const REAL8 timeorb, /**<< Time of the peak of the orbital frequency */
86  const REAL8 m1, /**<< Component mass 1 */
87  const REAL8 m2, /**<< Component mass 2 */
88  const REAL8 UNUSED deltaT /**<< Sampling interval */
89  );
90 
91 
92 
93 
94 
95 /*------------------------------------------------------------------------------------------
96  *
97  * Defintions of functions.
98  *
99  *------------------------------------------------------------------------------------------
100  */
101 
102 
103 /** ######################################################### */
104 /** FOR PRECESSING EOB
105  * This function calculates hlm mode factorized-resummed waveform
106  * for given dynamical variables. This is optimized for flux calculation,
107  * by ignoring complex arguments and keeping only absolute values.
108  * Changes:
109  * (i) Complex Argument of Tlm not exponentiated.
110  * (ii) exp(i deltalm) set to 1.
111  * Eq. 17 and the entire Appendix of PRD 86, 024011 (2012) + changes
112  * described in the section "Factorized waveforms" of https://dcc.ligo.org/T1400476
113  */
114 static INT4
115 XLALSimIMRSpinEOBFluxGetPrecSpinFactorizedWaveform(
116  COMPLEX16 * restrict hlmTab, /**< OUTPUT, hlm waveforms */
117  REAL8Vector * restrict values, /**< dyanmical variables: \f$(r,\phi,p_r,p_\phi)\f$ */
118  REAL8Vector * restrict cartvalues, /**< dyanmical variables */
119  const REAL8 v, /**< velocity */
120  const REAL8 Hreal, /**< real Hamiltonian */
121  const INT4 lMax, /**< maximum l mode to compute, compute 0 < m <= lMax */
122  SpinEOBParams * restrict params /**< Spin EOB parameters */
123 )
124 {
125  int debugPK = 0;
126  const REAL8 vPhiKepler = params->alignedSpins ?
127  XLALSimIMRSpinAlignedEOBNonKeplerCoeff(values->data, params) :
128  XLALSimIMRSpinPrecEOBNonKeplerCoeff(cartvalues->data, params);
129  if (XLAL_IS_REAL8_FAIL_NAN(vPhiKepler)) {
130  XLAL_ERROR(XLAL_EFUNC);
131  }
132 
133  for (INT4 l = 2; l <= lMax; l++)
134  {
135  for (INT4 m = 1; m <= l; m++)
136  {
137  COMPLEX16 *hlm = &hlmTab[l*(lMax+1)+m];
138 
139 
140  /* Status of function calls */
141  INT4 status;
142  INT4 i;
143 
144  REAL8 eta;
145  REAL8 UNUSED r , pp, Omega, v2, /* vh, vh3, */ k, hathatk, eulerlogxabs;
146  //pr
147  REAL8 UNUSED rcrossp_x, rcrossp_y, rcrossp_z;
148  REAL8 Slm , rholm, rholmPwrl;
149  REAL8 auxflm = 0.0;
150  REAL8 hathatksq4, hathatk4pi, Tlmprefac, Tlmprodfac;
151  REAL8 Tlm;
152  COMPLEX16 hNewton;
153  gsl_sf_result z2;
154 
155  /* Non-Keplerian velocity */
156  REAL8 vPhi , vPhi2;
157 
158  /* Pre-computed coefficients */
159  FacWaveformCoeffs *hCoeffs = params->eobParams->hCoeffs;
160 
161  if (abs(m) > (INT4) l) {
162  XLAL_ERROR(XLAL_EINVAL);
163  }
164  if (m == 0) {
165  XLAL_ERROR(XLAL_EINVAL);
166  }
167  eta = params->eobParams->eta;
168 
169  /* Check our eta was sensible */
170  if (eta > 0.25) {
171  XLALPrintError("XLAL Error - %s: Eta seems to be > 0.25 - this isn't allowed!\n", __func__);
172  XLAL_ERROR(XLAL_EINVAL);
173  }
174  /*
175  * else if ( eta == 0.25 && m % 2 ) { // If m is odd and dM = 0, hLM
176  * will be zero memset( hlm, 0, sizeof( COMPLEX16 ) ); return
177  * XLAL_SUCCESS; }
178  */
179 
180  //r = sqrt(values->data[0] * values->data[0] + values->data[1] * values->data[1] + values->data[2] * values->data[2]);
181  //pr = values->data[2];
182  r = values->data[0];
183  pp = values->data[3];
184 
185  rcrossp_x = cartvalues->data[1] * cartvalues->data[5] - cartvalues->data[2] * cartvalues->data[4];
186  rcrossp_y = cartvalues->data[2] * cartvalues->data[3] - cartvalues->data[0] * cartvalues->data[5];
187  rcrossp_z = cartvalues->data[0] * cartvalues->data[4] - cartvalues->data[1] * cartvalues->data[3];
188 
189  //pp = sqrt(rcrossp_x * rcrossp_x + rcrossp_y * rcrossp_y + rcrossp_z * rcrossp_z);
190 
191  v2 = v * v;
192  Omega = v2 * v;
193  //vh3 = Hreal * Omega;
194  //vh = cbrt(vh3);
195  eulerlogxabs = LAL_GAMMA + log(2.0 * (REAL8) m * v);
196 
197  /* Calculate the non-Keplerian velocity */
198  vPhi = vPhiKepler;
199 
200  vPhi = r * cbrt(vPhi);
201 
202  if (debugPK)
203  XLAL_PRINT_INFO("In XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole, getting rW = %.12e\n",
204  vPhi);
205  vPhi *= Omega;
206  vPhi2 = vPhi * vPhi;
207 
208  /*
209  * Calculate the newtonian multipole, 1st term in Eq. 17, given by
210  * Eq. A1
211  */
212  //debugPK
213  if (debugPK) {
214 // XLAL_PRINT_INFO("\nValues inside XLALSimIMRSpinEOBFluxGetPrecSpinFactorizedWaveform:\n");
215 // for (i = 0; i < 14; i++)
216 // XLAL_PRINT_INFO("values[%d] = %.12e\n", i, values->data[i]);
217 
218  XLAL_PRINT_INFO("Calculating hNewton, with v = %.12e, vPhi = %.12e, r = %.12e, Phi = %.12e, l = %d, m = %d\n",
219  v, vPhi, r, values->data[1], (UINT4) l, (UINT4) m);
220  }
221  status = XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole(&hNewton,
222  vPhi2, r, values->data[1], (UINT4) l, m, params->eobParams);
223  if (status == XLAL_FAILURE) {
224  XLAL_ERROR(XLAL_EFUNC);
225  }
226  /* Calculate the source term, 2nd term in Eq. 17, given by Eq. A5 */
227  if (((l + m) % 2) == 0) {
228  Slm = (Hreal * Hreal - 1.) / (2. * eta) + 1.;
229  } else {
230  Slm = v * pp;
231  //Slm = v * sqrt(rcrossp_x * rcrossp_x + rcrossp_y * rcrossp_y + rcrossp_z * rcrossp_z);
232  }
233  if (debugPK)
234  XLAL_PRINT_INFO("In XLALSimIMRSpinEOBFluxGetSpinFactorizedWaveform: Hreal = %e, Slm = %e, eta = %e\n", Hreal, Slm, eta);
235 
236  /*
237  * Calculate the absolute value of the Tail term, 3rd term in Eq. 17,
238  * given by Eq. A6, and Eq. (42) of
239  * http://arxiv.org/pdf/1212.4357.pdf
240  */
241  k = m * Omega;
242  hathatk = Hreal * k;
243  hathatksq4 = 4. * hathatk * hathatk;
244  hathatk4pi = 4. * LAL_PI * hathatk;
245  /*
246  * gsl_sf_result lnr1, arg1; XLAL_CALLGSL( status =
247  * gsl_sf_lngamma_complex_e( l+1.0, -2.0*hathatk, &lnr1, &arg1 ) );
248  * if (status != GSL_SUCCESS) { XLALPrintError("XLAL Error - %s:
249  * Error in GSL function\n", __func__ ); XLAL_ERROR( XLAL_EFUNC ); }
250  */
251  XLAL_CALLGSL(status = gsl_sf_fact_e(l, &z2));
252  if (status != GSL_SUCCESS) {
253  XLALPrintError("XLAL Error - %s: Error in GSL function\n", __func__);
254  XLAL_ERROR(XLAL_EFUNC);
255  }
256  /*
257  * COMPLEX16 Tlmold; Tlmold = cexp( ( lnr1.val + LAL_PI * hathatk ) +
258  * I * ( arg1.val + 2.0 * hathatk * log(4.0*k/sqrt(LAL_E)) ) );
259  * Tlmold /= z2.val;
260  */
261  /* Calculating the prefactor of Tlm, outside the multiple product */
262  Tlmprefac = sqrt(hathatk4pi / (1. - exp(-hathatk4pi))) / z2.val;
263 
264  /* Calculating the multiple product factor */
265  for (Tlmprodfac = 1., i = 1; i <= l; i++) {
266  Tlmprodfac *= (hathatksq4 + (REAL8) i * i);
267  }
268 
269  Tlm = Tlmprefac * sqrt(Tlmprodfac);
270 
271  /* Calculate the residue phase and amplitude terms */
272  /*
273  * deltalm is the 4th term in Eq. 17, delta 22 given by Eq. A15,
274  * others
275  */
276  /* rholm is the 5th term in Eq. 17, given by Eqs. A8 - A14 */
277  /*
278  * auxflm is a special part of the 5th term in Eq. 17, given by Eq.
279  * A15
280  */
281  /*
282  * Actual values of the coefficients are defined in the next function
283  * of this file
284  */
285  switch (l) {
286  case 2:
287  switch (abs(m)) {
288  case 2:
289  rholm = 1. + v2 * (hCoeffs->rho22v2 + v * (hCoeffs->rho22v3
290  + v * (hCoeffs->rho22v4
291  + v * (hCoeffs->rho22v5 + v * (hCoeffs->rho22v6
292  + hCoeffs->rho22v6l * eulerlogxabs + v * (hCoeffs->rho22v7
293  + v * (hCoeffs->rho22v8 + hCoeffs->rho22v8l * eulerlogxabs
294  + (hCoeffs->rho22v10 + hCoeffs->rho22v10l * eulerlogxabs) * v2)))))));
295  //FIXME
296  if (debugPK){
297  XLAL_PRINT_INFO("PK:: rho22v2 = %.12e, rho22v3 = %.12e, rho22v4 = %.12e,\n rho22v5 = %.16e, rho22v6 = %.16e, rho22v6LOG = %.16e, \n rho22v7 = %.12e, rho22v8 = %.16e, rho22v8LOG = %.16e, \n rho22v10 = %.16e, rho22v10LOG = %.16e\n, rho22v6 = %.12e, rho22v8 = %.12e, rho22v10 = %.12e\n",
298  hCoeffs->rho22v2, hCoeffs->rho22v3, hCoeffs->rho22v4,
299  hCoeffs->rho22v5, hCoeffs->rho22v6, hCoeffs->rho22v6l,
300  hCoeffs->rho22v7, hCoeffs->rho22v8, hCoeffs->rho22v8l,
301  hCoeffs->rho22v10, hCoeffs->rho22v10l,
302  hCoeffs->rho22v6 + hCoeffs->rho22v6l * eulerlogxabs,
303  hCoeffs->rho22v8 + hCoeffs->rho22v8l * eulerlogxabs,
304  hCoeffs->rho22v10 + hCoeffs->rho22v10l * eulerlogxabs);}
305  break;
306  case 1:
307  {
308  rholm = 1. + v * (hCoeffs->rho21v1
309  + v * (hCoeffs->rho21v2 + v * (hCoeffs->rho21v3 + v * (hCoeffs->rho21v4
310  + v * (hCoeffs->rho21v5 + v * (hCoeffs->rho21v6 + hCoeffs->rho21v6l * eulerlogxabs
311  + v * (hCoeffs->rho21v7 + hCoeffs->rho21v7l * eulerlogxabs
312  + v * (hCoeffs->rho21v8 + hCoeffs->rho21v8l * eulerlogxabs
313  + (hCoeffs->rho21v10 + hCoeffs->rho21v10l * eulerlogxabs) * v2))))))));
314  auxflm = v * hCoeffs->f21v1 + v2 * v * hCoeffs->f21v3;
315  }
316  break;
317  default:
318  XLAL_ERROR(XLAL_EINVAL);
319  break;
320  }
321  break;
322  case 3:
323  switch (m) {
324  case 3:
325  rholm = 1. + v2 * (hCoeffs->rho33v2 + v * (hCoeffs->rho33v3 + v * (hCoeffs->rho33v4
326  + v * (hCoeffs->rho33v5 + v * (hCoeffs->rho33v6 + hCoeffs->rho33v6l * eulerlogxabs
327  + v * (hCoeffs->rho33v7 + (hCoeffs->rho33v8 + hCoeffs->rho33v8l * eulerlogxabs) * v))))));
328  auxflm = v * v2 * hCoeffs->f33v3;
329  break;
330  case 2:
331  rholm = 1. + v * (hCoeffs->rho32v
332  + v * (hCoeffs->rho32v2 + v * (hCoeffs->rho32v3 + v * (hCoeffs->rho32v4 + v * (hCoeffs->rho32v5
333  + v * (hCoeffs->rho32v6 + hCoeffs->rho32v6l * eulerlogxabs
334  + (hCoeffs->rho32v8 + hCoeffs->rho32v8l * eulerlogxabs) * v2))))));
335  break;
336  case 1:
337  rholm = 1. + v2 * (hCoeffs->rho31v2 + v * (hCoeffs->rho31v3 + v * (hCoeffs->rho31v4
338  + v * (hCoeffs->rho31v5 + v * (hCoeffs->rho31v6 + hCoeffs->rho31v6l * eulerlogxabs
339  + v * (hCoeffs->rho31v7 + (hCoeffs->rho31v8 + hCoeffs->rho31v8l * eulerlogxabs) * v))))));
340  auxflm = v * v2 * hCoeffs->f31v3;
341  break;
342  default:
343  XLAL_ERROR(XLAL_EINVAL);
344  break;
345  }
346  break;
347  case 4:
348  switch (m) {
349  case 4:
350 
351  rholm = 1. + v2 * (hCoeffs->rho44v2
352  + v * (hCoeffs->rho44v3 + v * (hCoeffs->rho44v4
353  + v * (hCoeffs->rho44v5 + (hCoeffs->rho44v6
354  + hCoeffs->rho44v6l * eulerlogxabs) * v))));
355  break;
356  case 3:
357  rholm = 1. + v * (hCoeffs->rho43v
358  + v * (hCoeffs->rho43v2
359  + v2 * (hCoeffs->rho43v4 + v * (hCoeffs->rho43v5
360  + (hCoeffs->rho43v6 + hCoeffs->rho43v6l * eulerlogxabs) * v))));
361  auxflm = v * hCoeffs->f43v;
362  break;
363  case 2:
364  rholm = 1. + v2 * (hCoeffs->rho42v2
365  + v * (hCoeffs->rho42v3 + v * (hCoeffs->rho42v4 + v * (hCoeffs->rho42v5
366  + (hCoeffs->rho42v6 + hCoeffs->rho42v6l * eulerlogxabs) * v))));
367  break;
368  case 1:
369  rholm = 1. + v * (hCoeffs->rho41v
370  + v * (hCoeffs->rho41v2
371  + v2 * (hCoeffs->rho41v4 + v * (hCoeffs->rho41v5
372  + (hCoeffs->rho41v6 + hCoeffs->rho41v6l * eulerlogxabs) * v))));
373  auxflm = v * hCoeffs->f41v;
374  break;
375  default:
376  XLAL_ERROR(XLAL_EINVAL);
377  break;
378  }
379  break;
380  case 5:
381  switch (m) {
382  case 5:
383  rholm = 1. + v2 * (hCoeffs->rho55v2
384  + v * (hCoeffs->rho55v3 + v * (hCoeffs->rho55v4
385  + v * (hCoeffs->rho55v5 + hCoeffs->rho55v6 * v))));
386  break;
387  case 4:
388  rholm = 1. + v2 * (hCoeffs->rho54v2 + v * (hCoeffs->rho54v3
389  + hCoeffs->rho54v4 * v));
390  break;
391  case 3:
392  rholm = 1. + v2 * (hCoeffs->rho53v2
393  + v * (hCoeffs->rho53v3 + v * (hCoeffs->rho53v4 + hCoeffs->rho53v5 * v)));
394  break;
395  case 2:
396  rholm = 1. + v2 * (hCoeffs->rho52v2 + v * (hCoeffs->rho52v3
397  + hCoeffs->rho52v4 * v));
398  break;
399  case 1:
400  rholm = 1. + v2 * (hCoeffs->rho51v2
401  + v * (hCoeffs->rho51v3 + v * (hCoeffs->rho51v4 + hCoeffs->rho51v5 * v)));
402  break;
403  default:
404  XLAL_ERROR(XLAL_EINVAL);
405  break;
406  }
407  break;
408  case 6:
409  switch (m) {
410  case 6:
411  rholm = 1. + v2 * (hCoeffs->rho66v2 + v * (hCoeffs->rho66v3
412  + hCoeffs->rho66v4 * v));
413  break;
414  case 5:
415  rholm = 1. + v2 * (hCoeffs->rho65v2 + hCoeffs->rho65v3 * v);
416  break;
417  case 4:
418  rholm = 1. + v2 * (hCoeffs->rho64v2 + v * (hCoeffs->rho64v3
419  + hCoeffs->rho64v4 * v));
420  break;
421  case 3:
422  rholm = 1. + v2 * (hCoeffs->rho63v2 + hCoeffs->rho63v3 * v);
423  break;
424  case 2:
425  rholm = 1. + v2 * (hCoeffs->rho62v2 + v * (hCoeffs->rho62v3
426  + hCoeffs->rho62v4 * v));
427  break;
428  case 1:
429  rholm = 1. + v2 * (hCoeffs->rho61v2 + hCoeffs->rho61v3 * v);
430  break;
431  default:
432  XLAL_ERROR(XLAL_EINVAL);
433  break;
434  }
435  break;
436  case 7:
437  switch (m) {
438  case 7:
439  rholm = 1. + v2 * (hCoeffs->rho77v2 + hCoeffs->rho77v3 * v);
440  break;
441  case 6:
442  rholm = 1. + hCoeffs->rho76v2 * v2;
443  break;
444  case 5:
445  rholm = 1. + v2 * (hCoeffs->rho75v2 + hCoeffs->rho75v3 * v);
446  break;
447  case 4:
448  rholm = 1. + hCoeffs->rho74v2 * v2;
449  break;
450  case 3:
451  rholm = 1. + v2 * (hCoeffs->rho73v2 + hCoeffs->rho73v3 * v);
452  break;
453  case 2:
454  rholm = 1. + hCoeffs->rho72v2 * v2;
455  break;
456  case 1:
457  rholm = 1. + v2 * (hCoeffs->rho71v2 + hCoeffs->rho71v3 * v);
458  break;
459  default:
460  XLAL_ERROR(XLAL_EINVAL);
461  break;
462  }
463  break;
464  case 8:
465  switch (m) {
466  case 8:
467  rholm = 1. + hCoeffs->rho88v2 * v2;
468  break;
469  case 7:
470  rholm = 1. + hCoeffs->rho87v2 * v2;
471  break;
472  case 6:
473  rholm = 1. + hCoeffs->rho86v2 * v2;
474  break;
475  case 5:
476  rholm = 1. + hCoeffs->rho85v2 * v2;
477  break;
478  case 4:
479  rholm = 1. + hCoeffs->rho84v2 * v2;
480  break;
481  case 3:
482  rholm = 1. + hCoeffs->rho83v2 * v2;
483  break;
484  case 2:
485  rholm = 1. + hCoeffs->rho82v2 * v2;
486  break;
487  case 1:
488  rholm = 1. + hCoeffs->rho81v2 * v2;
489  break;
490  default:
491  XLAL_ERROR(XLAL_EINVAL);
492  break;
493  }
494  break;
495  default:
496  XLAL_ERROR(XLAL_EINVAL);
497  break;
498  }
499 
500  //debugPK
501  if (debugPK)
502  XLAL_PRINT_INFO("rho_%d_%d = %.12e \n", l, m, rholm);
503  /* Raise rholm to the lth power */
504  rholmPwrl = 1.0;
505  i = l;
506  while (i--) {
507  rholmPwrl *= rholm;
508  }
509  /*
510  * In the equal-mass odd m case, there is no contribution from
511  * nonspin terms, and the only contribution comes from the auxflm
512  * term that is proportional to chiA (asymmetric spins). In this
513  * case, we must ignore the nonspin terms directly, since the leading
514  * term defined by CalculateThisMultipolePrefix in
515  * LALSimIMREOBNewtonianMultipole.c is not zero (see comments there).
516  */
517  if (eta == 0.25 && m % 2) {
518  rholmPwrl = auxflm;
519  } else {
520  rholmPwrl += auxflm;
521  }
522 
523  if (r > 0.0 && debugPK) {
524  XLAL_PRINT_INFO("YP::dynamics variables in waveform: %i, %i, r = %.12e, v = %.12e\n", l, m, r, v);
525  XLAL_PRINT_INFO("rholm^l = %.16e, Tlm = %.16e + i %.16e, \nSlm = %.16e, hNewton = %.16e + i %.16e, delta = %.16e\n", rholmPwrl, Tlm, 0.0, Slm, creal(hNewton), cimag(hNewton), 0.0);
526  }
527  /* Put all factors in Eq. 17 together */
528  *hlm = Tlm * Slm * rholmPwrl;
529  *hlm *= hNewton;
530  /*
531  * if (r > 8.5) { XLAL_PRINT_INFO("YP::FullWave: %.16e,%.16e,
532  * %.16e\n",hlm->re,hlm->im,sqrt(hlm->re*hlm->re+hlm->im*hlm->im)); }
533  */
534  }
535  }
536  return XLAL_SUCCESS;
537 }
538 
539 /**
540  * This function calculates hlm mode factorized-resummed waveform
541  * for given dynamical variables.
542  * Eq. 17 and the entire Appendix of PRD 86, 024011 (2012) + changes
543  * described in the section "Factorized waveforms" of https://dcc.ligo.org/T1400476
544  */
545 static INT4
546 XLALSimIMRSpinEOBGetPrecSpinFactorizedWaveform(
547  COMPLEX16 * restrict hlm, /**< OUTPUT, hlm waveforms */
548  REAL8Vector * restrict values, /**< dyanmical variables: \f$(r,\phi,p_r,p_\phi)\f$ */
549  REAL8Vector * restrict cartvalues, /**< dyanmical variables */
550  const REAL8 v, /**< velocity */
551  const REAL8 Hreal, /**< real Hamiltonian */
552  const INT4 l, /**< l mode index */
553  const INT4 m, /**< m mode index */
554  SpinEOBParams * restrict params /**< Spin EOB parameters */
555 )
556 {
557  int debugPK = 0;
558  /* Status of function calls */
559  INT4 status;
560  INT4 i;
561 
562  REAL8 eta;
563  REAL8 UNUSED r , pp, Omega, v2, Omegav2, vh, vh3, k, hathatk, eulerlogxabs;
564  //pr
565  REAL8 UNUSED rcrossp_x, rcrossp_y, rcrossp_z;
566  REAL8 Slm , deltalm;
567  COMPLEX16 auxflm = 0.0;
568  REAL8 hathatksq4, hathatk4pi, Tlmprefac, Tlmprodfac;
569  COMPLEX16 Tlm, rholmPwrl,rholm;
570  COMPLEX16 hNewton;
571  gsl_sf_result z2;
572 
573  /* Non-Keplerian velocity */
574  REAL8 vPhi , vPhi2;
575 
576  /* Pre-computed coefficients */
577 
578  FacWaveformCoeffs *hCoeffs = params->eobParams->hCoeffs;
579 
580  if (abs(m) > (INT4) l) {
581  XLAL_ERROR(XLAL_EINVAL);
582  }
583  if (m == 0) {
584  XLAL_ERROR(XLAL_EINVAL);
585  }
586  eta = params->eobParams->eta;
587 
588  /* Check our eta was sensible */
589  if (eta > 0.25) {
590  XLALPrintError("XLAL Error - %s: Eta seems to be > 0.25 - this isn't allowed!\n", __func__);
591  XLAL_ERROR(XLAL_EINVAL);
592  }
593  /*
594  * else if ( eta == 0.25 && m % 2 ) { // If m is odd and dM = 0, hLM
595  * will be zero memset( hlm, 0, sizeof( COMPLEX16 ) ); return
596  * XLAL_SUCCESS; }
597  */
598 
599  r = values->data[0];
600  pp = values->data[3];
601 
602  rcrossp_x = cartvalues->data[1] * cartvalues->data[5] - cartvalues->data[2] * cartvalues->data[4];
603  rcrossp_y = cartvalues->data[2] * cartvalues->data[3] - cartvalues->data[0] * cartvalues->data[5];
604  rcrossp_z = cartvalues->data[0] * cartvalues->data[4] - cartvalues->data[1] * cartvalues->data[3];
605 
606  //pp = sqrt(rcrossp_x * rcrossp_x + rcrossp_y * rcrossp_y + rcrossp_z * rcrossp_z);
607 
608  v2 = v * v;
609  Omega = v2 * v;
610  Omegav2 = Omega*v2;
611  vh3 = Hreal * Omega;
612  vh = cbrt(vh3);
613  eulerlogxabs = LAL_GAMMA + log(2.0 * (REAL8) m * v);
614 
615  FacWaveformCoeffs oldCoeffs = *(params->eobParams->hCoeffs); //RC: the function XLALSimIMRSpinPrecEOBNonKeplerCoeff is calculating again the coefficients hCoeffs, for the omega. These are different because
616  // for the dynamics we are using different coefficients for the 21 mode. I store the old coefficients and the I recover them ofter vPhi is computed.
617 
618  /* Calculate the non-Keplerian velocity */
619  if (params->alignedSpins) {
620 
621  vPhi = XLALSimIMRSpinAlignedEOBNonKeplerCoeff(values->data, params);
622 
623  if (XLAL_IS_REAL8_FAIL_NAN(vPhi)) {
624  XLAL_ERROR(XLAL_EFUNC);
625  }
626  vPhi = r * cbrt(vPhi);
627 
628  if (debugPK)
629  XLAL_PRINT_INFO("In XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole, getting rW = %.12e\n",
630  vPhi);
631  vPhi *= Omega;
632  vPhi2 = vPhi * vPhi;
633  } else {
634  vPhi = XLALSimIMRSpinPrecEOBNonKeplerCoeff(cartvalues->data, params);
635  if (XLAL_IS_REAL8_FAIL_NAN(vPhi)) {
636  XLAL_ERROR(XLAL_EFUNC);
637  }
638  vPhi = r * cbrt(vPhi);
639 
640  if (debugPK)
641  XLAL_PRINT_INFO("In XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole, getting rW = %.12e\n",
642  vPhi);
643  vPhi *= Omega;
644  vPhi2 = vPhi * vPhi;
645  }
646  *hCoeffs = oldCoeffs; //RC: Here I recover the old coefficients
647  /*
648  * Calculate the newtonian multipole, 1st term in Eq. 17, given by
649  * Eq. A1
650  */
651  if (debugPK) {
652  XLAL_PRINT_INFO("\nValues inside XLALSimIMRSpinEOBFluxGetSpinFactorizedWaveform:\n");
653  for (i = 0; i < 11; i++)
654  XLAL_PRINT_INFO("values[%d] = %.12e\n", i, cartvalues->data[i]);
655 
656  XLAL_PRINT_INFO("Calculating hNewton, with v = %.12e, vPhi = %.12e, r = %.12e, Phi = %.12e, l = %d, m = %d\n",
657  v, vPhi, r, values->data[1], (UINT4) l, (UINT4) m);
658  }
659  status = XLALSimIMRSpinEOBCalculateNewtonianMultipole(&hNewton,
660  vPhi2, r, cartvalues->data[12] + cartvalues->data[13], (UINT4) l, m, params->eobParams);
661 
662  if (status == XLAL_FAILURE) {
663  XLAL_ERROR(XLAL_EFUNC);
664  }
665  /* Calculate the source term, 2nd term in Eq. 17, given by Eq. A5, Hreal is given by Eq.5 and Heff is in Eq.2 */
666  if (((l + m) % 2) == 0) {
667  Slm = (Hreal * Hreal - 1.) / (2. * eta) + 1.;
668  } else {
669  Slm = v * pp;
670  //Slm = v * sqrt(rcrossp_x * rcrossp_x + rcrossp_y * rcrossp_y + rcrossp_z * rcrossp_z);
671  }
672  if (debugPK)
673  XLAL_PRINT_INFO("In XLALSimIMRSpinEOBGetPrecSpinFactorizedWaveform: Hreal = %e, Slm = %e, eta = %e\n", Hreal, Slm, eta);
674 
675  /*
676  * Calculate the Tail term, 3rd term in Eq. 17,
677  * given by Eq. A6, and Eq. (42) of
678  * http://arxiv.org/pdf/1212.4357.pdf (or PRD 87 084035 (2013))
679  */
680  k = m * Omega;
681  hathatk = Hreal * k;
682 
683  gsl_sf_result lnr1, arg1;
684  XLAL_CALLGSL(status = gsl_sf_lngamma_complex_e(l + 1.0, -2.0 * hathatk, &lnr1, &arg1));
685  if (status != GSL_SUCCESS) {
686  XLALPrintError("XLAL Error - %s: Error in GSL function\n", __func__);
687  XLAL_ERROR(XLAL_EFUNC);
688  }
689  XLAL_CALLGSL(status = gsl_sf_fact_e(l, &z2));
690  if (status != GSL_SUCCESS) {
691  XLALPrintError("XLAL Error - %s: Error in GSL function\n", __func__);
692  XLAL_ERROR(XLAL_EFUNC);
693  }
694  Tlm = cexp((lnr1.val + LAL_PI * hathatk) + I * (
695  arg1.val + 2.0 * hathatk * log(4.0 * k / sqrt(LAL_E))));
696  Tlm /= z2.val;
697 
698 
699  if (debugPK){
700  hathatksq4 = 4. * hathatk * hathatk;
701  hathatk4pi = 4. * LAL_PI * hathatk;
702  /* Calculating the prefactor of Tlm, outside the multiple product */
703  Tlmprefac = sqrt(hathatk4pi / (1. - exp(-hathatk4pi))) / z2.val;
704 
705  /* Calculating the multiple product factor */
706  for (Tlmprodfac = 1., i = 1; i <= l; i++)
707  Tlmprodfac *= (hathatksq4 + (REAL8) i * i);
708 
709  REAL8 Tlmold;
710  Tlmold = Tlmprefac * sqrt(Tlmprodfac);
711  XLAL_PRINT_INFO("Tlm = %e + i%e, |Tlm| = %.16e (should be %.16e)\n", creal(Tlm), cimag(Tlm), cabs(Tlm), Tlmold);
712  }
713 
714  /* Calculate the residue phase and amplitude terms */
715  /*
716  * deltalm is the 4th term in Eq. 17, delta 22 given by Eq. A15,
717  * others
718  */
719  /* rholm is the 5th term in Eq. 17, given by Eqs. A8 - A14 */
720  /*
721  * auxflm is a special part of the 5th term in Eq. 17, given by Eq.
722  * A15
723  */
724  /*
725  * Actual values of the coefficients are defined in the another function
726  * see file LALSimIMRSpinEOBFactorizedWaveformCoefficientsPrec.c
727  */
728  switch (l) {
729  case 2:
730  switch (abs(m)) {
731  case 2:
732  deltalm = vh3 * (hCoeffs->delta22vh3 + vh3 * (hCoeffs->delta22vh6
733  + vh * vh * (hCoeffs->delta22vh9 * vh)))
734  + Omega*(hCoeffs->delta22v5 * v2 + Omega*(hCoeffs->delta22v6 + hCoeffs->delta22v8 *v2));
735  rholm = 1. + v2 * (hCoeffs->rho22v2 + v * (hCoeffs->rho22v3
736  + v * (hCoeffs->rho22v4
737  + v * (hCoeffs->rho22v5 + v * (hCoeffs->rho22v6
738  + hCoeffs->rho22v6l * eulerlogxabs + v * (hCoeffs->rho22v7
739  + v * (hCoeffs->rho22v8 + hCoeffs->rho22v8l * eulerlogxabs
740  + (hCoeffs->rho22v10 + hCoeffs->rho22v10l * eulerlogxabs) * v2)))))));
741  //FIXME
742  if (debugPK){
743  XLAL_PRINT_INFO("PK:: rho22v2 = %.12e, rho22v3 = %.12e, rho22v4 = %.12e,\n rho22v5 = %.16e, rho22v6 = %.16e, rho22v6LOG = %.16e, \n rho22v7 = %.12e, rho22v8 = %.16e, rho22v8LOG = %.16e, \n rho22v10 = %.16e, rho22v10LOG = %.16e\n, rho22v6 = %.12e, rho22v8 = %.12e, rho22v10 = %.12e\n",
744  hCoeffs->rho22v2, hCoeffs->rho22v3, hCoeffs->rho22v4,
745  hCoeffs->rho22v5, hCoeffs->rho22v6, hCoeffs->rho22v6l,
746  hCoeffs->rho22v7, hCoeffs->rho22v8, hCoeffs->rho22v8l,
747  hCoeffs->rho22v10, hCoeffs->rho22v10l,
748  hCoeffs->rho22v6 + hCoeffs->rho22v6l * eulerlogxabs,
749  hCoeffs->rho22v8 + hCoeffs->rho22v8l * eulerlogxabs,
750  hCoeffs->rho22v10 + hCoeffs->rho22v10l * eulerlogxabs);}
751  break;
752  case 1:
753  {
754  deltalm = vh3 * (hCoeffs->delta21vh3 + vh3 * (hCoeffs->delta21vh6
755  + vh * (hCoeffs->delta21vh7 + (hCoeffs->delta21vh9) * vh * vh)))
756  + Omegav2*(hCoeffs->delta21v5 + hCoeffs->delta21v7 * v2);
757  rholm = 1. + v * (hCoeffs->rho21v1
758  + v * (hCoeffs->rho21v2 + v * (hCoeffs->rho21v3 + v * (hCoeffs->rho21v4
759  + v * (hCoeffs->rho21v5 + v * (hCoeffs->rho21v6 + hCoeffs->rho21v6l * eulerlogxabs
760  + v * (hCoeffs->rho21v7 + hCoeffs->rho21v7l * eulerlogxabs
761  + v * (hCoeffs->rho21v8 + hCoeffs->rho21v8l * eulerlogxabs
762  + (hCoeffs->rho21v10 + hCoeffs->rho21v10l * eulerlogxabs) * v2))))))));
763  auxflm = v * (hCoeffs->f21v1 + v2 * (hCoeffs->f21v3 + v * hCoeffs->f21v4 + v2 * (hCoeffs->f21v5 + v * hCoeffs->f21v6 + v2 * hCoeffs->f21v7c)));
764  }
765  break;
766  default:
767  XLAL_ERROR(XLAL_EINVAL);
768  break;
769  }
770  break;
771  case 3:
772  switch (m) {
773  case 3:
774  deltalm = vh3 * (hCoeffs->delta33vh3 + vh3 * (hCoeffs->delta33vh6 + hCoeffs->delta33vh9 * vh3))
775  + hCoeffs->delta33v5 * v * v2 * v2 + hCoeffs->delta33v7 * v2 * v2 * v2 * v;
776  //R.C: delta33v7 is set to 0, whoever is adding it here as a coefficient is evil, TODO: double check that is zero and then remove it
777  //RC: This terms are in Eq.A6 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
778  rholm = 1. + v2 * (hCoeffs->rho33v2 + v * (hCoeffs->rho33v3 + v * (hCoeffs->rho33v4 + v * (hCoeffs->rho33v5 + v * (hCoeffs->rho33v6 +
779  hCoeffs->rho33v6l * eulerlogxabs + v * (hCoeffs->rho33v7 + v * (hCoeffs->rho33v8 + hCoeffs->rho33v8l * eulerlogxabs +
780  v2*(hCoeffs->rho33v10 + hCoeffs->rho33v10l*eulerlogxabs))))))));
781  //RC: This terms are in Eq.A10 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
782  auxflm = v * (v2 * (hCoeffs->f33v3 + v * (hCoeffs->f33v4 + v * (hCoeffs->f33v5 + v * hCoeffs->f33v6)))) + _Complex_I * vh3 * vh3 * hCoeffs->f33vh6;
783  break;
784  case 2:
785  deltalm = vh3 * (hCoeffs->delta32vh3 + vh * (hCoeffs->delta32vh4 + vh * vh * (hCoeffs->delta32vh6
786  + hCoeffs->delta32vh9 * vh3)));
787  rholm = 1. + v * (hCoeffs->rho32v
788  + v * (hCoeffs->rho32v2 + v * (hCoeffs->rho32v3 + v * (hCoeffs->rho32v4 + v * (hCoeffs->rho32v5
789  + v * (hCoeffs->rho32v6 + hCoeffs->rho32v6l * eulerlogxabs
790  + (hCoeffs->rho32v8 + hCoeffs->rho32v8l * eulerlogxabs) * v2))))));
791  break;
792  case 1:
793  deltalm = vh3 * (hCoeffs->delta31vh3 + vh3 * (hCoeffs->delta31vh6
794  + vh * (hCoeffs->delta31vh7 + hCoeffs->delta31vh9 * vh * vh)))
795  + hCoeffs->delta31v5 * v * v2 * v2;
796  rholm = 1. + v2 * (hCoeffs->rho31v2 + v * (hCoeffs->rho31v3 + v * (hCoeffs->rho31v4
797  + v * (hCoeffs->rho31v5 + v * (hCoeffs->rho31v6 + hCoeffs->rho31v6l * eulerlogxabs
798  + v * (hCoeffs->rho31v7 + (hCoeffs->rho31v8 + hCoeffs->rho31v8l * eulerlogxabs) * v))))));
799  auxflm = v * v2 * hCoeffs->f31v3;
800  break;
801  default:
802  XLAL_ERROR(XLAL_EINVAL);
803  break;
804  }
805  break;
806  case 4:
807  switch (m) {
808  case 4:
809  //RC: This terms are in Eq.A15 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
810  deltalm = vh3 * (hCoeffs->delta44vh3 + vh3 * (hCoeffs->delta44vh6 + vh3 * hCoeffs->delta44vh9))
811  + hCoeffs->delta44v5 * v2 * v2 * v;
812  //RC: This terms are in Eq.A8 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
813  rholm = 1. + v2 * (hCoeffs->rho44v2
814  + v * (hCoeffs->rho44v3 + v * (hCoeffs->rho44v4 + v * (hCoeffs->rho44v5 + v * (hCoeffs->rho44v6 + hCoeffs->rho44v6l *
815  eulerlogxabs + v2 *( hCoeffs->rho44v8 + hCoeffs->rho44v8l * eulerlogxabs +v2 * (hCoeffs->rho44v10 + hCoeffs->rho44v10l * eulerlogxabs) ) )))));
816  break;
817  case 3:
818  deltalm = vh3 * (hCoeffs->delta43vh3 + vh * (hCoeffs->delta43vh4
819  + hCoeffs->delta43vh6 * vh * vh));
820  rholm = 1. + v * (hCoeffs->rho43v
821  + v * (hCoeffs->rho43v2
822  + v2 * (hCoeffs->rho43v4 + v * (hCoeffs->rho43v5
823  + (hCoeffs->rho43v6 + hCoeffs->rho43v6l * eulerlogxabs) * v))));
824  auxflm = v * hCoeffs->f43v;
825  break;
826  case 2:
827  deltalm = vh3 * (hCoeffs->delta42vh3 + hCoeffs->delta42vh6 * vh3);
828  rholm = 1. + v2 * (hCoeffs->rho42v2
829  + v * (hCoeffs->rho42v3 + v * (hCoeffs->rho42v4 + v * (hCoeffs->rho42v5
830  + (hCoeffs->rho42v6 + hCoeffs->rho42v6l * eulerlogxabs) * v))));
831  break;
832  case 1:
833  deltalm = vh3 * (hCoeffs->delta41vh3 + vh * (hCoeffs->delta41vh4
834  + hCoeffs->delta41vh6 * vh * vh));
835  rholm = 1. + v * (hCoeffs->rho41v
836  + v * (hCoeffs->rho41v2
837  + v2 * (hCoeffs->rho41v4 + v * (hCoeffs->rho41v5
838  + (hCoeffs->rho41v6 + hCoeffs->rho41v6l * eulerlogxabs) * v))));
839  auxflm = v * hCoeffs->f41v;
840  break;
841  default:
842  XLAL_ERROR(XLAL_EINVAL);
843  break;
844  }
845  break;
846  case 5:
847  switch (m) {
848  case 5:
849  //RC: This terms are in Eq.A16 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
850  deltalm = vh3 *(hCoeffs->delta55vh3 +vh3*(hCoeffs->delta55vh6 +vh3 *(hCoeffs->delta55vh9)))
851  + hCoeffs->delta55v5 * v2 * v2 * v;
852  //RC: This terms are in Eq.A9 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
853  rholm = 1. + v2 * (hCoeffs->rho55v2 + v * (hCoeffs->rho55v3 + v * (hCoeffs->rho55v4 + v * (hCoeffs->rho55v5 +
854  v * (hCoeffs->rho55v6 + hCoeffs->rho55v6l * eulerlogxabs +
855  v2 * (hCoeffs->rho55v8 + hCoeffs->rho55v8l * eulerlogxabs +
856  v2 * (hCoeffs->rho55v10 + hCoeffs->rho55v10l * eulerlogxabs )))))));
857  //RC: This terms are in Eq.A12 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
858  auxflm = v2 * v * (hCoeffs->f55v3 + v * (hCoeffs->f55v4 + v * (hCoeffs->f55v5c) ));
859  break;
860  case 4:
861  deltalm = vh3 * (hCoeffs->delta54vh3 + hCoeffs->delta54vh4 * vh);
862  rholm = 1. + v2 * (hCoeffs->rho54v2 + v * (hCoeffs->rho54v3
863  + hCoeffs->rho54v4 * v));
864  break;
865  case 3:
866  deltalm = hCoeffs->delta53vh3 * vh3;
867  rholm = 1. + v2 * (hCoeffs->rho53v2
868  + v * (hCoeffs->rho53v3 + v * (hCoeffs->rho53v4 + hCoeffs->rho53v5 * v)));
869  break;
870  case 2:
871  deltalm = vh3 * (hCoeffs->delta52vh3 + hCoeffs->delta52vh4 * vh);
872  rholm = 1. + v2 * (hCoeffs->rho52v2 + v * (hCoeffs->rho52v3
873  + hCoeffs->rho52v4 * v));
874  break;
875  case 1:
876  deltalm = hCoeffs->delta51vh3 * vh3;
877  rholm = 1. + v2 * (hCoeffs->rho51v2
878  + v * (hCoeffs->rho51v3 + v * (hCoeffs->rho51v4 + hCoeffs->rho51v5 * v)));
879  break;
880  default:
881  XLAL_ERROR(XLAL_EINVAL);
882  break;
883  }
884  break;
885  case 6:
886  switch (m) {
887  case 6:
888  deltalm = hCoeffs->delta66vh3 * vh3;
889  rholm = 1. + v2 * (hCoeffs->rho66v2 + v * (hCoeffs->rho66v3
890  + hCoeffs->rho66v4 * v));
891  break;
892  case 5:
893  deltalm = hCoeffs->delta65vh3 * vh3;
894  rholm = 1. + v2 * (hCoeffs->rho65v2 + hCoeffs->rho65v3 * v);
895  break;
896  case 4:
897  deltalm = hCoeffs->delta64vh3 * vh3;
898  rholm = 1. + v2 * (hCoeffs->rho64v2 + v * (hCoeffs->rho64v3
899  + hCoeffs->rho64v4 * v));
900  break;
901  case 3:
902  deltalm = hCoeffs->delta63vh3 * vh3;
903  rholm = 1. + v2 * (hCoeffs->rho63v2 + hCoeffs->rho63v3 * v);
904  break;
905  case 2:
906  deltalm = hCoeffs->delta62vh3 * vh3;
907  rholm = 1. + v2 * (hCoeffs->rho62v2 + v * (hCoeffs->rho62v3
908  + hCoeffs->rho62v4 * v));
909  break;
910  case 1:
911  deltalm = hCoeffs->delta61vh3 * vh3;
912  rholm = 1. + v2 * (hCoeffs->rho61v2 + hCoeffs->rho61v3 * v);
913  break;
914  default:
915  XLAL_ERROR(XLAL_EINVAL);
916  break;
917  }
918  break;
919  case 7:
920  switch (m) {
921  case 7:
922  deltalm = hCoeffs->delta77vh3 * vh3;
923  rholm = 1. + v2 * (hCoeffs->rho77v2 + hCoeffs->rho77v3 * v);
924  break;
925  case 6:
926  deltalm = 0.0;
927  rholm = 1. + hCoeffs->rho76v2 * v2;
928  break;
929  case 5:
930  deltalm = hCoeffs->delta75vh3 * vh3;
931  rholm = 1. + v2 * (hCoeffs->rho75v2 + hCoeffs->rho75v3 * v);
932  break;
933  case 4:
934  deltalm = 0.0;
935  rholm = 1. + hCoeffs->rho74v2 * v2;
936  break;
937  case 3:
938  deltalm = hCoeffs->delta73vh3 * vh3;
939  rholm = 1. + v2 * (hCoeffs->rho73v2 + hCoeffs->rho73v3 * v);
940  break;
941  case 2:
942  deltalm = 0.0;
943  rholm = 1. + hCoeffs->rho72v2 * v2;
944  break;
945  case 1:
946  deltalm = hCoeffs->delta71vh3 * vh3;
947  rholm = 1. + v2 * (hCoeffs->rho71v2 + hCoeffs->rho71v3 * v);
948  break;
949  default:
950  XLAL_ERROR(XLAL_EINVAL);
951  break;
952  }
953  break;
954  case 8:
955  switch (m) {
956  case 8:
957  deltalm = 0.0;
958  rholm = 1. + hCoeffs->rho88v2 * v2;
959  break;
960  case 7:
961  deltalm = 0.0;
962  rholm = 1. + hCoeffs->rho87v2 * v2;
963  break;
964  case 6:
965  deltalm = 0.0;
966  rholm = 1. + hCoeffs->rho86v2 * v2;
967  break;
968  case 5:
969  deltalm = 0.0;
970  rholm = 1. + hCoeffs->rho85v2 * v2;
971  break;
972  case 4:
973  deltalm = 0.0;
974  rholm = 1. + hCoeffs->rho84v2 * v2;
975  break;
976  case 3:
977  deltalm = 0.0;
978  rholm = 1. + hCoeffs->rho83v2 * v2;
979  break;
980  case 2:
981  deltalm = 0.0;
982  rholm = 1. + hCoeffs->rho82v2 * v2;
983  break;
984  case 1:
985  deltalm = 0.0;
986  rholm = 1. + hCoeffs->rho81v2 * v2;
987  break;
988  default:
989  XLAL_ERROR(XLAL_EINVAL);
990  break;
991  }
992  break;
993  default:
994  XLAL_ERROR(XLAL_EINVAL);
995  break;
996  }
997 
998  //debugPK
999  if (debugPK)
1000  XLAL_PRINT_INFO("rho_%d_%d = %.12e + %.12e \n", l, m, creal(rholm),cimag(rholm));
1001  if (debugPK)
1002  XLAL_PRINT_INFO("exp(delta_%d_%d) = %.16e + i%.16e (abs=%e)\n", l, m, creal(cexp(I * deltalm)),
1003  cimag(cexp(I * deltalm)), cabs(cexp(I * deltalm)));
1004  /* Raise rholm to the lth power */
1005  rholmPwrl = 1.0;
1006  i = l;
1007  while (i--) {
1008  rholmPwrl *= rholm;
1009  }
1010  /*
1011  * In the equal-mass odd m case, there is no contribution from
1012  * nonspin terms, and the only contribution comes from the auxflm
1013  * term that is proportional to chiA (asymmetric spins). In this
1014  * case, we must ignore the nonspin terms directly, since the leading
1015  * term defined by CalculateThisMultipolePrefix in
1016  * LALSimIMREOBNewtonianMultipole.c is not zero (see comments there).
1017  */
1018  if (eta == 0.25 && m % 2) {
1019  rholmPwrl = auxflm;
1020  } else {
1021  rholmPwrl += auxflm;
1022  }
1023  // if (m==1) {
1024  // printf("f21v1 = %.16f f21v3 = %.16f f21v4 = %.16f f21v5 = %.16f\n", hCoeffs->f21v1, hCoeffs->f21v3, hCoeffs->f21v4, hCoeffs->f21v5);
1025  // }
1026 
1027  if (r > 0.0 && debugPK) {
1028  XLAL_PRINT_INFO("YP::dynamics variables in waveform: %i, %i, r = %.12e, v = %.12e\n", l, m, r, v);
1029  XLAL_PRINT_INFO("rholm^l = %.16e + %.16e, Tlm = %.16e + i %.16e, \nSlm = %.16e, hNewton = %.16e + i %.16e, delta = %.16e\n", creal(rholmPwrl), cimag(rholmPwrl), creal(Tlm), cimag(Tlm), Slm, creal(hNewton), cimag(hNewton), 0.0);
1030  }
1031  /* Put all factors in Eq. 17 together */
1032  *hlm = Tlm * cexp(I * deltalm) * Slm * rholmPwrl;
1033  *hlm *= hNewton;
1034  if (r > 8.5 && debugPK) {
1035  XLAL_PRINT_INFO("YP::FullWave: %.16e,%.16e, %.16e\n", creal(*hlm), cimag(*hlm), cabs(*hlm));
1036  }
1037  return XLAL_SUCCESS;
1038 }
1039 
1040 /**This function calculate the residue amplitude terms */
1041 /* rholm is the 4th term in Eq. 4.5 of https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701] given by Eqs. A6 - A9 */
1042 /* auxflm is a special part of the 4th term in Eq. 4.5 of https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701], given by Eq. A10-A12 */
1043 UNUSED static INT4 XLALSimIMRSpinEOBGetAmplitudeResidualPrec (COMPLEX16 * restrict rholmpwrl, const REAL8 v, const UNUSED REAL8 Hreal, const INT4 modeL, const INT4 modeM, SpinEOBParams * restrict params)
1044  {
1045  INT4 i = 0;
1046  REAL8 eta;
1047  REAL8 eulerlogxabs;
1048  REAL8 rholm;
1049  REAL8 v2 = v * v;
1050  COMPLEX16 auxflm = 0.0;
1051  COMPLEX16 rholmPwrl;
1052  FacWaveformCoeffs *hCoeffs = params->eobParams->hCoeffs;
1053  eulerlogxabs = LAL_GAMMA + log (2.0 * (REAL8) modeM * v);
1054 
1055  if (abs (modeM) > (INT4) modeL)
1056  {
1057  XLAL_ERROR (XLAL_EINVAL);
1058  }
1059  if (modeM == 0)
1060  {
1061  XLAL_ERROR (XLAL_EINVAL);
1062  }
1063  eta = params->eobParams->eta;
1064 
1065  /* Check our eta was sensible */
1066  if (eta > 0.25 && eta < 0.25 +1e-4) {
1067  eta = 0.25;
1068  }
1069  if (eta > 0.25)
1070  {
1071  XLALPrintError
1072  ("XLAL Error - %s: Eta seems to be > 0.25 - this isn't allowed!\n",
1073  __func__);
1074  XLAL_ERROR (XLAL_EINVAL);
1075  }
1076  switch (modeL)
1077  {
1078  case 2:
1079  switch (abs (modeM))
1080  {
1081  case 2:
1082  rholm =
1083  1. + v2 * (hCoeffs->rho22v2 +
1084  v * (hCoeffs->rho22v3 +
1085  v * (hCoeffs->rho22v4 +
1086  v * (hCoeffs->rho22v5 +
1087  v * (hCoeffs->rho22v6 +
1088  hCoeffs->rho22v6l * eulerlogxabs +
1089  v * (hCoeffs->rho22v7 +
1090  v * (hCoeffs->rho22v8 +
1091  hCoeffs->rho22v8l *
1092  eulerlogxabs +
1093  (hCoeffs->rho22v10 +
1094  hCoeffs->rho22v10l *
1095  eulerlogxabs) *
1096  v2)))))));
1097  break;
1098  case 1:
1099  {
1100  rholm =
1101  1. + v * (hCoeffs->rho21v1 +
1102  v * (hCoeffs->rho21v2 +
1103  v * (hCoeffs->rho21v3 +
1104  v * (hCoeffs->rho21v4 +
1105  v * (hCoeffs->rho21v5 +
1106  v * (hCoeffs->rho21v6 +
1107  hCoeffs->rho21v6l *
1108  eulerlogxabs +
1109  v * (hCoeffs->rho21v7 +
1110  hCoeffs->rho21v7l *
1111  eulerlogxabs +
1112  v * (hCoeffs->rho21v8 +
1113  hCoeffs->rho21v8l *
1114  eulerlogxabs +
1115  (hCoeffs->
1116  rho21v10 +
1117  hCoeffs->
1118  rho21v10l *
1119  eulerlogxabs) *
1120  v2))))))));
1121  auxflm = v * (hCoeffs->f21v1 + v2 * (hCoeffs->f21v3 + v * hCoeffs->f21v4 + v2 * (hCoeffs->f21v5 + v * hCoeffs->f21v6)));
1122  }
1123  break;
1124  default:
1125  XLAL_ERROR (XLAL_EINVAL);
1126  break;
1127  }
1128  break;
1129  case 3:
1130  switch (modeM)
1131  {
1132  case 3:
1133  rholm =
1134  1. + v2 * (hCoeffs->rho33v2 +
1135  v * (hCoeffs->rho33v3 +
1136  v * (hCoeffs->rho33v4 +
1137  v * (hCoeffs->rho33v5 +
1138  v * (hCoeffs->rho33v6 +
1139  hCoeffs->rho33v6l * eulerlogxabs +
1140  v * (hCoeffs->rho33v7 +
1141  v * (hCoeffs->rho33v8 +
1142  hCoeffs->rho33v8l *
1143  eulerlogxabs)))))));
1144  auxflm = v * v2 * hCoeffs->f33v3;
1145  break;
1146  case 2:
1147  rholm =
1148  1. + v * (hCoeffs->rho32v +
1149  v * (hCoeffs->rho32v2 +
1150  v * (hCoeffs->rho32v3 +
1151  v * (hCoeffs->rho32v4 +
1152  v * (hCoeffs->rho32v5 +
1153  v * (hCoeffs->rho32v6 +
1154  hCoeffs->rho32v6l *
1155  eulerlogxabs +
1156  (hCoeffs->rho32v8 +
1157  hCoeffs->rho32v8l *
1158  eulerlogxabs) * v2))))));
1159  break;
1160  case 1:
1161  rholm =
1162  1. + v2 * (hCoeffs->rho31v2 +
1163  v * (hCoeffs->rho31v3 +
1164  v * (hCoeffs->rho31v4 +
1165  v * (hCoeffs->rho31v5 +
1166  v * (hCoeffs->rho31v6 +
1167  hCoeffs->rho31v6l * eulerlogxabs +
1168  v * (hCoeffs->rho31v7 +
1169  (hCoeffs->rho31v8 +
1170  hCoeffs->rho31v8l *
1171  eulerlogxabs) * v))))));
1172  auxflm = v * v2 * hCoeffs->f31v3;
1173  break;
1174  default:
1175  XLAL_ERROR (XLAL_EINVAL);
1176  break;
1177  }
1178  break;
1179  case 4:
1180  switch (modeM)
1181  {
1182  case 4:
1183  rholm = 1. + v2 * (hCoeffs->rho44v2
1184  + v * (hCoeffs->rho44v3 + v * (hCoeffs->rho44v4
1185  +
1186  v *
1187  (hCoeffs->
1188  rho44v5 +
1189  (hCoeffs->
1190  rho44v6 +
1191  hCoeffs->
1192  rho44v6l *
1193  eulerlogxabs) *
1194  v))));
1195 
1196  break;
1197  case 3:
1198  rholm =
1199  1. + v * (hCoeffs->rho43v +
1200  v * (hCoeffs->rho43v2 +
1201  v2 * (hCoeffs->rho43v4 +
1202  v * (hCoeffs->rho43v5 +
1203  (hCoeffs->rho43v6 +
1204  hCoeffs->rho43v6l * eulerlogxabs) *
1205  v))));
1206  auxflm = v * hCoeffs->f43v;
1207  break;
1208  case 2:
1209  rholm = 1. + v2 * (hCoeffs->rho42v2
1210  + v * (hCoeffs->rho42v3 +
1211  v * (hCoeffs->rho42v4 +
1212  v * (hCoeffs->rho42v5 +
1213  (hCoeffs->rho42v6 +
1214  hCoeffs->rho42v6l *
1215  eulerlogxabs) * v))));
1216  break;
1217  case 1:
1218  rholm =
1219  1. + v * (hCoeffs->rho41v +
1220  v * (hCoeffs->rho41v2 +
1221  v2 * (hCoeffs->rho41v4 +
1222  v * (hCoeffs->rho41v5 +
1223  (hCoeffs->rho41v6 +
1224  hCoeffs->rho41v6l * eulerlogxabs) *
1225  v))));
1226  auxflm = v * hCoeffs->f41v;
1227  break;
1228  default:
1229  XLAL_ERROR (XLAL_EINVAL);
1230  break;
1231  }
1232  break;
1233  case 5:
1234  switch (modeM)
1235  {
1236  case 5:
1237  //RC: This terms are in Eq.A9 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
1238  rholm =
1239  1. + v2 * (hCoeffs->rho55v2 +
1240  v * (hCoeffs->rho55v3 +
1241  v * (hCoeffs->rho55v4 +
1242  v * (hCoeffs->rho55v5 +
1243  v * (hCoeffs->rho55v6 + hCoeffs->rho55v6l * eulerlogxabs +
1244  v2 * (hCoeffs->rho55v8 + hCoeffs->rho55v8l * eulerlogxabs +
1245  v2 * (hCoeffs->rho55v10 + hCoeffs->rho55v10l * eulerlogxabs )))))));
1246  //RC: This terms are in Eq.A12 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 [arXiv:1803.10701]
1247  auxflm = v2 * v * (hCoeffs->f55v3 + v * (hCoeffs->f55v4));
1248 
1249  break;
1250  case 4:
1251  rholm = 1. + v2 * (hCoeffs->rho54v2 + v * (hCoeffs->rho54v3
1252  + hCoeffs->rho54v4 * v));
1253  break;
1254  case 3:
1255  rholm = 1. + v2 * (hCoeffs->rho53v2
1256  + v * (hCoeffs->rho53v3 +
1257  v * (hCoeffs->rho53v4 +
1258  hCoeffs->rho53v5 * v)));
1259  break;
1260  case 2:
1261  rholm = 1. + v2 * (hCoeffs->rho52v2 + v * (hCoeffs->rho52v3
1262  + hCoeffs->rho52v4 * v));
1263  break;
1264  case 1:
1265  rholm = 1. + v2 * (hCoeffs->rho51v2
1266  + v * (hCoeffs->rho51v3 +
1267  v * (hCoeffs->rho51v4 +
1268  hCoeffs->rho51v5 * v)));
1269  break;
1270  default:
1271  XLAL_ERROR (XLAL_EINVAL);
1272  break;
1273  }
1274  break;
1275  case 6:
1276  switch (modeM)
1277  {
1278  case 6:
1279  rholm = 1. + v2 * (hCoeffs->rho66v2 + v * (hCoeffs->rho66v3
1280  + hCoeffs->rho66v4 * v));
1281  break;
1282  case 5:
1283  rholm = 1. + v2 * (hCoeffs->rho65v2 + hCoeffs->rho65v3 * v);
1284  break;
1285  case 4:
1286  rholm = 1. + v2 * (hCoeffs->rho64v2 + v * (hCoeffs->rho64v3
1287  + hCoeffs->rho64v4 * v));
1288  break;
1289  case 3:
1290  rholm = 1. + v2 * (hCoeffs->rho63v2 + hCoeffs->rho63v3 * v);
1291  break;
1292  case 2:
1293  rholm = 1. + v2 * (hCoeffs->rho62v2 + v * (hCoeffs->rho62v3
1294  + hCoeffs->rho62v4 * v));
1295  break;
1296  case 1:
1297  rholm = 1. + v2 * (hCoeffs->rho61v2 + hCoeffs->rho61v3 * v);
1298  break;
1299  default:
1300  XLAL_ERROR (XLAL_EINVAL);
1301  break;
1302  }
1303  break;
1304  case 7:
1305  switch (modeM)
1306  {
1307  case 7:
1308  rholm = 1. + v2 * (hCoeffs->rho77v2 + hCoeffs->rho77v3 * v);
1309  break;
1310  case 6:
1311  rholm = 1. + hCoeffs->rho76v2 * v2;
1312  break;
1313  case 5:
1314  rholm = 1. + v2 * (hCoeffs->rho75v2 + hCoeffs->rho75v3 * v);
1315  break;
1316  case 4:
1317  rholm = 1. + hCoeffs->rho74v2 * v2;
1318  break;
1319  case 3:
1320  rholm = 1. + v2 * (hCoeffs->rho73v2 + hCoeffs->rho73v3 * v);
1321  break;
1322  case 2:
1323  rholm = 1. + hCoeffs->rho72v2 * v2;
1324  break;
1325  case 1:
1326  rholm = 1. + v2 * (hCoeffs->rho71v2 + hCoeffs->rho71v3 * v);
1327  break;
1328  default:
1329  XLAL_ERROR (XLAL_EINVAL);
1330  break;
1331  }
1332  break;
1333  case 8:
1334  switch (modeM)
1335  {
1336  case 8:
1337  rholm = 1. + hCoeffs->rho88v2 * v2;
1338  break;
1339  case 7:
1340  rholm = 1. + hCoeffs->rho87v2 * v2;
1341  break;
1342  case 6:
1343  rholm = 1. + hCoeffs->rho86v2 * v2;
1344  break;
1345  case 5:
1346  rholm = 1. + hCoeffs->rho85v2 * v2;
1347  break;
1348  case 4:
1349  rholm = 1. + hCoeffs->rho84v2 * v2;
1350  break;
1351  case 3:
1352  rholm = 1. + hCoeffs->rho83v2 * v2;
1353  break;
1354  case 2:
1355  rholm = 1. + hCoeffs->rho82v2 * v2;
1356  break;
1357  case 1:
1358  rholm = 1. + hCoeffs->rho81v2 * v2;
1359  break;
1360  default:
1361  XLAL_ERROR (XLAL_EINVAL);
1362  break;
1363  }
1364  break;
1365  default:
1366  XLAL_ERROR (XLAL_EINVAL);
1367  break;
1368  }
1369 
1370  /* Raise rholm to the lth power */
1371  rholmPwrl = 1.0;
1372  i = modeL;
1373  while (i--)
1374  {
1375  rholmPwrl *= rholm;
1376  }
1377  /* In the equal-mass odd m case, there is no contribution from nonspin terms,
1378  * and the only contribution comes from the auxflm term that is proportional to chiA (asymmetric spins).
1379  * In this case, we must ignore the nonspin terms directly, since the leading term defined by
1380  * CalculateThisMultipolePrefix in LALSimIMREOBNewtonianMultipole.c is not zero (see comments there).
1381  */
1382  if (eta == 0.25 && modeM % 2)
1383  {
1384  rholmPwrl = auxflm;
1385  }
1386  else
1387  {
1388  rholmPwrl += auxflm;
1389  }
1390  *rholmpwrl = rholmPwrl;
1391 
1392 return XLAL_SUCCESS;
1393 }
1394 
1395 /**
1396 * This function calculate the calibration parameter for the amplitude of
1397 * the factorized-resummed waveforms in the case of the 21 and the 55 mode
1398 */
1399 UNUSED static INT4 XLALSimIMREOBCalcCalibCoefficientHigherModesPrec (
1400  SpinEOBParams * restrict UNUSED params, /**Output **/
1401  const UINT4 modeL, /*<< Mode index L */
1402  const UINT4 modeM, /*<< Mode index M */
1403  SEOBdynamics *seobdynamics, /*<< Dynamics vector */
1404  const REAL8 timeorb, /*<< Time of the peak of the orbital frequency */
1405  const REAL8 m1, /**Component mass 1 */
1406  const REAL8 m2, /**Component mass 2 */
1407  const REAL8 UNUSED deltaT /**<< Sampling interval */
1408  )
1409  {
1410  /* Status of function calls */
1411  UINT4 i;
1412  UINT4 debugRC = 0;
1413 
1414  /** Physical quantities */
1415  REAL8Vector *timeVec = NULL, *hLMdivrholmVec = NULL;
1416  REAL8Vector polarDynamics, values;
1417  REAL8 valuesdata[14] = {0.};
1418  REAL8 polarDynamicsdata[4] = {0.};
1419  polarDynamics.length = 4;
1420  values.length = 14;
1421  polarDynamics.data = polarDynamicsdata;
1422  values.data = valuesdata;
1423  REAL8 omegaAttachmentPoint, hLMdivrholmAttachmentPoint, rholmNRAttachmentPoint,rholmBeforeCalibAttachmentPoint;
1424  REAL8 v;
1425  COMPLEX16Vector *hLMVec = NULL, *rholmpwrlVec = NULL;
1426  REAL8Vector *rholmpwrlVecReal = NULL;
1427 
1428  REAL8Vector orbOmegaVec, HamVec;
1429  UINT4 retLen = seobdynamics->length;
1430 
1431  /** Find the vaulues of the final spins */
1432  REAL8 mtot = m1+m2;
1433  REAL8 eta = params->eobParams->eta;
1434  REAL8 s1dotZ = 0, s2dotZ = 0, chi1dotZ = 0, chi2dotZ = 0;
1435  REAL8 chiS = 0;
1436  REAL8 chiA = 0;
1437  REAL8 tplspin = 0;
1438  REAL8 spin1z_omegaPeak = params->spin1z_omegaPeak;
1439  REAL8 spin2z_omegaPeak = params->spin2z_omegaPeak;
1440  REAL8 chiS_omegaPeak = 0.5*(spin1z_omegaPeak+ spin2z_omegaPeak);
1441  REAL8 chiA_omegaPeak = 0.5*(spin1z_omegaPeak-spin2z_omegaPeak);
1442 
1443  /* Create dynamical arrays */
1444 
1445  orbOmegaVec.length = HamVec.length = retLen;
1446 
1447  hLMVec = XLALCreateCOMPLEX16Vector (retLen);
1448  rholmpwrlVec = XLALCreateCOMPLEX16Vector (retLen);
1449  timeVec = XLALCreateREAL8Vector (retLen);
1450  hLMdivrholmVec = XLALCreateREAL8Vector (retLen);
1451  rholmpwrlVecReal = XLALCreateREAL8Vector (retLen);
1452  if (!hLMVec || !rholmpwrlVec|| !timeVec|| !rholmpwrlVec || !rholmpwrlVecReal)
1453  XLAL_ERROR (XLAL_ENOMEM);
1454 
1455  orbOmegaVec.data = seobdynamics->omegaVec;
1456  HamVec.data = seobdynamics->hamVec;
1457 
1458 
1459  /* Stuff for interpolating function */
1460  gsl_spline *spline = NULL;
1461  gsl_interp_accel *acc = NULL;
1462 
1463  /* The calibration parameter is only used for 21 and 55 modes, if you try to use this function for other modes, you get an error */
1464 
1465  if (!((modeL == 2 && modeM == 1) || (modeL == 5 && modeM == 5)))
1466  {
1467  XLALPrintError ("Mode %d,%d is not supported by this function.\n", modeL,
1468  modeM);
1469  XLAL_ERROR (XLAL_EINVAL);
1470  }
1471  /* Populate time vector as necessary */
1472  for (i = 0; i < timeVec->length; i++)
1473  {
1474  timeVec->data[i] = i * deltaT;
1475  }
1476  /**Initializing stuff for interpolation */
1477  spline = gsl_spline_alloc (gsl_interp_cspline, orbOmegaVec.length);
1478  acc = gsl_interp_accel_alloc ();
1479  /* Calculation of the frequency at the attachment point */
1480  REAL8 timewavePeak = timeorb-XLALSimIMREOBGetNRSpinPeakDeltaTv4 (modeL, modeM, m1, m2, spin1z_omegaPeak, spin2z_omegaPeak);
1481  gsl_spline_init (spline, timeVec->data, orbOmegaVec.data, orbOmegaVec.length);
1482  gsl_interp_accel_reset (acc);
1483  omegaAttachmentPoint = gsl_spline_eval (spline, timewavePeak, acc);
1484 
1485  /** Calculation and interpolation at the matching point of rho_lm^l + f_lm */
1486  for(i=0; i<orbOmegaVec.length; i++){
1487  for (UINT4 j=0; j<14; j++) {
1488  values.data[j] = seobdynamics->array->data[i + (j+1)*retLen];
1489  }
1490  s1dotZ = seobdynamics->s1dotZVec[i];
1491  s2dotZ = seobdynamics->s2dotZVec[i];
1492  polarDynamics.data[0] = seobdynamics->polarrVec[i];
1493  polarDynamics.data[1] = seobdynamics->polarphiVec[i];
1494  polarDynamics.data[2] = seobdynamics->polarprVec[i];
1495  polarDynamics.data[3] = seobdynamics->polarpphiVec[i];
1496  chi1dotZ = s1dotZ * mtot*mtot / (m1*m1);
1497  chi2dotZ = s2dotZ * mtot*mtot / (m2*m2);
1498  chiS = 0.5*(chi1dotZ+chi2dotZ);
1499  chiA = 0.5*(chi1dotZ-chi2dotZ);
1500  tplspin = (1.-2.*eta) * chiS + (m1 - m2)/(m1 + m2) * chiA;
1501  v = cbrt (orbOmegaVec.data[i]);
1502  if ( XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients( params->eobParams->hCoeffs, m1, m2, eta, tplspin, chiS, chiA, v4Pwave ) == XLAL_FAILURE )
1503  {
1504  XLALPrintError("XLAL Error - %s: failure in XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients at step %d of the loop.\n", __func__, i);
1505  XLAL_ERROR( XLAL_EFUNC );
1506  }
1507  if ( XLALSimIMRSpinEOBGetPrecSpinFactorizedWaveform( &(hLMVec->data[i]), &polarDynamics, &values, v, HamVec.data[i], modeL, modeM, params ) == XLAL_FAILURE )
1508  {
1509  XLALPrintError("XLAL Error - %s: failure in XLALSimIMRSpinEOBGetPrecSpinFactorizedWaveform at step %d of the loop.\n", __func__, i);
1510  XLAL_ERROR( XLAL_EDOM );
1511  }
1512  if (XLALSimIMRSpinEOBGetAmplitudeResidualPrec (&(rholmpwrlVec->data[i]), v, HamVec.data[i], modeL, modeM, params) == XLAL_FAILURE)
1513  //RC: for the 21 and 55 mode rholmpwrlVec is always real. This is not true for the 33 mode. For this reason, in order to make this function general, we use a complex variable for it.
1514  {
1515  /* TODO: Clean-up */
1516  XLAL_ERROR (XLAL_EFUNC);
1517  }
1518  rholmpwrlVecReal->data[i] = (REAL8)creal(rholmpwrlVec->data[i]);
1519  hLMdivrholmVec->data[i] = ((REAL8)cabs(hLMVec->data[i]))/fabs(rholmpwrlVecReal->data[i]);
1520  }
1521  gsl_spline_init (spline, timeVec->data, hLMdivrholmVec->data, hLMdivrholmVec->length);
1522  gsl_interp_accel_reset (acc);
1523  hLMdivrholmAttachmentPoint = gsl_spline_eval (spline, timewavePeak, acc);
1524  REAL8 nra = 0;
1525  nra = XLALSimIMREOBGetNRSpinPeakAmplitudeV4 (modeL, modeM, m1, m2,chiS_omegaPeak, chiA_omegaPeak);
1526  if((fabs(nra/eta)< 3e-2) && ((modeL == 2) && (modeM == 1))){
1527  //R.C.: safeguard to avoid the 21 mode to go to 0
1528  nra = GSL_SIGN(nra)*eta*3e-2;
1529  }
1530  if((fabs(nra/eta)< 1e-4) && ((modeL == 5) && (modeM == 5))){
1531  //R.C.: safeguard to avoid the 55 mode to go to 0
1532  nra = GSL_SIGN(nra)*eta*1e-4;
1533  }
1534  rholmNRAttachmentPoint = nra/hLMdivrholmAttachmentPoint;
1535  gsl_spline_init (spline, timeVec->data, rholmpwrlVecReal->data, rholmpwrlVecReal->length);
1536  gsl_interp_accel_reset (acc);
1537  /* rho_lm^l + f_lm before the calibration parameter is set */
1538  rholmBeforeCalibAttachmentPoint = gsl_spline_eval (spline, timewavePeak, acc);
1539  if(debugRC==1){
1540  FILE *timeomegapeak = NULL;
1541  timeomegapeak = fopen ("timeomegapeak.dat", "w");
1542  fprintf(timeomegapeak, "%.16f\n", timewavePeak);
1543  fclose(timeomegapeak);
1544  }
1545  if ((modeL == 2) && (modeM == 1))
1546  {
1547  /* Here we compute ((rho_lm^l + f_lm + CalPar*omega^7/3)_NR - (rho_lm^l + f_lm)_EOB)/omega^7/3 to get CalPar.
1548  The factor rholmpwrlVecReal->data[0])/cabs(rholmpwrlVecReal->data[0]) is used to know the sign of the function (rho_lm^l + f_lm + CalPar*omega^7/3)_NR which is computed as absolute value */
1549  params->cal21 = (rholmNRAttachmentPoint - rholmBeforeCalibAttachmentPoint)/(pow(omegaAttachmentPoint,7./3.));
1550  //printf("cal21 = %.16f\n", params->cal21);
1551  }
1552  if ((modeL == 5) && (modeM == 5))
1553  {
1554  /* Here we compute ((rho_lm^l + f_lm + CalPar*omega^7/3)_NR - (rho_lm^l + f_lm)_EOB)/omega^7/3 to get CalPar.
1555  The factor rholmpwrlVecReal->data[0])/cabs(rholmpwrlVecReal->data[0]) is used to know the sign of the function (rho_lm^l + f_lm + CalPar*omega^7/3)_NR which is computed as absolute value */
1556  params->cal55 = (rholmNRAttachmentPoint - rholmBeforeCalibAttachmentPoint)/(pow(omegaAttachmentPoint,5./3.));
1557  //printf("params->cal55 = %.16f\n",params->cal55);
1558  }
1559 
1560 
1561  gsl_spline_free (spline);
1562  gsl_interp_accel_free (acc);
1563  XLALDestroyREAL8Vector (hLMdivrholmVec);
1564  XLALDestroyCOMPLEX16Vector (hLMVec);
1565  XLALDestroyREAL8Vector (timeVec);
1566  XLALDestroyCOMPLEX16Vector (rholmpwrlVec);
1567  XLALDestroyREAL8Vector (rholmpwrlVecReal);
1568 
1569 
1570 
1571  return XLAL_SUCCESS;}
1572 
1573 #endif /* _LALSIMIMRSPINPRECEOBFACTORIZEDWAVEFORM */
XLALSimIMREOBGetNRSpinPeakAmplitudeV4
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakAmplitudeV4(INT4 modeL, INT4 modeM, REAL8 m1, REAL8 m2, REAL8 chiS, REAL8 chiA)
Peak amplitude predicted by fitting NR results The coefficients for the mode (2,2) are in Eq.
Definition: LALSimIMREOBNQCCorrection.c:929
XLALSimIMREOBGetNRSpinPeakDeltaTv4
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakDeltaTv4(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED m1, REAL8 UNUSED m2, REAL8 UNUSED chi1, REAL8 UNUSED chi2)
The time difference between the orbital peak and the peak amplitude of the mode in question (currentl...
Definition: LALSimIMREOBNQCCorrection.c:1699
XLALSimIMRSpinEOBCalculateNewtonianMultipole
static UNUSED int XLALSimIMRSpinEOBCalculateNewtonianMultipole(COMPLEX16 *multipole, REAL8 x, UNUSED REAL8 r, REAL8 phi, UINT4 l, INT4 m, EOBParams *params)
This function calculates the Newtonian multipole part of the factorized waveform for the SEOBNRv1 mod...
Definition: LALSimIMREOBNewtonianMultipole.c:153
XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole
static UNUSED int XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole(COMPLEX16 *multipole, REAL8 x, UNUSED REAL8 r, UNUSED REAL8 phi, UINT4 l, INT4 m, EOBParams *params)
This function calculates the Newtonian multipole part of the factorized waveform for the SEOBNRv1 mod...
Definition: LALSimIMREOBNewtonianMultipole.c:191
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
XLALSimIMRSpinEOBGetAmplitudeResidualPrec
static UNUSED INT4 XLALSimIMRSpinEOBGetAmplitudeResidualPrec(COMPLEX16 *restrict rholmpwrl, const REAL8 v, const REAL8 Hreal, const INT4 modeL, const INT4 modeM, SpinEOBParams *restrict params)
XLALSimIMRSpinEOBGetPrecSpinFactorizedWaveform
static INT4 XLALSimIMRSpinEOBGetPrecSpinFactorizedWaveform(COMPLEX16 *restrict hlm, REAL8Vector *restrict values, REAL8Vector *restrict cartvalues, const REAL8 v, const REAL8 Hreal, const INT4 l, const INT4 m, SpinEOBParams *restrict params)
This function calculates hlm mode factorized-resummed waveform for given dynamical variables.
Definition: LALSimIMRSpinEOBFactorizedWaveformPrec.c:546
XLALSimIMREOBCalcCalibCoefficientHigherModesPrec
static UNUSED INT4 XLALSimIMREOBCalcCalibCoefficientHigherModesPrec(SpinEOBParams *restrict UNUSED params, const UINT4 modeL, const UINT4 modeM, SEOBdynamics *seobdynamics, const REAL8 timeorb, const REAL8 m1, const REAL8 m2, const REAL8 UNUSED deltaT)
This function calculate the calibration parameter for the amplitude of the factorized-resummed wavefo...
Definition: LALSimIMRSpinEOBFactorizedWaveformPrec.c:1399
XLALSimIMRSpinEOBFluxGetPrecSpinFactorizedWaveform
static INT4 XLALSimIMRSpinEOBFluxGetPrecSpinFactorizedWaveform(COMPLEX16 *restrict hlmTab, REAL8Vector *restrict values, REAL8Vector *restrict cartvalues, const REAL8 v, const REAL8 Hreal, const INT4 lMax, SpinEOBParams *restrict params)
FOR PRECESSING EOB This function calculates hlm mode factorized-resummed waveform for given dynamical...
Definition: LALSimIMRSpinEOBFactorizedWaveformPrec.c:115
XLALSimIMRSpinAlignedEOBNonKeplerCoeff
static REAL8 XLALSimIMRSpinAlignedEOBNonKeplerCoeff(const REAL8 values[], SpinEOBParams *funcParams)
Function to calculate the non-Keplerian coefficient for the spin-aligned EOB model.
Definition: LALSimIMRSpinEOBHamiltonian.c:1579
XLALSimIMRSpinPrecEOBNonKeplerCoeff
static UNUSED REAL8 XLALSimIMRSpinPrecEOBNonKeplerCoeff(const REAL8 values[], SpinEOBParams *funcParams)
Function to calculate the non-Keplerian coefficient for the PRECESSING EOB model.
Definition: LALSimIMRSpinEOBHamiltonianPrec.c:988
z2
static REAL8 UNUSED z2(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
Definition: LALSimInspiralOptimizedCoefficientsEccentricityFD.c:527
Hreal
REAL8 Hreal
Definition: SEOBNRv2_opt_.h:164
fprintf
#define fprintf
XLAL_CALLGSL
#define XLAL_CALLGSL(statement)
l
int l
Definition: bh_qnmode.c:135
i
double i
Definition: bh_ringdown.c:118
e
double e
Definition: bh_ringdown.c:117
pp
const double pp
Definition: exact_derivatives-Hreal.c:35
LAL_E
#define LAL_E
LAL_PI
#define LAL_PI
LAL_GAMMA
#define LAL_GAMMA
COMPLEX16
double complex COMPLEX16
REAL8
double REAL8
UINT4
uint32_t UINT4
INT4
int32_t INT4
r
static const INT4 r
m
static const INT4 m
XLALCreateCOMPLEX16Vector
COMPLEX16Vector * XLALCreateCOMPLEX16Vector(UINT4 length)
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLALDestroyCOMPLEX16Vector
void XLALDestroyCOMPLEX16Vector(COMPLEX16Vector *vector)
XLAL_PRINT_INFO
#define XLAL_PRINT_INFO(...)
XLAL_ERROR
#define XLAL_ERROR(...)
XLALPrintError
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_IS_REAL8_FAIL_NAN
#define XLAL_IS_REAL8_FAIL_NAN(val)
XLAL_ENOMEM
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EFUNC
XLAL_EDOM
XLAL_EDOM
XLAL_EINVAL
XLAL_EINVAL
XLAL_FAILURE
XLAL_FAILURE
status
string status
COMPLEX16Vector::data
COMPLEX16 * data
FacWaveformCoeffs
Structure containing the coefficients for calculating the factorized waveform.
Definition: LALSimIMREOBNRv2.h:97
REAL8Array::data
REAL8 * data
REAL8Vector::data
REAL8 * data
SEOBdynamics
Structure the EOB dynamics for precessing waveforms.
Definition: LALSimIMREOBNRv2.h:583
SEOBdynamics::array
REAL8Array * array
Definition: LALSimIMREOBNRv2.h:585
SEOBdynamics::polarpphiVec
REAL8 * polarpphiVec
Definition: LALSimIMREOBNRv2.h:615
SpinEOBParams
Parameters for the spinning EOB model.
Definition: LALSimIMRSpinEOB.h:116
params
Definition: burst.c:245
deltaT
double deltaT
Definition: unicorn.c:24