LALSimulation  5.4.0.1-fe68b98
LALSimIMRSpinEOBFactorizedWaveformPrec_v3opt.c
Go to the documentation of this file.
1 /**
2  * \author Craig Robinson, Yi Pan, Prayush Kumar, Stas Babak, Andrea Taracchini
3  *
4  * \brief Function to compute the factorized waveform as used in the SEOBNRv1 model.
5  * Waveform expressions are given by
6  * Taracchini et al. ( PRD 86, 024011 (2012), arXiv 1202.0790 ).
7  * All equation numbers in this file refer to equations of this paper,
8  * unless otherwise specified.
9  * Coefficients of the so-called "deltalm" terms are given by
10  * Damour et al. PRD 79, 064004 (2009) and Pan et al. PRD 83, 064003 (2011),
11  * henceforth DIN and PBFRT.
12  *
13  * Functions to compute the factorized waveform for the purpose of computing the
14  * flux, i.e. returning only the absolute value of the multipoles. The tail term
15  * Tlm is used in its resummed form, given by Eq. (42) of Damour, Nagar and
16  * Bernuzzi, PRD 87 084035 (2013), called DNB here.
17  *
18  */
19 
20 #ifndef _LALSIMIMRSPINPRECEOBFACTORIZEDWAVEFORM_V3OPT_C
21 #define _LALSIMIMRSPINPRECEOBFACTORIZEDWAVEFORM_V3OPT_C
22 #include <gsl/gsl_deriv.h>
23 #include <complex.h>
24 #include <lal/LALSimInspiral.h>
25 #include <lal/LALSimIMR.h>
26 
27 #include "LALSimIMREOBNRv2.h"
28 #include "LALSimIMRSpinEOB.h"
29 
35 
38 /* #include "LALSimIMRSpinEOBFactorizedWaveform.c" */
39 /*------------------------------------------------------------------------------------------
40  *
41  * Prototypes of functions defined in this code.
42  *
43  *------------------------------------------------------------------------------------------
44  */
45 static INT4
47  COMPLEX16 * restrict hlm, /**< OUTPUT, hlm waveforms */
48  REAL8Vector * restrict values, /**< dyanmical variables: \f$(r,\phi,p_r,p_\phi)\f$ */
49  REAL8Vector * restrict cartvalues, /**< dyanmical variables */
50  const REAL8 v, /**< velocity */
51  const REAL8 Hreal, /**< real Hamiltonian */
52  const INT4 l, /**< l mode index */
53  const INT4 m, /**< m mode index */
54  SpinEOBParams * restrict params /**< Spin EOB parameters */
55 );
56 
57 static INT4
59  COMPLEX16 * restrict hlmTab, /**< OUTPUT, hlm waveforms */
60  REAL8Vector * restrict values, /**< dyanmical variables: \f$(r,\phi,p_r,p_\phi)\f$ */
61  REAL8Vector * restrict cartvalues, /**< dyanmical variables */
62  const REAL8 v, /**< velocity */
63  const REAL8 Hreal, /**< real Hamiltonian */
64  const INT4 lMax, /**< maximum l mode to compute, compute 0 < m <= lMax */
65  SpinEOBParams * restrict params /**< Spin EOB parameters */
66 );
67 
68 
69 static REAL8
71  const REAL8 values[], /**<< Dynamical variables */
72  SpinEOBParams *funcParams /**<< EOB parameters */
73  );
74 
76  const REAL8 values[], /**<< Dynamical variables */
77  SpinEOBParams *funcParams /**<< EOB parameters */
78  );
79 
80 /*------------------------------------------------------------------------------------------
81  *
82  * Defintions of functions.
83  *
84  *------------------------------------------------------------------------------------------
85  */
86 
87 
88 /** ######################################################### */
89 /** FOR PRECESSING EOB
90  * This function calculates hlm mode factorized-resummed waveform
91  * for given dynamical variables. This is optimized for flux calculation,
92  * by ignoring complex arguments and keeping only absolute values.
93  * Changes:
94  * (i) Complex Argument of Tlm not exponentiated.
95  * (ii) exp(i deltalm) set to 1.
96  * Eq. 17 and the entire Appendix of PRD 86, 024011 (2012) + changes
97  * described in the section "Factorized waveforms" of https://dcc.ligo.org/T1400476
98  */
99 static INT4
101  COMPLEX16 * restrict hlmTab, /**< OUTPUT, hlm waveforms */
102  REAL8Vector * restrict values, /**< dyanmical variables: \f$(r,\phi,p_r,p_\phi)\f$ */
103  REAL8Vector * restrict cartvalues, /**< dyanmical variables */
104  const REAL8 v, /**< velocity */
105  const REAL8 Hreal, /**< real Hamiltonian */
106  const INT4 lMax, /**< maximum l mode to compute, compute 0 < m <= lMax */
107  SpinEOBParams * restrict params /**< Spin EOB parameters */
108 )
109 {
110  int debugPK = 0;
111  const REAL8 vPhiKepler = params->alignedSpins ?
114  if (XLAL_IS_REAL8_FAIL_NAN(vPhiKepler)) {
116  }
117 
118  for (INT4 l = 2; l <= lMax; l++)
119  {
120  for (INT4 m = 1; m <= l; m++)
121  {
122  COMPLEX16 *hlm = &hlmTab[l*(lMax+1)+m];
123 
124 
125  /* Status of function calls */
126  INT4 status;
127  INT4 i;
128 
129  REAL8 eta;
130  REAL8 UNUSED r , pp, Omega, v2, /* vh, vh3, */ k, hathatk, eulerlogxabs;
131  //pr
132  REAL8 UNUSED rcrossp_x, rcrossp_y, rcrossp_z;
133  REAL8 Slm , rholm, rholmPwrl;
134  REAL8 auxflm = 0.0;
135  REAL8 hathatksq4, hathatk4pi, Tlmprefac, Tlmprodfac;
136  REAL8 Tlm;
137  COMPLEX16 hNewton;
138  gsl_sf_result z2;
139 
140  /* Non-Keplerian velocity */
141  REAL8 vPhi , vPhi2;
142 
143  /* Pre-computed coefficients */
144  FacWaveformCoeffs *hCoeffs = params->eobParams->hCoeffs;
145 
146  if (abs(m) > (INT4) l) {
148  }
149  if (m == 0) {
151  }
152  eta = params->eobParams->eta;
153 
154  /* Check our eta was sensible */
155  if (eta > 0.25) {
156  XLALPrintError("XLAL Error - %s: Eta seems to be > 0.25 - this isn't allowed!\n", __func__);
158  }
159  /*
160  * else if ( eta == 0.25 && m % 2 ) { // If m is odd and dM = 0, hLM
161  * will be zero memset( hlm, 0, sizeof( COMPLEX16 ) ); return
162  * XLAL_SUCCESS; }
163  */
164 
165  //r = sqrt(values->data[0] * values->data[0] + values->data[1] * values->data[1] + values->data[2] * values->data[2]);
166  //pr = values->data[2];
167  r = values->data[0];
168  pp = values->data[3];
169 
170  rcrossp_x = cartvalues->data[1] * cartvalues->data[5] - cartvalues->data[2] * cartvalues->data[4];
171  rcrossp_y = cartvalues->data[2] * cartvalues->data[3] - cartvalues->data[0] * cartvalues->data[5];
172  rcrossp_z = cartvalues->data[0] * cartvalues->data[4] - cartvalues->data[1] * cartvalues->data[3];
173 
174  //pp = sqrt(rcrossp_x * rcrossp_x + rcrossp_y * rcrossp_y + rcrossp_z * rcrossp_z);
175 
176  v2 = v * v;
177  Omega = v2 * v;
178  //vh3 = Hreal * Omega;
179  //vh = cbrt(vh3);
180  eulerlogxabs = LAL_GAMMA + log(2.0 * (REAL8) m * v);
181 
182  /* Calculate the non-Keplerian velocity */
183  vPhi = vPhiKepler;
184 
185  vPhi = r * cbrt(vPhi);
186 
187  if (debugPK)
188  XLAL_PRINT_INFO("In XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole, getting rW = %.12e\n",
189  vPhi);
190  vPhi *= Omega;
191  vPhi2 = vPhi * vPhi;
192 
193  /*
194  * Calculate the newtonian multipole, 1st term in Eq. 17, given by
195  * Eq. A1
196  */
197  //debugPK
198  if (debugPK) {
199 // XLAL_PRINT_INFO("\nValues inside XLALSimIMRSpinEOBFluxGetPrecSpinFactorizedWaveform:\n");
200 // for (i = 0; i < 14; i++)
201 // XLAL_PRINT_INFO("values[%d] = %.12e\n", i, values->data[i]);
202 
203  XLAL_PRINT_INFO("Calculating hNewton, with v = %.12e, vPhi = %.12e, r = %.12e, Phi = %.12e, l = %d, m = %d\n",
204  v, vPhi, r, values->data[1], (UINT4) l, (UINT4) m);
205  }
207  vPhi2, r, values->data[1], (UINT4) l, m, params->eobParams);
208  if (status == XLAL_FAILURE) {
210  }
211  /* Calculate the source term, 2nd term in Eq. 17, given by Eq. A5 */
212  if (((l + m) % 2) == 0) {
213  Slm = (Hreal * Hreal - 1.) / (2. * eta) + 1.;
214  } else {
215  Slm = v * pp;
216  //Slm = v * sqrt(rcrossp_x * rcrossp_x + rcrossp_y * rcrossp_y + rcrossp_z * rcrossp_z);
217  }
218  if (debugPK)
219  XLAL_PRINT_INFO("In XLALSimIMRSpinEOBFluxGetSpinFactorizedWaveform: Hreal = %e, Slm = %e, eta = %e\n", Hreal, Slm, eta);
220 
221  /*
222  * Calculate the absolute value of the Tail term, 3rd term in Eq. 17,
223  * given by Eq. A6, and Eq. (42) of
224  * http://arxiv.org/pdf/1212.4357.pdf
225  */
226  k = m * Omega;
227  hathatk = Hreal * k;
228  hathatksq4 = 4. * hathatk * hathatk;
229  hathatk4pi = 4. * LAL_PI * hathatk;
230  /*
231  * gsl_sf_result lnr1, arg1; XLAL_CALLGSL( status =
232  * gsl_sf_lngamma_complex_e( l+1.0, -2.0*hathatk, &lnr1, &arg1 ) );
233  * if (status != GSL_SUCCESS) { XLALPrintError("XLAL Error - %s:
234  * Error in GSL function\n", __func__ ); XLAL_ERROR( XLAL_EFUNC ); }
235  */
236  XLAL_CALLGSL(status = gsl_sf_fact_e(l, &z2));
237  if (status != GSL_SUCCESS) {
238  XLALPrintError("XLAL Error - %s: Error in GSL function\n", __func__);
240  }
241  /*
242  * COMPLEX16 Tlmold; Tlmold = cexp( ( lnr1.val + LAL_PI * hathatk ) +
243  * I * ( arg1.val + 2.0 * hathatk * log(4.0*k/sqrt(LAL_E)) ) );
244  * Tlmold /= z2.val;
245  */
246  /* Calculating the prefactor of Tlm, outside the multiple product */
247  Tlmprefac = sqrt(hathatk4pi / (1. - exp(-hathatk4pi))) / z2.val;
248 
249  /* Calculating the multiple product factor */
250  for (Tlmprodfac = 1., i = 1; i <= l; i++) {
251  Tlmprodfac *= (hathatksq4 + (REAL8) i * i);
252  }
253 
254  Tlm = Tlmprefac * sqrt(Tlmprodfac);
255 
256  /* Calculate the residue phase and amplitude terms */
257  /*
258  * deltalm is the 4th term in Eq. 17, delta 22 given by Eq. A15,
259  * others
260  */
261  /* rholm is the 5th term in Eq. 17, given by Eqs. A8 - A14 */
262  /*
263  * auxflm is a special part of the 5th term in Eq. 17, given by Eq.
264  * A15
265  */
266  /*
267  * Actual values of the coefficients are defined in the next function
268  * of this file
269  */
270  switch (l) {
271  case 2:
272  switch (abs(m)) {
273  case 2:
274  rholm = 1. + v2 * (hCoeffs->rho22v2 + v * (hCoeffs->rho22v3
275  + v * (hCoeffs->rho22v4
276  + v * (hCoeffs->rho22v5 + v * (hCoeffs->rho22v6
277  + hCoeffs->rho22v6l * eulerlogxabs + v * (hCoeffs->rho22v7
278  + v * (hCoeffs->rho22v8 + hCoeffs->rho22v8l * eulerlogxabs
279  + (hCoeffs->rho22v10 + hCoeffs->rho22v10l * eulerlogxabs) * v2)))))));
280  //FIXME
281  if (debugPK){
282  XLAL_PRINT_INFO("PK:: rho22v2 = %.12e, rho22v3 = %.12e, rho22v4 = %.12e,\n rho22v5 = %.16e, rho22v6 = %.16e, rho22v6LOG = %.16e, \n rho22v7 = %.12e, rho22v8 = %.16e, rho22v8LOG = %.16e, \n rho22v10 = %.16e, rho22v10LOG = %.16e\n, rho22v6 = %.12e, rho22v8 = %.12e, rho22v10 = %.12e\n",
283  hCoeffs->rho22v2, hCoeffs->rho22v3, hCoeffs->rho22v4,
284  hCoeffs->rho22v5, hCoeffs->rho22v6, hCoeffs->rho22v6l,
285  hCoeffs->rho22v7, hCoeffs->rho22v8, hCoeffs->rho22v8l,
286  hCoeffs->rho22v10, hCoeffs->rho22v10l,
287  hCoeffs->rho22v6 + hCoeffs->rho22v6l * eulerlogxabs,
288  hCoeffs->rho22v8 + hCoeffs->rho22v8l * eulerlogxabs,
289  hCoeffs->rho22v10 + hCoeffs->rho22v10l * eulerlogxabs);}
290  break;
291  case 1:
292  {
293  rholm = 1. + v * (hCoeffs->rho21v1
294  + v * (hCoeffs->rho21v2 + v * (hCoeffs->rho21v3 + v * (hCoeffs->rho21v4
295  + v * (hCoeffs->rho21v5 + v * (hCoeffs->rho21v6 + hCoeffs->rho21v6l * eulerlogxabs
296  + v * (hCoeffs->rho21v7 + hCoeffs->rho21v7l * eulerlogxabs
297  + v * (hCoeffs->rho21v8 + hCoeffs->rho21v8l * eulerlogxabs
298  + (hCoeffs->rho21v10 + hCoeffs->rho21v10l * eulerlogxabs) * v2))))))));
299  auxflm = v * hCoeffs->f21v1 + v2 * v * hCoeffs->f21v3;
300  }
301  break;
302  default:
304  break;
305  }
306  break;
307  case 3:
308  switch (m) {
309  case 3:
310  rholm = 1. + v2 * (hCoeffs->rho33v2 + v * (hCoeffs->rho33v3 + v * (hCoeffs->rho33v4
311  + v * (hCoeffs->rho33v5 + v * (hCoeffs->rho33v6 + hCoeffs->rho33v6l * eulerlogxabs
312  + v * (hCoeffs->rho33v7 + (hCoeffs->rho33v8 + hCoeffs->rho33v8l * eulerlogxabs) * v))))));
313  auxflm = v * v2 * hCoeffs->f33v3;
314  break;
315  case 2:
316  rholm = 1. + v * (hCoeffs->rho32v
317  + v * (hCoeffs->rho32v2 + v * (hCoeffs->rho32v3 + v * (hCoeffs->rho32v4 + v * (hCoeffs->rho32v5
318  + v * (hCoeffs->rho32v6 + hCoeffs->rho32v6l * eulerlogxabs
319  + (hCoeffs->rho32v8 + hCoeffs->rho32v8l * eulerlogxabs) * v2))))));
320  break;
321  case 1:
322  rholm = 1. + v2 * (hCoeffs->rho31v2 + v * (hCoeffs->rho31v3 + v * (hCoeffs->rho31v4
323  + v * (hCoeffs->rho31v5 + v * (hCoeffs->rho31v6 + hCoeffs->rho31v6l * eulerlogxabs
324  + v * (hCoeffs->rho31v7 + (hCoeffs->rho31v8 + hCoeffs->rho31v8l * eulerlogxabs) * v))))));
325  auxflm = v * v2 * hCoeffs->f31v3;
326  break;
327  default:
329  break;
330  }
331  break;
332  case 4:
333  switch (m) {
334  case 4:
335 
336  rholm = 1. + v2 * (hCoeffs->rho44v2
337  + v * (hCoeffs->rho44v3 + v * (hCoeffs->rho44v4
338  + v * (hCoeffs->rho44v5 + (hCoeffs->rho44v6
339  + hCoeffs->rho44v6l * eulerlogxabs) * v))));
340  break;
341  case 3:
342  rholm = 1. + v * (hCoeffs->rho43v
343  + v * (hCoeffs->rho43v2
344  + v2 * (hCoeffs->rho43v4 + v * (hCoeffs->rho43v5
345  + (hCoeffs->rho43v6 + hCoeffs->rho43v6l * eulerlogxabs) * v))));
346  auxflm = v * hCoeffs->f43v;
347  break;
348  case 2:
349  rholm = 1. + v2 * (hCoeffs->rho42v2
350  + v * (hCoeffs->rho42v3 + v * (hCoeffs->rho42v4 + v * (hCoeffs->rho42v5
351  + (hCoeffs->rho42v6 + hCoeffs->rho42v6l * eulerlogxabs) * v))));
352  break;
353  case 1:
354  rholm = 1. + v * (hCoeffs->rho41v
355  + v * (hCoeffs->rho41v2
356  + v2 * (hCoeffs->rho41v4 + v * (hCoeffs->rho41v5
357  + (hCoeffs->rho41v6 + hCoeffs->rho41v6l * eulerlogxabs) * v))));
358  auxflm = v * hCoeffs->f41v;
359  break;
360  default:
362  break;
363  }
364  break;
365  case 5:
366  switch (m) {
367  case 5:
368  rholm = 1. + v2 * (hCoeffs->rho55v2
369  + v * (hCoeffs->rho55v3 + v * (hCoeffs->rho55v4
370  + v * (hCoeffs->rho55v5 + hCoeffs->rho55v6 * v))));
371  break;
372  case 4:
373  rholm = 1. + v2 * (hCoeffs->rho54v2 + v * (hCoeffs->rho54v3
374  + hCoeffs->rho54v4 * v));
375  break;
376  case 3:
377  rholm = 1. + v2 * (hCoeffs->rho53v2
378  + v * (hCoeffs->rho53v3 + v * (hCoeffs->rho53v4 + hCoeffs->rho53v5 * v)));
379  break;
380  case 2:
381  rholm = 1. + v2 * (hCoeffs->rho52v2 + v * (hCoeffs->rho52v3
382  + hCoeffs->rho52v4 * v));
383  break;
384  case 1:
385  rholm = 1. + v2 * (hCoeffs->rho51v2
386  + v * (hCoeffs->rho51v3 + v * (hCoeffs->rho51v4 + hCoeffs->rho51v5 * v)));
387  break;
388  default:
390  break;
391  }
392  break;
393  case 6:
394  switch (m) {
395  case 6:
396  rholm = 1. + v2 * (hCoeffs->rho66v2 + v * (hCoeffs->rho66v3
397  + hCoeffs->rho66v4 * v));
398  break;
399  case 5:
400  rholm = 1. + v2 * (hCoeffs->rho65v2 + hCoeffs->rho65v3 * v);
401  break;
402  case 4:
403  rholm = 1. + v2 * (hCoeffs->rho64v2 + v * (hCoeffs->rho64v3
404  + hCoeffs->rho64v4 * v));
405  break;
406  case 3:
407  rholm = 1. + v2 * (hCoeffs->rho63v2 + hCoeffs->rho63v3 * v);
408  break;
409  case 2:
410  rholm = 1. + v2 * (hCoeffs->rho62v2 + v * (hCoeffs->rho62v3
411  + hCoeffs->rho62v4 * v));
412  break;
413  case 1:
414  rholm = 1. + v2 * (hCoeffs->rho61v2 + hCoeffs->rho61v3 * v);
415  break;
416  default:
418  break;
419  }
420  break;
421  case 7:
422  switch (m) {
423  case 7:
424  rholm = 1. + v2 * (hCoeffs->rho77v2 + hCoeffs->rho77v3 * v);
425  break;
426  case 6:
427  rholm = 1. + hCoeffs->rho76v2 * v2;
428  break;
429  case 5:
430  rholm = 1. + v2 * (hCoeffs->rho75v2 + hCoeffs->rho75v3 * v);
431  break;
432  case 4:
433  rholm = 1. + hCoeffs->rho74v2 * v2;
434  break;
435  case 3:
436  rholm = 1. + v2 * (hCoeffs->rho73v2 + hCoeffs->rho73v3 * v);
437  break;
438  case 2:
439  rholm = 1. + hCoeffs->rho72v2 * v2;
440  break;
441  case 1:
442  rholm = 1. + v2 * (hCoeffs->rho71v2 + hCoeffs->rho71v3 * v);
443  break;
444  default:
446  break;
447  }
448  break;
449  case 8:
450  switch (m) {
451  case 8:
452  rholm = 1. + hCoeffs->rho88v2 * v2;
453  break;
454  case 7:
455  rholm = 1. + hCoeffs->rho87v2 * v2;
456  break;
457  case 6:
458  rholm = 1. + hCoeffs->rho86v2 * v2;
459  break;
460  case 5:
461  rholm = 1. + hCoeffs->rho85v2 * v2;
462  break;
463  case 4:
464  rholm = 1. + hCoeffs->rho84v2 * v2;
465  break;
466  case 3:
467  rholm = 1. + hCoeffs->rho83v2 * v2;
468  break;
469  case 2:
470  rholm = 1. + hCoeffs->rho82v2 * v2;
471  break;
472  case 1:
473  rholm = 1. + hCoeffs->rho81v2 * v2;
474  break;
475  default:
477  break;
478  }
479  break;
480  default:
482  break;
483  }
484 
485  //debugPK
486  if (debugPK)
487  XLAL_PRINT_INFO("rho_%d_%d = %.12e \n", l, m, rholm);
488  /* Raise rholm to the lth power */
489  rholmPwrl = 1.0;
490  i = l;
491  while (i--) {
492  rholmPwrl *= rholm;
493  }
494  /*
495  * In the equal-mass odd m case, there is no contribution from
496  * nonspin terms, and the only contribution comes from the auxflm
497  * term that is proportional to chiA (asymmetric spins). In this
498  * case, we must ignore the nonspin terms directly, since the leading
499  * term defined by CalculateThisMultipolePrefix in
500  * LALSimIMREOBNewtonianMultipole.c is not zero (see comments there).
501  */
502  if (eta == 0.25 && m % 2) {
503  rholmPwrl = auxflm;
504  } else {
505  rholmPwrl += auxflm;
506  }
507 
508  if (r > 0.0 && debugPK) {
509  XLAL_PRINT_INFO("YP::dynamics variables in waveform: %i, %i, r = %.12e, v = %.12e\n", l, m, r, v);
510  XLAL_PRINT_INFO("rholm^l = %.16e, Tlm = %.16e + i %.16e, \nSlm = %.16e, hNewton = %.16e + i %.16e, delta = %.16e\n", rholmPwrl, Tlm, 0.0, Slm, creal(hNewton), cimag(hNewton), 0.0);
511  }
512  /* Put all factors in Eq. 17 together */
513  *hlm = Tlm * Slm * rholmPwrl;
514  *hlm *= hNewton;
515  /*
516  * if (r > 8.5) { XLAL_PRINT_INFO("YP::FullWave: %.16e,%.16e,
517  * %.16e\n",hlm->re,hlm->im,sqrt(hlm->re*hlm->re+hlm->im*hlm->im)); }
518  */
519  }
520  }
521  return XLAL_SUCCESS;
522 }
523 
524 /**
525  * This function calculates hlm mode factorized-resummed waveform
526  * for given dynamical variables.
527  * Eq. 17 and the entire Appendix of PRD 86, 024011 (2012) + changes
528  * described in the section "Factorized waveforms" of https://dcc.ligo.org/T1400476
529  */
530 UNUSED
531 static INT4
533  COMPLEX16 * restrict hlm, /**< OUTPUT, hlm waveforms */
534  REAL8Vector * restrict values, /**< dyanmical variables: \f$(r,\phi,p_r,p_\phi)\f$ */
535  REAL8Vector * restrict cartvalues, /**< dyanmical variables */
536  const REAL8 v, /**< velocity */
537  const REAL8 Hreal, /**< real Hamiltonian */
538  const INT4 l, /**< l mode index */
539  const INT4 m, /**< m mode index */
540  SpinEOBParams * restrict params /**< Spin EOB parameters */
541 )
542 {
543  int debugPK = 0;
544  /* Status of function calls */
545  INT4 status;
546  INT4 i;
547 
548  REAL8 eta;
549  REAL8 UNUSED r , pp, Omega, v2, Omegav2, vh, vh3, k, hathatk, eulerlogxabs;
550  //pr
551  REAL8 UNUSED rcrossp_x, rcrossp_y, rcrossp_z;
552  REAL8 Slm , deltalm, rholm, rholmPwrl;
553  REAL8 auxflm = 0.0;
554  REAL8 hathatksq4, hathatk4pi, Tlmprefac, Tlmprodfac;
555  COMPLEX16 Tlm;
556  COMPLEX16 hNewton;
557  gsl_sf_result z2;
558 
559  /* Non-Keplerian velocity */
560  REAL8 vPhi , vPhi2;
561 
562  /* Pre-computed coefficients */
563 
564  FacWaveformCoeffs *hCoeffs = params->eobParams->hCoeffs;
565 
566  if (abs(m) > (INT4) l) {
568  }
569  if (m == 0) {
571  }
572  eta = params->eobParams->eta;
573 
574  /* Check our eta was sensible */
575  if (eta > 0.25) {
576  XLALPrintError("XLAL Error - %s: Eta seems to be > 0.25 - this isn't allowed!\n", __func__);
578  }
579  /*
580  * else if ( eta == 0.25 && m % 2 ) { // If m is odd and dM = 0, hLM
581  * will be zero memset( hlm, 0, sizeof( COMPLEX16 ) ); return
582  * XLAL_SUCCESS; }
583  */
584 
585  r = values->data[0];
586  pp = values->data[3];
587 
588  rcrossp_x = cartvalues->data[1] * cartvalues->data[5] - cartvalues->data[2] * cartvalues->data[4];
589  rcrossp_y = cartvalues->data[2] * cartvalues->data[3] - cartvalues->data[0] * cartvalues->data[5];
590  rcrossp_z = cartvalues->data[0] * cartvalues->data[4] - cartvalues->data[1] * cartvalues->data[3];
591 
592  //pp = sqrt(rcrossp_x * rcrossp_x + rcrossp_y * rcrossp_y + rcrossp_z * rcrossp_z);
593 
594  v2 = v * v;
595  Omega = v2 * v;
596  Omegav2 = Omega*v2;
597  vh3 = Hreal * Omega;
598  vh = cbrt(vh3);
599  eulerlogxabs = LAL_GAMMA + log(2.0 * (REAL8) m * v);
600 
601 
602  /* Calculate the non-Keplerian velocity */
603  if (params->alignedSpins) {
604 
605  vPhi = XLALSimIMRSpinAlignedEOBNonKeplerCoeff(values->data, params);
606 
607  if (XLAL_IS_REAL8_FAIL_NAN(vPhi)) {
609  }
610  vPhi = r * cbrt(vPhi);
611 
612  if (debugPK)
613  XLAL_PRINT_INFO("In XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole, getting rW = %.12e\n",
614  vPhi);
615  vPhi *= Omega;
616  vPhi2 = vPhi * vPhi;
617  } else {
618  vPhi = XLALSimIMRSpinPrecEOBNonKeplerCoeff_exact(cartvalues->data, params);
619 
620 
621  if (XLAL_IS_REAL8_FAIL_NAN(vPhi)) {
623  }
624  vPhi = r * cbrt(vPhi);
625 
626  if (debugPK)
627  XLAL_PRINT_INFO("In XLALSimIMRSpinEOBFluxCalculateNewtonianMultipole, getting rW = %.12e\n",
628  vPhi);
629  vPhi *= Omega;
630  vPhi2 = vPhi * vPhi;
631  }
632 
633  /*
634  * Calculate the newtonian multipole, 1st term in Eq. 17, given by
635  * Eq. A1
636  */
637  if (debugPK) {
638  XLAL_PRINT_INFO("\nValues inside XLALSimIMRSpinEOBFluxGetSpinFactorizedWaveform:\n");
639  for (i = 0; i < 11; i++)
640  XLAL_PRINT_INFO("values[%d] = %.12e\n", i, cartvalues->data[i]);
641 
642  XLAL_PRINT_INFO("Calculating hNewton, with v = %.12e, vPhi = %.12e, r = %.12e, Phi = %.12e, l = %d, m = %d\n",
643  v, vPhi, r, values->data[1], (UINT4) l, (UINT4) m);
644  }
646  vPhi2, r, cartvalues->data[12] + cartvalues->data[13], (UINT4) l, m, params->eobParams);
647  if (status == XLAL_FAILURE) {
649  }
650  /* Calculate the source term, 2nd term in Eq. 17, given by Eq. A5, Hreal is given by Eq.5 and Heff is in Eq.2 */
651  if (((l + m) % 2) == 0) {
652  Slm = (Hreal * Hreal - 1.) / (2. * eta) + 1.;
653  } else {
654  Slm = v * pp;
655  //Slm = v * sqrt(rcrossp_x * rcrossp_x + rcrossp_y * rcrossp_y + rcrossp_z * rcrossp_z);
656  }
657  if (debugPK)
658  XLAL_PRINT_INFO("In XLALSimIMRSpinEOBGetPrecSpinFactorizedWaveform: Hreal = %e, Slm = %e, eta = %e\n", Hreal, Slm, eta);
659 
660  /*
661  * Calculate the Tail term, 3rd term in Eq. 17,
662  * given by Eq. A6, and Eq. (42) of
663  * http://arxiv.org/pdf/1212.4357.pdf (or PRD 87 084035 (2013))
664  */
665  k = m * Omega;
666  hathatk = Hreal * k;
667 
668  gsl_sf_result lnr1, arg1;
669  XLAL_CALLGSL(status = gsl_sf_lngamma_complex_e(l + 1.0, -2.0 * hathatk, &lnr1, &arg1));
670  if (status != GSL_SUCCESS) {
671  XLALPrintError("XLAL Error - %s: Error in GSL function\n", __func__);
673  }
674  XLAL_CALLGSL(status = gsl_sf_fact_e(l, &z2));
675  if (status != GSL_SUCCESS) {
676  XLALPrintError("XLAL Error - %s: Error in GSL function\n", __func__);
678  }
679  Tlm = cexp((lnr1.val + LAL_PI * hathatk) + I * (
680  arg1.val + 2.0 * hathatk * log(4.0 * k / sqrt(LAL_E))));
681  Tlm /= z2.val;
682 
683 
684 
685  if (debugPK){
686  hathatksq4 = 4. * hathatk * hathatk;
687  hathatk4pi = 4. * LAL_PI * hathatk;
688  /* Calculating the prefactor of Tlm, outside the multiple product */
689  Tlmprefac = sqrt(hathatk4pi / (1. - exp(-hathatk4pi))) / z2.val;
690 
691  /* Calculating the multiple product factor */
692  for (Tlmprodfac = 1., i = 1; i <= l; i++)
693  Tlmprodfac *= (hathatksq4 + (REAL8) i * i);
694 
695  REAL8 Tlmold;
696  Tlmold = Tlmprefac * sqrt(Tlmprodfac);
697  XLAL_PRINT_INFO("Tlm = %e + i%e, |Tlm| = %.16e (should be %.16e)\n", creal(Tlm), cimag(Tlm), cabs(Tlm), Tlmold);
698  }
699 
700  /* Calculate the residue phase and amplitude terms */
701  /*
702  * deltalm is the 4th term in Eq. 17, delta 22 given by Eq. A15,
703  * others
704  */
705  /* rholm is the 5th term in Eq. 17, given by Eqs. A8 - A14 */
706  /*
707  * auxflm is a special part of the 5th term in Eq. 17, given by Eq.
708  * A15
709  */
710  /*
711  * Actual values of the coefficients are defined in the another function
712  * see file LALSimIMRSpinEOBFactorizedWaveformCoefficientsPrec.c
713  */
714  switch (l) {
715  case 2:
716  switch (abs(m)) {
717  case 2:
718  deltalm = vh3 * (hCoeffs->delta22vh3 + vh3 * (hCoeffs->delta22vh6
719  + vh * vh * (hCoeffs->delta22vh9 * vh)))
720  + Omega*(hCoeffs->delta22v5 * v2 + Omega*(hCoeffs->delta22v6 + hCoeffs->delta22v8 *v2));
721  rholm = 1. + v2 * (hCoeffs->rho22v2 + v * (hCoeffs->rho22v3
722  + v * (hCoeffs->rho22v4
723  + v * (hCoeffs->rho22v5 + v * (hCoeffs->rho22v6
724  + hCoeffs->rho22v6l * eulerlogxabs + v * (hCoeffs->rho22v7
725  + v * (hCoeffs->rho22v8 + hCoeffs->rho22v8l * eulerlogxabs
726  + (hCoeffs->rho22v10 + hCoeffs->rho22v10l * eulerlogxabs) * v2)))))));
727  //FIXME
728  if (debugPK){
729  XLAL_PRINT_INFO("PK:: rho22v2 = %.12e, rho22v3 = %.12e, rho22v4 = %.12e,\n rho22v5 = %.16e, rho22v6 = %.16e, rho22v6LOG = %.16e, \n rho22v7 = %.12e, rho22v8 = %.16e, rho22v8LOG = %.16e, \n rho22v10 = %.16e, rho22v10LOG = %.16e\n, rho22v6 = %.12e, rho22v8 = %.12e, rho22v10 = %.12e\n",
730  hCoeffs->rho22v2, hCoeffs->rho22v3, hCoeffs->rho22v4,
731  hCoeffs->rho22v5, hCoeffs->rho22v6, hCoeffs->rho22v6l,
732  hCoeffs->rho22v7, hCoeffs->rho22v8, hCoeffs->rho22v8l,
733  hCoeffs->rho22v10, hCoeffs->rho22v10l,
734  hCoeffs->rho22v6 + hCoeffs->rho22v6l * eulerlogxabs,
735  hCoeffs->rho22v8 + hCoeffs->rho22v8l * eulerlogxabs,
736  hCoeffs->rho22v10 + hCoeffs->rho22v10l * eulerlogxabs);}
737  break;
738  case 1:
739  {
740  deltalm = vh3 * (hCoeffs->delta21vh3 + vh3 * (hCoeffs->delta21vh6
741  + vh * (hCoeffs->delta21vh7 + (hCoeffs->delta21vh9) * vh * vh)))
742  + Omegav2*(hCoeffs->delta21v5 + hCoeffs->delta21v7 * v2);
743  rholm = 1. + v * (hCoeffs->rho21v1
744  + v * (hCoeffs->rho21v2 + v * (hCoeffs->rho21v3 + v * (hCoeffs->rho21v4
745  + v * (hCoeffs->rho21v5 + v * (hCoeffs->rho21v6 + hCoeffs->rho21v6l * eulerlogxabs
746  + v * (hCoeffs->rho21v7 + hCoeffs->rho21v7l * eulerlogxabs
747  + v * (hCoeffs->rho21v8 + hCoeffs->rho21v8l * eulerlogxabs
748  + (hCoeffs->rho21v10 + hCoeffs->rho21v10l * eulerlogxabs) * v2))))))));
749  auxflm = v * hCoeffs->f21v1 + v2 * v * hCoeffs->f21v3;
750  }
751  break;
752  default:
754  break;
755  }
756  break;
757  case 3:
758  switch (m) {
759  case 3:
760  deltalm = vh3 * (hCoeffs->delta33vh3 + vh3 * (hCoeffs->delta33vh6 + hCoeffs->delta33vh9 * vh3))
761  + hCoeffs->delta33v5 * v * v2 * v2 + hCoeffs->delta33v7 * v2 * v2 * v2 * v;
762  rholm = 1. + v2 * (hCoeffs->rho33v2 + v * (hCoeffs->rho33v3 + v * (hCoeffs->rho33v4
763  + v * (hCoeffs->rho33v5 + v * (hCoeffs->rho33v6 + hCoeffs->rho33v6l * eulerlogxabs
764  + v * (hCoeffs->rho33v7 + (hCoeffs->rho33v8 + hCoeffs->rho33v8l * eulerlogxabs) * v))))));
765  auxflm = v * v2 * hCoeffs->f33v3;
766  break;
767  case 2:
768  deltalm = vh3 * (hCoeffs->delta32vh3 + vh * (hCoeffs->delta32vh4 + vh * vh * (hCoeffs->delta32vh6
769  + hCoeffs->delta32vh9 * vh3)));
770  rholm = 1. + v * (hCoeffs->rho32v
771  + v * (hCoeffs->rho32v2 + v * (hCoeffs->rho32v3 + v * (hCoeffs->rho32v4 + v * (hCoeffs->rho32v5
772  + v * (hCoeffs->rho32v6 + hCoeffs->rho32v6l * eulerlogxabs
773  + (hCoeffs->rho32v8 + hCoeffs->rho32v8l * eulerlogxabs) * v2))))));
774  break;
775  case 1:
776  deltalm = vh3 * (hCoeffs->delta31vh3 + vh3 * (hCoeffs->delta31vh6
777  + vh * (hCoeffs->delta31vh7 + hCoeffs->delta31vh9 * vh * vh)))
778  + hCoeffs->delta31v5 * v * v2 * v2;
779  rholm = 1. + v2 * (hCoeffs->rho31v2 + v * (hCoeffs->rho31v3 + v * (hCoeffs->rho31v4
780  + v * (hCoeffs->rho31v5 + v * (hCoeffs->rho31v6 + hCoeffs->rho31v6l * eulerlogxabs
781  + v * (hCoeffs->rho31v7 + (hCoeffs->rho31v8 + hCoeffs->rho31v8l * eulerlogxabs) * v))))));
782  auxflm = v * v2 * hCoeffs->f31v3;
783  break;
784  default:
786  break;
787  }
788  break;
789  case 4:
790  switch (m) {
791  case 4:
792 
793  deltalm = vh3 * (hCoeffs->delta44vh3 + hCoeffs->delta44vh6 * vh3)
794  + hCoeffs->delta44v5 * v2 * v2 * v;
795  rholm = 1. + v2 * (hCoeffs->rho44v2
796  + v * (hCoeffs->rho44v3 + v * (hCoeffs->rho44v4
797  + v * (hCoeffs->rho44v5 + (hCoeffs->rho44v6
798  + hCoeffs->rho44v6l * eulerlogxabs) * v))));
799  break;
800  case 3:
801  deltalm = vh3 * (hCoeffs->delta43vh3 + vh * (hCoeffs->delta43vh4
802  + hCoeffs->delta43vh6 * vh * vh));
803  rholm = 1. + v * (hCoeffs->rho43v
804  + v * (hCoeffs->rho43v2
805  + v2 * (hCoeffs->rho43v4 + v * (hCoeffs->rho43v5
806  + (hCoeffs->rho43v6 + hCoeffs->rho43v6l * eulerlogxabs) * v))));
807  auxflm = v * hCoeffs->f43v;
808  break;
809  case 2:
810  deltalm = vh3 * (hCoeffs->delta42vh3 + hCoeffs->delta42vh6 * vh3);
811  rholm = 1. + v2 * (hCoeffs->rho42v2
812  + v * (hCoeffs->rho42v3 + v * (hCoeffs->rho42v4 + v * (hCoeffs->rho42v5
813  + (hCoeffs->rho42v6 + hCoeffs->rho42v6l * eulerlogxabs) * v))));
814  break;
815  case 1:
816  deltalm = vh3 * (hCoeffs->delta41vh3 + vh * (hCoeffs->delta41vh4
817  + hCoeffs->delta41vh6 * vh * vh));
818  rholm = 1. + v * (hCoeffs->rho41v
819  + v * (hCoeffs->rho41v2
820  + v2 * (hCoeffs->rho41v4 + v * (hCoeffs->rho41v5
821  + (hCoeffs->rho41v6 + hCoeffs->rho41v6l * eulerlogxabs) * v))));
822  auxflm = v * hCoeffs->f41v;
823  break;
824  default:
826  break;
827  }
828  break;
829  case 5:
830  switch (m) {
831  case 5:
832  deltalm = hCoeffs->delta55vh3 * vh3 + hCoeffs->delta55v5 * v2 * v2 * v;
833  rholm = 1. + v2 * (hCoeffs->rho55v2
834  + v * (hCoeffs->rho55v3 + v * (hCoeffs->rho55v4
835  + v * (hCoeffs->rho55v5 + hCoeffs->rho55v6 * v))));
836  break;
837  case 4:
838  deltalm = vh3 * (hCoeffs->delta54vh3 + hCoeffs->delta54vh4 * vh);
839  rholm = 1. + v2 * (hCoeffs->rho54v2 + v * (hCoeffs->rho54v3
840  + hCoeffs->rho54v4 * v));
841  break;
842  case 3:
843  deltalm = hCoeffs->delta53vh3 * vh3;
844  rholm = 1. + v2 * (hCoeffs->rho53v2
845  + v * (hCoeffs->rho53v3 + v * (hCoeffs->rho53v4 + hCoeffs->rho53v5 * v)));
846  break;
847  case 2:
848  deltalm = vh3 * (hCoeffs->delta52vh3 + hCoeffs->delta52vh4 * vh);
849  rholm = 1. + v2 * (hCoeffs->rho52v2 + v * (hCoeffs->rho52v3
850  + hCoeffs->rho52v4 * v));
851  break;
852  case 1:
853  deltalm = hCoeffs->delta51vh3 * vh3;
854  rholm = 1. + v2 * (hCoeffs->rho51v2
855  + v * (hCoeffs->rho51v3 + v * (hCoeffs->rho51v4 + hCoeffs->rho51v5 * v)));
856  break;
857  default:
859  break;
860  }
861  break;
862  case 6:
863  switch (m) {
864  case 6:
865  deltalm = hCoeffs->delta66vh3 * vh3;
866  rholm = 1. + v2 * (hCoeffs->rho66v2 + v * (hCoeffs->rho66v3
867  + hCoeffs->rho66v4 * v));
868  break;
869  case 5:
870  deltalm = hCoeffs->delta65vh3 * vh3;
871  rholm = 1. + v2 * (hCoeffs->rho65v2 + hCoeffs->rho65v3 * v);
872  break;
873  case 4:
874  deltalm = hCoeffs->delta64vh3 * vh3;
875  rholm = 1. + v2 * (hCoeffs->rho64v2 + v * (hCoeffs->rho64v3
876  + hCoeffs->rho64v4 * v));
877  break;
878  case 3:
879  deltalm = hCoeffs->delta63vh3 * vh3;
880  rholm = 1. + v2 * (hCoeffs->rho63v2 + hCoeffs->rho63v3 * v);
881  break;
882  case 2:
883  deltalm = hCoeffs->delta62vh3 * vh3;
884  rholm = 1. + v2 * (hCoeffs->rho62v2 + v * (hCoeffs->rho62v3
885  + hCoeffs->rho62v4 * v));
886  break;
887  case 1:
888  deltalm = hCoeffs->delta61vh3 * vh3;
889  rholm = 1. + v2 * (hCoeffs->rho61v2 + hCoeffs->rho61v3 * v);
890  break;
891  default:
893  break;
894  }
895  break;
896  case 7:
897  switch (m) {
898  case 7:
899  deltalm = hCoeffs->delta77vh3 * vh3;
900  rholm = 1. + v2 * (hCoeffs->rho77v2 + hCoeffs->rho77v3 * v);
901  break;
902  case 6:
903  deltalm = 0.0;
904  rholm = 1. + hCoeffs->rho76v2 * v2;
905  break;
906  case 5:
907  deltalm = hCoeffs->delta75vh3 * vh3;
908  rholm = 1. + v2 * (hCoeffs->rho75v2 + hCoeffs->rho75v3 * v);
909  break;
910  case 4:
911  deltalm = 0.0;
912  rholm = 1. + hCoeffs->rho74v2 * v2;
913  break;
914  case 3:
915  deltalm = hCoeffs->delta73vh3 * vh3;
916  rholm = 1. + v2 * (hCoeffs->rho73v2 + hCoeffs->rho73v3 * v);
917  break;
918  case 2:
919  deltalm = 0.0;
920  rholm = 1. + hCoeffs->rho72v2 * v2;
921  break;
922  case 1:
923  deltalm = hCoeffs->delta71vh3 * vh3;
924  rholm = 1. + v2 * (hCoeffs->rho71v2 + hCoeffs->rho71v3 * v);
925  break;
926  default:
928  break;
929  }
930  break;
931  case 8:
932  switch (m) {
933  case 8:
934  deltalm = 0.0;
935  rholm = 1. + hCoeffs->rho88v2 * v2;
936  break;
937  case 7:
938  deltalm = 0.0;
939  rholm = 1. + hCoeffs->rho87v2 * v2;
940  break;
941  case 6:
942  deltalm = 0.0;
943  rholm = 1. + hCoeffs->rho86v2 * v2;
944  break;
945  case 5:
946  deltalm = 0.0;
947  rholm = 1. + hCoeffs->rho85v2 * v2;
948  break;
949  case 4:
950  deltalm = 0.0;
951  rholm = 1. + hCoeffs->rho84v2 * v2;
952  break;
953  case 3:
954  deltalm = 0.0;
955  rholm = 1. + hCoeffs->rho83v2 * v2;
956  break;
957  case 2:
958  deltalm = 0.0;
959  rholm = 1. + hCoeffs->rho82v2 * v2;
960  break;
961  case 1:
962  deltalm = 0.0;
963  rholm = 1. + hCoeffs->rho81v2 * v2;
964  break;
965  default:
967  break;
968  }
969  break;
970  default:
972  break;
973  }
974 
975  //debugPK
976  if (debugPK)
977  XLAL_PRINT_INFO("rho_%d_%d = %.12e \n", l, m, rholm);
978  if (debugPK)
979  XLAL_PRINT_INFO("exp(delta_%d_%d) = %.16e + i%.16e (abs=%e)\n", l, m, creal(cexp(I * deltalm)),
980  cimag(cexp(I * deltalm)), cabs(cexp(I * deltalm)));
981  /* Raise rholm to the lth power */
982  rholmPwrl = 1.0;
983  i = l;
984  while (i--) {
985  rholmPwrl *= rholm;
986  }
987  /*
988  * In the equal-mass odd m case, there is no contribution from
989  * nonspin terms, and the only contribution comes from the auxflm
990  * term that is proportional to chiA (asymmetric spins). In this
991  * case, we must ignore the nonspin terms directly, since the leading
992  * term defined by CalculateThisMultipolePrefix in
993  * LALSimIMREOBNewtonianMultipole.c is not zero (see comments there).
994  */
995  if (eta == 0.25 && m % 2) {
996  rholmPwrl = auxflm;
997  } else {
998  rholmPwrl += auxflm;
999  }
1000 
1001  if (r > 0.0 && debugPK) {
1002  XLAL_PRINT_INFO("YP::dynamics variables in waveform: %i, %i, r = %.12e, v = %.12e\n", l, m, r, v);
1003  XLAL_PRINT_INFO("rholm^l = %.16e, Tlm = %.16e + i %.16e, \nSlm = %.16e, hNewton = %.16e + i %.16e, delta = %.16e\n", rholmPwrl, creal(Tlm), cimag(Tlm), Slm, creal(hNewton), cimag(hNewton), 0.0);
1004  }
1005  /* Put all factors in Eq. 17 together */
1006  *hlm = Tlm * cexp(I * deltalm) * Slm * rholmPwrl;
1007  *hlm *= hNewton;
1008  if (r > 8.5 && debugPK) {
1009  XLAL_PRINT_INFO("YP::FullWave: %.16e,%.16e, %.16e\n", creal(*hlm), cimag(*hlm), cabs(*hlm));
1010  }
1011  return XLAL_SUCCESS;
1012 }
1013 
1014 
1015 /*
1016  * Function to calculate the value of omega for the PRECESSING EOB waveform.
1017  * Needs the dynamics in Cartesian coordinates.
1018  *
1019  * First, the frame is rotated so that L is along the y-axis.
1020  * this rotation includes the spins.
1021  *
1022  * Second, \f$\vec{r}\f$ and \f$\vec{p}\f$ are converted to polar coordinates
1023  * (and not the spins). As L is along the y-axis, \f$\theta\f$ (defined as the
1024  * angle between L and the y-axis) is 0, which is a cyclic coordinate now and
1025  * that fixes
1026  * \f$p_\theta = 0\f$.
1027  *
1028  * Third, \f$p_r\f$ is set to 0.
1029  *
1030  * Fourth, the polar \f$(r,\phi,p_r=0,p_\phi)\f$ and the Cartesian spin vectors
1031  * are used to compute the derivative
1032  * \f$\partial Hreal/\partial p_\phi |p_r=0\f$.
1033  */
1034 
1036  const REAL8 values[], /**<< Dynamical variables */
1037  SpinEOBParams *funcParams /**<< EOB parameters */
1038  )
1039 {
1040  int debugPK = 1;
1041  if (debugPK){
1042  for(int i =0; i < 12; i++)
1043  if( isnan(values[i]) ) {
1044  XLAL_PRINT_INFO("XLALSimIMRSpinPrecEOBCalcOmega_exact::values %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f\n", values[0], values[1], values[2], values[3], values[4], values[5], values[6], values[7], values[8], values[9], values[10], values[11]);
1045  XLALPrintError( "XLAL Error - %s: nan in input values \n", __func__);
1047  }
1048  }
1049 
1050  /* ********************************************************************* */
1051  /* ************ Memory Allocation ************************************** */
1052  /* ********************************************************************* */
1053  static const REAL8 STEP_SIZE = 1.0e-4;
1054  REAL8 tmpvar = 0;
1055 
1057 
1058  /* Cartesian values for calculating the Hamiltonian */
1059  REAL8 cartValues[14] = {0.}, dvalues[14] = {0.};
1060  REAL8 cartvalues[14] = {0.}, polarvalues[6] = {0.}; /* The rotated cartesian/polar values */
1061  REAL8 polarRPcartSvalues[14] = {0.};
1062  memcpy( cartValues, values, 14 * sizeof(REAL8) );
1063 
1064  INT4 i, j;
1065 
1066  REAL8 rvec[3] = {0.,0,0}, pvec[3] = {0.,0,0};
1067  REAL8 s1vec[3] = {0.,0,0}, s2vec[3] = {0.,0,0};
1068 
1069  REAL8 rdotvec[3] = {0.,0,0};
1070  REAL8 rvecprime[3] = {0.,0,0}, pvecprime[3] = {0.,0,0},
1071  s1vecprime[3]= {0.,0,0}, s2vecprime[3]= {0.,0,0};
1072  REAL8 rvectmp[3] = {0.,0,0}, pvectmp[3] = {0.,0,0},
1073  s1vectmp[3] = {0.,0,0}, s2vectmp[3]= {0.,0,0};
1074  REAL8 LNhatprime[3]= {0.,0,0}, LNhatTmp[3]= {0.,0,0};
1075  REAL8 rcrossrdot[3] = {0.,0,0};
1076 
1077  REAL8 Rot1[3][3] ={{0.}}; // Rotation matrix for prevention of blowing up
1078  REAL8 Rot2[3][3] ={{0.}} ;
1079  REAL8 LNhat[3] = {0.,0,0};
1080 
1081  REAL8 Xhat[3] = {1, 0, 0};
1082  UNUSED REAL8 Yhat[3] = {0, 1, 0};
1083  UNUSED REAL8 Zhat[3] = {0, 0, 1};
1084 
1085  REAL8 Xprime[3] = {0.,0,0}, Yprime[3] = {0.,0,0}, Zprime[3] = {0.,0,0};
1086 
1087  gsl_function F;
1088  INT4 gslStatus;
1089 
1090  REAL8 omega;
1091 
1092  /* The error in a derivative as measured by GSL */
1093  REAL8 absErr;
1094 
1095  /* ********************************************************************* */
1096  /* ************ Main Logic begins ************************************ */
1097  /* ********************************************************************* */
1098 
1099  /* Copy over the coordinates and spins */
1100  memcpy( rvec, values, 3*sizeof(REAL8) );
1101  memcpy( pvec, values+3, 3*sizeof(REAL8) );
1102  memcpy( s1vec, values+6, 3*sizeof(REAL8) );
1103  memcpy( s2vec, values+9, 3*sizeof(REAL8) );
1104 
1105  /* Calculate rDot = \f$\partial Hreal / \partial p_r\f$ */
1106  memset( dvalues, 0, 14 * sizeof(REAL8) );
1107  if( XLALSpinPrecHcapRvecDerivative_exact( 0, values, dvalues,
1108  (void*) funcParams) == XLAL_FAILURE )
1109  {
1111  }
1112  memcpy( rdotvec, dvalues, 3*sizeof(REAL8) );
1113 
1114  if (debugPK){
1115  for(int ii =0; ii < 12; ii++)
1116  if( isnan(dvalues[ii]) ) {
1117  XLAL_PRINT_INFO("XLALSimIMRSpinPrecEOBCalcOmega::dvalues %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f\n", dvalues[0], dvalues[1], dvalues[2], dvalues[3], dvalues[4], dvalues[5], dvalues[6], dvalues[7], dvalues[8], dvalues[9], dvalues[10], dvalues[11]);
1118  XLALPrintError( "XLAL Error - %s: nan in dvalues \n", __func__);
1120  }
1121  }
1122 
1123  /* Calculate LN by computing r cross rDot and normalizing */
1124  cross_product( rvec, rdotvec, rcrossrdot );
1125  REAL8 rcrossrdotNorm = sqrt(inner_product( rcrossrdot, rcrossrdot ));
1126  for( i = 0; i < 3; i++ )
1127  rcrossrdot[i] /= rcrossrdotNorm;
1128  memcpy( LNhat, rcrossrdot, 3 * sizeof(REAL8) );
1129 
1130 
1131  /* ********************************************************************* */
1132  /* First, the frame is rotated so that L is along the y-axis. */
1133  /* this rotation includes the spins. */
1134  /* ********************************************************************* */
1135 
1136  /* For now, set first rotation matrix to identity. If LNhat and Xhat are
1137  too closely aligned, rotate LNhat by pi/4 in the x,y plane. Then
1138  LNhat determines Xprime, Yprime is the normalized cross product of
1139  Xprime and Xhat, and Zprime is the normalized cross product of Xprime
1140  and Yprime.
1141  */
1142  if( inner_product(LNhat, Xhat) < 0.9 )
1143  {
1144  Rot1[0][0] = 1.; Rot1[0][1] = 0; Rot1[0][2] = 0;
1145  Rot1[1][0] = 0.; Rot1[1][1] = 1; Rot1[1][2] = 0;
1146  Rot1[2][0] = 0.; Rot1[2][1] = 0; Rot1[2][2] = 1;
1147 
1148  memcpy(Xprime, LNhat, 3 * sizeof(REAL8));
1149  }
1150  else
1151  {
1152  Rot1[0][0] = 1./sqrt(2); Rot1[0][1] = -1/sqrt(2); Rot1[0][2] = 0;
1153  Rot1[1][0] = 1./sqrt(2); Rot1[1][1] = 1./sqrt(2); Rot1[1][2] = 0;
1154  Rot1[2][0] = 0.; Rot1[2][1] = 0; Rot1[2][2] = 1;
1155  LNhatTmp[0] = LNhatTmp[1] = LNhatTmp[2] = 0.;
1156 
1157  for(i=0; i<3; i++)
1158  for(j=0; j<3; j++)
1159  LNhatTmp[i] += Rot1[i][j]*LNhat[j];
1160 
1161  memcpy(Xprime, LNhatTmp, 3*sizeof(REAL8));
1162  }
1163 
1164  cross_product( Xprime, Xhat, Yprime );
1165  tmpvar = sqrt(inner_product(Yprime, Yprime));
1166 
1167  for( i=0; i<3; i++) Yprime[i] /= tmpvar;
1168 
1169  cross_product(Xprime, Yprime, Zprime);
1170  tmpvar = sqrt(inner_product(Zprime, Zprime));
1171  for( i=0; i<3; i++) Zprime[i] /= tmpvar;
1172 
1173  Rot2[0][0] = Xprime[0]; Rot2[0][1] = Xprime[1]; Rot2[0][2] = Xprime[2];
1174  Rot2[1][0] = Yprime[0]; Rot2[1][1] = Yprime[1]; Rot2[1][2] = Yprime[2];
1175  Rot2[2][0] = Zprime[0]; Rot2[2][1] = Zprime[1]; Rot2[2][2] = Zprime[2];
1176 
1177  memset( rvectmp, 0, 3 * sizeof(REAL8) );
1178  memset( pvectmp, 0, 3 * sizeof(REAL8) );
1179  memset( s1vectmp, 0, 3 * sizeof(REAL8) );
1180  memset( s2vectmp, 0, 3 * sizeof(REAL8) );
1181  memset( rvecprime, 0, 3 * sizeof(REAL8) );
1182  memset( pvecprime, 0, 3 * sizeof(REAL8) );
1183  memset( s1vecprime, 0, 3 * sizeof(REAL8) );
1184  memset( s2vecprime, 0, 3 * sizeof(REAL8) );
1185  memset( LNhatprime, 0, 3 * sizeof(REAL8) );
1186  memset( LNhatTmp, 0, 3 * sizeof(REAL8) );
1187 
1188  /* Perform the actual rotation */
1189  for (i=0; i<3; i++)
1190  for(j=0; j<3; j++)
1191  {
1192  rvectmp[i] += Rot1[i][j]*rvec[j];
1193  pvectmp[i] += Rot1[i][j]*pvec[j];
1194  s1vectmp[i] += Rot1[i][j]*s1vec[j];
1195  s2vectmp[i] += Rot1[i][j]*s2vec[j];
1196  LNhatTmp[i] += Rot1[i][j]*LNhat[j];
1197  }
1198  for (i=0; i<3; i++)
1199  for(j=0; j<3; j++)
1200  {
1201  rvecprime[i] += Rot2[i][j]*rvectmp[j];
1202  pvecprime[i] += Rot2[i][j]*pvectmp[j];
1203  s1vecprime[i] += Rot2[i][j]*s1vectmp[j];
1204  s2vecprime[i] += Rot2[i][j]*s2vectmp[j];
1205  LNhatprime[i] += Rot2[i][j]*LNhatTmp[j];
1206  }
1207 
1208  memcpy(cartvalues, rvecprime, 3*sizeof(REAL8));
1209  memcpy(cartvalues+3, pvecprime, 3*sizeof(REAL8));
1210  memcpy(cartvalues+6, s1vecprime, 3*sizeof(REAL8));
1211  memcpy(cartvalues+9, s2vecprime, 3*sizeof(REAL8));
1212 
1213  /* ********************************************************************* */
1214  /* Second, \f$\vec{r}\f$ and \f$\vec{p}\f$ are converted to polar
1215  * coordinates (and not the spins).
1216  * As L is along the y-axis, \f$\theta\f$ (defined as the angle between
1217  * L and the y-axis) is 0, which is a cyclic coordinate now and that fixes
1218  * \f$p_\theta = 0\f$. */
1219  /* ********************************************************************* */
1220 
1221  /** The polarvalues, respectively, are
1222  * \f${r, \theta, \phi, p_r, p_\theta, p_\phi}\f$.
1223  * Below we transform \f$\vec{r}\f$ and \f$\vec{p}\f$ to the six
1224  * corresponding polar coordinates. */
1225  polarvalues[0] = sqrt(inner_product(rvecprime,rvecprime));
1226  polarvalues[1] = acos(rvecprime[0] / polarvalues[0]);
1227  polarvalues[2] = atan2(-rvecprime[1], rvecprime[2]);
1228 
1229  /* FIX p_r = 0 */
1230  polarvalues[3] = 0;
1231 
1232  REAL8 rvecprime_x_xhat[3] = {0.}, rvecprime_x_xhat_x_rvecprime[3] = {0.};
1233  cross_product(rvecprime, Xhat, rvecprime_x_xhat);
1234  cross_product(rvecprime_x_xhat, rvecprime, rvecprime_x_xhat_x_rvecprime);
1235 
1236  polarvalues[4] = -inner_product(rvecprime_x_xhat_x_rvecprime, pvecprime)
1237  / polarvalues[0] / sin(polarvalues[1]);
1238  polarvalues[5] = -inner_product(rvecprime_x_xhat, pvecprime);
1239 
1240 
1241  /* ********************************************************************* */
1242  /* Finally, Differentiate Hamiltonian w.r.t. p_\phi, keeping p_r = 0 */
1243  /* ********************************************************************* */
1244 
1245  /* Populate the vector specifying the dynamical variables in mixed frames */
1246  memcpy( polarRPcartSvalues, cartvalues, 12*sizeof(REAL8));
1247  memcpy( polarRPcartSvalues, polarvalues, 6*sizeof(REAL8));
1248 
1249  /* Set up pointers for GSL */
1250  params.values = polarRPcartSvalues;
1251  params.params = funcParams;
1252 
1254  F.params = &params;
1255 
1256  /* Now calculate omega. In the chosen co-ordinate system, */
1257  /* we need dH/dpphi to calculate this, i.e. varyParam = 5 */
1258  params.varyParam = 5;
1259  XLAL_CALLGSL( gslStatus = gsl_deriv_central( &F, polarvalues[5],
1260  STEP_SIZE, &omega, &absErr ) );
1261 
1262  if ( gslStatus != GSL_SUCCESS )
1263  {
1264  XLALPrintError( "XLAL Error - %s: Failure in GSL function\n", __func__ );
1266  }
1267 
1268  return omega;
1269 }
1270 
1271 /**
1272  * Function to calculate the non-Keplerian coefficient for the PRECESSING EOB
1273  * model.
1274  *
1275  * radius \f$r\f$ times the cuberoot of the returned number is \f$r_\Omega\f$
1276  * defined in Eq. A2, i.e. the function returns
1277  * \f$(r_{\Omega} / r)^3\f$
1278  * = \f$1/(r^3 (\partial Hreal/\partial p_\phi |p_r=0)^2)\f$.
1279  */
1280 static REAL8
1282  const REAL8 values[], /**<< Dynamical variables */
1283  SpinEOBParams *funcParams /**<< EOB parameters */
1284  )
1285 {
1286  int debugPK = 1;
1287  if (debugPK){
1288  for(int i =0; i < 12; i++)
1289  if( isnan(values[i]) ) {
1290  XLAL_PRINT_INFO("XLALSimIMRSpinPrecEOBNonKeplerCoeff_exact::values %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f %3.10f\n",
1291  values[0], values[1], values[2], values[3], values[4], values[5],
1292  values[6], values[7], values[8], values[9], values[10], values[11]);
1293  XLALPrintError( "XLAL Error - %s: nan in values \n", __func__);
1295  }
1296  }
1297 
1298  REAL8 omegaCirc = 0;
1299  REAL8 tmpValues[14]= {0.};
1300  REAL8 r3;
1301 
1302  /* We need to find the values of omega assuming pr = 0 */
1303  memcpy( tmpValues, values, sizeof(tmpValues) );
1304  omegaCirc = XLALSimIMRSpinPrecEOBCalcOmega_exact( tmpValues, funcParams );
1305 
1306  if ( XLAL_IS_REAL8_FAIL_NAN( omegaCirc ) )
1307  {
1309  }
1310 
1311  /* inner_product operates on the first three entries of the vector
1312  * passed to it, and therefore inner_product(values,values) = r^2 */
1313  r3 = pow(inner_product(values, values), 3./2.);
1314  return 1.0/(omegaCirc*omegaCirc*r3);
1315 }
1316 
1317 #endif /* _LALSIMIMRSPINPRECEOBFACTORIZEDWAVEFORM_V3OPT */
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...
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...
static REAL8 XLALSimIMRSpinPrecEOBCalcOmega_exact(const REAL8 values[], SpinEOBParams *funcParams)
static REAL8 XLALSimIMRSpinPrecEOBNonKeplerCoeff_exact(const REAL8 values[], SpinEOBParams *funcParams)
Function to calculate the non-Keplerian coefficient for the PRECESSING EOB model.
static INT4 XLALSimIMRSpinEOBFluxGetPrecSpinFactorizedWaveform_exact(COMPLEX16 *restrict hlmTab, REAL8Vector *restrict values, REAL8Vector *restrict cartvalues, const REAL8 v, const REAL8 Hreal, const INT4 lMax, SpinEOBParams *restrict params)
FOR PRECESSING EOB This function calculates hlm mode factorized-resummed waveform for given dynamical...
static INT4 XLALSimIMRSpinEOBGetPrecSpinFactorizedWaveform_exact(COMPLEX16 *restrict hlm, REAL8Vector *restrict values, REAL8Vector *restrict cartvalues, const REAL8 v, const REAL8 Hreal, const INT4 l, const INT4 m, SpinEOBParams *restrict params)
This function calculates hlm mode factorized-resummed waveform for given dynamical variables.
static REAL8 XLALSimIMRSpinAlignedEOBNonKeplerCoeff(const REAL8 values[], SpinEOBParams *funcParams)
Function to calculate the non-Keplerian coefficient for the spin-aligned EOB model.
static REAL8 * cross_product(const REAL8 values1[], const REAL8 values2[], REAL8 result[])
static double GSLSpinPrecHamiltonianWrapperFordHdpphi(double x, void *params)
Wrapper for GSL to call the Hamiltonian function.
static REAL8 inner_product(const REAL8 values1[], const REAL8 values2[])
Functions to compute the inner product and cross products between vectors.
static REAL8 UNUSED z2(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
static int XLALSpinPrecHcapRvecDerivative_exact(double UNUSED t, const REAL8 values[], REAL8 dvalues[], void *funcParams)
Function to calculate exact derivatives of the spin EOB Hamiltonian, to get , as decribed in Eqs.
REAL8 Hreal
#define XLAL_CALLGSL(statement)
int l
Definition: bh_qnmode.c:135
double i
Definition: bh_ringdown.c:118
const double pp
#define LAL_E
#define LAL_PI
#define LAL_GAMMA
double complex COMPLEX16
double REAL8
uint32_t UINT4
int32_t INT4
static const INT4 r
static const INT4 m
#define XLAL_ERROR_REAL8(...)
#define XLAL_PRINT_INFO(...)
#define XLAL_ERROR(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_IS_REAL8_FAIL_NAN(val)
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EINVAL
XLAL_FAILURE
string status
Structure containing the coefficients for calculating the factorized waveform.
Parameters for the spinning EOB model.
Definition: burst.c:245
LALDict * params
Definition: inspiral.c:248