LALSimulation  5.4.0.1-fe68b98
LALSimIMRSpinEOBFactorizedWaveform.c
Go to the documentation of this file.
1 /*
2 * Copyright (C) 2010 Craig Robinson, Yi Pan, Prayush Kumar, Andrea Taracchini
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with with program; see the file COPYING. If not, write to the
16 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17 * MA 02110-1301 USA
18 */
19 
20 
21 /**
22  * \author Craig Robinson, Yi Pan, Prayush Kumar, Andrea Taracchini
23  *
24  * \brief Function to compute the factorized waveform as uses in the SEOBNRv1 model.
25  * Waveform expressions are given by
26  * Taracchini et al. ( PRD 86, 024011 (2012), arXiv 1202.0790 ).
27  * All equation numbers in this file refer to equations of this paper,
28  * unless otherwise specified.
29  * Coefficients of the so-called "deltalm" terms are given by
30  * Damour et al. PRD 79, 064004 (2009) [arXiv:0811.2069] and Pan et al. PRD 83, 064003 (2011) [arXiv:1006.0431],
31  * henceforth DIN and PBFRT.
32  *
33  * Functions to compute the factorized waveform for the purpose of computing the
34  * flux, i.e. returning only the absolute value of the multipoles. The tail term
35  * Tlm is used in its resummed form, given by Eq. (42) of Damour, Nagar and
36  * Bernuzzi, PRD 87 084035 (2013) [arxiv:1212.4357], called DNB here.
37  *
38  */
39 
40 #ifndef _LALSIMIMRSPINEOBFACTORIZEDWAVEFORM_C
41 #define _LALSIMIMRSPINEOBFACTORIZEDWAVEFORM_C
42 
43 #include <complex.h>
44 
45 #include <gsl/gsl_vector.h>
46 #include <gsl/gsl_matrix.h>
47 #include <gsl/gsl_spline.h>
48 #include <gsl/gsl_linalg.h>
49 
50 #include <lal/LALSimInspiral.h>
51 #include <lal/LALSimIMR.h>
52 
53 #include "LALSimIMREOBNRv2.h"
54 #include "LALSimIMRSpinEOB.h"
55 
56 #include "LALSimIMREOBNewtonianMultipole.c"
57 #include "LALSimIMRSpinEOBAuxFuncs.c"
58 #include "LALSimIMREOBNQCCorrection.c"
59 #include "LALSimIMRSpinEOBHamiltonian.c"
60 
61 /* OPTIMIZED */
62 #include "LALSimIMRSpinEOBHamiltonianOptimized.c"
63 /* END OPTIMIZED */
64 
65 /*------------------------------------------------------------------------------------------
66  *
67  * Prototypes of functions defined in this code.
68  *
69  *------------------------------------------------------------------------------------------
70  */
71 
72 
73 UNUSED static INT4 XLALSimIMRSpinEOBGetSpinFactorizedWaveform (COMPLEX16 * restrict hlm, REAL8Vector * restrict values, const REAL8 v, const REAL8 Hreal, const INT4 l, const INT4 m, SpinEOBParams * restrict params, INT4 use_optimized_v2
74  /**< Use optimized v2? */
75  );
76 
77 UNUSED static INT4 XLALSimIMRSpinEOBGetAmplitudeResidual (COMPLEX16 * restrict rholmpwrl, const REAL8 v, const REAL8 Hreal, const INT4 modeL, const INT4 modeM, SpinEOBParams * restrict params);
78 
79 static INT4 XLALSimIMRSpinEOBFluxGetSpinFactorizedWaveform (COMPLEX16 * restrict hlm, REAL8Vector * restrict values, const REAL8 v, const REAL8 Hreal, const INT4 l, const INT4 m, SpinEOBParams * restrict params, INT4 use_optimized_v2,
80  /**< Use optimized v2? */
81  REAL8 * vPhil2m2);
82 
83 static int XLALSimIMREOBCalcSpinFacWaveformCoefficients (FacWaveformCoeffs *
84  const coeffs,
85  SpinEOBParams * restrict params,
86  const REAL8 m1,
87  const REAL8 m2,
88  const REAL8 eta,
89  const REAL8 a,
90  const REAL8 chiS,
91  const REAL8 chiA,
92  const UINT4 SpinAlignedEOBversion);
93 
94 UNUSED static int XLALSimIMREOBCalcCalibCoefficientHigherModes(FacWaveformCoeffs *
95  const coeffs,/**< Output */
96  SpinEOBParams * restrict UNUSED params, /**< Additional structure which might be useful in the future */
97  const UINT4 modeL, /**< Mode index L */
98  const UINT4 modeM, /**< Mode index M */
99  REAL8Vector * restrict phaseVec, /**< Orbital phase, function of time */
100  REAL8Vector * restrict rVec, /**< Position-vector, function of time */
101  REAL8Vector * restrict prVec, /**< Momentum vector, function of time */
102  REAL8Vector * restrict orbOmegaVec, /**< Orbital frequency, func of time */
103  REAL8Vector * restrict HamVec, /**< Hamiltonian vector, func of time */
104  REAL8Vector * restrict pPhiVec, /**< Angular momentum vector, func of time */
105  const REAL8 timePeak, /**< Time of peak orbital frequency */
106  const REAL8 m1, /**< Component mass 1 */
107  const REAL8 m2, /**< Component mass 2 */
108  const REAL8 UNUSED chiS, /**< Symmetric dimensionless spin combination */
109  const REAL8 UNUSED chiA, /**< Asymmetric dimensionless spin combination */
110  const REAL8 UNUSED deltaT /**< Sampling interval */
111  );
112 
113 
114 
115 /*------------------------------------------------------------------------------------------
116  *
117  * Defintions of functions.
118  *
119  *------------------------------------------------------------------------------------------
120  */
121 /**
122  * Function to compute PN corrections to the tidal term that enters the (2,2) amplitude
123  */
124 static REAL8 XLALTEOBbeta221( REAL8 X /**<< NS mass */ ) {
125  return (-202. + 560.*X - 340.*X*X + 45.*X*X*X) / (42.*(3. - 2.*X));
126 }
127 
128 /**
129  * Function to compute the dynamical enhancement of the quadrupolar Love number
130  * that enters the tidal corrections to the mode amplitudes
131  */
132 static REAL8 XLALSimIMRTEOBk2effMode (
133  REAL8 Omega, /**<< Orbital frequency */
134  REAL8 k2TidaleffHam, /**<< Dynamical enhancement of k2Tidal that enters the Hamiltonian */
135  REAL8 omega02Tidal, /**<< f-mode frequency */
136  REAL8 XCompanion /**<< Mass of NS companion divided by M */
137 )
138 {
139  return ((-1. + k2TidaleffHam)*omega02Tidal*omega02Tidal + 6.*k2TidaleffHam*XCompanion*Omega*Omega) / (1. + 2*XCompanion) / (3.*Omega*Omega);
140 }
141 
142 /**
143  * Function to compute the tidal correction to the mode amplitudes
144  */
145 static COMPLEX16 XLALhTidal(
146  INT4 l, /**<< Mode index */
147  INT4 m, /**<< Mode index */
148  REAL8 phase, /**<< Orbutal phase */
149  REAL8 v, /**<< Orbital speed^2 */
150  REAL8 eta, /**<< Symmetric mass ratio */
151  TidalEOBParams *tidal1, /**<< Tidal parameters of body 1 */
152  TidalEOBParams *tidal2 /**<< Tidal parameters of body 2 */
153 )
154 {
155  COMPLEX16 hNewtonTidal = 0;
156  COMPLEX16 hhatTidal = 0;
157  REAL8 v2 = v*v;
158  REAL8 Omega = v*v*v;
159  REAL8 v10 = v2*v2*v2*v2*v2;
160  REAL8 X1 = tidal1->mByM;
161  REAL8 X2 = tidal2->mByM;
162  REAL8 lambda1 = tidal1->lambda2Tidal;
163  REAL8 lambda2 = tidal2->lambda2Tidal;
164  REAL8 omega02Tidal1 = tidal1->omega02Tidal;
165  REAL8 omega02Tidal2 = tidal2->omega02Tidal;
166  REAL8 k2Tidal1effHam = 0.;
167  REAL8 k2Tidal2effHam = 0.;
168  REAL8 k2Tidal1eff = 0.;
169  REAL8 k2Tidal2eff = 0.;
170  REAL8 u = 1./pow(Omega,-2./3);
171  if ( lambda1 != 0.) {
172  k2Tidal1effHam = XLALSimIMRTEOBkleff(2, u, eta, tidal1);
173  k2Tidal1eff = XLALSimIMRTEOBk2effMode (Omega, k2Tidal1effHam, omega02Tidal1, X2);
174  }
175  if ( lambda2 != 0.) {
176  k2Tidal2effHam = XLALSimIMRTEOBkleff(2, u, eta, tidal2);
177  k2Tidal2eff = XLALSimIMRTEOBk2effMode (Omega, k2Tidal2effHam, omega02Tidal2, X1);
178  }
179  REAL8 q = X2/X1;
180  switch (l) {
181  case 2:
182  switch (m) {
183  case 2:
184  hNewtonTidal = -8. * v2 * sqrt(LAL_PI/5.);
185  hhatTidal = (3.*q*lambda1*(X1/X2 + 3.)*(1. + XLALTEOBbeta221(X1)*v2)*k2Tidal1eff
186  + 3./q*lambda2*(X2/X1 + 3.)*(1. + XLALTEOBbeta221(X2)*v2)*k2Tidal2eff) * v10;
187  break;
188  case 1:
189  hNewtonTidal = 8./3. * I * v*v2 * sqrt(LAL_PI/5.);
190  hhatTidal = (-3*q*lambda1*(4.5 - 6.*X1) - (-3./q*lambda2*(4.5 - 6.*X2))) * v10;
191  break;
192  default:
193  return 0.;
194  break;
195  }
196  break;
197  case 3:
198  hhatTidal = -18.*(X2*q*lambda1 - X1/q*lambda2) * v10;
199  switch (m) {
200  case 3:
201  hNewtonTidal = -3. * I * v*v2 * sqrt(6.*LAL_PI/7.);
202  break;
203  case 2:
204  return 0.;
205  break;
206  case 1:
207  hNewtonTidal = 1./3. * I * v*v2 * sqrt(2.*LAL_PI/35.);
208  break;
209  case 0:
210  return 0.;
211  break;
212  default:
213  return 0.;
214  break;
215  }
216  break;
217  default:
218  return 0.;
219  break;
220  }
221  return eta*hNewtonTidal*hhatTidal*( cos(-m*phase) + I * sin(-m*phase) );
222 }
223 
224 /**This function calculate the residue amplitude terms */
225 /* rholm is the 4th term in Eq. 4.5 of https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028 given by Eqs. A6 - A9 */
226 /* 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, given by Eq. A10-A12 */
227 UNUSED static INT4 XLALSimIMRSpinEOBGetAmplitudeResidual (COMPLEX16 * restrict rholmpwrl, const REAL8 v, const REAL8 Hreal, const INT4 modeL, const INT4 modeM, SpinEOBParams * restrict params)
228  {
229  INT4 use_hm = 0;
230  INT4 i = 0;
231  REAL8 eta, vh3, Omega;
232  REAL8 eulerlogxabs;
233  REAL8 rholm;
234  REAL8 v2 = v * v;
235  COMPLEX16 auxflm = 0.0;
236  COMPLEX16 rholmPwrl;
237  FacWaveformCoeffs *hCoeffs = params->eobParams->hCoeffs;
238  Omega = v2 * v;
239  vh3 = Hreal * Omega;
240  eulerlogxabs = LAL_GAMMA + log (2.0 * (REAL8) modeM * v);
241 
242  if (abs (modeM) > modeL)
243  {
244  XLAL_ERROR (XLAL_EINVAL);
245  }
246  if (modeM == 0)
247  {
248  XLAL_ERROR (XLAL_EINVAL);
249  }
250  use_hm = params->use_hm;
251  eta = params->eobParams->eta;
252 
253  /* Check our eta was sensible */
254  if (eta > 0.25 && eta < 0.25 +1e-4) {
255  eta = 0.25;
256  }
257  if (eta > 0.25)
258  {
259  XLALPrintError
260  ("XLAL Error - %s: Eta = %.16f seems to be > 0.25 - this isn't allowed!\n",
261  __func__,eta);
262  XLAL_ERROR (XLAL_EINVAL);
263  }
264  switch (modeL)
265  {
266  case 2:
267  switch (abs (modeM))
268  {
269  case 2:
270  rholm =
271  1. + v2 * (hCoeffs->rho22v2 +
272  v * (hCoeffs->rho22v3 +
273  v * (hCoeffs->rho22v4 +
274  v * (hCoeffs->rho22v5 +
275  v * (hCoeffs->rho22v6 +
276  hCoeffs->rho22v6l * eulerlogxabs +
277  v * (hCoeffs->rho22v7 +
278  v * (hCoeffs->rho22v8 +
279  hCoeffs->rho22v8l *
280  eulerlogxabs +
281  (hCoeffs->rho22v10 +
282  hCoeffs->rho22v10l *
283  eulerlogxabs) *
284  v2)))))));
285  break;
286  case 1:
287  {
288  rholm =
289  1. + v * (hCoeffs->rho21v1 +
290  v * (hCoeffs->rho21v2 +
291  v * (hCoeffs->rho21v3 +
292  v * (hCoeffs->rho21v4 +
293  v * (hCoeffs->rho21v5 +
294  v * (hCoeffs->rho21v6 +
295  hCoeffs->rho21v6l *
296  eulerlogxabs +
297  v * (hCoeffs->rho21v7 +
298  hCoeffs->rho21v7l *
299  eulerlogxabs +
300  v * (hCoeffs->rho21v8 +
301  hCoeffs->rho21v8l *
302  eulerlogxabs +
303  v2*(hCoeffs->
304  rho21v10 +
305  hCoeffs->
306  rho21v10l *
307  eulerlogxabs)))))))));
308  auxflm = v * hCoeffs->f21v1;
309  if (use_hm){
310  //RC: This terms are in Eq.A11 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
311  auxflm = v * (hCoeffs->f21v1 + v2 * (hCoeffs->f21v3 + v * hCoeffs->f21v4 + v2 * (hCoeffs->f21v5 + v * hCoeffs->f21v6)));
312  }
313  }
314  break;
315  default:
316  XLAL_ERROR (XLAL_EINVAL);
317  break;
318  }
319  break;
320  case 3:
321  switch (modeM)
322  {
323  case 3:
324  if(use_hm){
325  //RC: This terms are in Eq.A6 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
326  rholm =
327  1. + v2 * (hCoeffs->rho33v2 +
328  v * (hCoeffs->rho33v3 +
329  v * (hCoeffs->rho33v4 +
330  v * (hCoeffs->rho33v5 +
331  v * (hCoeffs->rho33v6 +
332  hCoeffs->rho33v6l * eulerlogxabs +
333  v * (hCoeffs->rho33v7 +
334  v * (hCoeffs->rho33v8 +
335  hCoeffs->rho33v8l *
336  eulerlogxabs +
337  v2*(hCoeffs->rho33v10 + hCoeffs->rho33v10l*eulerlogxabs))))))));
338  //RC: This terms are in Eq.A10 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
339  auxflm = v * (v2 * (hCoeffs->f33v3 + v * (hCoeffs->f33v4 + v * (hCoeffs->f33v5 + v * hCoeffs->f33v6)))) + _Complex_I * vh3 * vh3 * hCoeffs->f33vh6;
340  }
341  else
342  {
343  rholm =
344  1. + v2 * (hCoeffs->rho33v2 +
345  v * (hCoeffs->rho33v3 +
346  v * (hCoeffs->rho33v4 +
347  v * (hCoeffs->rho33v5 +
348  v * (hCoeffs->rho33v6 +
349  hCoeffs->rho33v6l * eulerlogxabs +
350  v * (hCoeffs->rho33v7 +
351  v * (hCoeffs->rho33v8 +
352  hCoeffs->rho33v8l *
353  eulerlogxabs)))))));
354  auxflm = v * v2 * hCoeffs->f33v3;
355  }
356  break;
357  case 2:
358  rholm =
359  1. + v * (hCoeffs->rho32v +
360  v * (hCoeffs->rho32v2 +
361  v * (hCoeffs->rho32v3 +
362  v * (hCoeffs->rho32v4 +
363  v * (hCoeffs->rho32v5 +
364  v * (hCoeffs->rho32v6 +
365  hCoeffs->rho32v6l *
366  eulerlogxabs +
367  (hCoeffs->rho32v8 +
368  hCoeffs->rho32v8l *
369  eulerlogxabs) * v2))))));
370  break;
371  case 1:
372  rholm =
373  1. + v2 * (hCoeffs->rho31v2 +
374  v * (hCoeffs->rho31v3 +
375  v * (hCoeffs->rho31v4 +
376  v * (hCoeffs->rho31v5 +
377  v * (hCoeffs->rho31v6 +
378  hCoeffs->rho31v6l * eulerlogxabs +
379  v * (hCoeffs->rho31v7 +
380  (hCoeffs->rho31v8 +
381  hCoeffs->rho31v8l *
382  eulerlogxabs) * v))))));
383  auxflm = v * v2 * hCoeffs->f31v3;
384  break;
385  default:
386  XLAL_ERROR (XLAL_EINVAL);
387  break;
388  }
389  break;
390  case 4:
391  switch (modeM)
392  {
393  case 4:
394  if(use_hm){
395  //RC: This terms are in Eq.A8 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
396  rholm = 1. + v2 * (hCoeffs->rho44v2
397  + v * (hCoeffs->rho44v3 + v * (hCoeffs->rho44v4
398  +
399  v *
400  (hCoeffs->
401  rho44v5 +
402  v * (hCoeffs->
403  rho44v6 +
404  hCoeffs->
405  rho44v6l *
406  eulerlogxabs +
407  v2 *( hCoeffs->rho44v8 + hCoeffs->rho44v8l * eulerlogxabs
408  +v2 * (hCoeffs->rho44v10 + hCoeffs->rho44v10l * eulerlogxabs) ) )))));
409  }
410  else{
411  rholm = 1. + v2 * (hCoeffs->rho44v2
412  + v * (hCoeffs->rho44v3 + v * (hCoeffs->rho44v4
413  +
414  v *
415  (hCoeffs->
416  rho44v5 +
417  (hCoeffs->
418  rho44v6 +
419  hCoeffs->
420  rho44v6l *
421  eulerlogxabs) *
422  v))));
423  }
424 
425  break;
426  case 3:
427  rholm =
428  1. + v * (hCoeffs->rho43v +
429  v * (hCoeffs->rho43v2 +
430  v2 * (hCoeffs->rho43v4 +
431  v * (hCoeffs->rho43v5 +
432  (hCoeffs->rho43v6 +
433  hCoeffs->rho43v6l * eulerlogxabs) *
434  v))));
435  auxflm = v * hCoeffs->f43v;
436  break;
437  case 2:
438  rholm = 1. + v2 * (hCoeffs->rho42v2
439  + v * (hCoeffs->rho42v3 +
440  v * (hCoeffs->rho42v4 +
441  v * (hCoeffs->rho42v5 +
442  (hCoeffs->rho42v6 +
443  hCoeffs->rho42v6l *
444  eulerlogxabs) * v))));
445  break;
446  case 1:
447  rholm =
448  1. + v * (hCoeffs->rho41v +
449  v * (hCoeffs->rho41v2 +
450  v2 * (hCoeffs->rho41v4 +
451  v * (hCoeffs->rho41v5 +
452  (hCoeffs->rho41v6 +
453  hCoeffs->rho41v6l * eulerlogxabs) *
454  v))));
455  auxflm = v * hCoeffs->f41v;
456  break;
457  default:
458  XLAL_ERROR (XLAL_EINVAL);
459  break;
460  }
461  break;
462  case 5:
463  switch (modeM)
464  {
465  case 5:
466  if(use_hm){
467  //RC: This terms are in Eq.A9 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
468  rholm =
469  1. + v2 * (hCoeffs->rho55v2 +
470  v * (hCoeffs->rho55v3 +
471  v * (hCoeffs->rho55v4 +
472  v * (hCoeffs->rho55v5 +
473  v * (hCoeffs->rho55v6 + hCoeffs->rho55v6l * eulerlogxabs +
474  v2 * (hCoeffs->rho55v8 + hCoeffs->rho55v8l * eulerlogxabs +
475  v2 * (hCoeffs->rho55v10 + hCoeffs->rho55v10l * eulerlogxabs )))))));
476  //RC: This terms are in Eq.A12 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
477  auxflm = v2 * v * (hCoeffs->f55v3 + v * (hCoeffs->f55v4));
478  }
479  else
480  {
481  rholm =
482  1. + v2 * (hCoeffs->rho55v2 +
483  v * (hCoeffs->rho55v3 +
484  v * (hCoeffs->rho55v4 +
485  v * (hCoeffs->rho55v5 +
486  hCoeffs->rho55v6 * v))));
487  }
488  break;
489  case 4:
490  rholm = 1. + v2 * (hCoeffs->rho54v2 + v * (hCoeffs->rho54v3
491  + hCoeffs->rho54v4 * v));
492  break;
493  case 3:
494  rholm = 1. + v2 * (hCoeffs->rho53v2
495  + v * (hCoeffs->rho53v3 +
496  v * (hCoeffs->rho53v4 +
497  hCoeffs->rho53v5 * v)));
498  break;
499  case 2:
500  rholm = 1. + v2 * (hCoeffs->rho52v2 + v * (hCoeffs->rho52v3
501  + hCoeffs->rho52v4 * v));
502  break;
503  case 1:
504  rholm = 1. + v2 * (hCoeffs->rho51v2
505  + v * (hCoeffs->rho51v3 +
506  v * (hCoeffs->rho51v4 +
507  hCoeffs->rho51v5 * v)));
508  break;
509  default:
510  XLAL_ERROR (XLAL_EINVAL);
511  break;
512  }
513  break;
514  case 6:
515  switch (modeM)
516  {
517  case 6:
518  rholm = 1. + v2 * (hCoeffs->rho66v2 + v * (hCoeffs->rho66v3
519  + hCoeffs->rho66v4 * v));
520  break;
521  case 5:
522  rholm = 1. + v2 * (hCoeffs->rho65v2 + hCoeffs->rho65v3 * v);
523  break;
524  case 4:
525  rholm = 1. + v2 * (hCoeffs->rho64v2 + v * (hCoeffs->rho64v3
526  + hCoeffs->rho64v4 * v));
527  break;
528  case 3:
529  rholm = 1. + v2 * (hCoeffs->rho63v2 + hCoeffs->rho63v3 * v);
530  break;
531  case 2:
532  rholm = 1. + v2 * (hCoeffs->rho62v2 + v * (hCoeffs->rho62v3
533  + hCoeffs->rho62v4 * v));
534  break;
535  case 1:
536  rholm = 1. + v2 * (hCoeffs->rho61v2 + hCoeffs->rho61v3 * v);
537  break;
538  default:
539  XLAL_ERROR (XLAL_EINVAL);
540  break;
541  }
542  break;
543  case 7:
544  switch (modeM)
545  {
546  case 7:
547  rholm = 1. + v2 * (hCoeffs->rho77v2 + hCoeffs->rho77v3 * v);
548  break;
549  case 6:
550  rholm = 1. + hCoeffs->rho76v2 * v2;
551  break;
552  case 5:
553  rholm = 1. + v2 * (hCoeffs->rho75v2 + hCoeffs->rho75v3 * v);
554  break;
555  case 4:
556  rholm = 1. + hCoeffs->rho74v2 * v2;
557  break;
558  case 3:
559  rholm = 1. + v2 * (hCoeffs->rho73v2 + hCoeffs->rho73v3 * v);
560  break;
561  case 2:
562  rholm = 1. + hCoeffs->rho72v2 * v2;
563  break;
564  case 1:
565  rholm = 1. + v2 * (hCoeffs->rho71v2 + hCoeffs->rho71v3 * v);
566  break;
567  default:
568  XLAL_ERROR (XLAL_EINVAL);
569  break;
570  }
571  break;
572  case 8:
573  switch (modeM)
574  {
575  case 8:
576  rholm = 1. + hCoeffs->rho88v2 * v2;
577  break;
578  case 7:
579  rholm = 1. + hCoeffs->rho87v2 * v2;
580  break;
581  case 6:
582  rholm = 1. + hCoeffs->rho86v2 * v2;
583  break;
584  case 5:
585  rholm = 1. + hCoeffs->rho85v2 * v2;
586  break;
587  case 4:
588  rholm = 1. + hCoeffs->rho84v2 * v2;
589  break;
590  case 3:
591  rholm = 1. + hCoeffs->rho83v2 * v2;
592  break;
593  case 2:
594  rholm = 1. + hCoeffs->rho82v2 * v2;
595  break;
596  case 1:
597  rholm = 1. + hCoeffs->rho81v2 * v2;
598  break;
599  default:
600  XLAL_ERROR (XLAL_EINVAL);
601  break;
602  }
603  break;
604  default:
605  XLAL_ERROR (XLAL_EINVAL);
606  break;
607  }
608 
609  /* Raise rholm to the lth power */
610  rholmPwrl = 1.0;
611  for(i = 0 ; i < modeL ; i++) rholmPwrl *= rholm;
612  /* In the equal-mass odd m case, there is no contribution from nonspin terms,
613  * and the only contribution comes from the auxflm term that is proportional to chiA (asymmetric spins).
614  * In this case, we must ignore the nonspin terms directly, since the leading term defined by
615  * CalculateThisMultipolePrefix in LALSimIMREOBNewtonianMultipole.c is not zero (see comments there).
616  */
617  if (eta == 0.25 && modeM % 2)
618  {
619  rholmPwrl = auxflm;
620  }
621  else
622  {
623  rholmPwrl += auxflm;
624  }
625  *rholmpwrl = rholmPwrl;
626 return XLAL_SUCCESS;
627 }
628 
629 
630 
631 /**
632 * This function calculate the calibration parameter for the amplitude of
633 * the factorized-resummed waveforms in the case of the 21 and the 55 mode
634 */
635 static int XLALSimIMREOBCalcCalibCoefficientHigherModes (FacWaveformCoeffs *
636  const coeffs, /** Output **/
637  SpinEOBParams * restrict UNUSED params, /**Additional structure which might be useful in the future **/
638  const UINT4 modeL, /** Mode index L **/
639  const UINT4 modeM, /**Mode index M **/
640  REAL8Vector * restrict phaseVec, /** Orbital phase, function of time **/
641  REAL8Vector * restrict rVec, /**<< Position-vector, function of time */
642  REAL8Vector * restrict prVec, /**<< Momentum vector, function of time */
643  REAL8Vector * restrict orbOmegaVec, /**<< Orbital frequency, func of time */
644  REAL8Vector * restrict HamVec, /**<< Hamiltonian vector, func of time */
645  REAL8Vector * restrict pPhiVec, /**<< Angular momentum vector, func of time */
646  const REAL8 timewavePeak, /** Time of the attachment between inspiral and merger-RD */
647  const REAL8 m1, /**Component mass 1 */
648  const REAL8 m2, /**Component mass 2 */
649  const REAL8 UNUSED chiS, /** Symmetric dimensionless spin combination */
650  const REAL8 UNUSED chiA, /**Asymmetric dimensionless spin combination */
651  const REAL8 UNUSED deltaT /**<< Sampling interval */
652  )
653  {
654  /* Status of function calls */
655  UINT4 i;
656  REAL8 eta = m1*m2/((m1+m2)*(m1+m2));
657 
658  /** Physical quantities */
659  REAL8Vector *timeVec = NULL, *hLMdivrholmVec;
660  REAL8 omegaAttachmentPoint, hLMdivrholmAttachmentPoint, rholmNRAttachmentPoint,rholmBeforeCalibAttachmentPoint;
661  REAL8 v;
662  COMPLEX16Vector *hLMVec = NULL, *rholmpwrlVec = NULL;
663  REAL8Vector *values = NULL, *rholmpwrlVecReal = NULL;
664 
665  /* 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 */
666 
667  if (!((modeL == 2 && modeM == 1) || (modeL == 5 && modeM == 5)))
668  {
669  XLALPrintError ("Mode %d,%d is not supported by this function.\n", modeL,
670  modeM);
671  XLAL_ERROR (XLAL_EINVAL);
672  }
673 
674 
675  /* Create dynamical arrays */
676 
677  values = XLALCreateREAL8Vector (4);
678  memset (values->data, 0, values->length * sizeof (REAL8));
679  hLMVec = XLALCreateCOMPLEX16Vector (phaseVec->length);
680  rholmpwrlVec = XLALCreateCOMPLEX16Vector (phaseVec->length);
681  timeVec = XLALCreateREAL8Vector (phaseVec->length);
682  hLMdivrholmVec = XLALCreateREAL8Vector (phaseVec->length);
683  rholmpwrlVecReal = XLALCreateREAL8Vector (phaseVec->length);
684  if (!hLMVec || !rholmpwrlVec|| !values|| !timeVec|| !hLMdivrholmVec || !rholmpwrlVecReal){
685  if(hLMVec) XLALDestroyCOMPLEX16Vector(hLMVec);
686  if(rholmpwrlVec) XLALDestroyCOMPLEX16Vector(rholmpwrlVec);
687  if(values) XLALDestroyREAL8Vector(values);
688  if(timeVec) XLALDestroyREAL8Vector(timeVec);
689  if(hLMdivrholmVec) XLALDestroyREAL8Vector(hLMdivrholmVec);
690  if(rholmpwrlVecReal) XLALDestroyREAL8Vector(rholmpwrlVecReal);
691  XLAL_ERROR (XLAL_ENOMEM);
692  }
693 
694 
695 
696 
697  /* Stuff for interpolating function */
698  gsl_spline *spline = NULL;
699  gsl_interp_accel *acc = NULL;
700 
701  /* Populate time vector as necessary */
702  for (i = 0; i < timeVec->length; i++)
703  {
704  timeVec->data[i] = i * deltaT;
705  }
706  /**Initializing stuff for interpolation */
707  spline = gsl_spline_alloc (gsl_interp_cspline, phaseVec->length);
708  acc = gsl_interp_accel_alloc ();
709  /* Calculation of the frequency at the attachment point */
710  gsl_spline_init (spline, timeVec->data, phaseVec->data, phaseVec->length);
711  gsl_interp_accel_reset (acc);
712  omegaAttachmentPoint = gsl_spline_eval_deriv (spline, timewavePeak, acc);
713 
714  /** Calculation and interpolation at the matching point of rho_lm^l + f_lm */
715  for(i=0; i<phaseVec->length; i++){
716  values->data[0] = rVec->data[i];
717  values->data[1] = phaseVec->data[i];
718  values->data[2] = prVec->data[i];
719  values->data[3] = pPhiVec->data[i];
720  v = cbrt (orbOmegaVec->data[i]);
721  if (XLALSimIMRSpinEOBGetSpinFactorizedWaveform
722  (&(hLMVec->data[i]), values, v, HamVec->data[i], modeL, modeM, params,
723  params->postAdiabaticFlag) == XLAL_FAILURE)
724  {
725  if(hLMVec) XLALDestroyCOMPLEX16Vector(hLMVec);
726  if(rholmpwrlVec) XLALDestroyCOMPLEX16Vector(rholmpwrlVec);
727  if(values) XLALDestroyREAL8Vector(values);
728  if(timeVec) XLALDestroyREAL8Vector(timeVec);
729  if(hLMdivrholmVec) XLALDestroyREAL8Vector(hLMdivrholmVec);
730  if(rholmpwrlVecReal) XLALDestroyREAL8Vector(rholmpwrlVecReal);
731  gsl_spline_free (spline);
732  gsl_interp_accel_free (acc);
733  XLAL_ERROR (XLAL_EFUNC);
734  }
735  if (XLALSimIMRSpinEOBGetAmplitudeResidual (&(rholmpwrlVec->data[i]), v, HamVec->data[i], modeL, modeM, params) == XLAL_FAILURE)
736  //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.
737  {
738  if(hLMVec) XLALDestroyCOMPLEX16Vector(hLMVec);
739  if(rholmpwrlVec) XLALDestroyCOMPLEX16Vector(rholmpwrlVec);
740  if(values) XLALDestroyREAL8Vector(values);
741  if(timeVec) XLALDestroyREAL8Vector(timeVec);
742  if(hLMdivrholmVec) XLALDestroyREAL8Vector(hLMdivrholmVec);
743  if(rholmpwrlVecReal) XLALDestroyREAL8Vector(rholmpwrlVecReal);
744  gsl_spline_free (spline);
745  gsl_interp_accel_free (acc);
746  XLAL_ERROR (XLAL_EFUNC);
747  }
748  rholmpwrlVecReal->data[i] = (REAL8)creal(rholmpwrlVec->data[i]);
749  hLMdivrholmVec->data[i] = ((REAL8)cabs(hLMVec->data[i]))/fabs(rholmpwrlVecReal->data[i]);
750  }
751 
752  gsl_spline_init (spline, timeVec->data, hLMdivrholmVec->data, hLMdivrholmVec->length);
753  gsl_interp_accel_reset (acc);
754  hLMdivrholmAttachmentPoint = gsl_spline_eval (spline, timewavePeak, acc);
755  REAL8 nra = 0;
756  nra = XLALSimIMREOBGetNRSpinPeakAmplitudeV4 (modeL, modeM, m1, m2,chiS, chiA);
757  if((fabs(nra/eta)< 3e-2) && ((modeL == 2) && (modeM == 1))){
758  //R.C.: safeguard to avoid the 21 mode to go to 0
759  nra = GSL_SIGN(nra)*eta*3e-2;
760  }
761  if((fabs(nra/eta)< 1e-4) && ((modeL == 5) && (modeM == 5))){
762  //R.C.: safeguard to avoid the 55 mode to go to 0
763  nra = GSL_SIGN(nra)*eta*1e-4;
764  }
765  rholmNRAttachmentPoint = nra/hLMdivrholmAttachmentPoint;
766  gsl_spline_init (spline, timeVec->data, rholmpwrlVecReal->data, rholmpwrlVecReal->length);
767  gsl_interp_accel_reset (acc);
768  /* rho_lm^l + f_lm before the calibration parameter is set */
769  rholmBeforeCalibAttachmentPoint = gsl_spline_eval (spline, timewavePeak, acc);
770  if ((modeL == 2) && (modeM == 1))
771  {
772  /* 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. */
773  coeffs->f21v7c = (rholmNRAttachmentPoint - rholmBeforeCalibAttachmentPoint)/(pow(omegaAttachmentPoint,7./3.));
774  }
775  if ((modeL == 5) && (modeM == 5)){
776 
777  coeffs->f55v5c = (rholmNRAttachmentPoint - rholmBeforeCalibAttachmentPoint)/(pow(omegaAttachmentPoint,5./3.));
778  }
779 
780 
781 
782  gsl_spline_free (spline);
783  gsl_interp_accel_free (acc);
784  XLALDestroyREAL8Vector (hLMdivrholmVec);
785  XLALDestroyCOMPLEX16Vector (hLMVec);
786  XLALDestroyREAL8Vector (values);
787  XLALDestroyREAL8Vector (timeVec);
788  XLALDestroyCOMPLEX16Vector (rholmpwrlVec);
789  XLALDestroyREAL8Vector (rholmpwrlVecReal);
790 
791 
792  return XLAL_SUCCESS;}
793 
794 /**
795  * This function calculates hlm mode factorized-resummed waveform
796  * for given dynamical variables. This is optimized for flux calculation,
797  * by ignoring complex arguments and keeping only absolute values.
798  * Changes:
799  * (i) Complex Argument of Tlm not exponentiated.
800  * (ii) \f$exp(\ii deltalm)\f$ set to 1.
801  * Eq. 17 and the entire Appendix of the paper https://journals.aps.org/prd/abstract/10.1103/PhysRevD.86.024011.
802  */
803 static INT4
804 XLALSimIMRSpinEOBFluxGetSpinFactorizedWaveform (COMPLEX16 * restrict hlm,
805  /**< OUTPUT, hlm waveforms */
806  REAL8Vector * restrict values,
807  /**< dyanmical variables */
808  const REAL8 v,
809  /**< velocity */
810  const REAL8 Hreal,
811  /**< real Hamiltonian */
812  const INT4 l,
813  /**< l mode index */
814  const INT4 m,
815  /**< m mode index */
816  SpinEOBParams * restrict params,
817  /**< Spin EOB parameters */
818  INT4 use_optimized_v2,
819  /**< Use optimized v2? */
820  REAL8 * vPhiInput
821  /**< Reuse vPhi computed outside for efficiency! */
822  )
823 {
824  /* Status of function calls */
825  INT4 status;
826  INT4 i;
827 
828  REAL8 eta;
829  REAL8 r, pp, Omega, v2, /*vh, vh3, */ k, hathatk, eulerlogxabs; //pr
830  REAL8 Slm, rholm, rholmPwrl;
831  REAL8 auxflm = 0.0;
832  REAL8 hathatksq4, hathatk4pi, Tlmprefac, Tlmprodfac;
833  REAL8 Tlm;
834  COMPLEX16 hNewton;
835  gsl_sf_result z2;
836 
837  /* Non-Keplerian velocity */
838  REAL8 vPhi, vPhi2;
839 
840  /* Pre-computed coefficients */
841  FacWaveformCoeffs *hCoeffs = params->eobParams->hCoeffs;
842 
843  if (abs (m) > (INT4) l)
844  {
845  XLAL_ERROR (XLAL_EINVAL);
846  }
847  if (m == 0)
848  {
849  XLAL_ERROR (XLAL_EINVAL);
850  }
851 
852  eta = params->eobParams->eta;
853 
854 
855  /* Check our eta was sensible */
856  if (eta > 0.25 && eta < 0.25 +1e-4) {
857  eta = 0.25;
858  }
859  if (eta > 0.25)
860  {
861  XLALPrintError
862  ("XLAL Error - %s: Eta = %.16f seems to be > 0.25 - this isn't allowed!\n",
863  __func__,eta);
864  XLAL_ERROR (XLAL_EINVAL);
865  }
866  /*else if ( eta == 0.25 && m % 2 )
867  {
868  // If m is odd and dM = 0, hLM will be zero
869  memset( hlm, 0, sizeof( COMPLEX16 ) );
870  return XLAL_SUCCESS;
871  } */
872 
873  r = values->data[0];
874  //pr = values->data[2];
875  pp = values->data[3];
876 
877  v2 = v * v;
878  Omega = v2 * v;
879  //vh3 = Hreal * Omega;
880  //vh = cbrt(vh3);
881  eulerlogxabs = LAL_GAMMA + log (2.0 * (REAL8) m * v);
882 
883  /* Calculate the non-Keplerian velocity */
884  if (params->alignedSpins)
885  {
886  if (!use_optimized_v2)
887  {
888  vPhi =
889  XLALSimIMRSpinAlignedEOBNonKeplerCoeff (values->data, params);
890  }
891  else
892  {
893  /* OPTIMIZED */
894  vPhi = *vPhiInput;
895  /* END OPTIMIZED */
896  }
897 
898  if (XLAL_IS_REAL8_FAIL_NAN (vPhi))
899  {
900  XLAL_ERROR (XLAL_EFUNC);
901  }
902 
903  vPhi = r * cbrt (vPhi);
904  vPhi *= Omega;
905  vPhi2 = vPhi * vPhi;
906  }
907  else
908  {
909  vPhi = v;
910  vPhi2 = v2;
911  }
912 
913  /* Calculate the newtonian multipole, 1st term in Eq. 17, given by Eq. A1 */
914  status = XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole (&hNewton,
915  vPhi2, r,
916  values->data[1],
917  (UINT4) l, m,
918  params->
919  eobParams);
920  if (status == XLAL_FAILURE)
921  {
922  XLAL_ERROR (XLAL_EFUNC);
923  }
924 
925  /* Calculate the source term, 2nd term in Eq. 17, given by Eq. A5 */
926  if (((l + m) % 2) == 0)
927  {
928  Slm = (Hreal * Hreal - 1.) / (2. * eta) + 1.;
929  }
930  else
931  {
932  Slm = v * pp;
933  }
934  //printf( "Hreal = %e, Slm = %e, eta = %e\n", Hreal, Slm, eta );
935 
936  /* Calculate the absolute value of the Tail term,
937  * 3rd term in Eq. 17, given by Eq. A6, and Eq. (42) of
938  * http://arxiv.org/pdf/1212.4357.pdf */
939  k = m * Omega;
940  hathatk = Hreal * k;
941  hathatksq4 = 4. * hathatk * hathatk;
942  hathatk4pi = 4. * LAL_PI * hathatk;
943  /*
944  gsl_sf_result lnr1, arg1;
945  XLAL_CALLGSL( status = gsl_sf_lngamma_complex_e( l+1.0, -2.0*hathatk, &lnr1, &arg1 ) );
946  if (status != GSL_SUCCESS)
947  {
948  XLALPrintError("XLAL Error - %s: Error in GSL function\n", __func__ );
949  XLAL_ERROR( XLAL_EFUNC );
950  }
951  */
952  XLAL_CALLGSL (status = gsl_sf_fact_e (l, &z2));
953  if (status != GSL_SUCCESS)
954  {
955  XLALPrintError ("XLAL Error - %s: Error in GSL function: %s\n",
956  __func__, gsl_strerror (status));
957  XLAL_ERROR (XLAL_EFUNC);
958  }
959  /*
960  COMPLEX16 Tlmold;
961  Tlmold = cexp( ( lnr1.val + LAL_PI * hathatk ) + I * (
962  arg1.val + 2.0 * hathatk * log(4.0*k/sqrt(LAL_E)) ) );
963  Tlmold /= z2.val;
964  */
965  /* Calculating the prefactor of Tlm, outside the multiple product */
966  Tlmprefac = sqrt (hathatk4pi / (1. - exp (-hathatk4pi))) / z2.val;
967 
968  /* Calculating the multiple product factor */
969  for (Tlmprodfac = 1., i = 1; i <= l; i++)
970  {
971  Tlmprodfac *= (hathatksq4 + (REAL8) i * i);
972  }
973 
974  Tlm = Tlmprefac * sqrt (Tlmprodfac);
975 
976  /* Calculate the residue phase and amplitude terms */
977  /* deltalm is the 4th term in Eq. 17, delta 22 given by Eq. A15, others */
978  /* rholm is the 5th term in Eq. 17, given by Eqs. A8 - A14 */
979  /* auxflm is a special part of the 5th term in Eq. 17, given by Eq. A15 */
980  /* Actual values of the coefficients are defined in the next function of this file */
981  switch (l)
982  {
983  case 2:
984  switch (abs (m))
985  {
986  case 2:
987  rholm = 1. + v2 * (hCoeffs->rho22v2 + v * (hCoeffs->rho22v3
988  + v * (hCoeffs->rho22v4
989  +
990  v *
991  (hCoeffs->
992  rho22v5 +
993  v *
994  (hCoeffs->
995  rho22v6 +
996  hCoeffs->
997  rho22v6l *
998  eulerlogxabs +
999  v *
1000  (hCoeffs->
1001  rho22v7 +
1002  v *
1003  (hCoeffs->
1004  rho22v8 +
1005  hCoeffs->
1006  rho22v8l *
1007  eulerlogxabs +
1008  (hCoeffs->
1009  rho22v10 +
1010  hCoeffs->
1011  rho22v10l *
1012  eulerlogxabs)
1013  * v2)))))));
1014  break;
1015  case 1:
1016  {
1017  rholm = 1. + v * (hCoeffs->rho21v1
1018  + v * (hCoeffs->rho21v2 +
1019  v * (hCoeffs->rho21v3 +
1020  v * (hCoeffs->rho21v4 +
1021  v * (hCoeffs->rho21v5 +
1022  v * (hCoeffs->rho21v6 +
1023  hCoeffs->rho21v6l *
1024  eulerlogxabs +
1025  v *
1026  (hCoeffs->rho21v7 +
1027  hCoeffs->rho21v7l *
1028  eulerlogxabs +
1029  v *
1030  (hCoeffs->rho21v8 +
1031  hCoeffs->rho21v8l *
1032  eulerlogxabs +
1033  (hCoeffs->
1034  rho21v10 +
1035  hCoeffs->
1036  rho21v10l *
1037  eulerlogxabs) *
1038  v2))))))));
1039  auxflm = v * hCoeffs->f21v1 + v2 * v * hCoeffs->f21v3;
1040  }
1041  break;
1042  default:
1043  XLAL_ERROR (XLAL_EINVAL);
1044  break;
1045  }
1046  break;
1047  case 3:
1048  switch (m)
1049  {
1050  case 3:
1051  rholm =
1052  1. + v2 * (hCoeffs->rho33v2 +
1053  v * (hCoeffs->rho33v3 +
1054  v * (hCoeffs->rho33v4 +
1055  v * (hCoeffs->rho33v5 +
1056  v * (hCoeffs->rho33v6 +
1057  hCoeffs->rho33v6l * eulerlogxabs +
1058  v * (hCoeffs->rho33v7 +
1059  (hCoeffs->rho33v8 +
1060  hCoeffs->rho33v8l *
1061  eulerlogxabs) * v))))));
1062  auxflm = v * v2 * hCoeffs->f33v3;
1063  break;
1064  case 2:
1065  rholm = 1. + v * (hCoeffs->rho32v
1066  + v * (hCoeffs->rho32v2 +
1067  v * (hCoeffs->rho32v3 +
1068  v * (hCoeffs->rho32v4 +
1069  v * (hCoeffs->rho32v5 +
1070  v * (hCoeffs->rho32v6 +
1071  hCoeffs->rho32v6l *
1072  eulerlogxabs +
1073  (hCoeffs->rho32v8 +
1074  hCoeffs->rho32v8l *
1075  eulerlogxabs) *
1076  v2))))));
1077  break;
1078  case 1:
1079  rholm =
1080  1. + v2 * (hCoeffs->rho31v2 +
1081  v * (hCoeffs->rho31v3 +
1082  v * (hCoeffs->rho31v4 +
1083  v * (hCoeffs->rho31v5 +
1084  v * (hCoeffs->rho31v6 +
1085  hCoeffs->rho31v6l * eulerlogxabs +
1086  v * (hCoeffs->rho31v7 +
1087  (hCoeffs->rho31v8 +
1088  hCoeffs->rho31v8l *
1089  eulerlogxabs) * v))))));
1090  auxflm = v * v2 * hCoeffs->f31v3;
1091  break;
1092  default:
1093  XLAL_ERROR (XLAL_EINVAL);
1094  break;
1095  }
1096  break;
1097  case 4:
1098  switch (m)
1099  {
1100  case 4:
1101 
1102  rholm = 1. + v2 * (hCoeffs->rho44v2
1103  + v * (hCoeffs->rho44v3 + v * (hCoeffs->rho44v4
1104  +
1105  v *
1106  (hCoeffs->
1107  rho44v5 +
1108  (hCoeffs->
1109  rho44v6 +
1110  hCoeffs->
1111  rho44v6l *
1112  eulerlogxabs) *
1113  v))));
1114  break;
1115  case 3:
1116  rholm = 1. + v * (hCoeffs->rho43v
1117  + v * (hCoeffs->rho43v2
1118  + v2 * (hCoeffs->rho43v4 +
1119  v * (hCoeffs->rho43v5 +
1120  (hCoeffs->rho43v6 +
1121  hCoeffs->rho43v6l *
1122  eulerlogxabs) * v))));
1123  auxflm = v * hCoeffs->f43v;
1124  break;
1125  case 2:
1126  rholm = 1. + v2 * (hCoeffs->rho42v2
1127  + v * (hCoeffs->rho42v3 +
1128  v * (hCoeffs->rho42v4 +
1129  v * (hCoeffs->rho42v5 +
1130  (hCoeffs->rho42v6 +
1131  hCoeffs->rho42v6l *
1132  eulerlogxabs) * v))));
1133  break;
1134  case 1:
1135  rholm = 1. + v * (hCoeffs->rho41v
1136  + v * (hCoeffs->rho41v2
1137  + v2 * (hCoeffs->rho41v4 +
1138  v * (hCoeffs->rho41v5 +
1139  (hCoeffs->rho41v6 +
1140  hCoeffs->rho41v6l *
1141  eulerlogxabs) * v))));
1142  auxflm = v * hCoeffs->f41v;
1143  break;
1144  default:
1145  XLAL_ERROR (XLAL_EINVAL);
1146  break;
1147  }
1148  break;
1149  case 5:
1150  switch (m)
1151  {
1152  case 5:
1153  rholm = 1. + v2 * (hCoeffs->rho55v2
1154  + v * (hCoeffs->rho55v3 + v * (hCoeffs->rho55v4
1155  +
1156  v *
1157  (hCoeffs->
1158  rho55v5 +
1159  hCoeffs->
1160  rho55v6 * v))));
1161  break;
1162  case 4:
1163  rholm = 1. + v2 * (hCoeffs->rho54v2 + v * (hCoeffs->rho54v3
1164  + hCoeffs->rho54v4 * v));
1165  break;
1166  case 3:
1167  rholm = 1. + v2 * (hCoeffs->rho53v2
1168  + v * (hCoeffs->rho53v3 +
1169  v * (hCoeffs->rho53v4 +
1170  hCoeffs->rho53v5 * v)));
1171  break;
1172  case 2:
1173  rholm = 1. + v2 * (hCoeffs->rho52v2 + v * (hCoeffs->rho52v3
1174  + hCoeffs->rho52v4 * v));
1175  break;
1176  case 1:
1177  rholm = 1. + v2 * (hCoeffs->rho51v2
1178  + v * (hCoeffs->rho51v3 +
1179  v * (hCoeffs->rho51v4 +
1180  hCoeffs->rho51v5 * v)));
1181  break;
1182  default:
1183  XLAL_ERROR (XLAL_EINVAL);
1184  break;
1185  }
1186  break;
1187  case 6:
1188  switch (m)
1189  {
1190  case 6:
1191  rholm = 1. + v2 * (hCoeffs->rho66v2 + v * (hCoeffs->rho66v3
1192  + hCoeffs->rho66v4 * v));
1193  break;
1194  case 5:
1195  rholm = 1. + v2 * (hCoeffs->rho65v2 + hCoeffs->rho65v3 * v);
1196  break;
1197  case 4:
1198  rholm = 1. + v2 * (hCoeffs->rho64v2 + v * (hCoeffs->rho64v3
1199  + hCoeffs->rho64v4 * v));
1200  break;
1201  case 3:
1202  rholm = 1. + v2 * (hCoeffs->rho63v2 + hCoeffs->rho63v3 * v);
1203  break;
1204  case 2:
1205  rholm = 1. + v2 * (hCoeffs->rho62v2 + v * (hCoeffs->rho62v3
1206  + hCoeffs->rho62v4 * v));
1207  break;
1208  case 1:
1209  rholm = 1. + v2 * (hCoeffs->rho61v2 + hCoeffs->rho61v3 * v);
1210  break;
1211  default:
1212  XLAL_ERROR (XLAL_EINVAL);
1213  break;
1214  }
1215  break;
1216  case 7:
1217  switch (m)
1218  {
1219  case 7:
1220  rholm = 1. + v2 * (hCoeffs->rho77v2 + hCoeffs->rho77v3 * v);
1221  break;
1222  case 6:
1223  rholm = 1. + hCoeffs->rho76v2 * v2;
1224  break;
1225  case 5:
1226  rholm = 1. + v2 * (hCoeffs->rho75v2 + hCoeffs->rho75v3 * v);
1227  break;
1228  case 4:
1229  rholm = 1. + hCoeffs->rho74v2 * v2;
1230  break;
1231  case 3:
1232  rholm = 1. + v2 * (hCoeffs->rho73v2 + hCoeffs->rho73v3 * v);
1233  break;
1234  case 2:
1235  rholm = 1. + hCoeffs->rho72v2 * v2;
1236  break;
1237  case 1:
1238  rholm = 1. + v2 * (hCoeffs->rho71v2 + hCoeffs->rho71v3 * v);
1239  break;
1240  default:
1241  XLAL_ERROR (XLAL_EINVAL);
1242  break;
1243  }
1244  break;
1245  case 8:
1246  switch (m)
1247  {
1248  case 8:
1249  rholm = 1. + hCoeffs->rho88v2 * v2;
1250  break;
1251  case 7:
1252  rholm = 1. + hCoeffs->rho87v2 * v2;
1253  break;
1254  case 6:
1255  rholm = 1. + hCoeffs->rho86v2 * v2;
1256  break;
1257  case 5:
1258  rholm = 1. + hCoeffs->rho85v2 * v2;
1259  break;
1260  case 4:
1261  rholm = 1. + hCoeffs->rho84v2 * v2;
1262  break;
1263  case 3:
1264  rholm = 1. + hCoeffs->rho83v2 * v2;
1265  break;
1266  case 2:
1267  rholm = 1. + hCoeffs->rho82v2 * v2;
1268  break;
1269  case 1:
1270  rholm = 1. + hCoeffs->rho81v2 * v2;
1271  break;
1272  default:
1273  XLAL_ERROR (XLAL_EINVAL);
1274  break;
1275  }
1276  break;
1277  default:
1278  XLAL_ERROR (XLAL_EINVAL);
1279  break;
1280  }
1281 
1282  /* Raise rholm to the lth power */
1283  rholmPwrl = 1.0;
1284  i = l;
1285  while (i--)
1286  {
1287  rholmPwrl *= rholm;
1288  }
1289  /* In the equal-mass odd m case, there is no contribution from nonspin terms,
1290  * and the only contribution comes from the auxflm term that is proportional to chiA (asymmetric spins).
1291  * In this case, we must ignore the nonspin terms directly, since the leading term defined by
1292  * CalculateThisMultipolePrefix in LALSimIMREOBNewtonianMultipole.c is not zero (see comments there).
1293  */
1294  if (eta == 0.25 && m % 2)
1295  {
1296  rholmPwrl = auxflm;
1297  }
1298  else
1299  {
1300  rholmPwrl += auxflm;
1301  }
1302 
1303  /*if (r > 8.5)
1304  {
1305  printf("YP::dynamics variables in waveform: %i, %i, %e, %e\n",l,m,r,pp);
1306  printf( "rholm^l = %.16e, Tlm = %.16e + i %.16e, \nSlm = %.16e, hNewton = %.16e + i %.16e\n", rholmPwrl, Tlm, 0., Slm, creal(hNewton), cimag(hNewton) );} */
1307  /* Put all factors in Eq. 17 together */
1308  *hlm = Tlm * Slm * rholmPwrl;
1309  *hlm *= hNewton;
1310  /*if (r > 8.5)
1311  {
1312  printf("YP::FullWave: %.16e,%.16e, %.16e\n",hlm->re,hlm->im,sqrt(hlm->re*hlm->re+hlm->im*hlm->im));
1313  } */
1314  return XLAL_SUCCESS;
1315 }
1316 
1317 
1318 /*--------------------------------------------------------------*/
1319 /**
1320  * Spin Factors
1321  */
1322 
1323 /**
1324  * This function calculates coefficients for hlm mode factorized-resummed waveform.
1325  * The coefficients are pre-computed and stored in the SpinEOBParams structure.
1326  * Appendix of the paper, and papers DIN (PRD 79, 064004 (2009) [arXiv:0811.2069]) and PBFRT (PRD 83, 064003 (2011)).
1327  * Concerning SEOBNRv4 see also https://dcc.ligo.org/T1600383
1328  */
1329 
1330 UNUSED static int XLALSimIMREOBCalcSpinFacWaveformCoefficients (FacWaveformCoeffs * const coeffs,
1331  /**< OUTPUT, pre-computed waveform coefficients */
1332  SpinEOBParams * restrict params,
1333  /**< Others parameters like use_hm */
1334  const REAL8 m1,
1335  /**< mass 1 */
1336  const REAL8 m2,
1337  /**< mass 2 */
1338  const REAL8 eta,
1339  /**< symmetric mass ratio */
1340  REAL8 a,
1341  /**< Kerr spin parameter for test-particle terms */
1342  const REAL8 chiS,
1343  /**< (chi1+chi2)/2 */
1344  const REAL8 chiA,
1345  /**< (chi1-chi2)/2 */
1346  const UINT4 SpinAlignedEOBversion
1347  /**< 1 for SEOBNRv1; 2 for SEOBNRv2; 4 for SEOBNRv4 */
1348  )
1349 {
1350  if ( SpinAlignedEOBversion != 1 && SpinAlignedEOBversion != 2 && SpinAlignedEOBversion != 4)
1351  {
1352  XLALPrintError
1353  ("XLAL Error - %s: SpinAlignedEOBversion = %u is wrong, must be 1 or 2 or 4!\n",
1354  __func__, SpinAlignedEOBversion);
1355  XLAL_ERROR (XLAL_EINVAL);
1356  }
1357 
1358  INT4 use_hm = params->use_hm;
1359 
1360 
1361  REAL8 eta2 = eta * eta;
1362  REAL8 eta3 = eta2 * eta;
1363 
1364  REAL8 dM, dM2, chiA2, chiS2; //dM3;
1365  REAL8 aDelta, a2, a3;
1366 
1367  /* Combination which appears a lot */
1368  REAL8 m1Plus3eta, m1Plus3eta2, m1Plus3eta3;
1369 
1370  dM2 = 1. - 4. * eta;
1371  chiA2 = chiA * chiA;
1372  chiS2 = chiS * chiS;
1373  REAL8 chiA3 = chiA2 * chiA;
1374  REAL8 chiS3 = chiS2 * chiS;
1375 
1376  /* Combinations of kappa_1,2, coefficients of the spin-induced quadrupole */
1377  REAL8 kappaS, kappaA;
1378 
1379  //printf( "****************************** a = %e *********************************\n", a );
1380 
1381  /* Check that deltaM has a reasonable value */
1382  if (dM2 < 0. && dM2 > -1e-4 ) {
1383  dM2 = 0.;
1384  }
1385  if (dM2 < 0)
1386  {
1387  XLALPrintError
1388  ("XLAL Error - %s: dM2 = %.16f < -1e-4 this isn't allowed!\n",
1389  __func__,dM2);
1390  XLAL_ERROR (XLAL_EINVAL);
1391  }
1392 
1393  dM = sqrt (dM2);
1394  if (m1 < m2)
1395  {
1396  dM = -dM;
1397  }
1398  //dM3 = dM2 * dM;
1399 
1400  aDelta = 0.; // a value in delta_lm is 0 in SEOBNRv1, SEOBNRv2, and SEOBNRv4
1401  a2 = a * a;
1402  a3 = a2 * a;
1403 
1404  m1Plus3eta = -1. + 3. * eta;
1405  m1Plus3eta2 = m1Plus3eta * m1Plus3eta;
1406  m1Plus3eta3 = m1Plus3eta * m1Plus3eta2;
1407 
1408  /* Initialize all coefficients to zero */
1409  /* This is important, as we will not set some if dM is zero */
1410  memset (coeffs, 0, sizeof (FacWaveformCoeffs));
1411 
1412 
1413  /* l = 2, Eqs. A8a and A8b for rho, Eq. A15a for f,
1414  Eqs. 20 and 21 of DIN and Eqs. 27a and 27b of PBFRT for delta */
1415 
1416  coeffs->delta22vh3 = 7. / 3.;
1417  coeffs->delta22vh6 = (-4. * aDelta) / 3. + (428. * LAL_PI) / 105.;
1418  /* See https://dcc.ligo.org/T1600383 */
1419  if (SpinAlignedEOBversion == 4)
1420  {
1421  coeffs->delta22vh6 =
1422  -4. / 3. * (dM * chiA + chiS * (1 - 2 * eta)) +
1423  (428. * LAL_PI) / 105.;
1424  }
1425  coeffs->delta22v8 = (20. * aDelta) / 63.;
1426  coeffs->delta22vh9 = -2203. / 81. + (1712. * LAL_PI * LAL_PI) / 315.;
1427  coeffs->delta22v5 = -24. * eta;
1428  coeffs->delta22v6 = 0.0;
1429  if (SpinAlignedEOBversion == 2 && chiS + chiA * dM / (1. - 2. * eta) > 0.)
1430  {
1431  coeffs->delta22v6 = -540. * eta * (chiS + chiA * dM / (1. - 2. * eta));
1432  }
1433 
1434  coeffs->rho22v2 = -43. / 42. + (55. * eta) / 84.;
1435  coeffs->rho22v3 = (-2. * (chiS + chiA * dM - chiS * eta)) / 3.;
1436  switch (SpinAlignedEOBversion)
1437  {
1438  case 1:
1439  coeffs->rho22v4 = -20555. / 10584. + 0.5 * a2
1440  - (33025. * eta) / 21168. + (19583. * eta2) / 42336.;
1441  break;
1442  case 2:
1443  case 4:
1444  coeffs->rho22v4 =
1445  -20555. / 10584. + 0.5 * (chiS + chiA * dM) * (chiS + chiA * dM) -
1446  (33025. * eta) / 21168. + (19583. * eta2) / 42336.;
1447  /* If spin-induced quadrupole coefficients kappa_1,2 are not 1, include the leading-order correction at 2PN in rho22 */
1448  /* Relevant for BNS - TEOBv4 */
1449  if((params->seobCoeffs->tidal1->quadparam - 1.) != 0. || (params->seobCoeffs->tidal2->quadparam - 1.) != 0.) {
1450  kappaS = 0.5 * (params->seobCoeffs->tidal1->quadparam + params->seobCoeffs->tidal2->quadparam);
1451  kappaA = 0.5 * (params->seobCoeffs->tidal1->quadparam - params->seobCoeffs->tidal2->quadparam);
1452  coeffs->rho22v4 += chiA*chiA*(0.5*dM*kappaA - kappaS*eta + 0.5*kappaS + eta - 0.5) + chiA*chiS*(dM*kappaS - dM - 2*kappaA*eta + kappaA) + chiS*chiS*(0.5*dM*kappaA - kappaS*eta + 0.5*kappaS + eta - 0.5);
1453  }
1454  break;
1455  default:
1456  XLALPrintError
1457  ("XLAL Error - %s: SpinAlignedEOBversion = %u is wrong, must be 1 or 2 or 4!\n",
1458  __func__,SpinAlignedEOBversion);
1459  XLAL_ERROR (XLAL_EINVAL);
1460  break;
1461  }
1462  coeffs->rho22v5 = (-34. * a) / 21.;
1463  /* See https://dcc.ligo.org/T1600383 */
1464  if (SpinAlignedEOBversion == 4)
1465  {
1466  coeffs->rho22v5 =
1467  (-34. / 21. + 49. * eta / 18. + 209. * eta2 / 126.) * chiS +
1468  (-34. / 21. - 19. * eta / 42.) * dM * chiA;
1469  }
1470  coeffs->rho22v6 =
1471  1556919113. / 122245200. + (89. * a2) / 252. -
1472  (48993925. * eta) / 9779616. - (6292061. * eta2) / 3259872. +
1473  (10620745. * eta3) / 39118464. + (41. * eta * LAL_PI * LAL_PI) / 192.;
1474  coeffs->rho22v6l = -428. / 105.;
1475  coeffs->rho22v7 = (18733. * a) / 15876. + a * a2 / 3.;
1476  /* See https://dcc.ligo.org/T1600383 */
1477  if (SpinAlignedEOBversion == 4)
1478  {
1479  coeffs->rho22v7 =
1480  a3 / 3. + chiA * dM * (18733. / 15876. + (50140. * eta) / 3969. +
1481  (97865. * eta2) / 63504.) +
1482  chiS * (18733. / 15876. + (74749. * eta) / 5292. -
1483  (245717. * eta2) / 63504. + (50803. * eta3) / 63504.);
1484  }
1485  switch (SpinAlignedEOBversion)
1486  {
1487  case 1:
1488  coeffs->rho22v8 =
1489  eta * (-5.6 - 117.6 * eta) - 387216563023. / 160190110080. +
1490  (18353. * a2) / 21168. - a2 * a2 / 8.;
1491  break;
1492  case 2:
1493  case 4:
1494  coeffs->rho22v8 =
1495  -387216563023. / 160190110080. + (18353. * a2) / 21168. -
1496  a2 * a2 / 8.;
1497  break;
1498  default:
1499  XLALPrintError
1500  ("XLAL Error - %s: SpinAlignedEOBversion = %u is wrong, must be 1 or 2 or 4!\n",
1501  __func__,SpinAlignedEOBversion);
1502  XLAL_ERROR (XLAL_EINVAL);
1503  break;
1504  }
1505  coeffs->rho22v8l = 9202. / 2205.;
1506  coeffs->rho22v10 = -16094530514677. / 533967033600.;
1507  coeffs->rho22v10l = 439877. / 55566.;
1508 
1509  /*printf( "v2 = %.16e, v3 = %.16e, v4 = %.16e, v5 = %.16e\n"
1510  "v6 = %.16e, v6l = %.16e v7 = %.16e v8 = %.16e, v8l = %.16e v10 = %.16e v10l = %.16e\n",
1511  coeffs->rho22v2, coeffs->rho22v3, coeffs->rho22v4, coeffs->rho22v5, coeffs->rho22v6,
1512  coeffs->rho22v6l, coeffs->rho22v7, coeffs->rho22v8, coeffs->rho22v8l, coeffs->rho22v10,
1513  coeffs->rho22v10l ); */
1514 
1515  /*RC the HM model does not include spinning test mass terms in the higher-order modes */
1516  if(use_hm){
1517  a = 0;
1518  a2 = 0;
1519  a3 = 0;
1520  }
1521 
1522  //RC: The delta coefficient before were put inside the if(dM2). This is wrong.
1523  //it didn't effect the models before because this function is used only to
1524  //calculate the 22 mode and the flux of the other modes, but for the flux
1525  //you only need |h_lm|. This error was present in all the modes, and now is fixed.
1526  //Is not a problem for the 22 mode because there is no if(dM2) for the 22
1527  coeffs->delta21vh3 = 2. / 3.;
1528  coeffs->delta21vh6 = (-17. * aDelta) / 35. + (107. * LAL_PI) / 105.;
1529  coeffs->delta21vh7 = (3. * aDelta * aDelta) / 140.;
1530  coeffs->delta21vh9 = -272. / 81. + (214. * LAL_PI * LAL_PI) / 315.;
1531  coeffs->delta21v5 = -493. * eta / 42.;
1532  if (dM2)
1533  {
1534 
1535  //coeffs->rho21v1 = (-3.*(chiS+chiA/dM))/(4.);
1536  coeffs->rho21v1 = 0.0;
1537  //coeffs->rho21v2 = -59./56 - (9.*chiAPlusChiSdM*chiAPlusChiSdM)/(32.*dM2) + (23.*eta)/84.;
1538  switch (SpinAlignedEOBversion)
1539  {
1540  case 1:
1541  coeffs->rho21v2 = -59. / 56 + (23. * eta) / 84. - 9. / 32. * a2;
1542  coeffs->rho21v3 = 1177. / 672. * a - 27. / 128. * a3;
1543  break;
1544  case 2:
1545  case 4:
1546  coeffs->rho21v2 = -59. / 56 + (23. * eta) / 84.;
1547  coeffs->rho21v3 = 0.0;
1548  break;
1549  default:
1550  XLALPrintError
1551  ("XLAL Error - %s: SpinAlignedEOBversion = %u is wrong, must be 1 or 2 or 4!\n",
1552  __func__,SpinAlignedEOBversion);
1553  XLAL_ERROR (XLAL_EINVAL);
1554  break;
1555  }
1556  /*coeffs->rho21v3 = (-567.*chiA*chiA*chiA - 1701.*chiA*chiA*chiS*dM
1557  + chiA*(-4708. + 1701.*chiS*chiS - 2648.*eta)*(-1. + 4.*eta)
1558  + chiS* dM3 *(4708. - 567.*chiS*chiS
1559  + 1816.*eta))/(2688.*dM3); */
1560  coeffs->rho21v4 =
1561  -47009. / 56448. - (865. * a2) / 1792. - (405. * a2 * a2) / 2048. -
1562  (10993. * eta) / 14112. + (617. * eta2) / 4704.;
1563  coeffs->rho21v5 =
1564  (-98635. * a) / 75264. + (2031. * a * a2) / 7168. -
1565  (1701. * a2 * a3) / 8192.;
1566  coeffs->rho21v6 =
1567  7613184941. / 2607897600. + (9032393. * a2) / 1806336. +
1568  (3897. * a2 * a2) / 16384. - (15309. * a3 * a3) / 65536.;
1569  coeffs->rho21v6l = -107. / 105.;
1570  coeffs->rho21v7 =
1571  (-3859374457. * a) / 1159065600. - (55169. * a3) / 16384. +
1572  (18603. * a2 * a3) / 65536. - (72171. * a2 * a2 * a3) / 262144.;
1573  coeffs->rho21v7l = 107. * a / 140.;
1574  coeffs->rho21v8 = -1168617463883. / 911303737344.;
1575  coeffs->rho21v8l = 6313. / 5880.;
1576  coeffs->rho21v10 = -63735873771463. / 16569158860800.;
1577  coeffs->rho21v10l = 5029963. / 5927040.;
1578 
1579  coeffs->f21v1 = (-3. * (chiS + chiA / dM)) / (2.);
1580  switch (SpinAlignedEOBversion)
1581  {
1582  case 1:
1583  coeffs->f21v3 = 0.0;
1584  break;
1585  case 2:
1586  case 4:
1587  coeffs->f21v3 =
1588  (chiS * dM * (427. + 79. * eta) +
1589  chiA * (147. + 280. * dM * dM + 1251. * eta)) / 84. / dM;
1590  /*RC: New terms for SEOBNRv4HM, they are put to zero if use_hm == 0 */
1591  coeffs->f21v4 = 0.0;
1592  coeffs->f21v5 = 0.0;
1593  coeffs->f21v6 = 0.0;
1594  coeffs->f21v7c = 0;
1595  if(use_hm){
1596  //RC: This terms are in Eq.A11 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
1597  coeffs->f21v4 = (-3.-2.*eta)*chiA2 + (-3.+ eta/2.)*chiS2 + (-6.+21.*eta/2.)*chiS*chiA/dM;
1598  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
1599  + (-81./16. - 703.*eta2/112. + 8797.*eta/1008.+(9./4. - 6.*eta)*chiS2)*chiA/dM;
1600  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;
1601  coeffs->f21v7c = 0; //RC: this is the calibration parameter which is initially set to 0
1602  }
1603  /* End new terms for SEOBNRv4HM */
1604  break;
1605  default:
1606  XLALPrintError
1607  ("XLAL Error - %s: SpinAlignedEOBversion = %u is wrong, must be 1 or 2 or 4!\n",
1608  __func__,SpinAlignedEOBversion);
1609  XLAL_ERROR (XLAL_EINVAL);
1610  break;
1611  }
1612  }
1613  else
1614  {
1615  coeffs->f21v1 = -3. * chiA / 2.;
1616  switch (SpinAlignedEOBversion)
1617  {
1618  case 1:
1619  coeffs->f21v3 = 0.0;
1620  break;
1621  case 2:
1622  case 4:
1623  coeffs->f21v3 =
1624  (chiS * dM * (427. + 79. * eta) +
1625  chiA * (147. + 280. * dM * dM + 1251. * eta)) / 84.;
1626  /* New terms for SEOBNRv4HM, they are put to zero if use_hm == 0 */
1627  coeffs->f21v4 = 0.0;
1628  coeffs->f21v5 = 0.0;
1629  coeffs->f21v6 = 0.0;
1630  if(use_hm){
1631  //RC: This terms are in Eq.A11 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
1632  coeffs->f21v4 = (-6+21*eta/2.)*chiS*chiA;
1633  coeffs->f21v5 = (3./4.-3.*eta)*chiA3 + (-81./16. - 703.*eta2/112. + 8797.*eta/1008. + (9./4. - 6.*eta)*chiS2)*chiA;
1634  coeffs->f21v6 = (4163./126.-1636.*eta/21. + 1088.*eta2/63.)*chiS*chiA;
1635  }
1636  /* End new terms for SEOBNRv4HM */
1637  break;
1638  default:
1639  XLALPrintError
1640  ("XLAL Error - %s: SpinAlignedEOBversion = %u is wrong, must be 1 or 2 or 4!\n",
1641  __func__,SpinAlignedEOBversion);
1642  XLAL_ERROR (XLAL_EINVAL);
1643  break;
1644  }
1645  }
1646 
1647 
1648  /* l = 3, Eqs. A9a - A9c for rho, Eqs. A15b and A15c for f,
1649  Eqs. 22 - 24 of DIN and Eqs. 27c - 27e of PBFRT for delta */
1650  coeffs->delta33vh3 = 13. / 10.;
1651  coeffs->delta33vh6 = (-81. * aDelta) / 20. + (39. * LAL_PI) / 7.;
1652  coeffs->delta33vh9 = -227827. / 3000. + (78. * LAL_PI * LAL_PI) / 7.;
1653  coeffs->delta33v5 = -80897. * eta / 2430.;
1654  if (dM2)
1655  {
1656  coeffs->rho33v2 = -7. / 6. + (2. * eta) / 3.;
1657  //coeffs->rho33v3 = (chiS*dM*(-4. + 5.*eta) + chiA*(-4. + 19.*eta))/(6.*dM);
1658  coeffs->rho33v3 = 0.0;
1659  coeffs->rho33v4 =
1660  -6719. / 3960. + a2 / 2. - (1861. * eta) / 990. +
1661  (149. * eta2) / 330.;
1662  coeffs->rho33v5 = (-4. * a) / 3.;
1663  coeffs->rho33v6 = 3203101567. / 227026800. + (5. * a2) / 36.;
1664  coeffs->rho33v6l = -26. / 7.;
1665  coeffs->rho33v7 = (5297. * a) / 2970. + a * a2 / 3.;
1666  coeffs->rho33v8 = -57566572157. / 8562153600.;
1667  coeffs->rho33v8l = 13. / 3.;
1668  coeffs->rho33v10 = 0;
1669  coeffs->rho33v10l = 0;
1670  if(use_hm){
1671  //RC: This terms are in Eq.A6 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
1672  coeffs->rho33v6 = 3203101567. / 227026800. + (5. * a2) / 36. + (-129509./25740. + 41./192. * LAL_PI*LAL_PI)*eta - 274621./154440.*eta2+ 12011./46332.*eta3;
1673  coeffs->rho33v10 = -903823148417327./30566888352000.;
1674  coeffs->rho33v10l = 87347./13860.;
1675  }
1676 
1677 
1678 
1679 
1680  coeffs->f33v3 =
1681  (chiS * dM * (-4. + 5. * eta) + chiA * (-4. + 19. * eta)) / (2. * dM);
1682  coeffs->f33v4 = 0;
1683  coeffs->f33v5 = 0;
1684  coeffs->f33v6 = 0;
1685  coeffs->f33vh6 = 0;
1686  if(use_hm){
1687  //RC: This terms are in Eq.A10 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
1688  coeffs->f33v4 = (3./2. * chiS2 * dM + (3. - 12 * eta) * chiA * chiS + dM * (3./2. -6. * eta) * chiA2)/(dM);
1689  coeffs->f33v5 = (dM * (241./30. * eta2 + 11./20. * eta + 2./3.) * chiS + (407./30. * eta2 - 593./60. * eta + 2./3.)* chiA)/(dM);
1690  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;
1691  coeffs->f33vh6 = (dM * (593. / 108. * eta - 81./20.) * chiS + (7339./540. * eta - 81./20.) * chiA)/(dM);
1692  }
1693  }
1694  else
1695  {
1696  coeffs->f33v3 = chiA * 3. / 8.;
1697  coeffs->f33v4 = 0;
1698  coeffs->f33v5 = 0;
1699  coeffs->f33v6 = 0;
1700  coeffs->f33vh6 = 0;
1701  if(use_hm){
1702  //RC: This terms are in Eq.A10 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
1703  coeffs->f33v4 = ((3. - 12 * eta) * chiA * chiS);
1704  coeffs->f33v5 = ((407./30. * eta2 - 593./60. * eta + 2./3.)* chiA);
1705  coeffs->f33v6 = ((44. * eta2 - 1. * eta - 7./2.) * chiA * chiS);
1706  coeffs->f33vh6 = ((7339./540. * eta - 81./20.) * chiA);
1707  }
1708  }
1709 
1710  coeffs->delta32vh3 = (10. + 33. * eta) / (-15. * m1Plus3eta);
1711  coeffs->delta32vh4 = 4. * aDelta;
1712  coeffs->delta32vh6 = (-136. * aDelta) / 45. + (52. * LAL_PI) / 21.;
1713  coeffs->delta32vh9 = -9112. / 405. + (208. * LAL_PI * LAL_PI) / 63.;
1714 
1715  coeffs->rho32v = (4. * chiS * eta) / (-3. * m1Plus3eta);
1716  coeffs->rho32v2 = (-4. * a2 * eta2) / (9. * m1Plus3eta2) + (328. - 1115. * eta +
1717  320. * eta2) / (270. *
1718  m1Plus3eta);
1719  if (SpinAlignedEOBversion == 4) {
1720  coeffs->rho32v2 = (328. - 1115. * eta +
1721  320. * eta2) / (270. *
1722  m1Plus3eta);
1723  }
1724  coeffs->rho32v3 =
1725  (2. *
1726  (45. * a * m1Plus3eta3 -
1727  a * eta * (328. - 2099. * eta + 5. * (733. + 20. * a2) * eta2 -
1728  960. * eta3))) / (405. * m1Plus3eta3);
1729  coeffs->rho32v3 = 2. / 9. * a;
1730  coeffs->rho32v4 = a2 / 3. + (-1444528.
1731  + 8050045. * eta - 4725605. * eta2 -
1732  20338960. * eta3 +
1733  3085640. * eta2 * eta2) / (1603800. *
1734  m1Plus3eta2);
1735  coeffs->rho32v5 = (-2788. * a) / 1215.;
1736  coeffs->rho32v6 = 5849948554. / 940355325. + (488. * a2) / 405.;
1737  coeffs->rho32v6l = -104. / 63.;
1738  coeffs->rho32v8 = -10607269449358. / 3072140846775.;
1739  coeffs->rho32v8l = 17056. / 8505.;
1740 
1741  coeffs->delta31vh3 = 13. / 30.;
1742  coeffs->delta31vh6 = (61. * aDelta) / 20. + (13. * LAL_PI) / 21.;
1743  coeffs->delta31vh7 = (-24. * aDelta * aDelta) / 5.;
1744  coeffs->delta31vh9 = -227827. / 81000. + (26. * LAL_PI * LAL_PI) / 63.;
1745  coeffs->delta31v5 = -17. * eta / 10.;
1746  if (dM2)
1747  {
1748 
1749  coeffs->rho31v2 = -13. / 18. - (2. * eta) / 9.;
1750  //coeffs->rho31v3 = (chiA*(-4. + 11.*eta) + chiS*dM*(-4. + 13.*eta))/(6.*dM);
1751  coeffs->rho31v3 = 0.0;
1752  coeffs->rho31v4 = 101. / 7128.
1753  - (5. * a2) / 6. - (1685. * eta) / 1782. - (829. * eta2) / 1782.;
1754  coeffs->rho31v5 = (4. * a) / 9.;
1755  coeffs->rho31v6 = 11706720301. / 6129723600. - (49. * a2) / 108.;
1756  coeffs->rho31v6l = -26. / 63.;
1757  coeffs->rho31v7 = (-2579. * a) / 5346. + a * a2 / 9.;
1758  coeffs->rho31v8 = 2606097992581. / 4854741091200.;
1759  coeffs->rho31v8l = 169. / 567.;
1760 
1761  coeffs->f31v3 =
1762  (chiA * (-4. + 11. * eta) +
1763  chiS * dM * (-4. + 13. * eta)) / (2. * dM);
1764  }
1765  else
1766  {
1767  coeffs->f31v3 = -chiA * 5. / 8.;
1768  }
1769 
1770  /* l = 4, Eqs. A10a - A10d for delta, Eq. A15d for f
1771  Eqs. 25 - 28 of DIN and Eqs. 27f - 27i of PBFRT for delta */
1772 
1773  coeffs->delta44vh3 = (112. + 219. * eta) / (-120. * m1Plus3eta);
1774  coeffs->delta44vh6 = (-464. * aDelta) / 75. + (25136. * LAL_PI) / 3465.;
1775  coeffs->delta44vh9 = 0.;
1776  if(use_hm){
1777  //RC: This terms are in Eq.A15 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
1778  coeffs->delta44vh9 = -55144./375. + 201088.*LAL_PI*LAL_PI/10395.;
1779  }
1780 
1781  coeffs->rho44v2 =
1782  (1614. - 5870. * eta + 2625. * eta2) / (1320. * m1Plus3eta);
1783  coeffs->rho44v3 =
1784  (chiA * (10. - 39. * eta) * dM +
1785  chiS * (10. - 41. * eta + 42. * eta2)) / (15. * m1Plus3eta);
1786  coeffs->rho44v4 =
1787  a2 / 2. + (-511573572. + 2338945704. * eta - 313857376. * eta2 -
1788  6733146000. * eta3 +
1789  1252563795. * eta2 * eta2) / (317116800. * m1Plus3eta2);
1790  coeffs->rho44v5 = (-69. * a) / 55.;
1791  coeffs->rho44v8 = 0.;
1792  coeffs->rho44v8l = 0.;
1793  coeffs->rho44v10 = 0.;
1794  coeffs->rho44v10l = 0;
1795  if(use_hm){
1796  //RC: This terms are in Eq.A8 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
1797  coeffs->rho44v4 =
1798  (-511573572. + 2338945704. * eta - 313857376. * eta2 -
1799  6733146000. * eta3 +
1800  1252563795. * eta2 * eta2) / (317116800. * m1Plus3eta2)
1801  + chiS2/2. + dM*chiS*chiA + dM2*chiA2/2.;
1802  coeffs->rho44v5 = chiA*dM*(-8280. + 42716.*eta - 57990.*eta2 + 8955*eta3)/(6600.*m1Plus3eta2)
1803  + chiS*(-8280. + 66284.*eta-176418.*eta2+128085.*eta3 + 88650*eta2*eta2)/(6600.*m1Plus3eta2);
1804  coeffs->rho44v8 = -172066910136202271./19426955708160000.;
1805  coeffs->rho44v8l = 845198./190575.;
1806  coeffs->rho44v10 = - 17154485653213713419357./568432724020761600000.;
1807  coeffs->rho44v10l = 22324502267./3815311500;
1808  }
1809  coeffs->rho44v6 = 16600939332793. / 1098809712000. + (217. * a2) / 3960.;
1810  coeffs->rho44v6l = -12568. / 3465.;
1811 
1812  coeffs->delta43vh3 = (486. + 4961. * eta) / (810. * (1. - 2. * eta));
1813  coeffs->delta43vh4 = (11. * aDelta) / 4.;
1814  coeffs->delta43vh6 = 1571. * LAL_PI / 385.;
1815  if (dM2)
1816  {
1817 
1818  //coeffs->rho43v = (5.*(chiA - chiS*dM)*eta)/(8.*dM*(-1. + 2.*eta));
1819  coeffs->rho43v = 0.0;
1820  coeffs->rho43v2 =
1821  (222. - 547. * eta + 160. * eta2) / (176. * (-1. + 2. * eta));
1822  coeffs->rho43v4 = -6894273. / 7047040. + (3. * a2) / 8.;
1823  coeffs->rho43v5 = (-12113. * a) / 6160.;
1824  coeffs->rho43v6 = 1664224207351. / 195343948800.;
1825  coeffs->rho43v6l = -1571. / 770.;
1826 
1827  coeffs->f43v =
1828  (5. * (chiA - chiS * dM) * eta) / (2. * dM * (-1. + 2. * eta));
1829  }
1830  else
1831  {
1832  coeffs->f43v = -5. * chiA / 4.;
1833  }
1834 
1835  coeffs->delta42vh3 = (7. * (1. + 6. * eta)) / (-15. * m1Plus3eta);
1836  coeffs->delta42vh6 = (212. * aDelta) / 75. + (6284. * LAL_PI) / 3465.;
1837 
1838  coeffs->rho42v2 =
1839  (1146. - 3530. * eta + 285. * eta2) / (1320. * m1Plus3eta);
1840  coeffs->rho42v3 =
1841  (chiA * (10. - 21. * eta) * dM +
1842  chiS * (10. - 59. * eta + 78. * eta2)) / (15. * m1Plus3eta);
1843  coeffs->rho42v4 =
1844  a2 / 2. + (-114859044. + 295834536. * eta + 1204388696. * eta2 -
1845  3047981160. * eta3 -
1846  379526805. * eta2 * eta2) / (317116800. * m1Plus3eta2);
1847  coeffs->rho42v5 = (-7. * a) / 110.;
1848  coeffs->rho42v6 = 848238724511. / 219761942400. + (2323. * a2) / 3960.;
1849  coeffs->rho42v6l = -3142. / 3465.;
1850 
1851  coeffs->delta41vh3 = (2. + 507. * eta) / (10. * (1. - 2. * eta));
1852  coeffs->delta41vh4 = (11. * aDelta) / 12.;
1853  coeffs->delta41vh6 = 1571. * LAL_PI / 3465.;
1854 
1855  if (dM2)
1856  {
1857 
1858  //coeffs->rho41v = (5.*(chiA - chiS*dM)*eta)/(8.*dM*(-1. + 2.*eta));
1859  coeffs->rho41v = 0.0;
1860  coeffs->rho41v2 =
1861  (602. - 1385. * eta + 288. * eta2) / (528. * (-1. + 2. * eta));
1862  coeffs->rho41v4 = -7775491. / 21141120. + (3. * a2) / 8.;
1863  coeffs->rho41v5 = (-20033. * a) / 55440. - (5 * a * a2) / 6.;
1864  coeffs->rho41v6 = 1227423222031. / 1758095539200.;
1865  coeffs->rho41v6l = -1571. / 6930.;
1866 
1867  coeffs->f41v =
1868  (5. * (chiA - chiS * dM) * eta) / (2. * dM * (-1. + 2. * eta));
1869  }
1870  else
1871  {
1872  coeffs->f41v = -5. * chiA / 4.;
1873  }
1874 
1875  /* l = 5, Eqs. A11a - A11e for rho,
1876  Eq. 29 of DIN and Eqs. E1a and E1b of PBFRT for delta */
1877  coeffs->delta55vh3 =
1878 (96875. + 857528. * eta) / (131250. * (1 - 2 * eta));
1879  coeffs->delta55vh6 = 0;
1880  coeffs->delta55vh9 = 0;
1881  if(use_hm){
1882  //RC: This terms are in Eq.A16 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
1883  coeffs->delta55vh6 = 3865./429.*LAL_PI;
1884  coeffs->delta55vh9 = (-7686949127. + 954500400.*LAL_PI*LAL_PI)/31783752.;
1885  }
1886  if (dM2)
1887  {
1888 
1889  coeffs->rho55v2 =
1890  (487. - 1298. * eta + 512. * eta2) / (390. * (-1. + 2. * eta));
1891  coeffs->rho55v3 = (-2. * a) / 3.;
1892  coeffs->rho55v4 = -3353747. / 2129400. + a2 / 2.;
1893  coeffs->rho55v5 = -241. * a / 195.;
1894  coeffs->rho55v6 = 0.;
1895  coeffs->rho55v6l = 0.;
1896  coeffs->rho55v8 = 0.;
1897  coeffs->rho55v8l = 0.;
1898  coeffs->rho55v10 = 0.;
1899  coeffs->rho55v10l = 0.;
1900  coeffs->f55v3 = 0.;
1901  coeffs->f55v4 = 0.;
1902  coeffs->f55v5c = 0;
1903  if(use_hm){
1904  //RC: This terms are in Eq.A9 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
1905  coeffs->rho55v6 = 190606537999247./11957879934000.;
1906  coeffs->rho55v6l = - 1546./429.;
1907  coeffs->rho55v8 = - 1213641959949291437./118143853747920000.;
1908  coeffs->rho55v8l = 376451./83655.;
1909  coeffs->rho55v10 = -150082616449726042201261./4837990810977324000000.;
1910  coeffs->rho55v10l = 2592446431./456756300.;
1911 
1912  coeffs->f55v3 = chiA/dM *(10./(3.*(-1.+2.*eta)) - 70.*eta/(3.*(-1.+2.*eta)) + 110.*eta2/(3.*(-1.+2.*eta)) ) +
1913  chiS*(10./(3.*(-1.+2.*eta)) -10.*eta/(-1.+2.*eta) + 10*eta2/(-1.+2.*eta));
1914  coeffs->f55v4 = chiS2*(-5./(2.*(-1.+2.*eta)) + 5.*eta/(-1.+2.*eta)) +
1915  chiA*chiS/dM *(-5./(-1.+2.*eta) + 30.*eta/(-1.+2.*eta) - 40.*eta2/(-1.+2.*eta)) +
1916  chiA2*(-5./(2.*(-1.+2.*eta)) + 15.*eta/(-1.+2.*eta) - 20.*eta2/(-1.+2.*eta));
1917  coeffs->f55v5c = 0; //RC: this is the calibration parameter which is initially set to 0.
1918  }
1919 
1920  }
1921  else{
1922  coeffs->f55v3 = 0;
1923  coeffs->f55v4 = 0;
1924  coeffs->f55v5c = 0;
1925  if(use_hm){
1926  //RC: This terms are in Eq.A12 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
1927  coeffs->f55v3 = chiA *(10./(3.*(-1.+2.*eta)) - 70.*eta/(3.*(-1.+2.*eta)) + 110.*eta2/(3.*(-1.+2.*eta)) );
1928  coeffs->f55v4 = chiA*chiS *(-5./(-1.+2.*eta) + 30.*eta/(-1.+2.*eta) - 40.*eta2/(-1.+2.*eta));
1929  coeffs->f55v5c = 0;
1930  }
1931  }
1932 
1933 
1934 
1935  coeffs->delta54vh3 = 8. / 15.;
1936  coeffs->delta54vh4 = 12. * aDelta / 5.;
1937 
1938  coeffs->rho54v2 = (-17448. + 96019. * eta - 127610. * eta2
1939  + 33320. * eta3) / (13650. * (1. - 5. * eta +
1940  5. * eta2));
1941  coeffs->rho54v3 = (-2. * a) / 15.;
1942  coeffs->rho54v4 = -16213384. / 15526875. + (2. * a2) / 5.;
1943 
1944  coeffs->delta53vh3 = 31. / 70.;
1945  if (dM2)
1946  {
1947 
1948  coeffs->rho53v2 =
1949  (375. - 850. * eta + 176. * eta2) / (390. * (-1. + 2. * eta));
1950  coeffs->rho53v3 = (-2. * a) / 3.;
1951  coeffs->rho53v4 = -410833. / 709800. + a2 / 2.;
1952  coeffs->rho53v5 = -103. * a / 325.;
1953  }
1954 
1955  coeffs->delta52vh3 = 4. / 15.;
1956  coeffs->delta52vh4 = 6. * aDelta / 5.;
1957 
1958  coeffs->rho52v2 = (-15828. + 84679. * eta - 104930. * eta2
1959  + 21980. * eta3) / (13650. * (1. - 5. * eta +
1960  5. * eta2));
1961  coeffs->rho52v3 = (-2. * a) / 15.;
1962  coeffs->rho52v4 = -7187914. / 15526875. + (2. * a2) / 5.;
1963 
1964  coeffs->delta51vh3 = 31. / 210.;
1965  if (dM2)
1966  {
1967 
1968  coeffs->rho51v2 =
1969  (319. - 626. * eta + 8. * eta2) / (390. * (-1. + 2. * eta));
1970  coeffs->rho51v3 = (-2. * a) / 3.;
1971  coeffs->rho51v4 = -31877. / 304200. + a2 / 2.;
1972  coeffs->rho51v5 = 139. * a / 975.;
1973  }
1974 
1975  /* l = 6, Eqs. A12a - A12f for rho, Eqs. E1c and E1d of PBFRT for delta */
1976 
1977  coeffs->delta66vh3 = 43. / 70.;
1978 
1979  coeffs->rho66v2 = (-106. + 602. * eta - 861. * eta2
1980  + 273. * eta3) / (84. * (1. - 5. * eta + 5. * eta2));
1981  coeffs->rho66v3 = (-2. * a) / 3.;
1982  coeffs->rho66v4 = -1025435. / 659736. + a2 / 2.;
1983 
1984  coeffs->delta65vh3 = 10. / 21.;
1985  if (dM2)
1986  {
1987 
1988  coeffs->rho65v2 = (-185. + 838. * eta - 910. * eta2
1989  + 220. * eta3) / (144. * (dM2 + 3. * eta2));
1990  coeffs->rho65v3 = -2. * a / 9.;
1991  }
1992 
1993  coeffs->delta64vh3 = 43. / 105.;
1994 
1995  coeffs->rho64v2 = (-86. + 462. * eta - 581. * eta2
1996  + 133. * eta3) / (84. * (1. - 5. * eta + 5. * eta2));
1997  coeffs->rho64v3 = (-2. * a) / 3.;
1998  coeffs->rho64v4 = -476887. / 659736. + a2 / 2.;
1999 
2000  coeffs->delta63vh3 = 2. / 7.;
2001  if (dM2)
2002  {
2003 
2004  coeffs->rho63v2 = (-169. + 742. * eta - 750. * eta2
2005  + 156. * eta3) / (144. * (dM2 + 3. * eta2));
2006  coeffs->rho63v3 = -2. * a / 9.;
2007  }
2008 
2009  coeffs->delta62vh3 = 43. / 210.;
2010 
2011  coeffs->rho62v2 = (-74. + 378. * eta - 413. * eta2
2012  + 49. * eta3) / (84. * (1. - 5. * eta + 5. * eta2));
2013  coeffs->rho62v3 = (-2. * a) / 3.;
2014  coeffs->rho62v4 = -817991. / 3298680. + a2 / 2.;
2015 
2016  coeffs->delta61vh3 = 2. / 21.;
2017  if (dM2)
2018  {
2019 
2020  coeffs->rho61v2 = (-161. + 694. * eta - 670. * eta2
2021  + 124. * eta3) / (144. * (dM2 + 3. * eta2));
2022  coeffs->rho61v3 = -2. * a / 9.;
2023  }
2024 
2025  /* l = 7, Eqs. A13a - A13g for rho, Eqs. E1e and E1f of PBFRT for delta */
2026  coeffs->delta77vh3 = 19. / 36.;
2027  if (dM2)
2028  {
2029 
2030  coeffs->rho77v2 = (-906. + 4246. * eta - 4963. * eta2
2031  + 1380. * eta3) / (714. * (dM2 + 3. * eta2));
2032  coeffs->rho77v3 = -2. * a / 3.;
2033  }
2034 
2035  coeffs->rho76v2 = (2144. - 16185. * eta + 37828. * eta2 - 29351. * eta3
2036  + 6104. * eta2 * eta2) / (1666. * (-1 + 7 * eta -
2037  14 * eta2 +
2038  7 * eta3));
2039 
2040  coeffs->delta75vh3 = 95. / 252.;
2041  if (dM2)
2042  {
2043 
2044  coeffs->rho75v2 = (-762. + 3382. * eta - 3523. * eta2
2045  + 804. * eta3) / (714. * (dM2 + 3. * eta2));
2046  coeffs->rho75v3 = -2. * a / 3.;
2047  }
2048 
2049  coeffs->rho74v2 = (17756. - 131805. * eta + 298872. * eta2 - 217959. * eta3
2050  + 41076. * eta2 * eta2) / (14994. * (-1. + 7. * eta -
2051  14. * eta2 +
2052  7. * eta3));
2053 
2054  coeffs->delta73vh3 = 19. / 84.;
2055  if (dM2)
2056  {
2057 
2058  coeffs->rho73v2 = (-666. + 2806. * eta - 2563. * eta2
2059  + 420. * eta3) / (714. * (dM2 + 3. * eta2));
2060  coeffs->rho73v3 = -2. * a / 3.;
2061  }
2062 
2063  coeffs->rho72v2 = (16832. - 123489. * eta + 273924. * eta2 - 190239. * eta3
2064  + 32760. * eta2 * eta2) / (14994. * (-1. + 7. * eta -
2065  14. * eta2 +
2066  7. * eta3));
2067 
2068  coeffs->delta71vh3 = 19. / 252.;
2069  if (dM2)
2070  {
2071 
2072  coeffs->rho71v2 = (-618. + 2518. * eta - 2083. * eta2
2073  + 228. * eta3) / (714. * (dM2 + 3. * eta2));
2074  coeffs->rho71v3 = -2. * a / 3.;
2075  }
2076 
2077  /* l = 8, Eqs. A14a - A14h */
2078 
2079  coeffs->rho88v2 = (3482. - 26778. * eta + 64659. * eta2 - 53445. * eta3
2080  + 12243. * eta2 * eta2) / (2736. * (-1. + 7. * eta -
2081  14. * eta2 +
2082  7. * eta3));
2083 
2084  if (dM2)
2085  {
2086  coeffs->rho87v2 =
2087  (23478. - 154099. * eta + 309498. * eta2 - 207550. * eta3 +
2088  38920 * eta2 * eta2) / (18240. * (-1 + 6 * eta - 10 * eta2 +
2089  4 * eta3));
2090  }
2091 
2092  coeffs->rho86v2 = (1002. - 7498. * eta + 17269. * eta2 - 13055. * eta3
2093  + 2653. * eta2 * eta2) / (912. * (-1. + 7. * eta -
2094  14. * eta2 +
2095  7. * eta3));
2096 
2097  if (dM2)
2098  {
2099  coeffs->rho85v2 = (4350. - 28055. * eta + 54642. * eta2 - 34598. * eta3
2100  + 6056. * eta2 * eta2) / (3648. * (-1. + 6. * eta -
2101  10. * eta2 +
2102  4. * eta3));
2103  }
2104 
2105  coeffs->rho84v2 = (2666. - 19434. * eta + 42627. * eta2 - 28965. * eta3
2106  + 4899. * eta2 * eta2) / (2736. * (-1. + 7. * eta -
2107  14. * eta2 +
2108  7. * eta3));
2109 
2110  if (dM2)
2111  {
2112  coeffs->rho83v2 =
2113  (20598. - 131059. * eta + 249018. * eta2 - 149950. * eta3 +
2114  24520. * eta2 * eta2) / (18240. * (-1. + 6. * eta - 10. * eta2 +
2115  4. * eta3));
2116  }
2117 
2118  coeffs->rho82v2 = (2462. - 17598. * eta + 37119. * eta2 - 22845. * eta3
2119  + 3063. * eta2 * eta2) / (2736. * (-1. + 7. * eta -
2120  14. * eta2 +
2121  7. * eta3));
2122 
2123  if (dM2)
2124  {
2125  coeffs->rho81v2 =
2126  (20022. - 126451. * eta + 236922. * eta2 - 138430. * eta3 +
2127  21640. * eta2 * eta2) / (18240. * (-1. + 6. * eta - 10. * eta2 +
2128  4. * eta3));
2129  }
2130 
2131  /* All relevant coefficients should be set, so we return */
2132  return XLAL_SUCCESS;
2133 }
2134 
2135 
2136 
2137 
2138 /**
2139  * This function calculates hlm mode factorized-resummed waveform
2140  * for given dynamical variables.
2141  * Eq. 17 and the entire Appendix of the paper.
2142  */
2143 static INT4
2144 XLALSimIMRSpinEOBGetSpinFactorizedWaveform (COMPLEX16 * restrict hlm,
2145  /**< OUTPUT, hlm waveforms */
2146  REAL8Vector * restrict values,
2147  /**< dyanmical variables */
2148  const REAL8 v,
2149  /**< velocity */
2150  const REAL8 Hreal,
2151  /**< real Hamiltonian */
2152  const INT4 l,
2153  /**< l mode index */
2154  const INT4 m,
2155  /**< m mode index */
2156  SpinEOBParams * restrict params,
2157  /**< Spin EOB parameters */
2158  INT4 use_optimized_v2
2159  /**< Spin EOB parameters */
2160  )
2161 {
2162  /* Status of function calls */
2163  INT4 status;
2164  INT4 i;
2165  INT4 use_hm = 0;
2166 
2167  REAL8 eta;
2168  REAL8 r, pp, Omega, v2, vh, vh3, k, hathatk, eulerlogxabs; //pr
2169  REAL8 Slm, deltalm, rholm;
2170  COMPLEX16 auxflm = 0.0;
2171  COMPLEX16 Tlm, rholmPwrl;
2172  COMPLEX16 hNewton;
2173  gsl_sf_result lnr1, arg1, z2;
2174 
2175  /* Non-Keplerian velocity */
2176  REAL8 vPhi, vPhi2;
2177 
2178  /* Pre-computed coefficients */
2179  FacWaveformCoeffs *hCoeffs = params->eobParams->hCoeffs;
2180 
2181  if (abs (m) > l)
2182  {
2183  XLAL_ERROR (XLAL_EINVAL);
2184  }
2185  if (m == 0)
2186  {
2187  XLAL_ERROR (XLAL_EINVAL);
2188  }
2189  use_hm = params->use_hm;
2190  eta = params->eobParams->eta;
2191 
2192  /* Check our eta was sensible */
2193  if (eta > 0.25 && eta < 0.25 +1e-4) {
2194  eta = 0.25;
2195  }
2196  if (eta > 0.25)
2197  {
2198  XLALPrintError
2199  ("XLAL Error - %s: Eta = %.16f seems to be > 0.25 - this isn't allowed!\n",
2200  __func__,eta);
2201  XLAL_ERROR (XLAL_EINVAL);
2202  }
2203  /*else if ( eta == 0.25 && m % 2 )
2204  {
2205  // If m is odd and dM = 0, hLM will be zero
2206  memset( hlm, 0, sizeof( COMPLEX16 ) );
2207  return XLAL_SUCCESS;
2208  } */
2209 
2210  r = values->data[0];
2211  //pr = values->data[2];
2212  pp = values->data[3];
2213 
2214  v2 = v * v;
2215  Omega = v2 * v;
2216  vh3 = Hreal * Omega;
2217  vh = cbrt (vh3);
2218  eulerlogxabs = LAL_GAMMA + log (2.0 * (REAL8) m * v);
2219 
2220  /* Calculate the non-Keplerian velocity */
2221  if (params->alignedSpins)
2222  {
2223  // YP: !!!!! SEOBNRv3devel temporary change !!!!!
2224  if (use_optimized_v2)
2225  {
2226  /* OPTIMIZED */
2227  vPhi =
2228  XLALSimIMRSpinAlignedEOBNonKeplerCoeffOptimized (values->data,
2229  params);
2230  /* END OPTIMIZED */
2231  }
2232  else
2233  {
2234  vPhi =
2235  XLALSimIMRSpinAlignedEOBNonKeplerCoeff (values->data, params);
2236  }
2237 #if YPPrecWave
2238  vPhi = XLALSimIMRSpinEOBNonKeplerCoeff (values->data, params);
2239 #endif
2240  // YP: !!!!! SEOBNRv3devel temporary change !!!!!
2241 
2242  if (XLAL_IS_REAL8_FAIL_NAN (vPhi))
2243  {
2244  XLAL_ERROR (XLAL_EFUNC);
2245  }
2246 
2247  vPhi = r * cbrt (vPhi);
2248  vPhi *= Omega;
2249  vPhi2 = vPhi * vPhi;
2250  }
2251  else
2252  {
2253  vPhi = v;
2254  vPhi2 = v2;
2255  }
2256  /* Calculate the newtonian multipole, 1st term in Eq. 17, given by Eq. A1 */
2257  // YP: !!!!! SEOBNRv3devel temporary change !!!!!
2258  status = XLALSimIMRSpinEOBCalculateNewtonianMultipole (&hNewton, vPhi2, r,
2259  values->data[1],
2260  (UINT4) l, m,
2261  params->eobParams);
2262 #if YPPrecWave
2263  status = XLALSimIMRSpinEOBCalculateNewtonianMultipole (&hNewton, vPhi2, r,
2264  values->data[12] +
2265  values->data[13],
2266  (UINT4) l, m,
2267  params->eobParams);
2268 #endif
2269  // YP: !!!!! SEOBNRv3devel temporary change !!!!!
2270  if (status == XLAL_FAILURE)
2271  {
2272  XLAL_ERROR (XLAL_EFUNC);
2273  }
2274  /* Calculate the source term, 2nd term in Eq. 17, given by Eq. A5 */
2275  if (((l + m) % 2) == 0)
2276  {
2277  Slm = (Hreal * Hreal - 1.) / (2. * eta) + 1.;
2278  }
2279  else
2280  {
2281  Slm = v * pp;
2282  }
2283  //printf( "Hreal = %e, Slm = %e, eta = %e\n", Hreal, Slm, eta );
2284 
2285  /* Calculate the Tail term, 3rd term in Eq. 17, given by Eq. A6 */
2286  k = m * Omega;
2287  hathatk = Hreal * k;
2288  XLAL_CALLGSL (status =
2289  gsl_sf_lngamma_complex_e (l + 1.0, -2.0 * hathatk, &lnr1,
2290  &arg1));
2291  if (status != GSL_SUCCESS)
2292  {
2293  XLALPrintError ("XLAL Error - %s: Error in GSL function: %s\n",
2294  __func__, gsl_strerror (status));
2295  XLAL_ERROR (XLAL_EFUNC);
2296  }
2297  XLAL_CALLGSL (status = gsl_sf_fact_e (l, &z2));
2298  if (status != GSL_SUCCESS)
2299  {
2300  XLALPrintError ("XLAL Error - %s: Error in GSL function: %s\n",
2301  __func__, gsl_strerror (status));
2302  XLAL_ERROR (XLAL_EFUNC);
2303  }
2304  Tlm =
2305  cexp ((lnr1.val + LAL_PI * hathatk) +
2306  I * (arg1.val + 2.0 * hathatk * log (4.0 * k / sqrt (LAL_E))));
2307  Tlm /= z2.val;
2308 
2309 
2310  /* Calculate the residue phase and amplitude terms */
2311  /* deltalm is the 4th term in Eq. 17, delta 22 given by Eq. A15, others */
2312  /* rholm is the 5th term in Eq. 17, given by Eqs. A8 - A14 */
2313  /* auxflm is a special part of the 5th term in Eq. 17, given by Eq. A15 */
2314  /* Actual values of the coefficients are defined in the next function of this file */
2315  switch (l)
2316  {
2317  case 2:
2318  switch (abs (m))
2319  {
2320  case 2:
2321  deltalm = vh3 * (hCoeffs->delta22vh3 + vh3 * (hCoeffs->delta22vh6
2322  +
2323  vh * vh *
2324  (hCoeffs->delta22vh9 *
2325  vh))) +
2326  hCoeffs->delta22v5 * v * v2 * v2 +
2327  hCoeffs->delta22v6 * v2 * v2 * v2 +
2328  hCoeffs->delta22v8 * v2 * v2 * v2 * v2;
2329  rholm =
2330  1. + v2 * (hCoeffs->rho22v2 +
2331  v * (hCoeffs->rho22v3 +
2332  v * (hCoeffs->rho22v4 +
2333  v * (hCoeffs->rho22v5 +
2334  v * (hCoeffs->rho22v6 +
2335  hCoeffs->rho22v6l * eulerlogxabs +
2336  v * (hCoeffs->rho22v7 +
2337  v * (hCoeffs->rho22v8 +
2338  hCoeffs->rho22v8l *
2339  eulerlogxabs +
2340  (hCoeffs->rho22v10 +
2341  hCoeffs->rho22v10l *
2342  eulerlogxabs) *
2343  v2)))))));
2344  break;
2345  case 1:
2346  {
2347  deltalm = vh3 * (hCoeffs->delta21vh3 + vh3 * (hCoeffs->delta21vh6
2348  +
2349  vh *
2350  (hCoeffs->
2351  delta21vh7 +
2352  (hCoeffs->
2353  delta21vh9) * vh *
2354  vh))) +
2355  hCoeffs->delta21v5 * v * v2 * v2 +
2356  hCoeffs->delta21v7 * v2 * v2 * v2 * v;
2357  rholm =
2358  1. + v * (hCoeffs->rho21v1 +
2359  v * (hCoeffs->rho21v2 +
2360  v * (hCoeffs->rho21v3 +
2361  v * (hCoeffs->rho21v4 +
2362  v * (hCoeffs->rho21v5 +
2363  v * (hCoeffs->rho21v6 +
2364  hCoeffs->rho21v6l *
2365  eulerlogxabs +
2366  v * (hCoeffs->rho21v7 +
2367  hCoeffs->rho21v7l *
2368  eulerlogxabs +
2369  v * (hCoeffs->rho21v8 +
2370  hCoeffs->rho21v8l *
2371  eulerlogxabs +
2372  (hCoeffs->
2373  rho21v10 +
2374  hCoeffs->
2375  rho21v10l *
2376  eulerlogxabs) *
2377  v2))))))));
2378  if (use_hm){
2379  //RC: This terms are in Eq.A11 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
2380  auxflm = v * (hCoeffs->f21v1 + v2 * (hCoeffs->f21v3 + v * hCoeffs->f21v4 + v2 * (hCoeffs->f21v5 + v * hCoeffs->f21v6 + v2*hCoeffs->f21v7c)));
2381  }
2382  else{
2383  auxflm = v * hCoeffs->f21v1;
2384  }
2385  }
2386  break;
2387  default:
2388  XLAL_ERROR (XLAL_EINVAL);
2389  break;
2390  }
2391  break;
2392  case 3:
2393  switch (m)
2394  {
2395  case 3:
2396  deltalm =
2397  vh3 * (hCoeffs->delta33vh3 +
2398  vh3 * (hCoeffs->delta33vh6 + hCoeffs->delta33vh9 * vh3)) +
2399  hCoeffs->delta33v5 * v * v2 * v2;
2400  if(use_hm){
2401  //RC: This terms are in Eq.A6 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
2402  rholm =
2403  1. + v2 * (hCoeffs->rho33v2 +
2404  v * (hCoeffs->rho33v3 +
2405  v * (hCoeffs->rho33v4 +
2406  v * (hCoeffs->rho33v5 +
2407  v * (hCoeffs->rho33v6 +
2408  hCoeffs->rho33v6l * eulerlogxabs +
2409  v * (hCoeffs->rho33v7 +
2410  v * (hCoeffs->rho33v8 +
2411  hCoeffs->rho33v8l *
2412  eulerlogxabs +
2413  v2*(hCoeffs->rho33v10 + hCoeffs->rho33v10l*eulerlogxabs))))))));
2414  //RC: This terms are in Eq.A10 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
2415  auxflm = v * (v2 * (hCoeffs->f33v3 + v * (hCoeffs->f33v4 + v * (hCoeffs->f33v5 + v * hCoeffs->f33v6)))) + _Complex_I * vh3 * vh3 * hCoeffs->f33vh6;
2416  }
2417  else
2418  {
2419  rholm =
2420  1. + v2 * (hCoeffs->rho33v2 +
2421  v * (hCoeffs->rho33v3 +
2422  v * (hCoeffs->rho33v4 +
2423  v * (hCoeffs->rho33v5 +
2424  v * (hCoeffs->rho33v6 +
2425  hCoeffs->rho33v6l * eulerlogxabs +
2426  v * (hCoeffs->rho33v7 +
2427  v * (hCoeffs->rho33v8 +
2428  hCoeffs->rho33v8l *
2429  eulerlogxabs)))))));
2430  auxflm = v * v2 * hCoeffs->f33v3;
2431  }
2432  break;
2433  case 2:
2434  deltalm =
2435  vh3 * (hCoeffs->delta32vh3 +
2436  vh * (hCoeffs->delta32vh4 +
2437  vh * vh * (hCoeffs->delta32vh6 +
2438  hCoeffs->delta32vh9 * vh3)));
2439  rholm =
2440  1. + v * (hCoeffs->rho32v +
2441  v * (hCoeffs->rho32v2 +
2442  v * (hCoeffs->rho32v3 +
2443  v * (hCoeffs->rho32v4 +
2444  v * (hCoeffs->rho32v5 +
2445  v * (hCoeffs->rho32v6 +
2446  hCoeffs->rho32v6l *
2447  eulerlogxabs +
2448  (hCoeffs->rho32v8 +
2449  hCoeffs->rho32v8l *
2450  eulerlogxabs) * v2))))));
2451  break;
2452  case 1:
2453  deltalm = vh3 * (hCoeffs->delta31vh3 + vh3 * (hCoeffs->delta31vh6
2454  +
2455  vh *
2456  (hCoeffs->delta31vh7 +
2457  hCoeffs->delta31vh9 *
2458  vh * vh))) +
2459  hCoeffs->delta31v5 * v * v2 * v2;
2460  rholm =
2461  1. + v2 * (hCoeffs->rho31v2 +
2462  v * (hCoeffs->rho31v3 +
2463  v * (hCoeffs->rho31v4 +
2464  v * (hCoeffs->rho31v5 +
2465  v * (hCoeffs->rho31v6 +
2466  hCoeffs->rho31v6l * eulerlogxabs +
2467  v * (hCoeffs->rho31v7 +
2468  (hCoeffs->rho31v8 +
2469  hCoeffs->rho31v8l *
2470  eulerlogxabs) * v))))));
2471  auxflm = v * v2 * hCoeffs->f31v3;
2472  break;
2473  default:
2474  XLAL_ERROR (XLAL_EINVAL);
2475  break;
2476  }
2477  break;
2478  case 4:
2479  switch (m)
2480  {
2481  case 4:
2482  if(use_hm){
2483  //RC: This terms are in Eq.A15 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
2484  deltalm = vh3 * (hCoeffs->delta44vh3 + vh3 * (hCoeffs->delta44vh6 + vh3 * hCoeffs->delta44vh9))
2485  + hCoeffs->delta44v5 * v2 * v2 * v;
2486  //RC: This terms are in Eq.A8 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
2487  rholm = 1. + v2 * (hCoeffs->rho44v2
2488  + v * (hCoeffs->rho44v3 + v * (hCoeffs->rho44v4
2489  +
2490  v *
2491  (hCoeffs->
2492  rho44v5 +
2493  v * (hCoeffs->
2494  rho44v6 +
2495  hCoeffs->
2496  rho44v6l *
2497  eulerlogxabs +
2498  v2 *( hCoeffs->rho44v8 + hCoeffs->rho44v8l * eulerlogxabs
2499  +v2 * (hCoeffs->rho44v10 + hCoeffs->rho44v10l * eulerlogxabs) ) )))));
2500  }
2501  else{
2502  deltalm = vh3 * (hCoeffs->delta44vh3 + hCoeffs->delta44vh6 * vh3)
2503  + hCoeffs->delta44v5 * v2 * v2 * v;
2504 
2505  rholm = 1. + v2 * (hCoeffs->rho44v2
2506  + v * (hCoeffs->rho44v3 + v * (hCoeffs->rho44v4
2507  +
2508  v *
2509  (hCoeffs->
2510  rho44v5 +
2511  (hCoeffs->
2512  rho44v6 +
2513  hCoeffs->
2514  rho44v6l *
2515  eulerlogxabs) *
2516  v))));
2517  }
2518 
2519  break;
2520  case 3:
2521  deltalm = vh3 * (hCoeffs->delta43vh3 + vh * (hCoeffs->delta43vh4
2522  +
2523  hCoeffs->delta43vh6 *
2524  vh * vh));
2525  rholm =
2526  1. + v * (hCoeffs->rho43v +
2527  v * (hCoeffs->rho43v2 +
2528  v2 * (hCoeffs->rho43v4 +
2529  v * (hCoeffs->rho43v5 +
2530  (hCoeffs->rho43v6 +
2531  hCoeffs->rho43v6l * eulerlogxabs) *
2532  v))));
2533  auxflm = v * hCoeffs->f43v;
2534  break;
2535  case 2:
2536  deltalm = vh3 * (hCoeffs->delta42vh3 + hCoeffs->delta42vh6 * vh3);
2537  rholm = 1. + v2 * (hCoeffs->rho42v2
2538  + v * (hCoeffs->rho42v3 +
2539  v * (hCoeffs->rho42v4 +
2540  v * (hCoeffs->rho42v5 +
2541  (hCoeffs->rho42v6 +
2542  hCoeffs->rho42v6l *
2543  eulerlogxabs) * v))));
2544  break;
2545  case 1:
2546  deltalm = vh3 * (hCoeffs->delta41vh3 + vh * (hCoeffs->delta41vh4
2547  +
2548  hCoeffs->delta41vh6 *
2549  vh * vh));
2550  rholm =
2551  1. + v * (hCoeffs->rho41v +
2552  v * (hCoeffs->rho41v2 +
2553  v2 * (hCoeffs->rho41v4 +
2554  v * (hCoeffs->rho41v5 +
2555  (hCoeffs->rho41v6 +
2556  hCoeffs->rho41v6l * eulerlogxabs) *
2557  v))));
2558  auxflm = v * hCoeffs->f41v;
2559  break;
2560  default:
2561  XLAL_ERROR (XLAL_EINVAL);
2562  break;
2563  }
2564  break;
2565  case 5:
2566  switch (m)
2567  {
2568  case 5:
2569  if(use_hm){
2570  //RC: This terms are in Eq.A16 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
2571  deltalm =
2572  vh3 *(hCoeffs->delta55vh3 +vh3*(hCoeffs->delta55vh6 +vh3 *(hCoeffs->delta55vh9)))
2573  + hCoeffs->delta55v5 * v2 * v2 * v;
2574  //RC: This terms are in Eq.A9 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
2575  rholm =
2576  1. + v2 * (hCoeffs->rho55v2 +
2577  v * (hCoeffs->rho55v3 +
2578  v * (hCoeffs->rho55v4 +
2579  v * (hCoeffs->rho55v5 +
2580  v * (hCoeffs->rho55v6 + hCoeffs->rho55v6l * eulerlogxabs +
2581  v2 * (hCoeffs->rho55v8 + hCoeffs->rho55v8l * eulerlogxabs +
2582  v2 * (hCoeffs->rho55v10 + hCoeffs->rho55v10l * eulerlogxabs )))))));
2583  //RC: This terms are in Eq.A12 in https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084028
2584  auxflm = v2 * v * (hCoeffs->f55v3 + v * (hCoeffs->f55v4 + v * (hCoeffs->f55v5c) ));
2585  }
2586  else
2587  {
2588  deltalm =
2589  hCoeffs->delta55vh3 * vh3 + hCoeffs->delta55v5 * v2 * v2 * v;
2590  rholm =
2591  1. + v2 * (hCoeffs->rho55v2 +
2592  v * (hCoeffs->rho55v3 +
2593  v * (hCoeffs->rho55v4 +
2594  v * (hCoeffs->rho55v5 +
2595  hCoeffs->rho55v6 * v))));
2596  }
2597  break;
2598  case 4:
2599  deltalm = vh3 * (hCoeffs->delta54vh3 + hCoeffs->delta54vh4 * vh);
2600  rholm = 1. + v2 * (hCoeffs->rho54v2 + v * (hCoeffs->rho54v3
2601  + hCoeffs->rho54v4 * v));
2602  break;
2603  case 3:
2604  deltalm = hCoeffs->delta53vh3 * vh3;
2605  rholm = 1. + v2 * (hCoeffs->rho53v2
2606  + v * (hCoeffs->rho53v3 +
2607  v * (hCoeffs->rho53v4 +
2608  hCoeffs->rho53v5 * v)));
2609  break;
2610  case 2:
2611  deltalm = vh3 * (hCoeffs->delta52vh3 + hCoeffs->delta52vh4 * vh);
2612  rholm = 1. + v2 * (hCoeffs->rho52v2 + v * (hCoeffs->rho52v3
2613  + hCoeffs->rho52v4 * v));
2614  break;
2615  case 1:
2616  deltalm = hCoeffs->delta51vh3 * vh3;
2617  rholm = 1. + v2 * (hCoeffs->rho51v2
2618  + v * (hCoeffs->rho51v3 +
2619  v * (hCoeffs->rho51v4 +
2620  hCoeffs->rho51v5 * v)));
2621  break;
2622  default:
2623  XLAL_ERROR (XLAL_EINVAL);
2624  break;
2625  }
2626  break;
2627  case 6:
2628  switch (m)
2629  {
2630  case 6:
2631  deltalm = hCoeffs->delta66vh3 * vh3;
2632  rholm = 1. + v2 * (hCoeffs->rho66v2 + v * (hCoeffs->rho66v3
2633  + hCoeffs->rho66v4 * v));
2634  break;
2635  case 5:
2636  deltalm = hCoeffs->delta65vh3 * vh3;
2637  rholm = 1. + v2 * (hCoeffs->rho65v2 + hCoeffs->rho65v3 * v);
2638  break;
2639  case 4:
2640  deltalm = hCoeffs->delta64vh3 * vh3;
2641  rholm = 1. + v2 * (hCoeffs->rho64v2 + v * (hCoeffs->rho64v3
2642  + hCoeffs->rho64v4 * v));
2643  break;
2644  case 3:
2645  deltalm = hCoeffs->delta63vh3 * vh3;
2646  rholm = 1. + v2 * (hCoeffs->rho63v2 + hCoeffs->rho63v3 * v);
2647  break;
2648  case 2:
2649  deltalm = hCoeffs->delta62vh3 * vh3;
2650  rholm = 1. + v2 * (hCoeffs->rho62v2 + v * (hCoeffs->rho62v3
2651  + hCoeffs->rho62v4 * v));
2652  break;
2653  case 1:
2654  deltalm = hCoeffs->delta61vh3 * vh3;
2655  rholm = 1. + v2 * (hCoeffs->rho61v2 + hCoeffs->rho61v3 * v);
2656  break;
2657  default:
2658  XLAL_ERROR (XLAL_EINVAL);
2659  break;
2660  }
2661  break;
2662  case 7:
2663  switch (m)
2664  {
2665  case 7:
2666  deltalm = hCoeffs->delta77vh3 * vh3;
2667  rholm = 1. + v2 * (hCoeffs->rho77v2 + hCoeffs->rho77v3 * v);
2668  break;
2669  case 6:
2670  deltalm = 0.0;
2671  rholm = 1. + hCoeffs->rho76v2 * v2;
2672  break;
2673  case 5:
2674  deltalm = hCoeffs->delta75vh3 * vh3;
2675  rholm = 1. + v2 * (hCoeffs->rho75v2 + hCoeffs->rho75v3 * v);
2676  break;
2677  case 4:
2678  deltalm = 0.0;
2679  rholm = 1. + hCoeffs->rho74v2 * v2;
2680  break;
2681  case 3:
2682  deltalm = hCoeffs->delta73vh3 * vh3;
2683  rholm = 1. + v2 * (hCoeffs->rho73v2 + hCoeffs->rho73v3 * v);
2684  break;
2685  case 2:
2686  deltalm = 0.0;
2687  rholm = 1. + hCoeffs->rho72v2 * v2;
2688  break;
2689  case 1:
2690  deltalm = hCoeffs->delta71vh3 * vh3;
2691  rholm = 1. + v2 * (hCoeffs->rho71v2 + hCoeffs->rho71v3 * v);
2692  break;
2693  default:
2694  XLAL_ERROR (XLAL_EINVAL);
2695  break;
2696  }
2697  break;
2698  case 8:
2699  switch (m)
2700  {
2701  case 8:
2702  deltalm = 0.0;
2703  rholm = 1. + hCoeffs->rho88v2 * v2;
2704  break;
2705  case 7:
2706  deltalm = 0.0;
2707  rholm = 1. + hCoeffs->rho87v2 * v2;
2708  break;
2709  case 6:
2710  deltalm = 0.0;
2711  rholm = 1. + hCoeffs->rho86v2 * v2;
2712  break;
2713  case 5:
2714  deltalm = 0.0;
2715  rholm = 1. + hCoeffs->rho85v2 * v2;
2716  break;
2717  case 4:
2718  deltalm = 0.0;
2719  rholm = 1. + hCoeffs->rho84v2 * v2;
2720  break;
2721  case 3:
2722  deltalm = 0.0;
2723  rholm = 1. + hCoeffs->rho83v2 * v2;
2724  break;
2725  case 2:
2726  deltalm = 0.0;
2727  rholm = 1. + hCoeffs->rho82v2 * v2;
2728  break;
2729  case 1:
2730  deltalm = 0.0;
2731  rholm = 1. + hCoeffs->rho81v2 * v2;
2732  break;
2733  default:
2734  XLAL_ERROR (XLAL_EINVAL);
2735  break;
2736  }
2737  break;
2738  default:
2739  XLAL_ERROR (XLAL_EINVAL);
2740  break;
2741  }
2742 
2743  /* Raise rholm to the lth power */
2744  rholmPwrl = 1.0;
2745  for(i = 0 ; i < l ; i++) rholmPwrl *= rholm;
2746  /* In the equal-mass odd m case, there is no contribution from nonspin terms,
2747  * and the only contribution comes from the auxflm term that is proportional to chiA (asymmetric spins).
2748  * In this case, we must ignore the nonspin terms directly, since the leading term defined by
2749  * CalculateThisMultipolePrefix in LALSimIMREOBNewtonianMultipole.c is not zero (see comments there).
2750  */
2751  if (eta == 0.25 && m % 2)
2752  {
2753  rholmPwrl = auxflm;
2754  }
2755  else
2756  {
2757  rholmPwrl += auxflm;
2758  }
2759  /* Put all factors in Eq. 17 together */
2760  *hlm = Tlm * cexp (I * deltalm) * Slm * rholmPwrl;
2761  *hlm *= hNewton;
2762  /*if (r > 8.5)
2763  {
2764  printf("YP::FullWave: Reh = %.16e, Imh = %.16e, hAmp = %.16e, hPhi = %.16e\n",creal(*hlm),cimag(*hlm),cabs(*hlm),carg(*hlm));
2765  } */
2766  return XLAL_SUCCESS;
2767 }
2768 
2769 /**
2770  * This function calculates tidal correction to the hlm mode factorized-resummed waveform
2771  * for given dynamical variables.
2772  */
2773 static INT4
2774 XLALSimIMRSpinEOBWaveformTidal (COMPLEX16 * restrict hlm,
2775  /**< OUTPUT, hlm waveforms */
2776  REAL8Vector * restrict values,
2777  /**< dyanmical variables */
2778  const REAL8 v,
2779  /**< velocity */
2780  const INT4 l,
2781  /**< l mode index */
2782  const INT4 m,
2783  /**< m mode index */
2784  SpinEOBParams * restrict params
2785  /**< Spin EOB parameters */
2786 )
2787 {
2788  REAL8 eta;
2789  if (abs (m) > l)
2790  {
2791  XLAL_ERROR (XLAL_EINVAL);
2792  }
2793  if (m == 0)
2794  {
2795  XLAL_ERROR (XLAL_EINVAL);
2796  }
2797 
2798  eta = params->eobParams->eta;
2799 
2800  /* Check our eta was sensible */
2801  if (eta > 0.25 && eta < 0.25 +1e-4) {
2802  eta = 0.25;
2803  }
2804  if (eta > 0.25)
2805  {
2806  XLALPrintError
2807  ("XLAL Error - %s: Eta = %.16f seems to be > 0.25 - this isn't allowed!\n",
2808  __func__,eta);
2809  XLAL_ERROR (XLAL_EINVAL);
2810  }
2811 
2812  COMPLEX16 hTidal = XLALhTidal( l, m, values->data[1], v, eta, params->seobCoeffs->tidal1, params->seobCoeffs->tidal2 );
2813  *hlm = hTidal;
2814  return XLAL_SUCCESS;
2815 }
2816 
2817 #endif /* _LALSIMIMRSPINEOBFACTORIZEDWAVEFORM */
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
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
XLALSimIMRSpinEOBWaveformTidal
static INT4 XLALSimIMRSpinEOBWaveformTidal(COMPLEX16 *restrict hlm, REAL8Vector *restrict values, const REAL8 v, const INT4 l, const INT4 m, SpinEOBParams *restrict params)
This function calculates tidal correction to the hlm mode factorized-resummed waveform for given dyna...
Definition: LALSimIMRSpinEOBFactorizedWaveform.c:2774
XLALSimIMRTEOBk2effMode
static REAL8 XLALSimIMRTEOBk2effMode(REAL8 Omega, REAL8 k2TidaleffHam, REAL8 omega02Tidal, REAL8 XCompanion)
Function to compute the dynamical enhancement of the quadrupolar Love number that enters the tidal co...
Definition: LALSimIMRSpinEOBFactorizedWaveform.c:132
XLALSimIMREOBCalcCalibCoefficientHigherModes
static UNUSED int XLALSimIMREOBCalcCalibCoefficientHigherModes(FacWaveformCoeffs *const coeffs, SpinEOBParams *restrict UNUSED params, const UINT4 modeL, const UINT4 modeM, REAL8Vector *restrict phaseVec, REAL8Vector *restrict rVec, REAL8Vector *restrict prVec, REAL8Vector *restrict orbOmegaVec, REAL8Vector *restrict HamVec, REAL8Vector *restrict pPhiVec, const REAL8 timePeak, const REAL8 m1, const REAL8 m2, const REAL8 UNUSED chiS, const REAL8 UNUSED chiA, const REAL8 UNUSED deltaT)
This function calculate the calibration parameter for the amplitude of the factorized-resummed wavefo...
Definition: LALSimIMRSpinEOBFactorizedWaveform.c:635
XLALSimIMRSpinEOBFluxGetSpinFactorizedWaveform
static INT4 XLALSimIMRSpinEOBFluxGetSpinFactorizedWaveform(COMPLEX16 *restrict hlm, REAL8Vector *restrict values, const REAL8 v, const REAL8 Hreal, const INT4 l, const INT4 m, SpinEOBParams *restrict params, INT4 use_optimized_v2, REAL8 *vPhil2m2)
This function calculates hlm mode factorized-resummed waveform for given dynamical variables.
Definition: LALSimIMRSpinEOBFactorizedWaveform.c:804
XLALSimIMRSpinEOBGetAmplitudeResidual
static UNUSED INT4 XLALSimIMRSpinEOBGetAmplitudeResidual(COMPLEX16 *restrict rholmpwrl, const REAL8 v, const REAL8 Hreal, const INT4 modeL, const INT4 modeM, SpinEOBParams *restrict params)
This function calculate the residue amplitude terms.
Definition: LALSimIMRSpinEOBFactorizedWaveform.c:227
XLALTEOBbeta221
static REAL8 XLALTEOBbeta221(REAL8 X)
Function to compute PN corrections to the tidal term that enters the (2,2) amplitude.
Definition: LALSimIMRSpinEOBFactorizedWaveform.c:124
XLALSimIMREOBCalcSpinFacWaveformCoefficients
static int XLALSimIMREOBCalcSpinFacWaveformCoefficients(FacWaveformCoeffs *const coeffs, SpinEOBParams *restrict params, const REAL8 m1, const REAL8 m2, const REAL8 eta, const REAL8 a, const REAL8 chiS, const REAL8 chiA, const UINT4 SpinAlignedEOBversion)
Spin Factors.
Definition: LALSimIMRSpinEOBFactorizedWaveform.c:1330
XLALhTidal
static COMPLEX16 XLALhTidal(INT4 l, INT4 m, REAL8 phase, REAL8 v, REAL8 eta, TidalEOBParams *tidal1, TidalEOBParams *tidal2)
Function to compute the tidal correction to the mode amplitudes.
Definition: LALSimIMRSpinEOBFactorizedWaveform.c:145
XLALSimIMRSpinEOBGetSpinFactorizedWaveform
static UNUSED INT4 XLALSimIMRSpinEOBGetSpinFactorizedWaveform(COMPLEX16 *restrict hlm, REAL8Vector *restrict values, const REAL8 v, const REAL8 Hreal, const INT4 l, const INT4 m, SpinEOBParams *restrict params, INT4 use_optimized_v2)
This function calculates hlm mode factorized-resummed waveform for given dynamical variables.
Definition: LALSimIMRSpinEOBFactorizedWaveform.c:2144
XLALSimIMRSpinEOBNonKeplerCoeff
static REAL8 UNUSED XLALSimIMRSpinEOBNonKeplerCoeff(const REAL8 values[], SpinEOBParams *funcParams)
Function to calculate the non-Keplerian coefficient for the spin-aligned EOB model.
Definition: LALSimIMRSpinEOBHamiltonian.c:1617
XLALSimIMRTEOBkleff
static REAL8 XLALSimIMRTEOBkleff(INT4 l, REAL8 u, REAL8 eta, TidalEOBParams *tidal)
Function to compute the enhancement of kltidal due to the presence of f-mode resonance Implements (11...
Definition: LALSimIMRSpinEOBHamiltonian.c:392
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
XLALSimIMRSpinAlignedEOBNonKeplerCoeffOptimized
static REAL8 XLALSimIMRSpinAlignedEOBNonKeplerCoeffOptimized(const REAL8 values[], SpinEOBParams *funcParams)
Function to calculate the non-Keplerian coefficient for the spin-aligned EOB model.
Definition: LALSimIMRSpinEOBHamiltonianOptimized.c:179
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
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
u
const double u
Definition: exact_derivatives-HrealHeader.c:17
a2
const double a2
Definition: exact_derivatives-HrealHeader.c:33
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
q
static const INT4 q
a
static const INT4 a
XLALCreateCOMPLEX16Vector
COMPLEX16Vector * XLALCreateCOMPLEX16Vector(UINT4 length)
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLALDestroyCOMPLEX16Vector
void XLALDestroyCOMPLEX16Vector(COMPLEX16Vector *vector)
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_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
REAL8Vector::data
REAL8 * data
SpinEOBParams
Parameters for the spinning EOB model.
Definition: LALSimIMRSpinEOB.h:116
TidalEOBParams
Tidal parameters for EOB model of NS: mByM - dimensionless ratio m_{NS}/M lambda2Tidal - dimensionles...
Definition: LALSimIMREOBNRv2.h:60
params
Definition: burst.c:245
deltaT
double deltaT
Definition: unicorn.c:24