LALSimulation  5.4.0.1-fe68b98
LALSimIMREOBNQCCorrection.c
Go to the documentation of this file.
1 /*
2 * Copyright (C) 2010 Craig Robinson, Yi Pan
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, Andrea Taracchini
23  *
24  * \brief More recent versions of the EOB models, such as EOBNRv2 and SEOBNRv1, utilise
25  * a non-quasicircular correction (NQC) to bring the peak of the EOB frequency
26  * into agreement with that of NR simulations. This file contains the functions
27  * used to calculate these NQC corrections, described in DCC document T1100433.
28  * The fits to NR peak amplitude, frequency, and their derivatives, are taken
29  * from Pan et al. PRD 84 124052 (2011) [arXiv:1106.1021], for EOBNRv2, and
30  * from Taracchini et al. PRD 86, 024011 (2012) [arXiv:1202.0790], for SEOBNRv1.
31  */
32 
33 #include <complex.h>
34 #include <math.h>
35 
36 #include <gsl/gsl_vector.h>
37 #include <gsl/gsl_matrix.h>
38 #include <gsl/gsl_spline.h>
39 #include <gsl/gsl_linalg.h>
40 
41 #include "LALSimIMREOBNRv2.h"
42 
43 #ifndef _LALSIMIMRNQCCORRECTION_C
44 #define _LALSIMIMRNQCCORRECTION_C
45 
46 #ifdef __GNUC__
47 #define UNUSED __attribute__ ((unused))
48 #else
49 #define UNUSED
50 #endif
51 
52 #include "LALSimIMREOBNQCTables.c"
53 
54 /* ------------------------------------------------
55  * Non-spin (EOBNRv2)
56  * ------------------------------------------------*/
57 
58 /**
59  * Compute the time offset which should be used in computing the
60  * non-quasicircular correction and performing the ringdown attachment.
61  * These numbers were tuned to numerical relativity simulations, and
62  * are taken from Pan et al, PRD84, 124052(2011) [arXiv:1106.1021], lines 1-5 of Table III.
63  */
64 static REAL8
65 XLALSimIMREOBGetNRPeakDeltaT (INT4 l,
66  /**<< Mode l */
67  INT4 m,
68  /**<< Mode m */
69  REAL8 eta
70  /**<< Symmetric mass ratio */
71  )
72 {
73  switch (l)
74  {
75  case 2:
76  if (m == 2)
77  {
78  return 0.0;
79  }
80  else if (m == 1)
81  {
82  return 10.67 - 2.5 + 9.0 * eta - 41.41 * eta + 76.1 * eta * eta;
83  }
84  else
85  {
86  XLAL_ERROR_REAL8 (XLAL_EINVAL);
87  }
88  break;
89  case 3:
90  if (m == 3)
91  {
92  return 3.383 + 3.847 * eta + 8.979 * eta * eta;
93  }
94  else
95  {
96  XLAL_ERROR_REAL8 (XLAL_EINVAL);
97  }
98  break;
99  case 4:
100  if (m == 4)
101  {
102  return 5.57 - 49.86 * eta + 154.3 * eta * eta;
103  }
104  else
105  {
106  XLAL_ERROR_REAL8 (XLAL_EINVAL);
107  }
108  break;
109  case 5:
110  if (m == 5)
111  {
112  return 6.693 - 34.47 * eta + 102.7 * eta * eta;
113  }
114  else
115  {
116  XLAL_ERROR_REAL8 (XLAL_EINVAL);
117  }
118  break;
119  default:
120  XLAL_ERROR_REAL8 (XLAL_EINVAL);
121  break;
122  }
123 }
124 
125 /**
126  * Function which returns a value of the expected peak amplitude
127  * taken from a fit to numerical relativity simulations. The functions
128  * are taken from Pan et al, PRD84, 124052(2011) [arXiv:1106.1021], lines 1-5 of Table II.
129  */
130 static inline REAL8
131 GetNRPeakAmplitude (INT4 l, /**<< Mode l */
132  INT4 m, /**<< Mode m */
133  REAL8 eta /**<< Symmetric mass ratio */
134  )
135 {
136  switch (l)
137  {
138  case 2:
139  if (m == 2)
140  {
141  return eta * (1.422 + 0.3013 * eta + 1.246 * eta * eta);
142  }
143  else if (m == 1)
144  {
145  return eta * sqrt (1.0 - 4. * eta) * (0.4832 - 0.01032 * eta);
146  }
147  else
148  {
149  XLAL_ERROR_REAL8 (XLAL_EINVAL);
150  }
151  break;
152  case 3:
153  if (m == 3)
154  {
155  return eta * sqrt (1. - 4. * eta) * (0.5761 - 0.09638 * eta +
156  2.715 * eta * eta);
157  }
158  else
159  {
160  XLAL_ERROR_REAL8 (XLAL_EINVAL);
161  }
162  break;
163  case 4:
164  if (m == 4)
165  {
166  return eta * (0.354 - 1.779 * eta + 2.834 * eta * eta);
167  }
168  else
169  {
170  XLAL_ERROR_REAL8 (XLAL_EINVAL);
171  }
172  break;
173  case 5:
174  if (m == 5)
175  {
176  return eta * sqrt (1. - 4. * eta) * (0.1353 - 0.1485 * eta);
177  }
178  else
179  {
180  XLAL_ERROR_REAL8 (XLAL_EINVAL);
181  }
182  break;
183  default:
184  XLAL_ERROR_REAL8 (XLAL_EINVAL);
185  break;
186  }
187 }
188 
189 /**
190  * Function which returns second derivative of the amplitude at the peak
191  * taken from a fit to numerical relativity simulations. The functions
192  * are taken from Pan et al, PRD84, 124052(2011) [arXiv:1106.1021], lines 1-5 of Table II.
193  */
194 static inline REAL8
195 GetNRPeakADDot (INT4 l, /**<< Mode l */
196  INT4 m, /**<< Mode m */
197  REAL8 eta /**<< Symmetric mass ratio */
198  )
199 {
200  switch (l)
201  {
202  case 2:
203  if (m == 2)
204  {
205  return -0.01 * eta * (0.1679 + 1.44 * eta - 2.001 * eta * eta);
206  }
207  else if (m == 1)
208  {
209  return -0.01 * eta * sqrt (1. - 4. * eta) * (0.1867 + 0.6094 * eta);
210  }
211  else
212  {
213  XLAL_ERROR_REAL8 (XLAL_EINVAL);
214  }
215  break;
216  case 3:
217  if (m == 3)
218  {
219  return -0.01 * eta * sqrt (1. - 4. * eta) * (0.2518 - 0.8145 * eta +
220  5.731 * eta * eta);
221  }
222  else
223  {
224  XLAL_ERROR_REAL8 (XLAL_EINVAL);
225  }
226  break;
227  case 4:
228  if (m == 4)
229  {
230  return -0.01 * eta * (0.1813 - 0.9935 * eta + 1.858 * eta * eta);
231  }
232  else
233  {
234  XLAL_ERROR_REAL8 (XLAL_EINVAL);
235  }
236  break;
237  case 5:
238  if (m == 5)
239  {
240  return -0.01 * eta * sqrt (1. - 4. * eta) * (0.09051 -
241  0.1604 * eta);
242  }
243  else
244  {
245  XLAL_ERROR_REAL8 (XLAL_EINVAL);
246  }
247  break;
248  default:
249  XLAL_ERROR_REAL8 (XLAL_EINVAL);
250  break;
251  }
252 }
253 
254 
255 /**
256  * Function which returns a value of the expected peak frequency
257  * taken from a fit to numerical relativity simulations. The functions
258  * are taken from Pan et al, PRD84, 124052(2011) [arXiv:1106.1021], lines 1-5 of Table II.
259  */
260 static inline REAL8
261 GetNRPeakOmega (INT4 l, /**<< Mode l */
262  INT4 m, /**<< Mode m */
263  REAL8 eta /**<< Symmetric mass ratio */
264  )
265 {
266  switch (l)
267  {
268  case 2:
269  if (m == 2)
270  {
271  return 0.2733 + 0.2316 * eta + 0.4463 * eta * eta;
272  }
273  else if (m == 1)
274  {
275  return 0.2907 - 0.08338 * eta + 0.587 * eta * eta;
276  }
277  else
278  {
279  XLAL_ERROR_REAL8 (XLAL_EINVAL);
280  }
281  break;
282  case 3:
283  if (m == 3)
284  {
285  return 0.4539 + 0.5376 * eta + 1.042 * eta * eta;
286  }
287  else
288  {
289  XLAL_ERROR_REAL8 (XLAL_EINVAL);
290  }
291  break;
292  case 4:
293  if (m == 4)
294  {
295  return 0.6435 - 0.05103 * eta + 2.216 * eta * eta;
296  }
297  else
298  {
299  XLAL_ERROR_REAL8 (XLAL_EINVAL);
300  }
301  break;
302  case 5:
303  if (m == 5)
304  {
305  return 0.8217 + 0.2346 * eta + 2.599 * eta * eta;
306  }
307  else
308  {
309  XLAL_ERROR_REAL8 (XLAL_EINVAL);
310  }
311  break;
312  default:
313  XLAL_ERROR_REAL8 (XLAL_EINVAL);
314  break;
315  }
316 }
317 
318 /**
319  * Function which returns the derivative of the expected peak frequency
320  * taken from a fit to numerical relativity simulations. The functions
321  * are taken from Pan et al, PRD84, 124052(2011) [arXiv:1106.1021], lines 1-5 of Table II.
322  */
323 static inline REAL8
324 GetNRPeakOmegaDot (INT4 l, /**<< Mode l */
325  INT4 m, /**<< Mode m */
326  REAL8 eta /**<< Symmetric mass ratio */
327  )
328 {
329  switch (l)
330  {
331  case 2:
332  if (m == 2)
333  {
334  return 0.005862 + 0.01506 * eta + 0.02625 * eta * eta;
335  }
336  else if (m == 1)
337  {
338  return 0.00149 + 0.09197 * eta - 0.1909 * eta * eta;
339  }
340  else
341  {
342  XLAL_ERROR_REAL8 (XLAL_EINVAL);
343  }
344  break;
345  case 3:
346  if (m == 3)
347  {
348  return 0.01074 + 0.0293 * eta + 0.02066 * eta * eta;
349  }
350  else
351  {
352  XLAL_ERROR_REAL8 (XLAL_EINVAL);
353  }
354  break;
355  case 4:
356  if (m == 4)
357  {
358  return 0.01486 + 0.08529 * eta - 0.2174 * eta * eta;
359  }
360  else
361  {
362  XLAL_ERROR_REAL8 (XLAL_EINVAL);
363  }
364  break;
365  case 5:
366  if (m == 5)
367  {
368  return 0.01775 + 0.09801 * eta - 0.1686 * eta * eta;
369  }
370  else
371  {
372  XLAL_ERROR_REAL8 (XLAL_EINVAL);
373  }
374  break;
375  default:
376  XLAL_ERROR_REAL8 (XLAL_EINVAL);
377  break;
378  }
379 }
380 
381 
382 /**
383  * For the 2,2 mode, there are fits available for the NQC coefficients,
384  * given in Eqs.(40a)-(40c) of Pan et al, PRD84, 124052(2011) [arXiv:1106.1021].
385  * This function provides the values of these coefficients, so the
386  * correction can be used in the dynamics prior to finding the more
387  * accurate NQC values later on.
388  */
389 UNUSED static int
390 XLALSimIMREOBGetCalibratedNQCCoeffs (EOBNonQCCoeffs * coeffs,
391  /**<< OUTPUT, Structure for NQC coeffs */
392  INT4 l, /**<< Mode l */
393  INT4 m, /**<< Mode m */
394  REAL8 eta /**<< Symmetric mass ratio */
395  )
396 {
397 
398  if (!coeffs)
399  {
400  XLAL_ERROR (XLAL_EINVAL);
401  }
402 
403  if (l != 2 || m != 2)
404  {
405  XLALPrintError ("Mode %d,%d is not supported by this function.\n", l,
406  m);
407  XLAL_ERROR (XLAL_EINVAL);
408  }
409 
410  /* All NQC coefficients are set to zero here */
411  /* including coeffs->a3S, coeffs->a4 and coeffs->a5 that are not used in EOBNRv2 */
412  memset (coeffs, 0, sizeof (*coeffs));
413 
414  coeffs->a1 = -4.55919 + 18.761 * eta - 24.226 * eta * eta;
415  coeffs->a2 = 37.683 - 201.468 * eta + 324.591 * eta * eta;
416  coeffs->a3 = -39.6024 + 228.899 * eta - 387.222 * eta * eta;
417 
418  return XLAL_SUCCESS;
419 }
420 
421 /**
422  * This function calculates the non-quasicircular correction to apply to
423  * the waveform. The form of this correction can be found in Pan et al,
424  * PRD84, 124052(2011) [arXiv:1106.1021], Eq.(22), and also in the DCC document T1100433. Note
425  * that when calling this function, the NQC coefficients should already
426  * have been pre-computed.
427  */
428 UNUSED static int
429 XLALSimIMREOBNonQCCorrection (COMPLEX16 * restrict nqc, /**<< OUTPUT, The NQC correction */
430  REAL8Vector * restrict values,
431  /**<< Dynamics r, phi, pr, pphi */
432  const REAL8 omega, /**<< Angular frequency */
433  EOBNonQCCoeffs * restrict coeffs
434  /**<< NQC coefficients */
435  )
436 {
437 
438  REAL8 rOmega, rOmegaSq;
439  REAL8 r, p, sqrtR;
440 
441  REAL8 mag, phase;
442 
443 
444  r = values->data[0];
445  p = values->data[2];
446 
447  sqrtR = sqrt (r);
448 
449  rOmega = r * omega;
450  rOmegaSq = rOmega * rOmega;
451 /*printf("a1 = %.16e, a2 = %.16e, a3 = %.16e, a3S = %.16e, a4 = %.16e, a5 = %.16e\n",coeffs->a1,coeffs->a2,coeffs->a3,coeffs->a3S, coeffs->a4,coeffs->a5);
452 printf("b1 = %.16e, b2 = %.16e, b3 = %.16e, b4 = %.16e\n",coeffs->b1,coeffs->b2,coeffs->b3,coeffs->b4);*/
453  /* In EOBNRv2, coeffs->a3S, coeffs->a4 and coeffs->a5 are set to zero */
454  /* through XLALSimIMREOBGetCalibratedNQCCoeffs() */
455  /* and XLALSimIMREOBCalculateNQCCoefficients() */
456  mag = 1. + (p * p / rOmegaSq) * (coeffs->a1
457  + coeffs->a2 / r + (coeffs->a3 +
458  coeffs->a3S) / (r *
459  sqrtR)
460  + coeffs->a4 / (r * r) +
461  coeffs->a5 / (r * r * sqrtR));
462 //printf("NQC INFO mag = %.16e, r = %.16e, p = %.16e\n",mag,r,p);
463  phase = coeffs->b1 * p / rOmega + p * p * p / rOmega * (coeffs->b2
464  +
465  coeffs->b3 / sqrtR +
466  coeffs->b4 / r);
467 
468  *nqc = mag * cos (phase);
469  *nqc += I * mag * sin (phase);
470 /*printf("r = %.16e, pr = %.16e, omega = %.16e\n",r,p,omega);
471 printf("NQC mag = %.16e, arg = %.16e\n",mag,phase);*/
472  return XLAL_SUCCESS;
473 
474 }
475 
476 
477 
478 /**
479  * This function calculates the non-quasicircular correction to apply to
480  * the waveform. The form of this correction can be found in Pan et al,
481  * PRD84, 124052(2011) [arXiv:1106.1021], Eq.(22), and also in the DCC document T1100433. Note
482  * that when calling this function, the NQC coefficients should already
483  * have been pre-computed.
484  * This version is for generic precesing case where the dynamics variable
485  * values are given in Catesean coordinates.
486  */
487 UNUSED static int
488 XLALSimIMRSpinEOBNonQCCorrection (COMPLEX16 * restrict nqc,
489  /**<< OUTPUT, The NQC correction */
490  REAL8Vector * restrict values,
491  /**<< Dynamics r, phi, pr, pphi */
492  const REAL8 omega, /**<< Angular frequency */
493  EOBNonQCCoeffs * restrict coeffs
494  /**<< NQC coefficients */
495  )
496 {
497 
498  REAL8 rOmega, rOmegaSq;
499  REAL8 r, p, sqrtR;
500 
501  REAL8 mag, phase;
502 
503 
504  r =
505  sqrt (values->data[0] * values->data[0] +
506  values->data[1] * values->data[1] +
507  values->data[2] * values->data[2]);
508  p =
509  (values->data[0] * values->data[3] + values->data[1] * values->data[4] +
510  values->data[2] * values->data[5]) / r;
511 
512  sqrtR = sqrt (r);
513 
514  rOmega = r * omega;
515  rOmegaSq = rOmega * rOmega;
516 
517  /* In EOBNRv2, coeffs->a3S, coeffs->a4 and coeffs->a5 are set to zero */
518  /* through XLALSimIMREOBGetCalibratedNQCCoeffs() */
519  /* and XLALSimIMREOBCalculateNQCCoefficients() */
520  mag = 1. + (p * p / rOmegaSq) * (coeffs->a1
521  + coeffs->a2 / r + (coeffs->a3 +
522  coeffs->a3S) / (r *
523  sqrtR)
524  + coeffs->a4 / (r * r) +
525  coeffs->a5 / (r * r * sqrtR));
526 
527  phase = coeffs->b1 * p / rOmega + p * p * p / rOmega * (coeffs->b2
528  +
529  coeffs->b3 / sqrtR +
530  coeffs->b4 / r);
531 
532  *nqc = mag * cos (phase);
533  *nqc += I * mag * sin (phase);
534 
535  return XLAL_SUCCESS;
536 
537 }
538 
539 
540 
541 /**
542  * This function computes the coefficients a1, a2, etc. used in the
543  * non-quasicircular correction. The details of the calculation of these
544  * coefficients are found in the DCC document T1100433.
545  */
546 UNUSED static int
547 XLALSimIMREOBCalculateNQCCoefficients (EOBNonQCCoeffs * restrict coeffs,
548  /**<< OUTPUT, NQC coefficients */
549  REAL8Vector * restrict amplitude,
550  /**<< Waveform amplitude, func of time */
551  REAL8Vector * restrict phase,
552  /**<< Waveform phase(rad), func of time */
553  REAL8Vector * restrict q1,
554  /**<< Function of dynamics (see DCC doc) */
555  REAL8Vector * restrict q2,
556  /**<< Function of dynamics (see DCC doc) */
557  REAL8Vector * restrict q3,
558  /**<< Function of dynamics (see DCC doc) */
559  REAL8Vector * restrict p1,
560  /**<< Function of dynamics (see DCC doc) */
561  REAL8Vector * restrict p2,
562  /**<< Function of dynamics (see DCC doc) */
563  INT4 l, /**<< Mode l */
564  INT4 m, /**<< Mode m */
565  REAL8 timePeak,/**<< Time of peak orbital frequency */
566  REAL8 deltaT, /**<< Sampling interval */
567  REAL8 eta /**<< Symmetric mass ratio */
568  )
569 {
570 
571  UINT4 i;
572 
573  int signum;
574 
575  REAL8Vector *restrict timeVec = NULL;
576 
577  /* Since the vectors we actually want are q etc * A, we will have to generate them here */
578  REAL8Vector *q1LM = NULL;
579  REAL8Vector *q2LM = NULL;
580  REAL8Vector *q3LM = NULL;
581 
582  REAL8 a, aDot, aDDot;
583  REAL8 omega, omegaDot;
584 
585  REAL8 nra, nraDDot;
586  REAL8 nromega, nromegaDot;
587 
588  REAL8 nrDeltaT, nrTimePeak;
589 
590  /* Stuff for finding numerical derivatives */
591  gsl_spline *spline = NULL;
592  gsl_interp_accel *acc = NULL;
593 
594  /* Matrix stuff for calculating coefficients */
595  gsl_matrix *qMatrix = NULL, *pMatrix = NULL;
596  gsl_vector *aCoeff = NULL, *bCoeff = NULL;
597 
598  gsl_vector *amps = NULL, *omegaVec = NULL;
599 
600  gsl_permutation *perm1 = NULL, *perm2 = NULL;
601 
602  /* All NQC coefficients are set to zero here */
603  /* including coeffs->a4 that is not used in EOBNRv2 */
604  memset (coeffs, 0, sizeof (EOBNonQCCoeffs));
605 
606  /* Populate the time vector */
607  /* It is okay to assume initial t = 0 */
608  timeVec = XLALCreateREAL8Vector (q1->length);
609  q1LM = XLALCreateREAL8Vector (q1->length);
610  q2LM = XLALCreateREAL8Vector (q2->length);
611  q3LM = XLALCreateREAL8Vector (q3->length);
612 
613  /* Populate vectors as necessary */
614  for (i = 0; i < timeVec->length; i++)
615  {
616  timeVec->data[i] = i * deltaT;
617  q1LM->data[i] = q1->data[i] * amplitude->data[i];
618  q2LM->data[i] = q2->data[i] * amplitude->data[i];
619  q3LM->data[i] = q3->data[i] * amplitude->data[i];
620  }
621 
622  /* Allocate all the memory we need */
623  XLAL_CALLGSL (
624  /* a stuff */
625  qMatrix = gsl_matrix_alloc (3, 3);
626  aCoeff = gsl_vector_alloc (3);
627  amps = gsl_vector_alloc (3);
628  perm1 = gsl_permutation_alloc (3);
629  /* b stuff */
630  pMatrix = gsl_matrix_alloc (2, 2);
631  bCoeff = gsl_vector_alloc (2);
632  omegaVec = gsl_vector_alloc (2);
633  perm2 = gsl_permutation_alloc (2););
634 
635  if (!qMatrix || !aCoeff || !amps || !pMatrix || !bCoeff || !omegaVec)
636  {
637  gsl_matrix_free (qMatrix);
638  gsl_vector_free (amps);
639  gsl_vector_free (aCoeff);
640  gsl_permutation_free (perm1);
641  gsl_matrix_free (pMatrix);
642  gsl_vector_free (omegaVec);
643  gsl_vector_free (bCoeff);
644  gsl_permutation_free (perm2);
645  XLALDestroyREAL8Vector (q1LM);
646  XLALDestroyREAL8Vector (q2LM);
647  XLALDestroyREAL8Vector (q3LM);
648  XLALDestroyREAL8Vector (timeVec);
649  XLAL_ERROR (XLAL_ENOMEM);
650  }
651 
652  /* The time we want to take as the peak time depends on l and m */
653  /* Calculate the adjustment we need to make here */
654  nrDeltaT = XLALSimIMREOBGetNRPeakDeltaT (l, m, eta);
655  if (XLAL_IS_REAL8_FAIL_NAN (nrDeltaT))
656  {
657  XLALDestroyREAL8Vector (q1LM);
658  XLALDestroyREAL8Vector (q2LM);
659  XLALDestroyREAL8Vector (q3LM);
660  XLALDestroyREAL8Vector (timeVec);
661  XLAL_ERROR (XLAL_EFUNC);
662  }
663 
664  nrTimePeak = timePeak + nrDeltaT;
665  /* We are now in a position to use the interp stuff to calculate the derivatives we need */
666  /* We will start with the quantities used in the calculation of the a coefficients */
667  spline = gsl_spline_alloc (gsl_interp_cspline, amplitude->length);
668  acc = gsl_interp_accel_alloc ();
669 
670  /* Q1 */
671  gsl_spline_init (spline, timeVec->data, q1LM->data, q1LM->length);
672  gsl_matrix_set (qMatrix, 0, 0, gsl_spline_eval (spline, nrTimePeak, acc));
673  gsl_matrix_set (qMatrix, 1, 0,
674  gsl_spline_eval_deriv (spline, nrTimePeak, acc));
675  gsl_matrix_set (qMatrix, 2, 0,
676  gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
677 
678  /* Q2 */
679  gsl_spline_init (spline, timeVec->data, q2LM->data, q2LM->length);
680  gsl_interp_accel_reset (acc);
681  gsl_matrix_set (qMatrix, 0, 1, gsl_spline_eval (spline, nrTimePeak, acc));
682  gsl_matrix_set (qMatrix, 1, 1,
683  gsl_spline_eval_deriv (spline, nrTimePeak, acc));
684  gsl_matrix_set (qMatrix, 2, 1,
685  gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
686 
687  /* Q3 */
688  gsl_spline_init (spline, timeVec->data, q3LM->data, q3LM->length);
689  gsl_interp_accel_reset (acc);
690  gsl_matrix_set (qMatrix, 0, 2, gsl_spline_eval (spline, nrTimePeak, acc));
691  gsl_matrix_set (qMatrix, 1, 2,
692  gsl_spline_eval_deriv (spline, nrTimePeak, acc));
693  gsl_matrix_set (qMatrix, 2, 2,
694  gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
695 
696  /* Amplitude */
697  gsl_spline_init (spline, timeVec->data, amplitude->data, amplitude->length);
698  gsl_interp_accel_reset (acc);
699  a = gsl_spline_eval (spline, nrTimePeak, acc);
700  aDot = gsl_spline_eval_deriv (spline, nrTimePeak, acc);
701  aDDot = gsl_spline_eval_deriv2 (spline, nrTimePeak, acc);
702 
703  nra = GetNRPeakAmplitude (l, m, eta);
704  nraDDot = GetNRPeakADDot (l, m, eta);
705 
706  if (XLAL_IS_REAL8_FAIL_NAN (nra) || XLAL_IS_REAL8_FAIL_NAN (nraDDot))
707  {
708  XLALDestroyREAL8Vector (q1LM);
709  XLALDestroyREAL8Vector (q2LM);
710  XLALDestroyREAL8Vector (q3LM);
711  XLALDestroyREAL8Vector (timeVec);
712  XLAL_ERROR (XLAL_EFUNC);
713  }
714 
715  gsl_vector_set (amps, 0, nra - a);
716  gsl_vector_set (amps, 1, -aDot);
717  gsl_vector_set (amps, 2, nraDDot - aDDot);
718 
719  /* We have now set up all the stuff to calculate the a coefficients */
720  /* So let us do it! */
721  gsl_linalg_LU_decomp (qMatrix, perm1, &signum);
722  gsl_linalg_LU_solve (qMatrix, perm1, amps, aCoeff);
723 
724  /* Now we (should) have calculated the a values. Now we can do the b values */
725 
726  /* P1 */
727  gsl_spline_init (spline, timeVec->data, p1->data, p1->length);
728  gsl_interp_accel_reset (acc);
729  gsl_matrix_set (pMatrix, 0, 0,
730  -gsl_spline_eval_deriv (spline, nrTimePeak, acc));
731  gsl_matrix_set (pMatrix, 1, 0,
732  -gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
733 
734  /* P2 */
735  gsl_spline_init (spline, timeVec->data, p2->data, p2->length);
736  gsl_interp_accel_reset (acc);
737  gsl_matrix_set (pMatrix, 0, 1,
738  -gsl_spline_eval_deriv (spline, nrTimePeak, acc));
739  gsl_matrix_set (pMatrix, 1, 1,
740  -gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
741 
742  /* Phase */
743  gsl_spline_init (spline, timeVec->data, phase->data, phase->length);
744  gsl_interp_accel_reset (acc);
745  omega = gsl_spline_eval_deriv (spline, nrTimePeak, acc);
746  omegaDot = gsl_spline_eval_deriv2 (spline, nrTimePeak, acc);
747 
748  /* Since the phase can be decreasing, we need to take care not to have a -ve frequency */
749  if (omega * omegaDot > 0.0)
750  {
751  omega = fabs (omega);
752  omegaDot = fabs (omegaDot);
753  }
754  else
755  {
756  omega = fabs (omega);
757  omegaDot = -fabs (omegaDot);
758  }
759 
760  nromega = GetNRPeakOmega (l, m, eta);
761  nromegaDot = GetNRPeakOmegaDot (l, m, eta);
762 
763  if (XLAL_IS_REAL8_FAIL_NAN (nromega) || XLAL_IS_REAL8_FAIL_NAN (nromegaDot))
764  {
765  XLALDestroyREAL8Vector (q1LM);
766  XLALDestroyREAL8Vector (q2LM);
767  XLALDestroyREAL8Vector (q3LM);
768  XLALDestroyREAL8Vector (timeVec);
769  XLAL_ERROR (XLAL_EFUNC);
770  }
771 
772  gsl_vector_set (omegaVec, 0, nromega - omega);
773  gsl_vector_set (omegaVec, 1, nromegaDot - omegaDot);
774 
775  /* And now solve for the b coefficients */
776  gsl_linalg_LU_decomp (pMatrix, perm2, &signum);
777  gsl_linalg_LU_solve (pMatrix, perm2, omegaVec, bCoeff);
778 
779  /* We can now populate the coefficients structure */
780  coeffs->a1 = gsl_vector_get (aCoeff, 0);
781  coeffs->a2 = gsl_vector_get (aCoeff, 1);
782  coeffs->a3 = gsl_vector_get (aCoeff, 2);
783  coeffs->b1 = gsl_vector_get (bCoeff, 0);
784  coeffs->b2 = gsl_vector_get (bCoeff, 1);
785 
786  /* Free memory and exit */
787  gsl_matrix_free (qMatrix);
788  gsl_vector_free (amps);
789  gsl_vector_free (aCoeff);
790  gsl_permutation_free (perm1);
791 
792  gsl_matrix_free (pMatrix);
793  gsl_vector_free (omegaVec);
794  gsl_vector_free (bCoeff);
795  gsl_permutation_free (perm2);
796 
797  gsl_spline_free (spline);
798  gsl_interp_accel_free (acc);
799 
800  XLALDestroyREAL8Vector (q1LM);
801  XLALDestroyREAL8Vector (q2LM);
802  XLALDestroyREAL8Vector (q3LM);
803  XLALDestroyREAL8Vector (timeVec);
804 
805  return XLAL_SUCCESS;
806 }
807 
808 /* ------------------------------------------------
809  * Spin (SEOBNRv1 and SEOBNRv2)
810  * ------------------------------------------------*/
811 
812 /**
813  * The time difference between the orbital peak and the peak amplitude
814  * of the mode in question (currently only 2,2 implemented ).
815  * Eq. 33 of Taracchini et al. PRD 86, 024011 (2012) [arXiv:1202.0790].
816  */
817 UNUSED static inline REAL8
818 XLALSimIMREOBGetNRSpinPeakDeltaT (INT4 l,
819  /**<< Mode l */
820  INT4 m,
821  /**<< Mode m */
822  REAL8 UNUSED eta,
823  /**<< Symmetric mass ratio */
824  REAL8 a
825  /**<< Dimensionless spin */
826  )
827 {
828 
829  switch (l)
830  {
831  case 2:
832  switch (m)
833  {
834  case 2:
835  /* DeltaT22 defined here is a minus sign different from Eq. (33) of Taracchini et al. */
836  if (a <= 0.0)
837  {
838  return 2.5;
839  }
840  else
841  {
842  return (2.5 +
843  1.77 * a * a * a * a / (0.43655 * 0.43655 * 0.43655 *
844  0.43655) / (1.0 -
845  2.0 * eta) / (1.0 -
846  2.0 *
847  eta) /
848  (1.0 - 2.0 * eta) / (1.0 - 2.0 * eta));
849  }
850  break;
851  default:
852  XLAL_ERROR_REAL8 (XLAL_EINVAL);
853  }
854  break;
855  default:
856  XLAL_ERROR_REAL8 (XLAL_EINVAL);
857  }
858 
859  /* We should never get here, but I expect a compiler whinge without it... */
860  XLALPrintError ("XLAL Error %s - We should never get here!!\n", __func__);
861  XLAL_ERROR_REAL8 (XLAL_EINVAL);
862 }
863 
864 /**
865  * Peak amplitude predicted by fitting NR results (currently only 2,2 available).
866  * Unpublished. Used in building SEOBNRv2 tables.
867  */
868 UNUSED static inline REAL8
869 XLALSimIMREOBGetNRSpinPeakAmplitudeV2 (INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta,
870  REAL8 UNUSED a)
871 {
872  REAL8 chi = a, chi2 = chi * chi, chi3 = chi * chi2;
873  REAL8 eta2 = eta * eta;
874  REAL8 res;
875  if (chi > 0.8)
876  {
877  res = eta * (56.28859370276537 * (-0.1858184673895021 +
878  eta) * (0.18616944529501114 + eta) -
879  155.11365222671776 * chi3 * (-0.25025223669804486 +
880  eta) * (0.23614334403451426 +
881  eta) +
882  322.4309641674941 * chi2 * (-0.24986765309607953 +
883  eta) * (0.24802475468124208 +
884  eta) -
885  226.09242469439047 * chi * (-0.24993985462384588 +
886  eta) * (0.2573225045218015 +
887  eta));
888  }
889  else
890  {
891  res = eta * (1.449934273310975 +
892  3.8867234144877933 * chi * (-0.26967339454732164 +
893  eta) * (-0.15977862405445659 +
894  eta) +
895  2.2705573440821687 * chi2 * (-0.20039719578235954 +
896  eta) * (-0.06130397389190033 +
897  eta) -
898  8.094119513915285 * chi3 * (-0.2598144292071539 +
899  eta) * (-0.010564809220517786 +
900  eta) +
901  0.019756052721845246 * eta + 1.933934833691488 * eta2);
902  }
903  return res;
904 }
905 
906 /**
907  * Combine two spin-dependent fits of NR input values in a quardatic-in-eta polynomial
908  * with the linear-in-eta coefficient A1 being provided by a polyinamial-in-spin global
909  * fit to NR in the whole eta-chi plane. See Eqs. (20)-(21) in https://dcc.ligo.org/T1600383
910  */
911 static inline REAL8 CombineTPLEQMFits (REAL8 eta, REAL8 A1, REAL8 fEQ, REAL8 fTPL)
912 {
913  REAL8 A0, A2;
914  REAL8 eta2 = eta * eta;
915  // Impose that TPL and equal-mass limit are exactly recovered
916  A0 = -0.00099601593625498 * A1 - 0.00001600025600409607 * fEQ + 1.000016000256004 * fTPL;
917  A2 = -3.984063745019967 * A1 + 16.00025600409607 * fEQ - 16.0002560041612 * fTPL;
918  // Final formula
919  return A0 + A1 * eta + A2 * eta2;
920 }
921 
922 
923 /**
924  * Peak amplitude predicted by fitting NR results
925  * The coefficients for the mode (2,2) are in Eq.(A1) of https://journals.aps.org/prd/pdf/10.1103/PhysRevD.95.044028
926  * The coefficients for the modes (2,1), (3,3), (4,4), (5,5) are in Eqs. (B7-B10) of https://arxiv.org/pdf/1803.10701.pdf
927  */
928 UNUSED static inline REAL8
929 XLALSimIMREOBGetNRSpinPeakAmplitudeV4 (INT4 modeL, INT4 modeM, REAL8 m1, REAL8 m2,
930  REAL8 chiS, REAL8 chiA)
931 {
932  REAL8 eta = (m1 * m2) / ((m1 + m2) * (m1 + m2));
933  REAL8 dM = sqrt(1.-4.*eta);
934  REAL8 tempm1;
935  if (m1 < m2)
936  {
937  //RC: The fits for the HMs are done under the assumption m1>m2, so if m2>m1 we just swap the two bodies
938  tempm1 = m1;
939  m1 = m2;
940  m2 = tempm1;
941  chiA = -chiA;
942  }
943  REAL8 eta2 = eta*eta;
944  REAL8 chi = chiS + chiA * (m1 - m2) / (m1 + m2) / (1. -
945  2. *
946  eta);
947  REAL8 chi21 = chiS*dM/(1.-1.3*eta) + chiA;
948  REAL8 chi33 = chiS*dM + chiA;
949  REAL8 chi44 = chiS*(1-5*eta) + chiA*dM;
950  REAL8 chi2 = chi * chi, chi3 = chi * chi2;
951  REAL8 res;
952  REAL8 fTPL, fEQ, A1, e0, e1, e2, e3;
953  switch (modeL) {
954  case 2:
955  switch (modeM) {
956  case 2:
957  // TPL fit
958  fTPL = 1.4528573105413543 + 0.16613449160880395 * chi + 0.027355646661735258 * chi2 - 0.020072844926136438 * chi3;
959  // Equal-mass fit
960  fEQ = 1.577457498227 - 0.0076949474494639085 * chi + 0.02188705616693344 * chi2 + 0.023268366492696667 * chi3;
961  // Global fit coefficients
962  e0 = -0.03442402416125921;
963  e1 = -1.218066264419839;
964  e2 = -0.5683726304811634;
965  e3 = 0.4011143761465342;
966  A1 = e0 + e1 * chi + e2 * chi2 + e3 * chi3;
967  res = eta * CombineTPLEQMFits(eta, A1, fEQ, fTPL);
968  break;
969  case 1:
970 
971  res = -((0.29256703361640224-0.19710255145276584*eta)*eta*chi21 + dM*eta*(-0.42817941710649793 + 0.11378918021042442*eta-0.7736772957051212*eta2
972  +chi21*chi21*(0.047004057952214004 - eta*0.09326128322462478)) +dM*eta*chi21*(-0.010195081244587765 + 0.016876911550777807*chi21*chi21));
973 
974  break;
975 
976  default:
977  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
978  XLAL_ERROR (XLAL_EINVAL);
979  break;
980  }
981  break;
982 
983  case 3:
984  switch (modeM) {
985  case 3:
986  res = (0.10109183988848384*eta - 0.4704095462146807*eta2 + 1.0735457779890183*eta2*eta)*chi33 +
987  dM*(0.5636580081367962*eta - 0.054609013952480856*eta2 + 2.3093699480319234*eta2*eta + chi33*chi33*(0.029812986680919126*eta - 0.09688097244145283*eta2) );
988  break;
989 
990  default:
991  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
992  XLAL_ERROR (XLAL_EINVAL);
993  break;
994  }
995  break;
996 
997  case 4:
998  switch (modeM) {
999  case 4:
1000  res = eta*(0.2646580063832686 + 0.067584186955327*chi44 +0.02925102905737779*chi44*chi44) + eta2 *(-0.5658246076387973 -0.8667455348964268*chi44 +0.005234192027729502*chi44*chi44)
1001  + eta*eta2*(-2.5008294352355405 + 6.880772754797872*chi44 -1.0234651570264885*chi44*chi44) + eta2*eta2*(7.6974501716202735 -16.551524307203252*chi44);
1002  break;
1003  default:
1004  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1005  XLAL_ERROR (XLAL_EINVAL);
1006  break;
1007  }
1008  break;
1009  case 5:
1010  switch (modeM) {
1011  case 5:
1012 
1013  res = 0.128621*dM *eta -0.474201 *dM *eta*eta +1.0833 * dM *eta*eta*eta + 0.0322784 * eta * chi33 -0.134511 *chi33 *eta * eta +0.0990202 *chi33*eta*eta*eta;
1014 
1015  break;
1016 
1017  default:
1018  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1019  XLAL_ERROR (XLAL_EINVAL);
1020  break;
1021  }
1022  break;
1023 
1024  default:
1025  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1026  XLAL_ERROR (XLAL_EINVAL);
1027  break;
1028  }
1029 // printf("A %.16e\n", res);
1030  return res;
1031 }
1032 
1033 /**
1034  * Peak amplitude slope predicted by fitting NR results.
1035  * The coefficient for the (2,2) is 0 because the
1036  * attachment is done at the peak of the (2,2) mode see Eq.(2.7) in https://journals.aps.org/prd/pdf/10.1103/PhysRevD.95.044028
1037  * The coefficients for the modes (2,1), (3,3), (4,4) and (5,5) are in Eqs.(B11-B14) of https://arxiv.org/pdf/1803.10701.pdf
1038  */
1039 UNUSED static inline REAL8
1040 XLALSimIMREOBGetNRSpinPeakADotV4 (INT4 modeL, INT4 modeM, REAL8 UNUSED m1, REAL8 UNUSED m2,
1041  REAL8 UNUSED chiS, REAL8 UNUSED chiA)
1042 
1043 {
1044  REAL8 eta = (m1 * m2) / ((m1 + m2) * (m1 + m2));
1045  REAL8 dM = sqrt(1.-4.*eta);
1046  REAL8 dM2 = dM*dM;
1047  REAL8 tempm1;
1048  if (m1 < m2)
1049  {
1050  //RC: The fits for the HMs are done under the assumption m1>m2, so if m2>m1 we just swap the two bodies
1051  tempm1 = m1;
1052  m1 = m2;
1053  m2 = tempm1;
1054  chiA = -chiA;
1055  }
1056  REAL8 eta2 = eta*eta;
1057  REAL8 chi21 = chiS*dM/(1.-2.*eta) + chiA;
1058  REAL8 chi33 = chiS*dM + chiA;
1059  REAL8 chi44 = chiS*(1-7*eta) + chiA*dM;
1060  REAL8 res;
1061  switch (modeL) {
1062  case 2:
1063  switch (modeM) {
1064  case 2:
1065  res = 0.;
1066  break;
1067  case 1:
1068  res = dM*eta*(0.007147528020812309-eta*0.035644027582499495) + dM*eta*chi21*(-0.0087785131749995 + eta*0.03054672006241107) + eta*0.00801714459112299*fabs(-dM*(0.7875612917853588 + eta*1.161274164728927 + eta2*11.306060006923605)+chi21);
1069 
1070  break;
1071  default:
1072  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1073  XLAL_ERROR (XLAL_EINVAL);
1074  break;
1075  }
1076  break;
1077  case 3:
1078  switch (modeM) {
1079  case 3:
1080  res = dM*eta*(-0.00309943555972098 + eta*0.010076527264663805)*chi33*chi33 +
1081 
1082  eta*0.0016309606446766923*sqrt(dM2*(8.811660714437027 + 104.47752236009688*eta) + dM*chi33*(-5.352043503655119 + eta*49.68621807460999) + chi33*chi33);
1083  break;
1084  default:
1085  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1086  XLAL_ERROR (XLAL_EINVAL);
1087  break;
1088  }
1089  break;
1090  case 4:
1091  switch (modeM) {
1092  case 4:
1093  res = eta*(0.004347588211099233 -0.0014612210699052148*chi44 -0.002428047910361957*chi44*chi44) + eta2*(0.023320670701084355-0.02240684127113227*chi44+0.011427087840231389*chi44*chi44)+
1094  eta*eta2*(-0.46054477257132803 + 0.433526632115367*chi44) + eta2*eta2*(1.2796262150829425-1.2400051122897835*chi44);
1095  break;
1096  default:
1097  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1098  XLAL_ERROR (XLAL_EINVAL);
1099  break;
1100  }
1101  break;
1102 
1103  case 5:
1104  switch (modeM) {
1105  case 5:
1106  res = eta * (dM*(-0.008389798844109389 + 0.04678354680410954*eta) + dM*chi33*(-0.0013605616383929452 + 0.004302712487297126*eta) +dM*chi33*chi33*(-0.0011412109287400596 + 0.0018590391891716925*eta) +
1107 
1108  0.0002944221308683548*fabs(dM*(37.11125499129578 - 157.79906814398277*eta) + chi33));
1109  break;
1110  default:
1111  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1112  XLAL_ERROR (XLAL_EINVAL);
1113  break;
1114  }
1115  break;
1116  default:
1117  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1118  XLAL_ERROR (XLAL_EINVAL);
1119  break;
1120  }
1121  // printf("ddA %.16e\n", res);
1122  return res;
1123 }
1124 
1125 
1126 /**
1127  * Peak amplitude curvature predicted by fitting NR results (currently only 2,2 available).
1128  * Unpublished. Used in building SEOBNRv2 tables.
1129  */
1130 UNUSED static inline REAL8
1131 XLALSimIMREOBGetNRSpinPeakADDotV2 (INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta,
1132  REAL8 UNUSED a)
1133 {
1134  REAL8 chi = a;
1135  REAL8 res;
1136  if (chi > 0.8)
1137  {
1138  res = eta * (0.04818855743392375 * chi * (-0.15271017198256195 +
1139  eta) * (0.2154794478856639 +
1140  eta) -
1141  0.038199344157345716 * (-0.06621184971565616 +
1142  eta) * (0.31317077454081577 +
1143  eta));
1144  }
1145  else
1146  {
1147  res = eta * (0.010916757595083287 * (-1.0608229327701018 +
1148  eta) * (0.19667724848989968 +
1149  eta) -
1150  0.007331284524315633 * chi * (-0.753628708015681 +
1151  eta) * (0.341046049832081 +
1152  eta) +
1153  chi * chi * (0.0006958660609341137 -
1154  0.01113385697494582 * eta * eta) +
1155  chi * chi * chi * (-0.00029607425270115136 +
1156  0.004737188043218422 * eta * eta));
1157  }
1158  return res;
1159 }
1160 
1161 /**
1162  * Peak amplitude curvature predicted by fitting NR results.
1163  * The coefficients for the (2,2) are in Eq.(A3-A4) of https://journals.aps.org/prd/pdf/10.1103/PhysRevD.95.044028
1164  * The coefficients for the modes (2,1), (3,3), (4,4) and (5,5) are in Eqs.(B15-B18) of https://arxiv.org/pdf/1803.10701.pdf
1165  */
1166 UNUSED static inline REAL8
1167 XLALSimIMREOBGetNRSpinPeakADDotV4 (INT4 modeL, INT4 modeM, REAL8 m1, REAL8 m2,
1168  REAL8 chiS, REAL8 chiA)
1169 {
1170  REAL8 eta = (m1 * m2) / ((m1 + m2) * (m1 + m2));
1171  REAL8 dM = sqrt(1.-4.*eta);
1172  REAL8 tempm1;
1173  if (m1 < m2)
1174  {
1175  //RC: The fits for the HMs are done under the assumption m1>m2, so if m2>m1 we just swap the two bodies
1176  tempm1 = m1;
1177  m1 = m2;
1178  m2 = tempm1;
1179  chiA = -chiA;
1180  }
1181  REAL8 eta2 = eta*eta;
1182  REAL8 chi21 = chiS*dM/(1.-2.*eta) + chiA;
1183  REAL8 chi = chiS + chiA * (m1 - m2) / (m1 + m2) / (1. -
1184  2. *
1185  eta);
1186  REAL8 chi33 = chiS*dM + chiA;
1187  REAL8 chiMinus1 = -1. + chi;
1188  REAL8 res;
1189  REAL8 fTPL, fEQ, A1, e0, e1;
1190  switch (modeL) {
1191  case 2:
1192  switch (modeM) {
1193  case 2:
1194  // TPL fit
1195  fTPL = 0.002395610769995033 * chiMinus1 - 0.00019273850675004356 * chiMinus1 * chiMinus1 - 0.00029666193167435337 * chiMinus1 * chiMinus1 * chiMinus1;
1196  // Equal-mass fit
1197  fEQ = -0.004126509071377509 + 0.002223999138735809 * chi;
1198  // Global fit coefficients
1199  e0 = -0.005776537350356959;
1200  e1 = 0.001030857482885267;
1201  A1 = e0 + e1 * chi;
1202  res = eta * CombineTPLEQMFits(eta, A1, fEQ, fTPL);;
1203  break;
1204  case 1:
1205 
1206  res = eta*dM*0.00037132201959950333 -
1207  fabs(dM*eta*(-0.0003650874948532221 - eta*0.003054168419880019)
1208  +dM*eta*chi21*chi21*(-0.0006306232037821514-eta*0.000868047918883389 + eta2*0.022306229435339213)+eta*chi21*chi21*chi21*0.0003402427901204342+dM*eta*chi21*0.00028398490492743);
1209 
1210  break;
1211  default:
1212  XLALPrintError("XLAL Error - %s: At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__);
1213  XLAL_ERROR (XLAL_EINVAL);
1214  break;
1215  }
1216  break;
1217 
1218  case 3:
1219  switch (modeM) {
1220  case 3:
1221  res = dM*eta*(0.0009605689249339088 - 0.00019080678283595965*eta)*chi33 - 0.00015623760412359145*eta*fabs(dM*(4.676662024170895 + 79.20189790272218*eta - 1097.405480250759*eta2 + 6512.959044311574*eta*eta2 -13263.36920919937*eta2*eta2) + chi33);
1222  break;
1223  default:
1224  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1225  XLAL_ERROR (XLAL_EINVAL);
1226  break;
1227  }
1228  break;
1229  case 4:
1230  switch (modeM) {
1231  case 4:
1232  res = eta*(-0.000301722928925693 + 0.0003215952388023551*chi) + eta2*(0.006283048344165004 + 0.0011598784110553046*chi) + eta2*eta*(-0.08143521096050622 - 0.013819464720298994*chi)+
1233  eta2*eta2*(0.22684871200570564 + 0.03275749240408555*chi);
1234  break;
1235  default:
1236  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1237  XLAL_ERROR (XLAL_EINVAL);
1238  break;
1239  }
1240  break;
1241  case 5:
1242  switch (modeM) {
1243  case 5:
1244  res = eta * (dM *(0.00012727220842255978 + 0.0003211670856771251*eta) + dM*chi33*(-0.00006621677859895541 + 0.000328855327605536*eta) + chi33*chi33*(-0.00005824622885648688 + 0.00013944293760663706*eta));
1245 
1246  break;
1247  default:
1248  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1249  XLAL_ERROR (XLAL_EINVAL);
1250  break;
1251  }
1252  break;
1253 
1254  default:
1255  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1256  XLAL_ERROR (XLAL_EINVAL);
1257  break;
1258  }
1259  // printf("ddA %.16e\n", res);
1260  return res;
1261 }
1262 
1263 
1264 /**
1265  * Peak frequency predicted by fitting NR results (currently only 2,2 available).
1266  * Unpublished. Used in building SEOBNRv2 tables.
1267  */
1268 UNUSED static inline REAL8
1269 XLALSimIMREOBGetNRSpinPeakOmegaV2 (INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 a)
1270 {
1271  REAL8 chi = a;
1272  REAL8 eta2 = eta * eta, eta3 = eta * eta2;
1273  REAL8 res;
1274  if (eta > 50. / 51. / 51.)
1275  {
1276  res = 0.43747541927878864 + (-0.10933208665273314 -
1277  0.007325831113333813 * chi) *
1278  log (4.500844771420863 - 9.681916048928946 * eta +
1279  chi * (-4.254886879579986 + 11.513558950322647 * eta));
1280  }
1281  else
1282  {
1283  res = 1.5609526077704716 - 122.25721149839733 * eta +
1284  3586.2192688666914 * eta2 -
1285  13869.506144441548 * eta3 + (-0.25 +
1286  eta) * (1651.5823693445805 *
1287  (-0.019223375977400495 +
1288  eta) * (-0.01922337527211892 +
1289  eta) +
1290  66.87492814925524 * chi *
1291  (0.0003695381704106058 -
1292  0.03844675124951941 * eta +
1293  eta2)) * log (5600.67382718678 -
1294  5555.824895398546
1295  * chi) +
1296  (-1412.8186461833657 + 67.66455403259023 * chi) * (-0.001 +
1297  eta) *
1298  (0.0003695381704106056 - 0.038446751249519406 * eta +
1299  eta2) * log (0.5680439481719505 - 0.36813967358200156 * chi) +
1300  0.012328326527732041 * log (4.500844771420863 -
1301  9.681916048928946 * eta +
1302  chi * (-4.254886879579986 +
1303  11.513558950322647 * eta)) +
1304  0.0008260634258180991 * chi * log (4.500844771420863 -
1305  9.681916048928946 * eta +
1306  chi * (-4.254886879579986 +
1307  11.513558950322647 * eta)) -
1308  12.6575493872956 * eta * log (4.500844771420863 -
1309  9.681916048928946 * eta +
1310  chi * (-4.254886879579986 +
1311  11.513558950322647 * eta)) -
1312  0.8481231078533651 * chi * eta * log (4.500844771420863 -
1313  9.681916048928946 * eta +
1314  chi * (-4.254886879579986 +
1315  11.513558950322647 *
1316  eta)) +
1317  329.2228595635586 * eta2 * log (4.500844771420863 -
1318  9.681916048928946 * eta +
1319  chi * (-4.254886879579986 +
1320  11.513558950322647 * eta)) +
1321  22.05968203526603 * chi * eta2 * log (4.500844771420863 -
1322  9.681916048928946 * eta +
1323  chi * (-4.254886879579986 +
1324  11.513558950322647 *
1325  eta));
1326  }
1327  return res;
1328 }
1329 
1330 /**
1331  * Peak frequency predicted by fitting NR results
1332  * The coefficients for the (2,2) are in Eq.(A6) of https://journals.aps.org/prd/pdf/10.1103/PhysRevD.95.044028
1333  * The coefficients for the modes (2,1), (3,3), (4,4) and (5,5) are in Eqs.(B19-B22) of https://arxiv.org/pdf/1803.10701.pdf
1334  */
1335 UNUSED static inline REAL8
1336 XLALSimIMREOBGetNRSpinPeakOmegaV4 (INT4 modeL, INT4 modeM, REAL8 UNUSED eta, REAL8 a)
1337 {
1338  REAL8 chi = a, eta2 = eta*eta;
1339  REAL8 res;
1340  REAL8 c0, c1, c2, c3, c4, d2, d3, A3, A4;
1341  switch (modeL) {
1342  case 2:
1343  switch (modeM) {
1344  case 2:
1345  // From TPL fit
1346  c0 = 0.5626787200433265;
1347  c1 = -0.08706198756945482;
1348  c2 = 25.81979479453255;
1349  c3 = 25.85037751197443;
1350  // From equal-mass fit
1351  d2 = 7.629921628648589;
1352  d3 = 10.26207326082448;
1353  // Combine TPL and equal-mass
1354  A4 = d2 + 4 * (d2 - c2) * (eta - 0.25);
1355  A3 = d3 + 4 * (d3 - c3) * (eta - 0.25);
1356  c4 = 0.00174345193125868;
1357  // Final formula
1358  res = c0 + (c1 + c4 * chi) * log(A3 - A4 * chi);
1359  break;
1360  case 1:
1361  res = (0.1743194440996283 + eta*0.1938944514123048 + 0.1670063050527942*eta2 + 0.053508705425291826 *chi - eta*chi*0.18460213023455802 + eta2*chi*0.2187305149636044
1362  +chi*chi*0.030228846150378793 - eta*chi*chi*0.11222178038468673);
1363  break;
1364  default:
1365  XLALPrintError("XLAL Error - %s: At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__);
1366  XLAL_ERROR (XLAL_EINVAL);
1367  break;
1368  }
1369  break;
1370 
1371  case 3:
1372  switch (modeM) {
1373  case 3:
1374  res = 0.3973947703114506 + 0.16419332207671075*chi + 0.1635531186118689*chi*chi + 0.06140164491786984*chi*chi*chi+
1375  eta*(0.6995063984915486-0.3626744855912085*chi -0.9775469868881651*chi*chi)+ eta2*(-0.3455328417046369+0.31952307610699876*chi + 1.9334166149686984*chi*chi);
1376  break;
1377  default:
1378  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1379  XLAL_ERROR (XLAL_EINVAL);
1380  break;
1381  }
1382  break;
1383  case 4:
1384  switch (modeM) {
1385  case 4:
1386  res = 0.5389359134370971 + 0.16635177426821202*chi + 0.2075386047689103*chi*chi + 0.15268115749910835*chi*chi*chi +
1387  eta*(0.7617423831337586 + 0.009587856087825369*chi - 1.302303785053009*chi*chi - 0.5562751887042064*chi*chi*chi)
1388  + eta2*(0.9675153069365782 - 0.22059322127958586*chi + 2.678097398558074*chi*chi) - eta2*eta*4.895381222514275;
1389  break;
1390  default:
1391  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1392  XLAL_ERROR (XLAL_EINVAL);
1393  break;
1394  }
1395  break;
1396  case 5:
1397  switch (modeM) {
1398  case 5:
1399  res = 0.6437545281817488 + 0.22315530037902315*chi + 0.2956893357624277*chi*chi + 0.17327819169083758*chi*chi*chi +
1400  eta*(-0.47017798518175785 - 0.3929010618358481*chi - 2.2653368626130654*chi*chi - 0.5512998466154311*chi*chi*chi) +
1401  eta2*(2.311483807604238 + 0.8829339243493562*chi + 5.817595866020152*chi*chi);
1402  break;
1403  default:
1404  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1405  XLAL_ERROR (XLAL_EINVAL);
1406  break;
1407  }
1408  break;
1409 
1410  default:
1411  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1412  XLAL_ERROR (XLAL_EINVAL);
1413  break;
1414  }
1415 // printf("w %.16e\n", res);
1416  return res;
1417 }
1418 
1419 /**
1420  * Peak frequency predicted by fitting NR results used in SEOBNRv5
1421  * The coefficients are in Eq.(C29) of https://dcc.ligo.org/LIGO-P2300070
1422  */
1423 UNUSED static inline REAL8
1424 XLALSimIMREOBGetNRSpinPeakOmegaV5 (INT4 modeL, INT4 modeM, REAL8 UNUSED eta, REAL8 a)
1425 {
1426  REAL8 chi = a, chi2=chi*chi, chi3=chi2*chi, chi4=chi2*chi2;
1427  REAL8 eta2 = eta*eta, eta3 = eta2*eta, eta4 = eta2*eta2;
1428  REAL8 res;
1429  switch (modeL) {
1430  case 2:
1431  switch (modeM) {
1432  case 2:
1433  res = -1. * (5.89352329617707670906 * eta4
1434  + 3.75145580491965446868 * eta3 * chi
1435  - 3.34930536209472551334 * eta3
1436  - 0.97140932625194231775 * eta2 * chi2
1437  - 1.69734302394369973577 * eta2 * chi
1438  + 0.28539204856044564362 * eta2
1439  + 0.2419483723662931296 * eta * chi3
1440  + 0.51801427018052081941 * eta * chi2
1441  + 0.25096450064948544467 * eta * chi
1442  - 0.31709602351033533418 * eta
1443  - 0.01525897668158244028 * chi4
1444  - 0.06692658483513916345 * chi3
1445  - 0.08715176045684569495 * chi2
1446  - 0.09133931944098934441 * chi
1447  - 0.2685414392185025978);
1448  break;
1449  default:
1450  XLALPrintError("XLAL Error - %s: At present only fits for the (2,2) mode are available.\n", __func__);
1451  XLAL_ERROR (XLAL_EINVAL);
1452  break;
1453  }
1454  break;
1455  }
1456 // printf("w %.16e\n", res);
1457  return res;
1458 }
1459 
1460 /**
1461  * Peak frequency slope predicted by fitting NR results (currently only 2,2 available).
1462  * Unpublished. Used in building SEOBNRv2 tables.
1463  */
1464 UNUSED static inline REAL8
1465 XLALSimIMREOBGetNRSpinPeakOmegaDotV2 (INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta,
1466  REAL8 UNUSED a)
1467 {
1468  REAL8 chi = a;
1469  REAL8 res;
1470  if (chi > 0.8)
1471  {
1472  res =
1473  -0.10069512275335238 * (-0.46107388514323044 +
1474  eta) * (0.2832795481380979 + eta) +
1475  0.2614619716504706 * chi * (-0.24838163750494138 +
1476  eta) * (0.320112993649413 + eta) +
1477  chi * chi * (0.010000160002560042 - 0.16000256004096067 * eta * eta);
1478  }
1479  else
1480  {
1481  res = -0.07086074186161867 * chi * (-0.26367236731979804 +
1482  eta) * (-0.0010019969893089581 +
1483  eta) +
1484  0.2893863668183948 * (-0.16845695144529893 +
1485  eta) * (0.23032241797163952 + eta) +
1486  (0.004086861548547749 - 0.06538978477676398 * eta * eta +
1487  chi * (0.0006334026884930817 -
1488  0.010134443015889307 * eta * eta)) * log (68.47466578101876 -
1489  58.30148755701496 *
1490  chi);
1491  }
1492  return res;
1493 }
1494 
1495 /**
1496  * Peak frequency slope predicted by fitting NR results
1497  * The coefficients for the (2,2) are in Eq.(A10) of https://journals.aps.org/prd/pdf/10.1103/PhysRevD.95.044028
1498  * The coefficients for the modes (2,1), (3,3), (4,4) and (5,5) are in Eqs.(B23-B26) of https://arxiv.org/pdf/1803.10701.pdf
1499  */
1500 UNUSED static inline REAL8
1501 XLALSimIMREOBGetNRSpinPeakOmegaDotV4 (INT4 modeL, INT4 modeM, REAL8 UNUSED eta,
1502  REAL8 UNUSED a)
1503 {
1504  REAL8 chi = a, eta2 = eta*eta;
1505  REAL8 res;
1506  REAL8 fTPL, fEQ, A1, e0, e1;
1507  switch (modeL) {
1508  case 2:
1509  switch (modeM) {
1510  case 2:
1511  // TPL fit
1512  fTPL = -0.011209791668428353 + (0.0040867958978563915 + 0.0006333925136134493 * chi) * log(68.47466578100956 - 58.301487557007206 * chi);
1513  // Equal-mass fit
1514  fEQ = 0.01128156666995859 + 0.0002869276768158971* chi;
1515  // Global fit coefficients
1516  e0 = 0.01574321112717377;
1517  e1 = 0.02244178140869133;
1518  A1 = e0 + e1 * chi;
1519  res = CombineTPLEQMFits(eta, A1, fEQ, fTPL);;
1520  break;
1521  case 1:
1522  res = (0.0070987396362959514 + eta*0.024816844694685373 -eta2*0.050428973182277494 + eta*eta2*0.03442040062259341-chi*0.0017751850002442097+eta*chi*0.004244058872768811
1523  -eta2*chi*0.031996494883796855-chi*chi*0.0035627260615894584+eta*chi*chi*0.01471807973618255 - chi*chi*chi*0.0019020967877681962);
1524  break;
1525  default:
1526  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1527  XLAL_ERROR (XLAL_EINVAL);
1528  break;
1529  }
1530  break;
1531 
1532  case 3:
1533  switch (modeM) {
1534  case 3:
1535  res = 0.010337157192240338 - 0.0053067782526697764*chi*chi - 0.005087932726777773*chi*chi*chi+
1536  eta*(0.027735564986787684 + 0.018864151181629343*chi + 0.021754491131531044*chi*chi + 0.01785477515931398*chi*chi*chi)+
1537  eta2*(0.018084233854540898 - 0.08204268775495138*chi);
1538  break;
1539 
1540  default:
1541  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1542 
1543  XLAL_ERROR (XLAL_EINVAL);
1544  break;
1545  }
1546  break;
1547 
1548  case 4:
1549  switch (modeM) {
1550  case 4:
1551  res = 0.013997911323773867 - 0.0051178205260273574*chi - 0.0073874256262988*chi*chi +
1552  eta*(0.0528489379269367 + 0.01632304766334543*chi + 0.02539072293029433*chi*chi)
1553  +eta2*(-0.06529992724396189 + 0.05782894076431308*chi);
1554  break;
1555 
1556  default:
1557  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1558  XLAL_ERROR (XLAL_EINVAL);
1559  break;
1560  }
1561  break;
1562  case 5:
1563  switch (modeM) {
1564  case 5:
1565  res = 0.01763430670755021 - 0.00024925743340389135*chi - 0.009240404217656968*chi*chi - 0.007907831334704586*chi*chi*chi+
1566  eta*(-0.1366002854361568 + 0.0561378177186783*chi + 0.16406275673019852*chi*chi + 0.07736232247880881*chi*chi*chi)+
1567  eta2*(0.9875890632901151 - 0.31392112794887855*chi - 0.5926145463423832*chi*chi) - 1.6943356548192614*eta2*eta;
1568  break;
1569 
1570  default:
1571  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1572  XLAL_ERROR (XLAL_EINVAL);
1573  break;
1574  }
1575  break;
1576  default:
1577  XLALPrintError("XLAL Error - %s: (%d,%d) mode. At present only fits for the (2,2), (2,1), (3,3), (4,4) and (5,5) mode are available.\n", __func__, modeL, modeM);
1578  XLAL_ERROR (XLAL_EINVAL);
1579  break;
1580  }
1581 // printf("dw %.16e\n", res);
1582  return res;
1583 }
1584 
1585 /**
1586  * Peak amplitude predicted by fitting NR results (currently only 2,2 available).
1587  * Tables IV and V and Eq. 42 of Taracchini et al. PRD 86, 024011 (2012) [arXiv:1202.0790].
1588  */
1589 UNUSED static inline REAL8
1590 XLALSimIMREOBGetNRSpinPeakAmplitude (INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta,
1591  REAL8 UNUSED a)
1592 {
1593  /* Fit for HOMs missing */
1594  return 1.3547468629743946 * eta + 0.9187885481024214 * eta * eta;
1595 }
1596 
1597 /**
1598  * Peak amplitude curvature predicted by fitting NR results (currently only 2,2 available).
1599  * Tables IV and V and Eq. 42 of Taracchini et al. PRD 86, 024011 (2012) [arXiv:1202.0790].
1600  */
1601 UNUSED static inline REAL8
1602 XLALSimIMREOBGetNRSpinPeakADDot (INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta,
1603  REAL8 UNUSED a)
1604 {
1605  /* Fit for HOMs missing */
1606  return eta * (-0.0024971911410897156 +
1607  (-0.006128515435641139 +
1608  0.01732656 * a / (2.0 - 4.0 * eta)) * eta);
1609 }
1610 
1611 /**
1612  * Peak frequency predicted by fitting NR results (currently only 2,2 available).
1613  * Tables IV and V and Eq. 42 of Taracchini et al. PRD 86, 024011 (2012) [arXiv:1202.0790].
1614  */
1615 UNUSED static inline REAL8
1616 XLALSimIMREOBGetNRSpinPeakOmega (INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 a)
1617 {
1618  /* Fit for HOMs missing */
1619  return 0.27581190323955274 + 0.19347381066059993 * eta
1620  - 0.08898338208573725 * log (1.0 - a / (1.0 - 2.0 * eta))
1621  +
1622  eta * eta * (1.78832 * (0.2690779744133912 + a / (2.0 - 4.0 * eta)) *
1623  (1.2056469070395925 + a / (2.0 - 4.0 * eta)) +
1624  1.423734113371796 * log (1.0 - a / (1.0 - 2.0 * eta)));
1625 }
1626 
1627 /**
1628  * Peak frequency slope predicted by fitting NR results (currently only 2,2 available).
1629  * Tables IV and V and Eq. 42 of Taracchini et al. PRD 86, 024011 (2012) [arXiv:1202.0790].
1630  */
1631 UNUSED static inline REAL8
1632 XLALSimIMREOBGetNRSpinPeakOmegaDot (INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta,
1633  REAL8 UNUSED a)
1634 {
1635  /* Fit for HOMs missing */
1636  return 0.006075014646800278 + 0.012040017219351778 * eta
1637  + (0.0007353536801336875 +
1638  0.0015592659912461832 * a / (1.0 - 2.0 * eta)) * log (1.0 - a / (1.0 -
1639  2.0 *
1640  eta))
1641  + eta * eta * (0.03575969677378844 +
1642  (-0.011765658882139 -
1643  0.02494825585993893 * a / (1.0 - 2.0 * eta)) * log (1.0 -
1644  a /
1645  (1.0 -
1646  2.0 *
1647  eta)));
1648 }
1649 
1650 /**
1651  * The time difference between the orbital peak and the peak amplitude
1652  * of the mode in question (currently only 2,2 implemented ).
1653  */
1654 UNUSED static inline REAL8
1655 XLALSimIMREOBGetNRSpinPeakDeltaTv2 (INT4 UNUSED l,
1656  /**<< Mode l */
1657  INT4 UNUSED m,
1658  /**<< Mode m */
1659  REAL8 UNUSED m1,
1660  /**<< mass 1 */
1661  REAL8 UNUSED m2,
1662  /**<< mass 22 */
1663  REAL8 chi1,
1664  /**<< Dimensionless spin1 */
1665  REAL8 chi2
1666  /**<< Dimensionless spin2 */
1667  )
1668 {
1669  REAL8 chi, chichi;
1670  REAL8 eta = m1 * m2 / ((m1 + m2) * (m1 + m2));
1671  chi =
1672  (chi1 + chi2) / 2. + (chi1 - chi2) / 2. * ((m1 - m2) / (m1 + m2)) / (1 -
1673  2. *
1674  eta);
1675 
1676  chichi = chi * chi;
1677  if (chi > 0.8)
1678  {
1679  return (0.75 * eta * chi + sqrt (1. - 4. * eta)) * (57.1755 -
1680  48.0564 * chi);
1681  }
1682  else if (chi > 0.0)
1683  {
1684  return (0.75 * eta * chi + sqrt (1. - 4. * eta)) * (2.5 + 10. * chichi +
1685  24. * chichi *
1686  chichi);
1687  }
1688  else
1689  {
1690  return 2.5 + (1. + 2.5 * chi) * (-2.5 + 2.5 * sqrt (1. - 4. * eta));
1691  }
1692 }
1693 
1694 /**
1695  * The time difference between the orbital peak and the peak amplitude
1696  * of the mode in question (currently only 2,2 implemented ).
1697  */
1698 UNUSED static inline REAL8
1699 XLALSimIMREOBGetNRSpinPeakDeltaTv4 (INT4 UNUSED l, /**<< Mode l */
1700  INT4 UNUSED m, /**<< Mode m */
1701  REAL8 UNUSED m1, /**<< mass 1 */
1702  REAL8 UNUSED m2, /**<< mass 2 */
1703  REAL8 UNUSED chi1, /**<< Dimensionless spin1 */
1704  REAL8 UNUSED chi2 /**<< Dimensionless spin2 */
1705  )
1706 {
1707  REAL8 eta = m1 * m2 / (m1 + m2) / (m1 + m2);
1708  REAL8 chi =
1709  0.5 * (chi1 + chi2) + 0.5 * (chi1 - chi2) * (m1 - m2) / (m1 + m2) / (1. -
1710  2. *
1711  eta);
1712  REAL8 eta2 = eta * eta, eta3 = eta2 * eta;
1713  REAL8 chiTo2 = chi * chi, chiTo3 = chiTo2 * chi;
1714  REAL8 coeff00, coeff01, coeff02, coeff03;
1715  REAL8 coeff10, coeff11, coeff12, coeff13;
1716  REAL8 coeff20, coeff21, coeff22, coeff23;
1717  REAL8 coeff30, coeff31, coeff32, coeff33;
1718  REAL8 res;
1719  // Calibrationv21_Sep8a
1720  coeff00 = 2.50499;
1721  coeff01 = 13.0064;
1722  coeff02 = 11.5435;
1723  coeff03 = 0;
1724  coeff10 = 45.8838;
1725  coeff11 = -40.3183;
1726  coeff12 = 0;
1727  coeff13 = -19.0538;
1728  coeff20 = 13.0879;
1729  coeff21 = 0;
1730  coeff22 = 0;
1731  coeff23 = 0.192775;
1732  coeff30 = -716.044;
1733  coeff31 = 0;
1734  coeff32 = 0;
1735  coeff33 = 0;
1736  res = coeff00 + coeff01 * chi + coeff02 * chiTo2 + coeff03 * chiTo3 +
1737  coeff10 * eta + coeff11 * eta * chi + coeff12 * eta * chiTo2 +
1738  coeff13 * eta * chiTo3 + coeff20 * eta2 + coeff21 * eta2 * chi +
1739  coeff22 * eta2 * chiTo2 + coeff23 * eta2 * chiTo3 + coeff30 * eta3 +
1740  coeff31 * eta3 * chi + coeff32 * eta3 * chiTo2 + coeff33 * eta3 * chiTo3;
1741 
1742  //RC: for the 55 mode the attachment is done at tpeak22 -10M, note that since here deltat22 is defined as -deltat22 with respect
1743  //to SEOBNRv4 paper, here I need to add 10 M
1744 if((l == 5)&&(m == 5)){
1745  res = res + 10;
1746 }
1747 
1748  return res;
1749 }
1750 
1751 /**
1752  * Peak frequency predicted by fitting NR results (currently only 2,2 available).
1753  * Take from unpublished SEOBNRv2 results.
1754  */
1755 UNUSED static inline REAL8
1756 XLALSimIMREOBGetNRSpinPeakOmegav2 (INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 a)
1757 {
1758  REAL8 chi = a / (1.0 - 2.0 * eta);
1759  REAL8 eta2 = eta * eta;
1760  if (eta > 50. / 51. / 51.)
1761  {
1762  return 0.43747541927878864 + (-0.10933208665273314 -
1763  0.007325831113333813 * chi) *
1764  log (4.500844771420863 - 9.681916048928946 * eta +
1765  chi * (-4.254886879579986 + 11.513558950322647 * eta));
1766  }
1767  else
1768  {
1769  return 1.5609526077704716 - 122.25721149839733 * eta +
1770  3586.2192688666914 * eta2 - 13869.506144441548 * eta * eta2 +
1771  (eta - 0.25) * (1651.5823693445805 * (-0.01922337588094282 + eta) *
1772  (-0.01922337536857659 + eta) +
1773  66.87492814925524 * chi * (0.0003695381704106058 -
1774  0.03844675124951941 * eta +
1775  eta2)) *
1776  log (5600.67382718678 - 5555.824895398546 * chi) +
1777  (-1412.8186461833657 + 67.66455403259023 * chi) * (eta -
1778  0.001) *
1779  (0.0003695381704106056 - 0.038446751249519406 * eta +
1780  eta2) * log (0.5680439481719505 - 0.36813967358200156 * chi) +
1781  0.012328326527732041 * log (4.500844771420863 -
1782  9.681916048928946 * eta +
1783  chi * (-4.254886879579986 +
1784  11.513558950322647 * eta)) +
1785  0.0008260634258180991 * chi * log (4.500844771420863 -
1786  9.681916048928946 * eta +
1787  chi * (-4.254886879579986 +
1788  11.513558950322647 * eta)) -
1789  12.6575493872956 * eta * log (4.500844771420863 -
1790  9.681916048928946 * eta +
1791  chi * (-4.254886879579986 +
1792  11.513558950322647 * eta)) -
1793  0.8481231078533651 * chi * eta * log (4.500844771420863 -
1794  9.681916048928946 * eta +
1795  chi * (-4.254886879579986 +
1796  11.513558950322647 *
1797  eta)) +
1798  329.2228595635586 * eta2 * log (4.500844771420863 -
1799  9.681916048928946 * eta +
1800  chi * (-4.254886879579986 +
1801  11.513558950322647 * eta)) +
1802  22.05968203526603 * chi * eta2 * log (4.500844771420863 -
1803  9.681916048928946 * eta +
1804  chi * (-4.254886879579986 +
1805  11.513558950322647 *
1806  eta));
1807  }
1808 }
1809 
1810 /**
1811  * Peak frequency slope predicted by fitting NR results (currently only 2,2 available).
1812  * Take from unpublished SEOBNRv2 results.
1813  */
1814 UNUSED static inline REAL8
1815 XLALSimIMREOBGetNRSpinPeakOmegaDotv2 (INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta,
1816  REAL8 UNUSED a)
1817 {
1818  REAL8 chi = a / (1.0 - 2.0 * eta);
1819  REAL8 eta2 = eta * eta;
1820  /* Fit for HOMs missing */
1821  if (chi < 0.8)
1822  {
1823  return -0.07086074186161867 * chi * (-0.26367236731979804 + eta) *
1824  (-0.0010019969893089581 + eta) + 0.2893863668183948 *
1825  (-0.16845695144529893 + eta) * (0.23032241797163952 + eta) +
1826  (0.004086861548547749 - 0.06538978477676398 * eta2 +
1827  chi * (0.0006334026884930817 - 0.010134443015889307 * eta2)) *
1828  log (68.47466578101876 - 58.30148755701496 * chi);
1829  }
1830  else
1831  {
1832  return -0.10069512275335238 * (-0.46107388514323044 + eta) *
1833  (0.2832795481380979 + eta) + 0.2614619716504706 * chi *
1834  (-0.24838163750494138 + eta) * (0.320112993649413 + eta) +
1835  chi * chi * (0.010000160002560042 - 0.16000256004096067 * eta2);
1836  }
1837 }
1838 
1839 
1840 /**
1841  * This function computes the coefficients a3s, a4, etc. used in the
1842  * non-quasicircular correction. The details of the calculation of these
1843  * coefficients are found in the DCC document T1100433.
1844  * In brief, this function populates and solves the linear equations
1845  * Eq. 18 (for amplitude) and Eq. 19 (for phase) of the DCC document T1100433v2.
1846  */
1847 UNUSED static int
1848 XLALSimIMRSpinEOBCalculateNQCCoefficients (REAL8Vector * restrict amplitude,
1849  /**<< Waveform amplitude, func of time */
1850  REAL8Vector * restrict phase,
1851  /**<< Waveform phase(rad), func of time */
1852  REAL8Vector * restrict rVec,
1853  /**<< Position-vector, function of time */
1854  REAL8Vector * restrict prVec,
1855  /**<< Momentum vector, function of time */
1856  REAL8Vector * restrict orbOmegaVec,
1857  /**<< Orbital frequency, func of time */
1858  INT4 l, /**<< Mode index l */
1859  INT4 m, /**<< Mode index m */
1860  REAL8 timePeak,
1861  /**<< Time of peak orbital frequency */
1862  REAL8 deltaT,/**<< Sampling interval */
1863  REAL8 m1, /**<< Component mass 1 */
1864  REAL8 m2, /**<< Component mass 2 */
1865  REAL8 a, /**<< Normalized spin of deformed-Kerr */
1866  REAL8 chiA, /**<< Assymmetric dimensionless spin combination */
1867  REAL8 chiS, /**<< Symmetric dimensionless spin combination */
1868  EOBNonQCCoeffs * restrict coeffs,
1869  /**<< OUTPUT, NQC coefficients */
1870  UINT4 SpinAlignedEOBversion
1871  /**<< 1 for SEOBNRv1, 2 for SEOBNRv2 */
1872  )
1873 {
1874 
1875  /* For gsl permutation stuff */
1876 
1877  int signum;
1878 
1879  REAL8Vector *restrict timeVec = NULL;
1880 
1881  /* Vectors which are used in the computation of the NQC coefficients */
1882  REAL8Vector *q3 = NULL, *q4 = NULL, *q5 = NULL;
1883  REAL8Vector *p3 = NULL, *p4 = NULL;
1884 
1885  REAL8Vector *qNS = NULL, *pNS = NULL;
1886 
1887  /* Since the vectors we actually want are q etc * A, we will have to generate them here */
1888  REAL8Vector *q3LM = NULL;
1889  REAL8Vector *q4LM = NULL;
1890  REAL8Vector *q5LM = NULL;
1891  REAL8Vector *qNSLM = NULL;
1892 
1893  REAL8 eta = (m1 * m2) / ((m1 + m2) * (m1 + m2));
1894  REAL8 amp, aDot, aDDot;
1895  REAL8 omega, omegaDot;
1896 
1897  REAL8 qNSLMPeak, qNSLMDot, qNSLMDDot;
1898  REAL8 pNSLMDot, pNSLMDDot;
1899 
1900  REAL8 nra, nraDDot;
1901  REAL8 nromega, nromegaDot;
1902 
1903  REAL8 nrDeltaT, nrTimePeak;
1904  REAL8 chi1 = chiS + chiA;
1905  REAL8 chi2 = chiS - chiA;
1906 
1907  /* Stuff for finding numerical derivatives */
1908  gsl_spline *spline = NULL;
1909  gsl_interp_accel *acc = NULL;
1910 
1911  /* Matrix stuff for calculating coefficients */
1912  gsl_matrix *qMatrix = NULL, *pMatrix = NULL;
1913  gsl_vector *aCoeff = NULL, *bCoeff = NULL;
1914 
1915  gsl_vector *amps = NULL, *omegaVec = NULL;
1916 
1917  gsl_permutation *perm1 = NULL, *perm2 = NULL;
1918 
1919  memset (coeffs, 0, sizeof (EOBNonQCCoeffs));
1920 
1921  /* Populate the time vector */
1922  /* It is okay to assume initial t = 0 */
1923  timeVec = XLALCreateREAL8Vector (rVec->length);
1924  q3 = XLALCreateREAL8Vector (rVec->length);
1925  q4 = XLALCreateREAL8Vector (rVec->length);
1926  q5 = XLALCreateREAL8Vector (rVec->length);
1927  p3 = XLALCreateREAL8Vector (rVec->length);
1928  p4 = XLALCreateREAL8Vector (rVec->length);
1929  qNS = XLALCreateREAL8Vector (rVec->length);
1930  pNS = XLALCreateREAL8Vector (rVec->length);
1931  q3LM = XLALCreateREAL8Vector (rVec->length);
1932  q4LM = XLALCreateREAL8Vector (rVec->length);
1933  q5LM = XLALCreateREAL8Vector (rVec->length);
1934  qNSLM = XLALCreateREAL8Vector (rVec->length);
1935 
1936  if (!timeVec || !q3 || !q4 || !q5 || !p3 || !p4 || !qNS || !pNS || !q3LM
1937  || !q4LM || !q5LM || !qNSLM)
1938  {
1939  XLALDestroyREAL8Vector (timeVec);
1940  XLALDestroyREAL8Vector (q3);
1941  XLALDestroyREAL8Vector (q4);
1942  XLALDestroyREAL8Vector (q5);
1943  XLALDestroyREAL8Vector (p3);
1944  XLALDestroyREAL8Vector (p4);
1945  XLALDestroyREAL8Vector (qNS);
1946  XLALDestroyREAL8Vector (pNS);
1947  XLALDestroyREAL8Vector (q3LM);
1948  XLALDestroyREAL8Vector (q4LM);
1949  XLALDestroyREAL8Vector (q5LM);
1950  XLALDestroyREAL8Vector (qNSLM);
1951  XLAL_ERROR (XLAL_EFUNC);
1952  }
1953 
1954  /* We need the calibrated non-spinning NQC coefficients */
1955  switch (SpinAlignedEOBversion)
1956  {
1957  case 1:
1958  if (XLALSimIMRGetEOBCalibratedSpinNQC (coeffs, l, m, eta, a) ==
1959  XLAL_FAILURE)
1960  {
1961  XLALDestroyREAL8Vector (timeVec);
1962  XLALDestroyREAL8Vector (q3);
1963  XLALDestroyREAL8Vector (q4);
1964  XLALDestroyREAL8Vector (q5);
1965  XLALDestroyREAL8Vector (p3);
1966  XLALDestroyREAL8Vector (p4);
1967  XLALDestroyREAL8Vector (qNS);
1968  XLALDestroyREAL8Vector (pNS);
1969  XLALDestroyREAL8Vector (q3LM);
1970  XLALDestroyREAL8Vector (q4LM);
1971  XLALDestroyREAL8Vector (q5LM);
1972  XLALDestroyREAL8Vector (qNSLM);
1973  XLAL_ERROR (XLAL_EFUNC);
1974  }
1975  break;
1976  case 2:
1977  // if ( XLALSimIMRGetEOBCalibratedSpinNQCv2( coeffs, l, m, eta, a ) == XLAL_FAILURE )
1978  if (XLALSimIMRGetEOBCalibratedSpinNQC3D (coeffs, l, m, m1, m2, a, chiA)
1979  == XLAL_FAILURE)
1980  {
1981  XLALDestroyREAL8Vector (timeVec);
1982  XLALDestroyREAL8Vector (q3);
1983  XLALDestroyREAL8Vector (q4);
1984  XLALDestroyREAL8Vector (q5);
1985  XLALDestroyREAL8Vector (p3);
1986  XLALDestroyREAL8Vector (p4);
1987  XLALDestroyREAL8Vector (qNS);
1988  XLALDestroyREAL8Vector (pNS);
1989  XLALDestroyREAL8Vector (q3LM);
1990  XLALDestroyREAL8Vector (q4LM);
1991  XLALDestroyREAL8Vector (q5LM);
1992  XLALDestroyREAL8Vector (qNSLM);
1993  XLAL_ERROR (XLAL_EFUNC);
1994  }
1995  break;
1996  default:
1997  XLALPrintError
1998  ("XLAL Error - %s: Unknown SEOBNR version!\nAt present only v1 and v2 are available.\n",
1999  __func__);
2000  XLAL_ERROR (XLAL_EINVAL);
2001  break;
2002  }
2003 
2004  /* Populate vectors as necessary. Eqs. 14 - 17 of the LIGO DCC document T1100433v2 */
2005  for (unsigned int i = 0; i < timeVec->length; i++)
2006  {
2007 
2008  REAL8 rootR = sqrt (rVec->data[i]);
2009  REAL8 rOmega = rVec->data[i] * orbOmegaVec->data[i];
2010 
2011  /* We don't need these as vectors as their coefficients are calibrated */
2012  REAL8 q1, q2, p1, p2;
2013 
2014  timeVec->data[i] = i * deltaT;
2015  q1 = prVec->data[i] * prVec->data[i] / (rOmega * rOmega);
2016  q2 = q1 / rVec->data[i];
2017  q3->data[i] = q2 / rootR;
2018  q4->data[i] = q2 / rVec->data[i];
2019  q5->data[i] = q3->data[i] / rVec->data[i];
2020 
2021  p1 = prVec->data[i] / rOmega;
2022  p2 = p1 * prVec->data[i] * prVec->data[i];
2023  p3->data[i] = p2 / rootR;
2024  p4->data[i] = p2 / rVec->data[i];
2025 
2026  qNS->data[i] =
2027  coeffs->a1 * q1 + coeffs->a2 * q2 + coeffs->a3 * q3->data[i];
2028  pNS->data[i] = coeffs->b1 * p1 + coeffs->b2 * p2;
2029  q3LM->data[i] = q3->data[i] * amplitude->data[i];
2030  q4LM->data[i] = q4->data[i] * amplitude->data[i];
2031  q5LM->data[i] = q5->data[i] * amplitude->data[i];
2032 
2033  qNSLM->data[i] = qNS->data[i] * amplitude->data[i];
2034  }
2035  /* Allocate all the memory we need */
2036  XLAL_CALLGSL (
2037  /* a stuff */
2038  qMatrix = gsl_matrix_alloc (3, 3);
2039  aCoeff = gsl_vector_alloc (3);
2040  amps = gsl_vector_alloc (3);
2041  perm1 = gsl_permutation_alloc (3);
2042  /* b stuff */
2043  pMatrix = gsl_matrix_alloc (2, 2);
2044  bCoeff = gsl_vector_alloc (2);
2045  omegaVec = gsl_vector_alloc (2);
2046  perm2 = gsl_permutation_alloc (2););
2047 
2048  if (!qMatrix || !aCoeff || !amps || !pMatrix || !bCoeff || !omegaVec)
2049  {
2050  XLALDestroyREAL8Vector (timeVec);
2051  XLALDestroyREAL8Vector (q3);
2052  XLALDestroyREAL8Vector (q4);
2053  XLALDestroyREAL8Vector (q5);
2054  XLALDestroyREAL8Vector (p3);
2055  XLALDestroyREAL8Vector (p4);
2056  XLALDestroyREAL8Vector (qNS);
2057  XLALDestroyREAL8Vector (pNS);
2058  XLALDestroyREAL8Vector (q3LM);
2059  XLALDestroyREAL8Vector (q4LM);
2060  XLALDestroyREAL8Vector (q5LM);
2061  XLALDestroyREAL8Vector (qNSLM);
2062  XLAL_ERROR (XLAL_ENOMEM);
2063  }
2064 
2065  /* The time we want to take as the peak time depends on l and m */
2066  /* Calculate the adjustment we need to make here */
2067  switch (SpinAlignedEOBversion)
2068  {
2069  case 1:
2070  nrDeltaT = XLALSimIMREOBGetNRSpinPeakDeltaT (l, m, eta, a);
2071  break;
2072  case 2:
2073  nrDeltaT =
2074  XLALSimIMREOBGetNRSpinPeakDeltaTv2 (l, m, m1, m2, chi1, chi2);
2075  break;
2076  default:
2077  XLALPrintError
2078  ("XLAL Error - %s: Unknown SEOBNR version!\nAt present only v1 and v2 are available.\n",
2079  __func__);
2080  XLAL_ERROR (XLAL_EINVAL);
2081  break;
2082  }
2083 
2084 
2085  if (XLAL_IS_REAL8_FAIL_NAN (nrDeltaT))
2086  {
2087  XLALDestroyREAL8Vector (timeVec);
2088  XLALDestroyREAL8Vector (q3);
2089  XLALDestroyREAL8Vector (q4);
2090  XLALDestroyREAL8Vector (q5);
2091  XLALDestroyREAL8Vector (p3);
2092  XLALDestroyREAL8Vector (p4);
2093  XLALDestroyREAL8Vector (qNS);
2094  XLALDestroyREAL8Vector (pNS);
2095  XLALDestroyREAL8Vector (q3LM);
2096  XLALDestroyREAL8Vector (q4LM);
2097  XLALDestroyREAL8Vector (q5LM);
2098  XLALDestroyREAL8Vector (qNSLM);
2099  XLAL_ERROR (XLAL_EFUNC);
2100  }
2101 
2102  /* nrDeltaT defined in XLALSimIMREOBGetNRSpinPeakDeltaT is a minus sign different from Eq. (33) of Taracchini et al.
2103  * Therefore, the plus sign in Eq. (21) of Taracchini et al and Eq. (18) of DCC document T1100433v2 is
2104  * changed to a minus sign here.
2105  */
2106  nrTimePeak = timePeak - nrDeltaT;
2107 
2108  /* We are now in a position to use the interp stuff to calculate the derivatives we need */
2109  /* We will start with the quantities used in the calculation of the a coefficients */
2110  spline = gsl_spline_alloc (gsl_interp_cspline, amplitude->length);
2111  acc = gsl_interp_accel_alloc ();
2112 
2113  /* Populate the Q matrix in Eq. 18 of the LIGO DCC document T1100433v2 */
2114  /* Q3 */
2115  gsl_spline_init (spline, timeVec->data, q3LM->data, q3LM->length);
2116  gsl_matrix_set (qMatrix, 0, 0, gsl_spline_eval (spline, nrTimePeak, acc));
2117  gsl_matrix_set (qMatrix, 1, 0,
2118  gsl_spline_eval_deriv (spline, nrTimePeak, acc));
2119  gsl_matrix_set (qMatrix, 2, 0,
2120  gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
2121 
2122  /* Q4 */
2123  gsl_spline_init (spline, timeVec->data, q4LM->data, q4LM->length);
2124  gsl_interp_accel_reset (acc);
2125  gsl_matrix_set (qMatrix, 0, 1, gsl_spline_eval (spline, nrTimePeak, acc));
2126  gsl_matrix_set (qMatrix, 1, 1,
2127  gsl_spline_eval_deriv (spline, nrTimePeak, acc));
2128  gsl_matrix_set (qMatrix, 2, 1,
2129  gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
2130 
2131  /* Q5 */
2132  gsl_spline_init (spline, timeVec->data, q5LM->data, q5LM->length);
2133  gsl_interp_accel_reset (acc);
2134  gsl_matrix_set (qMatrix, 0, 2, gsl_spline_eval (spline, nrTimePeak, acc));
2135  gsl_matrix_set (qMatrix, 1, 2,
2136  gsl_spline_eval_deriv (spline, nrTimePeak, acc));
2137  gsl_matrix_set (qMatrix, 2, 2,
2138  gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
2139 
2140  /* Populate the r.h.s vector of Eq. 18 of the LIGO DCC document T1100433v2 */
2141  /* Amplitude */
2142  gsl_spline_init (spline, timeVec->data, amplitude->data, amplitude->length);
2143  gsl_interp_accel_reset (acc);
2144  amp = gsl_spline_eval (spline, nrTimePeak, acc);
2145  aDot = gsl_spline_eval_deriv (spline, nrTimePeak, acc);
2146  aDDot = gsl_spline_eval_deriv2 (spline, nrTimePeak, acc);
2147 
2148  /* qNSLM */
2149  gsl_spline_init (spline, timeVec->data, qNSLM->data, qNSLM->length);
2150  gsl_interp_accel_reset (acc);
2151  qNSLMPeak = gsl_spline_eval (spline, nrTimePeak, acc);
2152  qNSLMDot = gsl_spline_eval_deriv (spline, nrTimePeak, acc);
2153  qNSLMDDot = gsl_spline_eval_deriv2 (spline, nrTimePeak, acc);
2154 
2155  nra = XLALSimIMREOBGetNRSpinPeakAmplitude (l, m, eta, a);
2156  nraDDot = -XLALSimIMREOBGetNRSpinPeakADDot (l, m, eta, a);
2157 
2158  if (XLAL_IS_REAL8_FAIL_NAN (nra) || XLAL_IS_REAL8_FAIL_NAN (nraDDot))
2159  {
2160  XLALDestroyREAL8Vector (timeVec);
2161  XLALDestroyREAL8Vector (q3);
2162  XLALDestroyREAL8Vector (q4);
2163  XLALDestroyREAL8Vector (q5);
2164  XLALDestroyREAL8Vector (p3);
2165  XLALDestroyREAL8Vector (p4);
2166  XLALDestroyREAL8Vector (qNS);
2167  XLALDestroyREAL8Vector (pNS);
2168  XLALDestroyREAL8Vector (q3LM);
2169  XLALDestroyREAL8Vector (q4LM);
2170  XLALDestroyREAL8Vector (q5LM);
2171  XLALDestroyREAL8Vector (qNSLM);
2172  XLAL_ERROR (XLAL_EFUNC);
2173  }
2174 
2175  gsl_vector_set (amps, 0, nra - amp - qNSLMPeak);
2176  gsl_vector_set (amps, 1, -aDot - qNSLMDot);
2177  gsl_vector_set (amps, 2, nraDDot - aDDot - qNSLMDDot);
2178 
2179  /* We have now set up all the stuff to calculate the a coefficients */
2180  /* So let us do it! */
2181  gsl_linalg_LU_decomp (qMatrix, perm1, &signum);
2182  gsl_linalg_LU_solve (qMatrix, perm1, amps, aCoeff);
2183 
2184  /* Now we (should) have calculated the a values. Now we can do the b values */
2185 
2186  /* Populate the P matrix in Eq. 18 of the LIGO DCC document T1100433v2 */
2187  /* P3 */
2188  gsl_spline_init (spline, timeVec->data, p3->data, p3->length);
2189  gsl_interp_accel_reset (acc);
2190  gsl_matrix_set (pMatrix, 0, 0,
2191  -gsl_spline_eval_deriv (spline, nrTimePeak, acc));
2192  gsl_matrix_set (pMatrix, 1, 0,
2193  -gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
2194 
2195  /* P4 */
2196  gsl_spline_init (spline, timeVec->data, p4->data, p4->length);
2197  gsl_interp_accel_reset (acc);
2198  gsl_matrix_set (pMatrix, 0, 1,
2199  -gsl_spline_eval_deriv (spline, nrTimePeak, acc));
2200  gsl_matrix_set (pMatrix, 1, 1,
2201  -gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
2202 
2203  /* Populate the r.h.s vector of Eq. 18 of the LIGO DCC document T1100433v2 */
2204  /* Phase */
2205  gsl_spline_init (spline, timeVec->data, phase->data, phase->length);
2206  gsl_interp_accel_reset (acc);
2207  omega = gsl_spline_eval_deriv (spline, nrTimePeak, acc);
2208  omegaDot = gsl_spline_eval_deriv2 (spline, nrTimePeak, acc);
2209 
2210  /* pNSLM */
2211  gsl_spline_init (spline, timeVec->data, pNS->data, pNS->length);
2212  gsl_interp_accel_reset (acc);
2213  pNSLMDot = gsl_spline_eval_deriv (spline, nrTimePeak, acc);
2214  pNSLMDDot = gsl_spline_eval_deriv2 (spline, nrTimePeak, acc);
2215 
2216  /* Since the phase can be decreasing, we need to take care not to have a -ve frequency */
2217  if (omega * omegaDot > 0.0)
2218  {
2219  omega = fabs (omega);
2220  omegaDot = fabs (omegaDot);
2221  }
2222  else
2223  {
2224  omega = fabs (omega);
2225  omegaDot = -fabs (omegaDot);
2226  }
2227 
2228  //nromega = GetNRPeakOmega( l, m, eta );
2229  //nromegaDot = GetNRPeakOmegaDot( l, m, eta );
2230  switch (SpinAlignedEOBversion)
2231  {
2232  case 1:
2233  nromega = XLALSimIMREOBGetNRSpinPeakOmega (l, m, eta, a);
2234  nromegaDot = XLALSimIMREOBGetNRSpinPeakOmegaDot (l, m, eta, a);
2235  break;
2236  case 2:
2237  nromega = XLALSimIMREOBGetNRSpinPeakOmegav2 (l, m, eta, a);
2238  nromegaDot = XLALSimIMREOBGetNRSpinPeakOmegaDotv2 (l, m, eta, a);
2239  break;
2240  default:
2241  XLALPrintError
2242  ("XLAL Error - %s: Unknown SEOBNR version!\nAt present only v1 and v2 are available.\n",
2243  __func__);
2244  XLAL_ERROR (XLAL_EINVAL);
2245  break;
2246  }
2247 
2248  /*printf("NR inputs: %.16e, %.16e, %.16e, %.16e\n",nra,nraDDot,nromega,nromegaDot);
2249  printf("NR inputs: %.16e, %.16e, %.16e, %.16e\n",pNSLMDot, pNSLMDDot,omega,omegaDot); */
2250 
2251  if (XLAL_IS_REAL8_FAIL_NAN (nromega) || XLAL_IS_REAL8_FAIL_NAN (nromegaDot))
2252  {
2253  XLALDestroyREAL8Vector (timeVec);
2254  XLALDestroyREAL8Vector (q3);
2255  XLALDestroyREAL8Vector (q4);
2256  XLALDestroyREAL8Vector (q5);
2257  XLALDestroyREAL8Vector (p3);
2258  XLALDestroyREAL8Vector (p4);
2259  XLALDestroyREAL8Vector (qNS);
2260  XLALDestroyREAL8Vector (pNS);
2261  XLALDestroyREAL8Vector (q3LM);
2262  XLALDestroyREAL8Vector (q4LM);
2263  XLALDestroyREAL8Vector (q5LM);
2264  XLALDestroyREAL8Vector (qNSLM);
2265  XLAL_ERROR (XLAL_EFUNC);
2266  }
2267 
2268  gsl_vector_set (omegaVec, 0, nromega - omega + pNSLMDot);
2269  gsl_vector_set (omegaVec, 1, nromegaDot - omegaDot + pNSLMDDot);
2270 
2271  /*printf( "P MATRIX\n" );
2272  for (unsigned int i = 0; i < 2; i++ )
2273  {
2274  for (unsigned int j = 0; j < 2; j++ )
2275  {
2276  printf( "%.12e\t", gsl_matrix_get( pMatrix, i, j ));
2277  }
2278  printf( "= %.12e\n", gsl_vector_get( omegaVec, i ) );
2279  } */
2280 
2281  /* And now solve for the b coefficients */
2282  gsl_linalg_LU_decomp (pMatrix, perm2, &signum);
2283  gsl_linalg_LU_solve (pMatrix, perm2, omegaVec, bCoeff);
2284 
2285  /* We can now populate the coefficients structure */
2286 /* coeffs->a3S = gsl_vector_get( aCoeff, 0 );
2287  coeffs->a4 = gsl_vector_get( aCoeff, 1 );
2288  coeffs->a5 = gsl_vector_get( aCoeff, 2 );*/
2289  switch (SpinAlignedEOBversion)
2290  {
2291  case 1:
2292  coeffs->b3 = gsl_vector_get (bCoeff, 0);
2293  coeffs->b4 = gsl_vector_get (bCoeff, 1);
2294  break;
2295  case 2:
2296  //coeffs->b3 = gsl_vector_get( bCoeff, 0 );
2297  //coeffs->b4 = gsl_vector_get( bCoeff, 1 );
2298  break;
2299  default:
2300  XLALPrintError
2301  ("XLAL Error - %s: Unknown SEOBNR version!\nAt present only v1 and v2 are available.\n",
2302  __func__);
2303  XLAL_ERROR (XLAL_EINVAL);
2304  break;
2305  }
2306  coeffs->b3 *= 1.0;
2307  coeffs->b4 *= 1.0;
2308 // coeffs->b3 = -778.891568845;
2309 // coeffs->b4 = 1237.46952422;
2310 // coeffs->b3 = -876.669217307;
2311 // coeffs->b4 = 1386.13223658;
2312 // coeffs->b3 = 41583.9402122;
2313 // coeffs->b4 = 68359.70064;
2314 
2315  //printf ("NQC coefficients:\n");
2316  //printf ("{%f,%f,%f,%f,%f,%f}\n", coeffs->a1, coeffs->a2, coeffs->a3,
2317  // coeffs->a3S, coeffs->a4, coeffs->a5);
2318 
2319  //printf ("{%f,%f,%f,%f}\n", coeffs->b1, coeffs->b2, coeffs->b3, coeffs->b4);
2320 
2321  /* Free memory and exit */
2322  gsl_matrix_free (qMatrix);
2323  gsl_vector_free (amps);
2324  gsl_vector_free (aCoeff);
2325  gsl_permutation_free (perm1);
2326 
2327  gsl_matrix_free (pMatrix);
2328  gsl_vector_free (omegaVec);
2329  gsl_vector_free (bCoeff);
2330  gsl_permutation_free (perm2);
2331 
2332  gsl_spline_free (spline);
2333  gsl_interp_accel_free (acc);
2334 
2335  XLALDestroyREAL8Vector (timeVec);
2336  XLALDestroyREAL8Vector (q3);
2337  XLALDestroyREAL8Vector (q4);
2338  XLALDestroyREAL8Vector (q5);
2339  XLALDestroyREAL8Vector (p3);
2340  XLALDestroyREAL8Vector (p4);
2341  XLALDestroyREAL8Vector (qNS);
2342  XLALDestroyREAL8Vector (pNS);
2343  XLALDestroyREAL8Vector (q3LM);
2344  XLALDestroyREAL8Vector (q4LM);
2345  XLALDestroyREAL8Vector (q5LM);
2346  XLALDestroyREAL8Vector (qNSLM);
2347 
2348  return XLAL_SUCCESS;
2349 }
2350 
2351 /**
2352  * This function computes the coefficients a3s, a4, etc. used in the
2353  * non-quasicircular correction. The details of the calculation of these
2354  * coefficients are found in the DCC document T1100433.
2355  * In brief, this function populates and solves the linear equations
2356  * Eq. 18 (for amplitude) and Eq. 19 (for phase) of https://dcc.ligo.org/T1100433
2357  */
2358 UNUSED static int
2359 XLALSimIMRSpinEOBCalculateNQCCoefficientsV4 (REAL8Vector * restrict amplitude, /**<< Waveform amplitude, func of time */
2360  REAL8Vector * restrict phase, /**<< Waveform phase(rad), func of time */
2361  REAL8Vector * restrict rVec, /**<< Position-vector, function of time */
2362  REAL8Vector * restrict prVec, /**<< Momentum vector, function of time */
2363  REAL8Vector * restrict orbOmegaVec, /**<< Orbital frequency, func of time */
2364  INT4 modeL, /**<< Mode index l */
2365  INT4 modeM, /**<< Mode index m */
2366  REAL8 timePeak, /**<< Time of peak orbital frequency */
2367  REAL8 deltaT, /**<< Sampling interval */
2368  REAL8 m1, /**<< Component mass 1 */
2369  REAL8 m2, /**<< Component mass 2 */
2370  REAL8 chiA, /**<< Assymmetric dimensionless spin combination */
2371  REAL8 chiS, /**<< Symmetric dimensionless spin combination */
2372  EOBNonQCCoeffs * restrict coeffs /**<< OUTPUT, NQC coefficients */
2373  )
2374 {
2375  int debugAT = 0;
2376  /* For gsl permutation stuff */
2377 
2378  int signum;
2379 
2380  REAL8Vector *restrict timeVec = NULL;
2381 
2382  /* Vectors which are used in the computation of the NQC coefficients */
2383  REAL8Vector *q3 = NULL, *q4 = NULL, *q5 = NULL;
2384  REAL8Vector *p3 = NULL, *p4 = NULL;
2385 
2386  REAL8Vector *qNS = NULL, *pNS = NULL;
2387 
2388  /* Since the vectors we actually want are q etc * A, we will have to generate them here */
2389  REAL8Vector *q3LM = NULL;
2390  REAL8Vector *q4LM = NULL;
2391  REAL8Vector *q5LM = NULL;
2392  REAL8Vector *qNSLM = NULL;
2393 
2394  REAL8 eta = (m1 * m2) / ((m1 + m2) * (m1 + m2));
2395  REAL8 amp, aDot, aDDot;
2396  REAL8 omega, omegaDot;
2397 
2398  UNUSED REAL8 qNSLMPeak, qNSLMDot, qNSLMDDot;
2399  UNUSED REAL8 pNSLMDot, pNSLMDDot;
2400 
2401  REAL8 nra, nraDot, nraDDot;
2402  REAL8 nromega, nromegaDot;
2403 
2404  REAL8 nrDeltaT, nrTimePeak;
2405  REAL8 chi1 = chiS + chiA;
2406  REAL8 chi2 = chiS - chiA;
2407 
2408  /* Stuff for finding numerical derivatives */
2409  gsl_spline *spline = NULL;
2410  gsl_interp_accel *acc = NULL;
2411 
2412  /* Matrix stuff for calculating coefficients */
2413  gsl_matrix *qMatrix = NULL, *pMatrix = NULL;
2414  gsl_vector *aCoeff = NULL, *bCoeff = NULL;
2415 
2416  gsl_vector *amps = NULL, *omegaVec = NULL;
2417 
2418  gsl_permutation *perm1 = NULL, *perm2 = NULL;
2419 
2420  memset (coeffs, 0, sizeof (EOBNonQCCoeffs));
2421 
2422  /* Populate the time vector */
2423  /* It is okay to assume initial t = 0 */
2424  timeVec = XLALCreateREAL8Vector (rVec->length);
2425  q3 = XLALCreateREAL8Vector (rVec->length);
2426  q4 = XLALCreateREAL8Vector (rVec->length);
2427  q5 = XLALCreateREAL8Vector (rVec->length);
2428  p3 = XLALCreateREAL8Vector (rVec->length);
2429  p4 = XLALCreateREAL8Vector (rVec->length);
2430  qNS = XLALCreateREAL8Vector (rVec->length);
2431  pNS = XLALCreateREAL8Vector (rVec->length);
2432  q3LM = XLALCreateREAL8Vector (rVec->length);
2433  q4LM = XLALCreateREAL8Vector (rVec->length);
2434  q5LM = XLALCreateREAL8Vector (rVec->length);
2435  qNSLM = XLALCreateREAL8Vector (rVec->length);
2436 
2437  if (!timeVec || !q3 || !q4 || !q5 || !p3 || !p4 || !qNS || !pNS || !q3LM
2438  || !q4LM || !q5LM || !qNSLM)
2439  {
2440  XLALDestroyREAL8Vector (timeVec);
2441  XLALDestroyREAL8Vector (q3);
2442  XLALDestroyREAL8Vector (q4);
2443  XLALDestroyREAL8Vector (q5);
2444  XLALDestroyREAL8Vector (p3);
2445  XLALDestroyREAL8Vector (p4);
2446  XLALDestroyREAL8Vector (qNS);
2447  XLALDestroyREAL8Vector (pNS);
2448  XLALDestroyREAL8Vector (q3LM);
2449  XLALDestroyREAL8Vector (q4LM);
2450  XLALDestroyREAL8Vector (q5LM);
2451  XLALDestroyREAL8Vector (qNSLM);
2452  XLAL_ERROR (XLAL_EFUNC);
2453  }
2454 
2455  /* Populate vectors as necessary. Eqs. 14 - 17 of the LIGO DCC document T1100433v2 */
2456 // FILE *out = fopen( "out.dat","w");
2457  for (unsigned int i = 0; i < timeVec->length; i++)
2458  {
2459 
2460  REAL8 rootR = sqrt (rVec->data[i]);
2461  REAL8 rOmega = rVec->data[i] * orbOmegaVec->data[i];
2462 
2463  /* We don't need these as vectors as their coefficients are calibrated */
2464  REAL8 q1, q2, p1, p2;
2465 
2466  timeVec->data[i] = i * deltaT;
2467  q1 = prVec->data[i] * prVec->data[i] / (rOmega * rOmega);
2468  q2 = q1 / rVec->data[i];
2469  q3->data[i] = q1;
2470  q4->data[i] = q2;
2471  q5->data[i] = q2 / rootR;
2472 
2473  p1 = prVec->data[i] / rOmega;
2474  p2 = p1 * prVec->data[i] * prVec->data[i];
2475  p3->data[i] = p1;
2476  p4->data[i] = p2;
2477 
2478  qNS->data[i] = 0.;
2479  pNS->data[i] = 0.;
2480  q3LM->data[i] = q3->data[i] * amplitude->data[i];
2481  q4LM->data[i] = q4->data[i] * amplitude->data[i];
2482  q5LM->data[i] = q5->data[i] * amplitude->data[i];
2483 
2484 
2485  qNSLM->data[i] = 0.;
2486 // fprintf(out, "%.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n", amplitude->data[i], prVec->data[i],rVec->data[i],orbOmegaVec->data[i], q3->data[i],q4->data[i],q5->data[i], p3->data[i], p4->data[i],phase->data[i]);
2487  }
2488 // fclose(out);
2489 
2490  /* Allocate all the memory we need */
2491  XLAL_CALLGSL (
2492  /* a stuff */
2493  qMatrix = gsl_matrix_alloc (3, 3);
2494  aCoeff = gsl_vector_alloc (3);
2495  amps = gsl_vector_alloc (3);
2496  perm1 = gsl_permutation_alloc (3);
2497  /* b stuff */
2498  pMatrix = gsl_matrix_alloc (2, 2);
2499  bCoeff = gsl_vector_alloc (2);
2500  omegaVec = gsl_vector_alloc (2);
2501  perm2 = gsl_permutation_alloc (2););
2502 
2503  if (!qMatrix || !aCoeff || !amps || !pMatrix || !bCoeff || !omegaVec)
2504  {
2505  XLALDestroyREAL8Vector (timeVec);
2506  XLALDestroyREAL8Vector (q3);
2507  XLALDestroyREAL8Vector (q4);
2508  XLALDestroyREAL8Vector (q5);
2509  XLALDestroyREAL8Vector (p3);
2510  XLALDestroyREAL8Vector (p4);
2511  XLALDestroyREAL8Vector (qNS);
2512  XLALDestroyREAL8Vector (pNS);
2513  XLALDestroyREAL8Vector (q3LM);
2514  XLALDestroyREAL8Vector (q4LM);
2515  XLALDestroyREAL8Vector (q5LM);
2516  XLALDestroyREAL8Vector (qNSLM);
2517  XLAL_ERROR (XLAL_ENOMEM);
2518  }
2519 
2520  /* The time we want to take as the peak time depends on l and m */
2521  /* Calculate the adjustment we need to make here */
2522 
2523  nrDeltaT = XLALSimIMREOBGetNRSpinPeakDeltaTv4 (modeL, modeM, m1, m2, chi1, chi2);
2524 
2525 
2526  if (XLAL_IS_REAL8_FAIL_NAN (nrDeltaT))
2527  {
2528  XLALDestroyREAL8Vector (timeVec);
2529  XLALDestroyREAL8Vector (q3);
2530  XLALDestroyREAL8Vector (q4);
2531  XLALDestroyREAL8Vector (q5);
2532  XLALDestroyREAL8Vector (p3);
2533  XLALDestroyREAL8Vector (p4);
2534  XLALDestroyREAL8Vector (qNS);
2535  XLALDestroyREAL8Vector (pNS);
2536  XLALDestroyREAL8Vector (q3LM);
2537  XLALDestroyREAL8Vector (q4LM);
2538  XLALDestroyREAL8Vector (q5LM);
2539  XLALDestroyREAL8Vector (qNSLM);
2540  XLAL_ERROR (XLAL_EFUNC);
2541  }
2542 
2543  /* nrDeltaT defined in XLALSimIMREOBGetNRSpinPeakDeltaT is a minus sign different from Eq. (33) of Taracchini et al.
2544  * Therefore, the plus sign in Eq. (21) of Taracchini et al and Eq. (18) of DCC document T1100433v2 is
2545  * changed to a minus sign here.
2546  */
2547  nrTimePeak = timePeak - nrDeltaT;
2548  if (debugAT)
2549  printf ("nrTimePeak, timePeak %.16e %.16e\n", nrTimePeak, timePeak);
2550  /* We are now in a position to use the interp stuff to calculate the derivatives we need */
2551  /* We will start with the quantities used in the calculation of the a coefficients */
2552  spline = gsl_spline_alloc (gsl_interp_cspline, amplitude->length);
2553  acc = gsl_interp_accel_alloc ();
2554 
2555  /* Populate the Q matrix in Eq. 18 of the LIGO DCC document T1100433v2 */
2556  /* Q3 */
2557  gsl_spline_init (spline, timeVec->data, q3LM->data, q3LM->length);
2558  gsl_matrix_set (qMatrix, 0, 0, gsl_spline_eval (spline, nrTimePeak, acc));
2559  gsl_matrix_set (qMatrix, 1, 0,
2560  gsl_spline_eval_deriv (spline, nrTimePeak, acc));
2561  gsl_matrix_set (qMatrix, 2, 0,
2562  gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
2563 
2564  /* Q4 */
2565  gsl_spline_init (spline, timeVec->data, q4LM->data, q4LM->length);
2566  gsl_interp_accel_reset (acc);
2567  gsl_matrix_set (qMatrix, 0, 1, gsl_spline_eval (spline, nrTimePeak, acc));
2568  gsl_matrix_set (qMatrix, 1, 1,
2569  gsl_spline_eval_deriv (spline, nrTimePeak, acc));
2570  gsl_matrix_set (qMatrix, 2, 1,
2571  gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
2572 
2573  /* Q5 */
2574  gsl_spline_init (spline, timeVec->data, q5LM->data, q5LM->length);
2575  gsl_interp_accel_reset (acc);
2576  gsl_matrix_set (qMatrix, 0, 2, gsl_spline_eval (spline, nrTimePeak, acc));
2577  gsl_matrix_set (qMatrix, 1, 2,
2578  gsl_spline_eval_deriv (spline, nrTimePeak, acc));
2579  gsl_matrix_set (qMatrix, 2, 2,
2580  gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
2581 
2582  /* Populate the r.h.s vector of Eq. 18 of the LIGO DCC document T1100433v2 */
2583  /* Amplitude */
2584  gsl_spline_init (spline, timeVec->data, amplitude->data, amplitude->length);
2585  gsl_interp_accel_reset (acc);
2586  amp = gsl_spline_eval (spline, nrTimePeak, acc);
2587  aDot = gsl_spline_eval_deriv (spline, nrTimePeak, acc);
2588  aDDot = gsl_spline_eval_deriv2 (spline, nrTimePeak, acc);
2589 
2590  /* qNSLM */
2591  gsl_spline_init (spline, timeVec->data, qNSLM->data, qNSLM->length);
2592  gsl_interp_accel_reset (acc);
2593  qNSLMPeak = gsl_spline_eval (spline, nrTimePeak, acc);
2594  qNSLMDot = gsl_spline_eval_deriv (spline, nrTimePeak, acc);
2595  qNSLMDDot = gsl_spline_eval_deriv2 (spline, nrTimePeak, acc);
2596 
2597 
2598  nra =
2599  fabs(XLALSimIMREOBGetNRSpinPeakAmplitudeV4 (modeL, modeM, m1,m2,chiS,chiA));
2600 
2601  nraDot =
2602  XLALSimIMREOBGetNRSpinPeakADotV4 (modeL, modeM, m1,m2,chiS,chiA);
2603  //RC: In SEOBNRv4 nraDot is zero because the NQC are defining the peak of the 22 mode
2604  // which by definition has a first derivative equal to 0.
2605  //For SEOBNRv4HM we are not computing the NQC at the peak fo the modes (see Eq.(4.3))
2606  //of https://arxiv.org/pdf/1803.10701.pdf, so first the derivative
2607  //is entering as a fitted coefficient.
2608 
2609  nraDDot =
2610  XLALSimIMREOBGetNRSpinPeakADDotV4 (modeL, modeM, m1,m2,chiS,chiA);
2611 // printf("eta, chiS, chiA, dM/M, chi = %.16e %.16e %.16e %.16e %.16e\n",eta,chiS,chiA, (m1 - m2)/(m1 + m2),chiS + chiA*(m1 - m2)/(m1 + m2)/(1. - 2.*eta));
2612  if (XLAL_IS_REAL8_FAIL_NAN (nra) || XLAL_IS_REAL8_FAIL_NAN (nraDot) || XLAL_IS_REAL8_FAIL_NAN (nraDDot))
2613  {
2614  XLALDestroyREAL8Vector (timeVec);
2615  XLALDestroyREAL8Vector (q3);
2616  XLALDestroyREAL8Vector (q4);
2617  XLALDestroyREAL8Vector (q5);
2618  XLALDestroyREAL8Vector (p3);
2619  XLALDestroyREAL8Vector (p4);
2620  XLALDestroyREAL8Vector (qNS);
2621  XLALDestroyREAL8Vector (pNS);
2622  XLALDestroyREAL8Vector (q3LM);
2623  XLALDestroyREAL8Vector (q4LM);
2624  XLALDestroyREAL8Vector (q5LM);
2625  XLALDestroyREAL8Vector (qNSLM);
2626  XLAL_ERROR (XLAL_EFUNC);
2627  }
2628 
2629  gsl_vector_set (amps, 0, nra - amp);
2630  gsl_vector_set (amps, 1,nraDot -aDot);
2631 
2632  gsl_vector_set (amps, 2, nraDDot - aDDot);
2633 // printf("Amps %.16e %.16e %.16e\n", nra, amp, qNSLMPeak);
2634 // printf("dAmps %.16e %.16e\n", aDot, qNSLMDot);
2635 // printf("ddAmps %.16e %.16e %.16e\n", nraDDot, aDDot, qNSLMDDot);
2636 
2637 
2638  /* We have now set up all the stuff to calculate the a coefficients */
2639  /* So let us do it! */
2640  gsl_linalg_LU_decomp (qMatrix, perm1, &signum);
2641  gsl_linalg_LU_solve (qMatrix, perm1, amps, aCoeff);
2642 
2643  if (debugAT)
2644  {
2645  printf ("Q MATRIX\n");
2646  for (unsigned int i = 0; i < 3; i++)
2647  {
2648  for (unsigned int j = 0; j < 3; j++)
2649  {
2650  printf ("%.12e\t", gsl_matrix_get (qMatrix, i, j));
2651  }
2652  printf ("= %.12e\n", gsl_vector_get (amps, i));
2653  }
2654  }
2655 
2656 
2657 
2658  /* Now we (should) have calculated the a values. Now we can do the b values */
2659 
2660  /* Populate the P matrix in Eq. 18 of the LIGO DCC document T1100433v2 */
2661  /* P3 */
2662  gsl_spline_init (spline, timeVec->data, p3->data, p3->length);
2663  gsl_interp_accel_reset (acc);
2664  gsl_matrix_set (pMatrix, 0, 0,
2665  -gsl_spline_eval_deriv (spline, nrTimePeak, acc));
2666  gsl_matrix_set (pMatrix, 1, 0,
2667  -gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
2668 
2669  /* P4 */
2670  gsl_spline_init (spline, timeVec->data, p4->data, p4->length);
2671  gsl_interp_accel_reset (acc);
2672  gsl_matrix_set (pMatrix, 0, 1,
2673  -gsl_spline_eval_deriv (spline, nrTimePeak, acc));
2674  gsl_matrix_set (pMatrix, 1, 1,
2675  -gsl_spline_eval_deriv2 (spline, nrTimePeak, acc));
2676 
2677  /* Populate the r.h.s vector of Eq. 18 of the LIGO DCC document T1100433v2 */
2678  /* Phase */
2679  gsl_spline_init (spline, timeVec->data, phase->data, phase->length);
2680  gsl_interp_accel_reset (acc);
2681  omega = gsl_spline_eval_deriv (spline, nrTimePeak, acc);
2682  omegaDot = gsl_spline_eval_deriv2 (spline, nrTimePeak, acc);
2683 
2684  /* pNSLM */
2685  gsl_spline_init (spline, timeVec->data, pNS->data, pNS->length);
2686  gsl_interp_accel_reset (acc);
2687  pNSLMDot = gsl_spline_eval_deriv (spline, nrTimePeak, acc);
2688  pNSLMDDot = gsl_spline_eval_deriv2 (spline, nrTimePeak, acc);
2689 
2690  /* Since the phase can be decreasing, we need to take care not to have a -ve frequency */
2691  if (omega * omegaDot > 0.0)
2692  {
2693  omega = fabs (omega);
2694  omegaDot = fabs (omegaDot);
2695  }
2696  else
2697  {
2698  omega = fabs (omega);
2699  omegaDot = -fabs (omegaDot);
2700  }
2701 
2702  nromega =
2703  XLALSimIMREOBGetNRSpinPeakOmegaV4 (modeL, modeM, eta,
2704  chiS + chiA * (m1 - m2) / (m1 + m2) / (1. -
2705  2. * eta));
2706  nromegaDot =
2707  XLALSimIMREOBGetNRSpinPeakOmegaDotV4 (modeL, modeM, eta,
2708  chiS + chiA * (m1 - m2) / (m1 + m2) / (1. -
2709  2. *
2710  eta));
2711 
2712  if (debugAT)
2713  printf ("NR inputs: %.16e, %.16e, %.16e, %.16e\n", nra, nraDDot, nromega,
2714  nromegaDot);
2715 
2716 /* printf("NR inputs: %.16e, %.16e, %.16e, %.16e\n",pNSLMDot, pNSLMDDot,omega,omegaDot);*/
2717 
2718  if (XLAL_IS_REAL8_FAIL_NAN (nromega) || XLAL_IS_REAL8_FAIL_NAN (nromegaDot))
2719  {
2720  XLALDestroyREAL8Vector (timeVec);
2721  XLALDestroyREAL8Vector (q3);
2722  XLALDestroyREAL8Vector (q4);
2723  XLALDestroyREAL8Vector (q5);
2724  XLALDestroyREAL8Vector (p3);
2725  XLALDestroyREAL8Vector (p4);
2726  XLALDestroyREAL8Vector (qNS);
2727  XLALDestroyREAL8Vector (pNS);
2728  XLALDestroyREAL8Vector (q3LM);
2729  XLALDestroyREAL8Vector (q4LM);
2730  XLALDestroyREAL8Vector (q5LM);
2731  XLALDestroyREAL8Vector (qNSLM);
2732  XLAL_ERROR (XLAL_EFUNC);
2733  }
2734 
2735  gsl_vector_set (omegaVec, 0, nromega - omega);
2736  gsl_vector_set (omegaVec, 1, nromegaDot - omegaDot);
2737 
2738  if (debugAT)
2739  {
2740  printf ("P MATRIX\n");
2741  for (unsigned int i = 0; i < 2; i++)
2742  {
2743  for (unsigned int j = 0; j < 2; j++)
2744  {
2745  printf ("%.12e\t", gsl_matrix_get (pMatrix, i, j));
2746  }
2747  printf ("= %.12e\n", gsl_vector_get (omegaVec, i));
2748  }
2749  }
2750 
2751  /* And now solve for the b coefficients */
2752  gsl_linalg_LU_decomp (pMatrix, perm2, &signum);
2753  gsl_linalg_LU_solve (pMatrix, perm2, omegaVec, bCoeff);
2754 
2755  /* We can now populate the coefficients structure */
2756  coeffs->a3S = 0.;
2757  coeffs->a4 = 0.;
2758  coeffs->a5 = 0.;
2759  coeffs->b3 = 0.;
2760  coeffs->b4 = 0.;
2761 
2762 
2763  coeffs->a1 = gsl_vector_get (aCoeff, 0);
2764  coeffs->a2 = gsl_vector_get (aCoeff, 1);
2765  coeffs->a3 = gsl_vector_get (aCoeff, 2);
2766 
2767  coeffs->b1 = gsl_vector_get (bCoeff, 0);
2768  coeffs->b2 = gsl_vector_get (bCoeff, 1);
2769 
2770 
2771 
2772 // printf( "NQC coefficients:\n" );
2773 // printf( "{%f,%f,%f,%f,%f,%f}\n", coeffs->a1, coeffs->a2, coeffs->a3, coeffs->a3S, coeffs->a4, coeffs->a5 );
2774 //
2775 // printf( "{%f,%f,%f,%f}\n", coeffs->b1, coeffs->b2, coeffs->b3, coeffs->b4 );
2776 
2777 
2778  /* Free memory and exit */
2779  gsl_matrix_free (qMatrix);
2780  gsl_vector_free (amps);
2781  gsl_vector_free (aCoeff);
2782  gsl_permutation_free (perm1);
2783 
2784  gsl_matrix_free (pMatrix);
2785  gsl_vector_free (omegaVec);
2786  gsl_vector_free (bCoeff);
2787  gsl_permutation_free (perm2);
2788 
2789  gsl_spline_free (spline);
2790  gsl_interp_accel_free (acc);
2791 
2792  XLALDestroyREAL8Vector (timeVec);
2793  XLALDestroyREAL8Vector (q3);
2794  XLALDestroyREAL8Vector (q4);
2795  XLALDestroyREAL8Vector (q5);
2796  XLALDestroyREAL8Vector (p3);
2797  XLALDestroyREAL8Vector (p4);
2798  XLALDestroyREAL8Vector (qNS);
2799  XLALDestroyREAL8Vector (pNS);
2800  XLALDestroyREAL8Vector (q3LM);
2801  XLALDestroyREAL8Vector (q4LM);
2802  XLALDestroyREAL8Vector (q5LM);
2803  XLALDestroyREAL8Vector (qNSLM);
2804 
2805 
2806  return XLAL_SUCCESS;
2807 }
2808 
2809 #endif /*_LALSIMIMRNQCCORRECTION_C*/
XLALSimIMREOBGetNRSpinPeakOmegaDot
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakOmegaDot(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 UNUSED a)
Peak frequency slope predicted by fitting NR results (currently only 2,2 available).
Definition: LALSimIMREOBNQCCorrection.c:1632
XLALSimIMRSpinEOBNonQCCorrection
static UNUSED int XLALSimIMRSpinEOBNonQCCorrection(COMPLEX16 *restrict nqc, REAL8Vector *restrict values, const REAL8 omega, EOBNonQCCoeffs *restrict coeffs)
This function calculates the non-quasicircular correction to apply to the waveform.
Definition: LALSimIMREOBNQCCorrection.c:488
CombineTPLEQMFits
static REAL8 CombineTPLEQMFits(REAL8 eta, REAL8 A1, REAL8 fEQ, REAL8 fTPL)
Combine two spin-dependent fits of NR input values in a quardatic-in-eta polynomial with the linear-i...
Definition: LALSimIMREOBNQCCorrection.c:911
XLALSimIMREOBGetNRSpinPeakOmegaV2
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakOmegaV2(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 a)
Peak frequency predicted by fitting NR results (currently only 2,2 available).
Definition: LALSimIMREOBNQCCorrection.c:1269
XLALSimIMREOBCalculateNQCCoefficients
static UNUSED int XLALSimIMREOBCalculateNQCCoefficients(EOBNonQCCoeffs *restrict coeffs, REAL8Vector *restrict amplitude, REAL8Vector *restrict phase, REAL8Vector *restrict q1, REAL8Vector *restrict q2, REAL8Vector *restrict q3, REAL8Vector *restrict p1, REAL8Vector *restrict p2, INT4 l, INT4 m, REAL8 timePeak, REAL8 deltaT, REAL8 eta)
This function computes the coefficients a1, a2, etc.
Definition: LALSimIMREOBNQCCorrection.c:547
XLALSimIMREOBNonQCCorrection
static UNUSED int XLALSimIMREOBNonQCCorrection(COMPLEX16 *restrict nqc, REAL8Vector *restrict values, const REAL8 omega, EOBNonQCCoeffs *restrict coeffs)
This function calculates the non-quasicircular correction to apply to the waveform.
Definition: LALSimIMREOBNQCCorrection.c:429
XLALSimIMREOBGetNRSpinPeakADotV4
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakADotV4(INT4 modeL, INT4 modeM, REAL8 UNUSED m1, REAL8 UNUSED m2, REAL8 UNUSED chiS, REAL8 UNUSED chiA)
Peak amplitude slope predicted by fitting NR results.
Definition: LALSimIMREOBNQCCorrection.c:1040
XLALSimIMREOBGetNRSpinPeakAmplitude
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakAmplitude(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 UNUSED a)
Peak amplitude predicted by fitting NR results (currently only 2,2 available).
Definition: LALSimIMREOBNQCCorrection.c:1590
XLALSimIMREOBGetNRSpinPeakDeltaTv2
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakDeltaTv2(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED m1, REAL8 UNUSED m2, REAL8 chi1, REAL8 chi2)
The time difference between the orbital peak and the peak amplitude of the mode in question (currentl...
Definition: LALSimIMREOBNQCCorrection.c:1655
XLALSimIMREOBGetNRSpinPeakOmegaDotV2
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakOmegaDotV2(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 UNUSED a)
Peak frequency slope predicted by fitting NR results (currently only 2,2 available).
Definition: LALSimIMREOBNQCCorrection.c:1465
XLALSimIMREOBGetNRSpinPeakDeltaT
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakDeltaT(INT4 l, INT4 m, REAL8 UNUSED eta, REAL8 a)
The time difference between the orbital peak and the peak amplitude of the mode in question (currentl...
Definition: LALSimIMREOBNQCCorrection.c:818
XLALSimIMREOBGetNRSpinPeakOmegaDotV4
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakOmegaDotV4(INT4 modeL, INT4 modeM, REAL8 UNUSED eta, REAL8 UNUSED a)
Peak frequency slope predicted by fitting NR results The coefficients for the (2,2) are in Eq.
Definition: LALSimIMREOBNQCCorrection.c:1501
GetNRPeakADDot
static REAL8 GetNRPeakADDot(INT4 l, INT4 m, REAL8 eta)
Function which returns second derivative of the amplitude at the peak taken from a fit to numerical r...
Definition: LALSimIMREOBNQCCorrection.c:195
XLALSimIMREOBGetNRSpinPeakOmega
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakOmega(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 a)
Peak frequency predicted by fitting NR results (currently only 2,2 available).
Definition: LALSimIMREOBNQCCorrection.c:1616
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
XLALSimIMREOBGetNRSpinPeakADDot
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakADDot(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 UNUSED a)
Peak amplitude curvature predicted by fitting NR results (currently only 2,2 available).
Definition: LALSimIMREOBNQCCorrection.c:1602
GetNRPeakAmplitude
static REAL8 GetNRPeakAmplitude(INT4 l, INT4 m, REAL8 eta)
Function which returns a value of the expected peak amplitude taken from a fit to numerical relativit...
Definition: LALSimIMREOBNQCCorrection.c:131
XLALSimIMREOBGetNRSpinPeakOmegav2
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakOmegav2(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 a)
Peak frequency predicted by fitting NR results (currently only 2,2 available).
Definition: LALSimIMREOBNQCCorrection.c:1756
XLALSimIMRSpinEOBCalculateNQCCoefficients
static UNUSED int XLALSimIMRSpinEOBCalculateNQCCoefficients(REAL8Vector *restrict amplitude, REAL8Vector *restrict phase, REAL8Vector *restrict rVec, REAL8Vector *restrict prVec, REAL8Vector *restrict orbOmegaVec, INT4 l, INT4 m, REAL8 timePeak, REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 a, REAL8 chiA, REAL8 chiS, EOBNonQCCoeffs *restrict coeffs, UINT4 SpinAlignedEOBversion)
This function computes the coefficients a3s, a4, etc.
Definition: LALSimIMREOBNQCCorrection.c:1848
XLALSimIMREOBGetCalibratedNQCCoeffs
static UNUSED int XLALSimIMREOBGetCalibratedNQCCoeffs(EOBNonQCCoeffs *coeffs, INT4 l, INT4 m, REAL8 eta)
For the 2,2 mode, there are fits available for the NQC coefficients, given in Eqs.
Definition: LALSimIMREOBNQCCorrection.c:390
XLALSimIMREOBGetNRPeakDeltaT
static REAL8 XLALSimIMREOBGetNRPeakDeltaT(INT4 l, INT4 m, REAL8 eta)
Compute the time offset which should be used in computing the non-quasicircular correction and perfor...
Definition: LALSimIMREOBNQCCorrection.c:65
XLALSimIMREOBGetNRSpinPeakOmegaV5
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakOmegaV5(INT4 modeL, INT4 modeM, REAL8 UNUSED eta, REAL8 a)
Peak frequency predicted by fitting NR results used in SEOBNRv5 The coefficients are in Eq.
Definition: LALSimIMREOBNQCCorrection.c:1424
XLALSimIMREOBGetNRSpinPeakOmegaV4
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakOmegaV4(INT4 modeL, INT4 modeM, REAL8 UNUSED eta, REAL8 a)
Peak frequency predicted by fitting NR results The coefficients for the (2,2) are in Eq.
Definition: LALSimIMREOBNQCCorrection.c:1336
XLALSimIMREOBGetNRSpinPeakADDotV4
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakADDotV4(INT4 modeL, INT4 modeM, REAL8 m1, REAL8 m2, REAL8 chiS, REAL8 chiA)
Peak amplitude curvature predicted by fitting NR results.
Definition: LALSimIMREOBNQCCorrection.c:1167
XLALSimIMREOBGetNRSpinPeakDeltaTv4
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakDeltaTv4(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED m1, REAL8 UNUSED m2, REAL8 UNUSED chi1, REAL8 UNUSED chi2)
The time difference between the orbital peak and the peak amplitude of the mode in question (currentl...
Definition: LALSimIMREOBNQCCorrection.c:1699
XLALSimIMREOBGetNRSpinPeakADDotV2
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakADDotV2(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 UNUSED a)
Peak amplitude curvature predicted by fitting NR results (currently only 2,2 available).
Definition: LALSimIMREOBNQCCorrection.c:1131
GetNRPeakOmega
static REAL8 GetNRPeakOmega(INT4 l, INT4 m, REAL8 eta)
Function which returns a value of the expected peak frequency taken from a fit to numerical relativit...
Definition: LALSimIMREOBNQCCorrection.c:261
GetNRPeakOmegaDot
static REAL8 GetNRPeakOmegaDot(INT4 l, INT4 m, REAL8 eta)
Function which returns the derivative of the expected peak frequency taken from a fit to numerical re...
Definition: LALSimIMREOBNQCCorrection.c:324
XLALSimIMREOBGetNRSpinPeakOmegaDotv2
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakOmegaDotv2(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 UNUSED a)
Peak frequency slope predicted by fitting NR results (currently only 2,2 available).
Definition: LALSimIMREOBNQCCorrection.c:1815
XLALSimIMREOBGetNRSpinPeakAmplitudeV2
static UNUSED REAL8 XLALSimIMREOBGetNRSpinPeakAmplitudeV2(INT4 UNUSED l, INT4 UNUSED m, REAL8 UNUSED eta, REAL8 UNUSED a)
Peak amplitude predicted by fitting NR results (currently only 2,2 available).
Definition: LALSimIMREOBNQCCorrection.c:869
XLALSimIMRSpinEOBCalculateNQCCoefficientsV4
static UNUSED int XLALSimIMRSpinEOBCalculateNQCCoefficientsV4(REAL8Vector *restrict amplitude, REAL8Vector *restrict phase, REAL8Vector *restrict rVec, REAL8Vector *restrict prVec, REAL8Vector *restrict orbOmegaVec, INT4 modeL, INT4 modeM, REAL8 timePeak, REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 chiA, REAL8 chiS, EOBNonQCCoeffs *restrict coeffs)
This function computes the coefficients a3s, a4, etc.
Definition: LALSimIMREOBNQCCorrection.c:2359
XLALSimIMRGetEOBCalibratedSpinNQC
static UNUSED int XLALSimIMRGetEOBCalibratedSpinNQC(EOBNonQCCoeffs *restrict coeffs, INT4 UNUSED l, INT4 UNUSED m, REAL8 eta, REAL8 a)
More recent versions of the EOB models, such as EOBNRv2 and SEOBNRv1, utilise a non-quasicircular cor...
Definition: LALSimIMREOBNQCTables.c:52
XLALSimIMRGetEOBCalibratedSpinNQC3D
static UNUSED int XLALSimIMRGetEOBCalibratedSpinNQC3D(EOBNonQCCoeffs *restrict coeffs, INT4 UNUSED l, INT4 UNUSED m, REAL8 m1, REAL8 m2, REAL8 a, REAL8 chiAin)
Function to 3D-interpolate known amplitude NQC coeffcients of spin terms for SEOBNRv2,...
Definition: LALSimIMREOBNQCTables.c:2424
c1
const double c1
Definition: LALSimIMRLackeyTidal2013.c:51
c2
const double c2
Definition: LALSimIMRLackeyTidal2013.c:52
c0
const double c0
Definition: LALSimIMRLackeyTidal2013.c:50
q5
static REAL8 UNUSED q5(REAL8 UNUSED f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
Definition: LALSimInspiralOptimizedCoefficientsEccentricityFD.c:1237
q3
static REAL8 UNUSED q3(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
Definition: LALSimInspiralOptimizedCoefficientsEccentricityFD.c:1195
q2
static REAL8 UNUSED q2(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
Definition: LALSimInspiralOptimizedCoefficientsEccentricityFD.c:1175
q1
static REAL8 UNUSED q1(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
Definition: LALSimInspiralOptimizedCoefficientsEccentricityFD.c:1155
q4
static REAL8 UNUSED q4(REAL8 e0, REAL8 f0, REAL8 inc, REAL8 bet, REAL8 FPlus, REAL8 FCross)
Definition: LALSimInspiralOptimizedCoefficientsEccentricityFD.c:1215
XLAL_CALLGSL
#define XLAL_CALLGSL(statement)
l
int l
Definition: bh_qnmode.c:135
i
double i
Definition: bh_ringdown.c:118
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
a
static const INT4 a
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLAL_ERROR_REAL8
#define XLAL_ERROR_REAL8(...)
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
p
list p
EOBNonQCCoeffs
The coefficients which are used in calculating the non-quasicircular correction to the EOBNRv2 model.
Definition: LALSimIMREOBNRv2.h:467
REAL8Vector::data
REAL8 * data
deltaT
double deltaT
Definition: unicorn.c:24