LALSimulation  5.4.0.1-fe68b98
LALSimInspiralSpinDominatedWaveform.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2014 M. Tapai
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 #include <math.h>
21 #include <complex.h>
22 
23 #include <lal/LALStdlib.h>
24 #include <lal/LALConstants.h>
25 #include <lal/Date.h>
26 #include <lal/FrequencySeries.h>
27 #include <lal/TimeSeries.h>
28 #include <lal/TimeFreqFFT.h>
29 #include <lal/Units.h>
30 #include <lal/LALSimInspiral.h>
31 #include <lal/LALSimIMR.h>
32 #include <lal/LALAdaptiveRungeKuttaIntegrator.h>
33 
34 #define LAL_SDW_ABSOLUTE_TOLERANCE 1.e-12
35 #define LAL_SDW_RELATIVE_TOLERANCE 1.e-12
36 #define LAL_SDW_NUM_VARIABLES 3
37 #define LAL_SDW_MAX_PN_PARAM 0.8
38 
39 #ifdef __GNUC__
40 #define UNUSED __attribute__ ((unused))
41 #else
42 #define UNUSED
43 #endif
44 
45 
46 
47 static const REAL8 G_CP2 = LAL_G_SI / LAL_C_SI / LAL_C_SI;
48 
49 enum {
50  PNDEF = -1 ,PN00 = 0, PN05 = 1, PN10 = 2, PN15 =3, PN20 =4, PN25 =5, PN3 =6, PN_ORDER =7,
51 };
52 
53 typedef enum {
54  PHI, OMEGA, PSI, PHASE = PSI,
55 } TIME_DEPENDANT;
56 
57 typedef enum {
58  PLUS_ = 0, MINUS = 1, CROSS_ = 1, PLUS_MINUS_DIM = 2, PLUS_CROSS_DIM = 2,
59 } COMPONENTS;
60 
61 typedef enum {
62  PN00DIM = 2, PN05DIM = 11, PN10DIM = 15, PN15DIM = 17,
63 } COEFFICIENT_DIMENSIONS;
64 
65 typedef enum {
66  TRIGONOMETRIC_POWER = 5, AMPCOEFF_DIM = 11, OMEGA_POWER_DIM = 6, PHI_PSI_DIM = 6,
67 } CONSTANTS;
68 
69 #define vectorProd(lhs, rhs, denominator, result) \
70  result[X] = (lhs[Y] * rhs[Z] - lhs[Z] * rhs[Y]) / denominator; \
71  result[Y] = (lhs[Z] * rhs[X] - lhs[X] * rhs[Z]) / denominator; \
72  result[Z] = (lhs[X] * rhs[Y] - lhs[Y] * rhs[X]) / denominator;
73 
74 /**
75  * Structure containing the prefered variabloes for Spin-Dominated waveforms.
76  */
77 typedef struct tagLALSDWaveformParams {
78  REAL8 totalmass; //total mass of the binary
79  REAL8 nu; // mass ratio, <1
80  REAL8 chi1; // chi1 dimensionless spin parameter
81  REAL8 dist; // distance to the source
82  REAL8 kappa1; // angle between L and S1
83  REAL8 beta1; // angle between J and S1
84  REAL8 theta; // angle between J and N
85  REAL8 eps; // PN paramter
86  REAL8 xi; // second small parameter
87  REAL8 omega;
88  int pnamp;
89  int pnphase;
90  REAL8 ccoeff00pn[4];
91  REAL8 ccoeff05pn[22];
92  REAL8 ccoeff10pn[30];
93  REAL8 ccoeff15pn[34];
94  REAL8 prevdomega;
95  REAL8 polarizationangle;
96 } LALSDWaveformParams;
97 
98 static INT4 XLALSpinDominatedWaveformStoppingTest(UNUSED REAL8 t, const REAL8 values[], REAL8 dvalues[],
99 UNUSED void *mparams);
100 
101 static INT4 XLALSpinDominatedWaveformDerivatives(UNUSED REAL8 t, const REAL8 values[], REAL8 dvalues[], void *mparams);
102 
103 int XLALSpinDominatedWaveformBuild(LALSDWaveformParams *params, REAL8 expr[], REAL8TimeSeries **hplus,
104  REAL8TimeSeries **hcross, int idx);
105 
106 int XLALSpinDominatedWaveformConstantCoefficients(LALSDWaveformParams * params);
107 
108 /**
109  * Function for allocating memory for a matrix
110  */
111 static REAL8 *XLALDmatrix(INT8 nrh, INT8 nch) {
112  INT8 size = (nrh) * (nch) * sizeof(REAL8);
113  REAL8 *ptr = (REAL8 *) LALMalloc(size);
114  if (ptr != NULL) {
115  return ptr;
116  }
117  printf("malloc error");
118  return NULL;
119 }
120 
121 /**
122  * Function for freeing memory for a matrix
123  */
124 static void XLALFreeDmatrix(REAL8 *m) {
125  LALFree(m);
126 }
127 
128 /**
129  * Function for calculating the constant coefficients of Spin-Dominated waveforms
130  * See tables 1 to 5 in the appendix of Arxiv:1209.1722
131  */
132 int XLALSpinDominatedWaveformConstantCoefficients(LALSDWaveformParams * params) {
133 
134  int i, j;
135  REAL8 *acoeff00pn, *b0coeff0pn, *d0coeff0pn, *acoeff0_5pn, *b0coeff0_5pn, *d0coeff0_5pn, *acoeff1pn, *b0coeff1pn,
136  *d0coeff1pn, *b1coeff1pn, *d1coeff1pn, *acoeff1_5pn, *b0coeff1_5pn, *d0coeff1_5pn, *b1coeff1_5pn,
137  *d1coeff1_5pn;
138 
139  REAL8 skp[TRIGONOMETRIC_POWER], ckp[TRIGONOMETRIC_POWER], stp[TRIGONOMETRIC_POWER], ctp[TRIGONOMETRIC_POWER];
140  skp[0] = ckp[0] = stp[0] = ctp[0] = 1.;
141  skp[1] = sin(params->kappa1);
142  ckp[1] = cos(params->kappa1);
143  stp[1] = sin(params->theta);
144  ctp[1] = cos(params->theta);
145  for (i = 2; i < TRIGONOMETRIC_POWER; ++i) {
146  skp[i] = skp[1] * skp[i - 1];
147  ckp[i] = ckp[1] * ckp[i - 1];
148  stp[i] = stp[1] * stp[i - 1];
149  ctp[i] = ctp[1] * ctp[i - 1];
150  }
151  REAL8 k[PLUS_MINUS_DIM] = { ckp[1] - 1., ckp[1] + 1. };
152  REAL8 c2t = cos(2. * params->theta);
153 
154 
155  acoeff00pn = XLALDmatrix(PN00DIM, PLUS_MINUS_DIM);
156  b0coeff0pn = XLALDmatrix(PN00DIM, PLUS_MINUS_DIM);
157  d0coeff0pn = XLALDmatrix(PN00DIM, PLUS_MINUS_DIM);
158  acoeff0_5pn = XLALDmatrix(PN05DIM, PLUS_MINUS_DIM);
159  b0coeff0_5pn = XLALDmatrix(PN05DIM, PLUS_MINUS_DIM);
160  d0coeff0_5pn = XLALDmatrix(PN05DIM, PLUS_MINUS_DIM);
161  acoeff1pn = XLALDmatrix(PN10DIM, PLUS_MINUS_DIM);
162  b0coeff1pn = XLALDmatrix(PN10DIM, PLUS_MINUS_DIM);
163  d0coeff1pn = XLALDmatrix(PN10DIM, PLUS_MINUS_DIM);
164  b1coeff1pn = XLALDmatrix(PN10DIM, PLUS_MINUS_DIM);
165  d1coeff1pn = XLALDmatrix(PN10DIM, PLUS_MINUS_DIM);
166  acoeff1_5pn = XLALDmatrix(PN15DIM, PLUS_MINUS_DIM);
167  b0coeff1_5pn = XLALDmatrix(PN15DIM, PLUS_MINUS_DIM);
168  d0coeff1_5pn = XLALDmatrix(PN15DIM, PLUS_MINUS_DIM);
169  b1coeff1_5pn = XLALDmatrix(PN15DIM, PLUS_MINUS_DIM);
170  d1coeff1_5pn = XLALDmatrix(PN15DIM, PLUS_MINUS_DIM);
171  /*
172  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
173  %%%%%%%%%% 0 PN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
174  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
175  */
176  acoeff00pn[0 + PN00DIM * PLUS_] = -2. * k[PLUS_];
177  acoeff00pn[0 + PN00DIM * MINUS] = +2. * k[MINUS];
178  acoeff00pn[1 + PN00DIM * PLUS_] = -k[PLUS_];
179  acoeff00pn[1 + PN00DIM * MINUS] = +k[MINUS];
180  b0coeff0pn[0 + PN00DIM * PLUS_] = -1.;
181  b0coeff0pn[0 + PN00DIM * MINUS] = -1.;
182  b0coeff0pn[1 + PN00DIM * PLUS_] = -2.;
183  b0coeff0pn[1 + PN00DIM * MINUS] = -2.;
184  d0coeff0pn[0 + PN00DIM * PLUS_] = +0.;
185  d0coeff0pn[0 + PN00DIM * MINUS] = +0.;
186  d0coeff0pn[1 + PN00DIM * PLUS_] = +0.;
187  d0coeff0pn[1 + PN00DIM * MINUS] = +0.;
188  for (i = 0; i <= 1; i++) {
189  for (j = 0; j <= 1; j++) {
190  params->ccoeff00pn[i + PN00DIM * j] = acoeff00pn[i + PN00DIM * j]
191  + skp[2] * (b0coeff0pn[i + PN00DIM * j] + ckp[1] * d0coeff0pn[i + PN00DIM * j]);
192  }
193  }
194  /*
195  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
196  %%%%%%%%%% 0.5 PN %%%%%%%%%%%%%%%%%%%%%%%%%%%
197  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
198  */
199  acoeff0_5pn[0 + PN05DIM * PLUS_] = +4. * k[PLUS_] * (6. - stp[2]);
200  acoeff0_5pn[0 + PN05DIM * MINUS] = +4. * k[MINUS] * (6. - stp[2]);
201  acoeff0_5pn[1 + PN05DIM * PLUS_] = +4. * k[PLUS_];
202  acoeff0_5pn[1 + PN05DIM * MINUS] = +4. * k[MINUS];
203  acoeff0_5pn[2 + PN05DIM * PLUS_] = +2. * k[PLUS_] * (6. - stp[2]);
204  acoeff0_5pn[2 + PN05DIM * MINUS] = -2. * k[MINUS] * (6. - stp[2]);
205  acoeff0_5pn[3 + PN05DIM * PLUS_] = 12. * k[PLUS_];
206  acoeff0_5pn[3 + PN05DIM * MINUS] = 12. * k[MINUS];
207  acoeff0_5pn[4 + PN05DIM * PLUS_] = +2. * k[PLUS_] * (2. * stp[2] - 3.);
208  acoeff0_5pn[4 + PN05DIM * MINUS] = -2. * k[MINUS] * (2. * stp[2] - 3.);
209  acoeff0_5pn[5 + PN05DIM * PLUS_] = -2. * k[PLUS_];
210  acoeff0_5pn[5 + PN05DIM * MINUS] = -2. * k[MINUS];
211  acoeff0_5pn[6 + PN05DIM * PLUS_] = -2. * params->ccoeff00pn[1 + PN00DIM * PLUS_] * (6. - stp[2]);
212  acoeff0_5pn[6 + PN05DIM * MINUS] = +2. * params->ccoeff00pn[1 + PN00DIM * MINUS] * (6. - stp[2]);
213  acoeff0_5pn[7 + PN05DIM * PLUS_] = +2. * k[PLUS_];
214  acoeff0_5pn[7 + PN05DIM * MINUS] = -2. * k[MINUS];
215  acoeff0_5pn[8 + PN05DIM * PLUS_] = (44. - 34. * stp[2] + 2. * (5. * stp[2] - 46.) * ckp[1]);
216  acoeff0_5pn[8 + PN05DIM * MINUS] = (44. - 34. * stp[2] - 2. * (5. * stp[2] - 46.) * ckp[1]);
217  acoeff0_5pn[9 + PN05DIM * PLUS_] = 22. + 46. * ckp[1];
218  acoeff0_5pn[9 + PN05DIM * MINUS] = 22. - 46. * ckp[1];
219  acoeff0_5pn[10 + PN05DIM * PLUS_] = -2. * k[PLUS_] * (3. - 2. * stp[2]);
220  acoeff0_5pn[10 + PN05DIM * MINUS] = -2. * k[MINUS] * (3. - 2. * stp[2]);
221  b0coeff0_5pn[0 + PN05DIM * PLUS_] = +(46. - 5. * stp[2]);
222  b0coeff0_5pn[0 + PN05DIM * MINUS] = -(46. - 5. * stp[2]);
223  b0coeff0_5pn[1 + PN05DIM * PLUS_] = +3.;
224  b0coeff0_5pn[1 + PN05DIM * MINUS] = -3.;
225  b0coeff0_5pn[2 + PN05DIM * PLUS_] = (2. - 3. * stp[2]);
226  b0coeff0_5pn[2 + PN05DIM * MINUS] = (2. - 3. * stp[2]);
227  b0coeff0_5pn[3 + PN05DIM * PLUS_] = +23.;
228  b0coeff0_5pn[3 + PN05DIM * MINUS] = -23.;
229  b0coeff0_5pn[4 + PN05DIM * PLUS_] = -c2t;
230  b0coeff0_5pn[4 + PN05DIM * MINUS] = -c2t;
231  b0coeff0_5pn[5 + PN05DIM * PLUS_] = -4.;
232  b0coeff0_5pn[5 + PN05DIM * MINUS] = +4.;
233  b0coeff0_5pn[6 + PN05DIM * PLUS_] = +0.;
234  b0coeff0_5pn[6 + PN05DIM * MINUS] = +0.;
235  b0coeff0_5pn[7 + PN05DIM * PLUS_] = +3.;
236  b0coeff0_5pn[7 + PN05DIM * MINUS] = +3.;
237  b0coeff0_5pn[8 + PN05DIM * PLUS_] = -15. * (2. - 3. * stp[2]);
238  b0coeff0_5pn[8 + PN05DIM * MINUS] = -15. * (2. - 3. * stp[2]);
239  b0coeff0_5pn[9 + PN05DIM * PLUS_] = -15.;
240  b0coeff0_5pn[9 + PN05DIM * MINUS] = -15.;
241  b0coeff0_5pn[10 + PN05DIM * PLUS_] = -4. * (3. - 2. * stp[2]);
242  b0coeff0_5pn[10 + PN05DIM * MINUS] = +4. * (3. - 2. * stp[2]);
243  d0coeff0_5pn[0 + PN05DIM * PLUS_] = +5. * (3. * stp[2] - 2.);
244  d0coeff0_5pn[0 + PN05DIM * MINUS] = +5. * (3. * stp[2] - 2.);
245  d0coeff0_5pn[1 + PN05DIM * PLUS_] = -1.;
246  d0coeff0_5pn[1 + PN05DIM * MINUS] = -1.;
247  d0coeff0_5pn[2 + PN05DIM * PLUS_] = +0.;
248  d0coeff0_5pn[2 + PN05DIM * MINUS] = +0.;
249  d0coeff0_5pn[3 + PN05DIM * PLUS_] = -5.;
250  d0coeff0_5pn[3 + PN05DIM * MINUS] = -5.;
251  d0coeff0_5pn[4 + PN05DIM * PLUS_] = +0.;
252  d0coeff0_5pn[4 + PN05DIM * MINUS] = +0.;
253  d0coeff0_5pn[5 + PN05DIM * PLUS_] = +3.;
254  d0coeff0_5pn[5 + PN05DIM * MINUS] = +3.;
255  d0coeff0_5pn[6 + PN05DIM * PLUS_] = -3. * (2. - 3. * stp[2]);
256  d0coeff0_5pn[6 + PN05DIM * MINUS] = -3. * (2. - 3. * stp[2]);
257  d0coeff0_5pn[7 + PN05DIM * PLUS_] = +0.;
258  d0coeff0_5pn[7 + PN05DIM * MINUS] = +0.;
259  d0coeff0_5pn[8 + PN05DIM * PLUS_] = +0.;
260  d0coeff0_5pn[8 + PN05DIM * MINUS] = +0.;
261  d0coeff0_5pn[9 + PN05DIM * PLUS_] = +0.;
262  d0coeff0_5pn[9 + PN05DIM * MINUS] = +0.;
263  d0coeff0_5pn[10 + PN05DIM * PLUS_] = +3. * c2t;
264  d0coeff0_5pn[10 + PN05DIM * MINUS] = +3. * c2t;
265  for (i = 0; i <= 10; i++) {
266  for (j = 0; j <= 1; j++) {
267  params->ccoeff05pn[i + PN05DIM * j] = acoeff0_5pn[i + PN05DIM * j]
268  + skp[2] * (b0coeff0_5pn[i + PN05DIM * j] + ckp[1] * d0coeff0_5pn[i + PN05DIM * j]);
269  }
270  }
271  /*
272  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
273  %%%%%%%%%% 1 PN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
274  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
275  */
276  acoeff1pn[0 + PN10DIM * PLUS_] = +8. * k[PLUS_];
277  acoeff1pn[0 + PN10DIM * MINUS] = -8. * k[MINUS];
278  acoeff1pn[1 + PN10DIM * PLUS_] = +6. * k[PLUS_] * (stp[2] + 5.);
279  acoeff1pn[1 + PN10DIM * MINUS] = +6. * k[MINUS] * (stp[2] + 5.);
280  acoeff1pn[2 + PN10DIM * PLUS_] = +2. * k[PLUS_] * (4. - stp[2]);
281  acoeff1pn[2 + PN10DIM * MINUS] = +2. * k[MINUS] * (4. - stp[2]);
282  acoeff1pn[3 + PN10DIM * PLUS_] = +2. * k[PLUS_] * (2. * stp[4] + 11. * stp[2] - 38.);
283  acoeff1pn[3 + PN10DIM * MINUS] = -2. * k[MINUS] * (2. * stp[4] + 11. * stp[2] - 38.);
284  acoeff1pn[4 + PN10DIM * PLUS_] = +6. * k[PLUS_] * (3. * stp[2] + 5.);
285  acoeff1pn[4 + PN10DIM * MINUS] = +6. * k[MINUS] * (3. * stp[2] + 5.);
286  acoeff1pn[5 + PN10DIM * PLUS_] = +2. * k[PLUS_] * (4. * stp[2] + 19.);
287  acoeff1pn[5 + PN10DIM * MINUS] = -2. * k[MINUS] * (4. * stp[2] + 19.);
288  acoeff1pn[6 + PN10DIM * PLUS_] = -2. * k[PLUS_] * (3. * stp[2] - 4.);
289  acoeff1pn[6 + PN10DIM * MINUS] = -2. * k[MINUS] * (3. * stp[2] - 4.);
290  acoeff1pn[7 + PN10DIM * PLUS_] = +2. * k[PLUS_] * (4. - stp[2]);
291  acoeff1pn[7 + PN10DIM * MINUS] = -2. * k[MINUS] * (4. - stp[2]);
292  acoeff1pn[8 + PN10DIM * PLUS_] = +6. * k[PLUS_] * (5. + stp[2]);
293  acoeff1pn[8 + PN10DIM * MINUS] = -6. * k[MINUS] * (5. + stp[2]);
294  acoeff1pn[9 + PN10DIM * PLUS_] = -4. * k[PLUS_];
295  acoeff1pn[9 + PN10DIM * MINUS] = +4. * k[MINUS];
296  acoeff1pn[10 + PN10DIM * PLUS_] = k[PLUS_] * (22. + 29. * stp[2] - 16. * stp[4]);
297  acoeff1pn[10 + PN10DIM * MINUS] = k[MINUS] * (22. + 29. * stp[2] - 16. * stp[4]);
298  acoeff1pn[11 + PN10DIM * PLUS_] = +2. * k[PLUS_];
299  acoeff1pn[11 + PN10DIM * MINUS] = +2. * k[MINUS];
300  acoeff1pn[12 + PN10DIM * PLUS_] = +6. * k[PLUS_] * (3. * stp[2] + 5.);
301  acoeff1pn[12 + PN10DIM * MINUS] = -6. * k[MINUS] * (3. * stp[2] + 5.);
302  acoeff1pn[13 + PN10DIM * PLUS_] = -k[PLUS_] * (20. * stp[2] + 11.);
303  acoeff1pn[13 + PN10DIM * MINUS] = -k[MINUS] * (20. * stp[2] + 11.);
304  acoeff1pn[14 + PN10DIM * PLUS_] = -2. * k[PLUS_] * (3. * stp[2] - 4.);
305  acoeff1pn[14 + PN10DIM * MINUS] = +2. * k[MINUS] * (3. * stp[2] - 4.);
306  b0coeff1pn[0 + PN10DIM * PLUS_] = +8.;
307  b0coeff1pn[0 + PN10DIM * MINUS] = +8.;
308  b0coeff1pn[1 + PN10DIM * PLUS_] = -18. + 7. * stp[2];
309  b0coeff1pn[1 + PN10DIM * MINUS] = +18. - 7. * stp[2];
310  b0coeff1pn[2 + PN10DIM * PLUS_] = -3. * stp[2] + 6.;
311  b0coeff1pn[2 + PN10DIM * MINUS] = +3. * stp[2] - 6.;
312  b0coeff1pn[3 + PN10DIM * PLUS_] = (-22. - 29. * stp[2] + 16. * stp[4]);
313  b0coeff1pn[3 + PN10DIM * MINUS] = (-22. - 29. * stp[2] + 16. * stp[4]);
314  b0coeff1pn[4 + PN10DIM * PLUS_] = +26. * stp[2] - 18.;
315  b0coeff1pn[4 + PN10DIM * MINUS] = -26. * stp[2] + 18.;
316  b0coeff1pn[5 + PN10DIM * PLUS_] = +11. + 20. * stp[2];
317  b0coeff1pn[5 + PN10DIM * MINUS] = +11. + 20. * stp[2];
318  b0coeff1pn[6 + PN10DIM * PLUS_] = -6. * stp[2] + 6.;
319  b0coeff1pn[6 + PN10DIM * MINUS] = +6. * stp[2] - 6.;
320  b0coeff1pn[7 + PN10DIM * PLUS_] = +2. * (11. - 5. * stp[2]);
321  b0coeff1pn[7 + PN10DIM * MINUS] = +2. * (11. - 5. * stp[2]);
322  b0coeff1pn[8 + PN10DIM * PLUS_] = +6. * (7. + 9. * stp[2]);
323  b0coeff1pn[8 + PN10DIM * MINUS] = +6. * (7. + 9. * stp[2]);
324  b0coeff1pn[9 + PN10DIM * PLUS_] = -11.;
325  b0coeff1pn[9 + PN10DIM * MINUS] = -11.;
326  b0coeff1pn[10 + PN10DIM * PLUS_] = -3. * (8. - 20. * stp[2] + 7. * stp[4]);
327  b0coeff1pn[10 + PN10DIM * MINUS] = +3. * (8. - 20. * stp[2] + 7. * stp[4]);
328  b0coeff1pn[11 + PN10DIM * PLUS_] = +3.;
329  b0coeff1pn[11 + PN10DIM * MINUS] = -3.;
330  b0coeff1pn[12 + PN10DIM * PLUS_] = +3. * (19. * stp[2] + 14.);
331  b0coeff1pn[12 + PN10DIM * MINUS] = +3. * (19. * stp[2] + 14.);
332  b0coeff1pn[13 + PN10DIM * PLUS_] = +12. * c2t;
333  b0coeff1pn[13 + PN10DIM * MINUS] = -12. * c2t;
334  b0coeff1pn[14 + PN10DIM * PLUS_] = (22. - 21. * stp[2]);
335  b0coeff1pn[14 + PN10DIM * MINUS] = (22. - 21. * stp[2]);
336  d0coeff1pn[0 + PN10DIM * PLUS_] = -4.;
337  d0coeff1pn[0 + PN10DIM * MINUS] = +4.;
338  d0coeff1pn[1 + PN10DIM * PLUS_] = (6. - 14. * stp[2]);
339  d0coeff1pn[1 + PN10DIM * MINUS] = (6. - 14. * stp[2]);
340  d0coeff1pn[2 + PN10DIM * PLUS_] = +2. * (stp[2] - 1.);
341  d0coeff1pn[2 + PN10DIM * MINUS] = +2. * (stp[2] - 1.);
342  d0coeff1pn[3 + PN10DIM * PLUS_] = -2. * (8. - 20. * stp[2] + 7. * stp[4]);
343  d0coeff1pn[3 + PN10DIM * MINUS] = +2. * (8. - 20. * stp[2] + 7. * stp[4]);
344  d0coeff1pn[4 + PN10DIM * PLUS_] = (6. - 7. * stp[2]);
345  d0coeff1pn[4 + PN10DIM * MINUS] = (6. - 7. * stp[2]);
346  d0coeff1pn[5 + PN10DIM * PLUS_] = (-16. * stp[2] + 8.);
347  d0coeff1pn[5 + PN10DIM * MINUS] = (+16. * stp[2] - 8.);
348  d0coeff1pn[6 + PN10DIM * PLUS_] = (3. * stp[2] - 2.);
349  d0coeff1pn[6 + PN10DIM * MINUS] = (3. * stp[2] - 2.);
350  d0coeff1pn[7 + PN10DIM * PLUS_] = +9. * (stp[2] - 2.);
351  d0coeff1pn[7 + PN10DIM * MINUS] = -9. * (stp[2] - 2.);
352  d0coeff1pn[8 + PN10DIM * PLUS_] = +3. * (18. - 7. * stp[2]);
353  d0coeff1pn[8 + PN10DIM * MINUS] = -3. * (18. - 7. * stp[2]);
354  d0coeff1pn[9 + PN10DIM * PLUS_] = +9.;
355  d0coeff1pn[9 + PN10DIM * MINUS] = -9.;
356  d0coeff1pn[10 + PN10DIM * PLUS_] = +4. * (2. - 8. * stp[2] + 7. * stp[4]);
357  d0coeff1pn[10 + PN10DIM * MINUS] = +4. * (2. - 8. * stp[2] + 7. * stp[4]);
358  d0coeff1pn[11 + PN10DIM * PLUS_] = -2.;
359  d0coeff1pn[11 + PN10DIM * MINUS] = -2.;
360  d0coeff1pn[12 + PN10DIM * PLUS_] = +6. * (9. - 13. * stp[2]);
361  d0coeff1pn[12 + PN10DIM * MINUS] = -6. * (9. - 13. * stp[2]);
362  d0coeff1pn[13 + PN10DIM * PLUS_] = +2. * (7. * stp[2] - 2.);
363  d0coeff1pn[13 + PN10DIM * MINUS] = +2. * (7. * stp[2] - 2.);
364  d0coeff1pn[14 + PN10DIM * PLUS_] = -18. * ctp[2];
365  d0coeff1pn[14 + PN10DIM * MINUS] = +18. * ctp[2];
366  b1coeff1pn[0 + PN10DIM * PLUS_] = -1.;
367  b1coeff1pn[0 + PN10DIM * MINUS] = -1.;
368  b1coeff1pn[1 + PN10DIM * PLUS_] = +0.;
369  b1coeff1pn[1 + PN10DIM * MINUS] = +0.;
370  b1coeff1pn[2 + PN10DIM * PLUS_] = +0.;
371  b1coeff1pn[2 + PN10DIM * MINUS] = +0.;
372  b1coeff1pn[3 + PN10DIM * PLUS_] = -2. * (2. - 8. * stp[2] + 7. * stp[4]);
373  b1coeff1pn[3 + PN10DIM * MINUS] = -2. * (2. - 8. * stp[2] + 7. * stp[4]);
374  b1coeff1pn[4 + PN10DIM * PLUS_] = +0.;
375  b1coeff1pn[4 + PN10DIM * MINUS] = +0.;
376  b1coeff1pn[5 + PN10DIM * PLUS_] = (2. - 7. * stp[2]);
377  b1coeff1pn[5 + PN10DIM * MINUS] = (2. - 7. * stp[2]);
378  b1coeff1pn[6 + PN10DIM * PLUS_] = +0.;
379  b1coeff1pn[6 + PN10DIM * MINUS] = +0.;
380  b1coeff1pn[7 + PN10DIM * PLUS_] = -8. * ctp[2];
381  b1coeff1pn[7 + PN10DIM * MINUS] = -8. * ctp[2];
382  b1coeff1pn[8 + PN10DIM * PLUS_] = +8. * (3. - 7. * stp[2]);
383  b1coeff1pn[8 + PN10DIM * MINUS] = +8. * (3. - 7. * stp[2]);
384  b1coeff1pn[9 + PN10DIM * PLUS_] = +4.;
385  b1coeff1pn[9 + PN10DIM * MINUS] = +4.;
386  b1coeff1pn[10 + PN10DIM * PLUS_] = +0.;
387  b1coeff1pn[10 + PN10DIM * MINUS] = +0.;
388  b1coeff1pn[11 + PN10DIM * PLUS_] = +0.;
389  b1coeff1pn[11 + PN10DIM * MINUS] = +0.;
390  b1coeff1pn[12 + PN10DIM * PLUS_] = -4. * (7. * stp[2] - 6.);
391  b1coeff1pn[12 + PN10DIM * MINUS] = -4. * (7. * stp[2] - 6.);
392  b1coeff1pn[13 + PN10DIM * PLUS_] = +0.;
393  b1coeff1pn[13 + PN10DIM * MINUS] = +0.;
394  b1coeff1pn[14 + PN10DIM * PLUS_] = -4. * (2. - 3. * stp[2]);
395  b1coeff1pn[14 + PN10DIM * MINUS] = -4. * (2. - 3. * stp[2]);
396  for (i = 0; i < 15; i++) {
397  d1coeff1pn[i + PN10DIM * PLUS_] = +0.;
398  d1coeff1pn[i + PN10DIM * MINUS] = +0.;
399  }
400  for (i = 0; i <= 14; i++) {
401  for (j = 0; j <= 1; j++) {
402  params->ccoeff10pn[i + PN10DIM * j] = acoeff1pn[i + PN10DIM * j]
403  + skp[2] * (b0coeff1pn[i + PN10DIM * j] + ckp[1] * d0coeff1pn[i + PN10DIM * j])
404  + skp[4] * (b1coeff1pn[i + PN10DIM * j] + ckp[1] * d1coeff1pn[i + PN10DIM * j]);
405  }
406  }
407  /*
408  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
409  %%%%%%%%%% 1.5 PN %%%%%%%%%%%%%%%%%%%%%%%%%%%
410  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
411  */
412  acoeff1_5pn[0 + PN15DIM * PLUS_] = +4. * k[PLUS_] * (stp[2] - 6.);
413  acoeff1_5pn[0 + PN15DIM * MINUS] = -4. * k[MINUS] * (stp[2] - 6.);
414  acoeff1_5pn[1 + PN15DIM * PLUS_] = +4. * k[PLUS_] * (stp[4] + 42. * stp[2] - 166.);
415  acoeff1_5pn[1 + PN15DIM * MINUS] = -4. * k[MINUS] * (stp[4] + 42. * stp[2] - 166.);
416  acoeff1_5pn[2 + PN15DIM * PLUS_] = +16. * k[PLUS_];
417  acoeff1_5pn[2 + PN15DIM * MINUS] = +16. * k[MINUS];
418  acoeff1_5pn[3 + PN15DIM * PLUS_] = +8. * k[PLUS_] * (stp[4] + 8. * stp[2] - 28.);
419  acoeff1_5pn[3 + PN15DIM * MINUS] = +8. * k[MINUS] * (stp[4] + 8. * stp[2] - 28.);
420  acoeff1_5pn[4 + PN15DIM * PLUS_] = +8. * k[PLUS_] * (-332. + 94. * stp[2] + stp[4]);
421  acoeff1_5pn[4 + PN15DIM * MINUS] = +8. * k[MINUS] * (-332. + 94. * stp[2] + stp[4]);
422  acoeff1_5pn[5 + PN15DIM * PLUS_] = +8. * k[PLUS_] * (38. - 42. * stp[2] - 9. * stp[4]);
423  acoeff1_5pn[5 + PN15DIM * MINUS] = -8. * k[MINUS] * (38. - 42. * stp[2] - 9. * stp[4]);
424  acoeff1_5pn[6 + PN15DIM * PLUS_] = -16. * k[PLUS_] * (152. - 46. * stp[2] - 9. * stp[4]);
425  acoeff1_5pn[6 + PN15DIM * MINUS] = -16. * k[MINUS] * (152. - 46. * stp[2] - 9. * stp[4]);
426  acoeff1_5pn[7 + PN15DIM * PLUS_] = +24. * k[PLUS_] * (3. * stp[2] - 10.);
427  acoeff1_5pn[7 + PN15DIM * MINUS] = -24. * k[MINUS] * (3. * stp[2] - 10.);
428  acoeff1_5pn[8 + PN15DIM * PLUS_] = -8. * k[PLUS_] * (160. - 204. * stp[2] - 63. * stp[4]);
429  acoeff1_5pn[8 + PN15DIM * MINUS] = -8. * k[MINUS] * (160. - 204. * stp[2] - 63. * stp[4]);
430  acoeff1_5pn[9 + PN15DIM * PLUS_] = +4. * k[PLUS_] * (3. - 2. * stp[2]);
431  acoeff1_5pn[9 + PN15DIM * MINUS] = -4. * k[MINUS] * (3. - 2. * stp[2]);
432  acoeff1_5pn[10 + PN15DIM * PLUS_] = -8. * k[PLUS_] * (14. + 3. * stp[2]);
433  acoeff1_5pn[10 + PN15DIM * MINUS] = -8. * k[MINUS] * (14. + 3. * stp[2]);
434  acoeff1_5pn[11 + PN15DIM * PLUS_] = -16. * k[PLUS_] * (15. * stp[2] + 76.);
435  acoeff1_5pn[11 + PN15DIM * MINUS] = -16. * k[MINUS] * (15. * stp[2] + 76.);
436  acoeff1_5pn[12 + PN15DIM * PLUS_] = -8. * k[PLUS_] * (5. * stp[2] + 166.);
437  acoeff1_5pn[12 + PN15DIM * MINUS] = -8. * k[MINUS] * (5. * stp[2] + 166.);
438  acoeff1_5pn[13 + PN15DIM * PLUS_] = -8. * k[PLUS_] * (80. + 63. * stp[2]);
439  acoeff1_5pn[13 + PN15DIM * MINUS] = -8. * k[MINUS] * (80. + 63. * stp[2]);
440  acoeff1_5pn[14 + PN15DIM * PLUS_] = +4. * k[PLUS_] * (166. - 125. * stp[2] - 8. * stp[4]);
441  acoeff1_5pn[14 + PN15DIM * MINUS] = -4. * k[MINUS] * (166. - 125. * stp[2] - 8. * stp[4]);
442  acoeff1_5pn[15 + PN15DIM * PLUS_] = -8. * k[PLUS_] * (38. - 61. * stp[2] - 24. * stp[4]);
443  acoeff1_5pn[15 + PN15DIM * MINUS] = +8. * k[MINUS] * (38. - 61. * stp[2] - 24. * stp[4]);
444  acoeff1_5pn[16 + PN15DIM * PLUS_] = +8. * k[PLUS_] * (5. - 4. * stp[2]);
445  acoeff1_5pn[16 + PN15DIM * MINUS] = -8. * k[MINUS] * (5. - 4. * stp[2]);
446  b0coeff1_5pn[0 + PN15DIM * PLUS_] = (5. * stp[2] - 6.);
447  b0coeff1_5pn[0 + PN15DIM * MINUS] = (5. * stp[2] - 6.);
448  b0coeff1_5pn[1 + PN15DIM * PLUS_] = (18. * stp[4] + 252. * stp[2] - 188.);
449  b0coeff1_5pn[1 + PN15DIM * MINUS] = (18. * stp[4] + 252. * stp[2] - 188.);
450  b0coeff1_5pn[2 + PN15DIM * PLUS_] = +20.;
451  b0coeff1_5pn[2 + PN15DIM * MINUS] = -20.;
452  b0coeff1_5pn[3 + PN15DIM * PLUS_] = (+9. * stp[4] - 90. * stp[2] + 56.);
453  b0coeff1_5pn[3 + PN15DIM * MINUS] = (-9. * stp[4] + 90. * stp[2] - 56.);
454  b0coeff1_5pn[4 + PN15DIM * PLUS_] = -4. * (1184. - 172. * stp[2] - 7. * stp[4]);
455  b0coeff1_5pn[4 + PN15DIM * MINUS] = +4. * (1184. - 172. * stp[2] - 7. * stp[4]);
456  b0coeff1_5pn[5 + PN15DIM * PLUS_] = +2. * (46. + 48. * stp[2] - 99. * stp[4]);
457  b0coeff1_5pn[5 + PN15DIM * MINUS] = +2. * (46. + 48. * stp[2] - 99. * stp[4]);
458  b0coeff1_5pn[6 + PN15DIM * PLUS_] = -12. * (72. + 110. * stp[2] - 63. * stp[4]);
459  b0coeff1_5pn[6 + PN15DIM * MINUS] = +12. * (72. + 110. * stp[2] - 63. * stp[4]);
460  b0coeff1_5pn[7 + PN15DIM * PLUS_] = +144. * (stp[2] - 2.);
461  b0coeff1_5pn[7 + PN15DIM * MINUS] = +144. * (stp[2] - 2.);
462  b0coeff1_5pn[8 + PN15DIM * PLUS_] = -3. * (-204. + 406. * stp[2] - 189. * stp[4]);
463  b0coeff1_5pn[8 + PN15DIM * MINUS] = +3. * (-204. + 406. * stp[2] - 189. * stp[4]);
464  b0coeff1_5pn[9 + PN15DIM * PLUS_] = +3. - 4. * stp[2];
465  b0coeff1_5pn[9 + PN15DIM * MINUS] = +3. - 4. * stp[2];
466  b0coeff1_5pn[10 + PN15DIM * PLUS_] = +28. - 31. * stp[2];
467  b0coeff1_5pn[10 + PN15DIM * MINUS] = -28. + 31. * stp[2];
468  b0coeff1_5pn[11 + PN15DIM * PLUS_] = -432. - 876. * stp[2];
469  b0coeff1_5pn[11 + PN15DIM * MINUS] = +432. + 876. * stp[2];
470  b0coeff1_5pn[12 + PN15DIM * PLUS_] = -4. * (71. * stp[2] + 592.);
471  b0coeff1_5pn[12 + PN15DIM * MINUS] = +4. * (71. * stp[2] + 592.);
472  b0coeff1_5pn[13 + PN15DIM * PLUS_] = +306. - 651. * stp[2];
473  b0coeff1_5pn[13 + PN15DIM * MINUS] = -306. + 651. * stp[2];
474  b0coeff1_5pn[14 + PN15DIM * PLUS_] = +2. * (94. - 173. * stp[2] - 24. * stp[4]);
475  b0coeff1_5pn[14 + PN15DIM * MINUS] = +2. * (94. - 173. * stp[2] - 24. * stp[4]);
476  b0coeff1_5pn[15 + PN15DIM * PLUS_] = -2. * (46. + 25. * stp[2] - 180. * stp[4]);
477  b0coeff1_5pn[15 + PN15DIM * MINUS] = -2. * (46. + 25. * stp[2] - 180. * stp[4]);
478  b0coeff1_5pn[16 + PN15DIM * PLUS_] = +48. * ctp[2];
479  b0coeff1_5pn[16 + PN15DIM * MINUS] = +48. * ctp[2];
480  d0coeff1_5pn[0 + PN15DIM * PLUS_] = +0.;
481  d0coeff1_5pn[0 + PN15DIM * MINUS] = +0.;
482  d0coeff1_5pn[1 + PN15DIM * PLUS_] = (-6. * stp[4] + 72. * stp[2] - 20.);
483  d0coeff1_5pn[1 + PN15DIM * MINUS] = (+6. * stp[4] - 72. * stp[2] + 20.);
484  d0coeff1_5pn[2 + PN15DIM * PLUS_] = -12.;
485  d0coeff1_5pn[2 + PN15DIM * MINUS] = -12.;
486  d0coeff1_5pn[3 + PN15DIM * PLUS_] = (-15. * stp[4] + 22. * stp[2] - 8.);
487  d0coeff1_5pn[3 + PN15DIM * MINUS] = (-15. * stp[4] + 22. * stp[2] - 8.);
488  d0coeff1_5pn[4 + PN15DIM * PLUS_] = (1920. - 2832. * stp[2] - 84. * stp[4]);
489  d0coeff1_5pn[4 + PN15DIM * MINUS] = (1920. - 2832. * stp[2] - 84. * stp[4]);
490  d0coeff1_5pn[5 + PN15DIM * PLUS_] = +6. * (10. - 44. * stp[2] + 27. * stp[4]);
491  d0coeff1_5pn[5 + PN15DIM * MINUS] = -6. * (10. - 44. * stp[2] + 27. * stp[4]);
492  d0coeff1_5pn[6 + PN15DIM * PLUS_] = -4. * (88. - 422. * stp[2] + 171. * stp[4]);
493  d0coeff1_5pn[6 + PN15DIM * MINUS] = -4. * (88. - 422. * stp[2] + 171. * stp[4]);
494  d0coeff1_5pn[7 + PN15DIM * PLUS_] = +12. * (14. - 9. * stp[2]);
495  d0coeff1_5pn[7 + PN15DIM * MINUS] = -12. * (14. - 9. * stp[2]);
496  d0coeff1_5pn[8 + PN15DIM * PLUS_] = -9. * (28. - 126. * stp[2] + 105. * stp[4]);
497  d0coeff1_5pn[8 + PN15DIM * MINUS] = -9. * (28. - 126. * stp[2] + 105. * stp[4]);
498  d0coeff1_5pn[9 + PN15DIM * PLUS_] = +0.;
499  d0coeff1_5pn[9 + PN15DIM * MINUS] = +0.;
500  d0coeff1_5pn[10 + PN15DIM * PLUS_] = (9. * stp[2] - 4.);
501  d0coeff1_5pn[10 + PN15DIM * MINUS] = (9. * stp[2] - 4.);
502  d0coeff1_5pn[11 + PN15DIM * PLUS_] = (-176. + 756. * stp[2]);
503  d0coeff1_5pn[11 + PN15DIM * MINUS] = (-176. + 756. * stp[2]);
504  d0coeff1_5pn[12 + PN15DIM * PLUS_] = +12. * (7. * stp[2] + 80.);
505  d0coeff1_5pn[12 + PN15DIM * MINUS] = +12. * (7. * stp[2] + 80.);
506  d0coeff1_5pn[13 + PN15DIM * PLUS_] = (-126. + 189. * stp[2]);
507  d0coeff1_5pn[13 + PN15DIM * MINUS] = (-126. + 189. * stp[2]);
508  d0coeff1_5pn[14 + PN15DIM * PLUS_] = +2. * (10. - 41. * stp[2] + 36. * stp[4]);
509  d0coeff1_5pn[14 + PN15DIM * MINUS] = -2. * (10. - 41. * stp[2] + 36. * stp[4]);
510  d0coeff1_5pn[15 + PN15DIM * PLUS_] = -6. * (10. - 49. * stp[2] + 44. * stp[4]);
511  d0coeff1_5pn[15 + PN15DIM * MINUS] = +6. * (10. - 49. * stp[2] + 44. * stp[4]);
512  d0coeff1_5pn[16 + PN15DIM * PLUS_] = -4. * (7. - 8. * stp[2]);
513  d0coeff1_5pn[16 + PN15DIM * MINUS] = +4. * (7. - 8. * stp[2]);
514  b1coeff1_5pn[0 + PN15DIM * PLUS_] = +0.;
515  b1coeff1_5pn[0 + PN15DIM * MINUS] = +0.;
516  b1coeff1_5pn[1 + PN15DIM * PLUS_] = (-15. * stp[4] + 12. * stp[2] - 2.);
517  b1coeff1_5pn[1 + PN15DIM * MINUS] = (-15. * stp[4] + 12. * stp[2] - 2.);
518  b1coeff1_5pn[2 + PN15DIM * PLUS_] = -5.;
519  b1coeff1_5pn[2 + PN15DIM * MINUS] = +5.;
520  b1coeff1_5pn[3 + PN15DIM * PLUS_] = +0.;
521  b1coeff1_5pn[3 + PN15DIM * MINUS] = +0.;
522  b1coeff1_5pn[4 + PN15DIM * PLUS_] = -(236. - 294. * stp[2] + 21. * stp[4]);
523  b1coeff1_5pn[4 + PN15DIM * MINUS] = +(236. - 294. * stp[2] + 21. * stp[4]);
524  b1coeff1_5pn[5 + PN15DIM * PLUS_] = +3. * (6. - 36. * stp[2] + 45. * stp[4]);
525  b1coeff1_5pn[5 + PN15DIM * MINUS] = +3. * (6. - 36. * stp[2] + 45. * stp[4]);
526  b1coeff1_5pn[6 + PN15DIM * PLUS_] = -3. * (232. - 510. * stp[2] + 243. * stp[4]);
527  b1coeff1_5pn[6 + PN15DIM * MINUS] = +3. * (232. - 510. * stp[2] + 243. * stp[4]);
528  b1coeff1_5pn[7 + PN15DIM * PLUS_] = +9. * (6. - 5. * stp[2]);
529  b1coeff1_5pn[7 + PN15DIM * MINUS] = +9. * (6. - 5. * stp[2]);
530  b1coeff1_5pn[8 + PN15DIM * PLUS_] = +0.;
531  b1coeff1_5pn[8 + PN15DIM * MINUS] = +0.;
532  b1coeff1_5pn[9 + PN15DIM * PLUS_] = +0.;
533  b1coeff1_5pn[9 + PN15DIM * MINUS] = +0.;
534  b1coeff1_5pn[10 + PN15DIM * PLUS_] = +0.;
535  b1coeff1_5pn[10 + PN15DIM * MINUS] = +0.;
536  b1coeff1_5pn[11 + PN15DIM * PLUS_] = (-348. + 591. * stp[2]);
537  b1coeff1_5pn[11 + PN15DIM * MINUS] = (+348. - 591. * stp[2]);
538  b1coeff1_5pn[12 + PN15DIM * PLUS_] = (+273. * stp[2] - 118.);
539  b1coeff1_5pn[12 + PN15DIM * MINUS] = (-273. * stp[2] + 118.);
540  b1coeff1_5pn[13 + PN15DIM * PLUS_] = +0.;
541  b1coeff1_5pn[13 + PN15DIM * MINUS] = +0.;
542  b1coeff1_5pn[14 + PN15DIM * PLUS_] = (2. - 13. * stp[2] + 12. * stp[4]);
543  b1coeff1_5pn[14 + PN15DIM * MINUS] = (2. - 13. * stp[2] + 12. * stp[4]);
544  b1coeff1_5pn[15 + PN15DIM * PLUS_] = -9. * (2. - 13. * stp[2] + 12. * stp[4]);
545  b1coeff1_5pn[15 + PN15DIM * MINUS] = -9. * (2. - 13. * stp[2] + 12. * stp[4]);
546  b1coeff1_5pn[16 + PN15DIM * PLUS_] = -3. * (3. - 4. * stp[2]);
547  b1coeff1_5pn[16 + PN15DIM * MINUS] = -3. * (3. - 4. * stp[2]);
548  d1coeff1_5pn[0 + PN15DIM * PLUS_] = +0.;
549  d1coeff1_5pn[0 + PN15DIM * MINUS] = +0.;
550  d1coeff1_5pn[1 + PN15DIM * PLUS_] = +0.;
551  d1coeff1_5pn[1 + PN15DIM * MINUS] = +0.;
552  d1coeff1_5pn[2 + PN15DIM * PLUS_] = +1.;
553  d1coeff1_5pn[2 + PN15DIM * MINUS] = +1.;
554  d1coeff1_5pn[3 + PN15DIM * PLUS_] = +0.;
555  d1coeff1_5pn[3 + PN15DIM * MINUS] = +0.;
556  d1coeff1_5pn[4 + PN15DIM * PLUS_] = (28. - 126. * stp[2] + 105. * stp[4]);
557  d1coeff1_5pn[4 + PN15DIM * MINUS] = (28. - 126. * stp[2] + 105. * stp[4]);
558  d1coeff1_5pn[5 + PN15DIM * PLUS_] = +0.;
559  d1coeff1_5pn[5 + PN15DIM * MINUS] = +0.;
560  d1coeff1_5pn[6 + PN15DIM * PLUS_] = +27. * (8. - 22. * stp[2] + 15. * stp[4]);
561  d1coeff1_5pn[6 + PN15DIM * MINUS] = +27. * (8. - 22. * stp[2] + 15. * stp[4]);
562  d1coeff1_5pn[7 + PN15DIM * PLUS_] = +0.;
563  d1coeff1_5pn[7 + PN15DIM * MINUS] = +0.;
564  d1coeff1_5pn[8 + PN15DIM * PLUS_] = +0.;
565  d1coeff1_5pn[8 + PN15DIM * MINUS] = +0.;
566  d1coeff1_5pn[9 + PN15DIM * PLUS_] = +0.;
567  d1coeff1_5pn[9 + PN15DIM * MINUS] = +0.;
568  d1coeff1_5pn[10 + PN15DIM * PLUS_] = +0.;
569  d1coeff1_5pn[10 + PN15DIM * MINUS] = +0.;
570  d1coeff1_5pn[11 + PN15DIM * PLUS_] = (-243. * stp[2] + 108.);
571  d1coeff1_5pn[11 + PN15DIM * MINUS] = (-243. * stp[2] + 108.);
572  d1coeff1_5pn[12 + PN15DIM * PLUS_] = +7. * (2. - 3. * stp[2]);
573  d1coeff1_5pn[12 + PN15DIM * MINUS] = +7. * (2. - 3. * stp[2]);
574  d1coeff1_5pn[13 + PN15DIM * PLUS_] = +0.;
575  d1coeff1_5pn[13 + PN15DIM * MINUS] = +0.;
576  d1coeff1_5pn[14 + PN15DIM * PLUS_] = +0.;
577  d1coeff1_5pn[14 + PN15DIM * MINUS] = +0.;
578  d1coeff1_5pn[15 + PN15DIM * PLUS_] = +0.;
579  d1coeff1_5pn[15 + PN15DIM * MINUS] = +0.;
580  d1coeff1_5pn[16 + PN15DIM * PLUS_] = +0.;
581  d1coeff1_5pn[16 + PN15DIM * MINUS] = +0.;
582  for (i = 0; i <= 16; i++) {
583  for (j = 0; j <= 1; j++) {
584  params->ccoeff15pn[i + PN15DIM * j] = acoeff1_5pn[i + PN15DIM * j]
585  + skp[2] * (b0coeff1_5pn[i + PN15DIM * j] + ckp[1] * d0coeff1_5pn[i + PN15DIM * j])
586  + skp[4] * (b1coeff1_5pn[i + PN15DIM * j] + ckp[1] * d1coeff1_5pn[i + PN15DIM * j]);
587  }
588  }
589  XLALFreeDmatrix(acoeff00pn);
590  XLALFreeDmatrix(b0coeff0pn);
591  XLALFreeDmatrix(d0coeff0pn);
592  XLALFreeDmatrix(acoeff0_5pn);
593  XLALFreeDmatrix(b0coeff0_5pn);
594  XLALFreeDmatrix(d0coeff0_5pn);
595  XLALFreeDmatrix(acoeff1pn);
596  XLALFreeDmatrix(b0coeff1pn);
597  XLALFreeDmatrix(d0coeff1pn);
598  XLALFreeDmatrix(b1coeff1pn);
599  XLALFreeDmatrix(d1coeff1pn);
600  XLALFreeDmatrix(acoeff1_5pn);
601  XLALFreeDmatrix(b0coeff1_5pn);
602  XLALFreeDmatrix(d0coeff1_5pn);
603  XLALFreeDmatrix(b1coeff1_5pn);
604  XLALFreeDmatrix(d1coeff1_5pn);
605  return XLAL_SUCCESS;
606 }
607 
608 /**
609  * Function building the wavefrom from the calculated parameters at a given time
610  * For the formulae see the appendix of Arxiv:1209.1722
611  */
612 int XLALSpinDominatedWaveformBuild(LALSDWaveformParams *params, /**< The SDW parameters */
613 REAL8 expr[], /**< The 3 time dependent variables of the waveform at the time indexed by idx */
614 REAL8TimeSeries **hplus, /**< +-polarization waveform */
615 REAL8TimeSeries **hcross, /**< x-polarization waveform */
616 int idx /**< index of the time dependent variables */) {
617  enum {
618  PN00_A, PN00_B, PN05_A, PN05_B, PN10_A, PN10_B, PN10_C, PN10_D, PN15_A, PN15_B, PN15_C,
619  };
620  REAL8 omega0 = 1;
621  REAL8 chi1_1 = params->chi1;
622  REAL8 skp[TRIGONOMETRIC_POWER], ckp[TRIGONOMETRIC_POWER], stp[TRIGONOMETRIC_POWER], ctp[TRIGONOMETRIC_POWER];
623  skp[0] = ckp[0] = stp[0] = ctp[0] = 1.;
624  skp[1] = sin(params->kappa1);
625  ckp[1] = cos(params->kappa1);
626  stp[1] = sin(params->theta);
627  ctp[1] = cos(params->theta);
628  for (int i = 2; i < TRIGONOMETRIC_POWER; ++i) {
629  skp[i] = skp[1] * skp[i - 1];
630  ckp[i] = ckp[1] * ckp[i - 1];
631  stp[i] = stp[1] * stp[i - 1];
632  ctp[i] = ctp[1] * ctp[i - 1];
633  }
634  REAL8 s2t = sin(2. * params->theta);
635  REAL8 c2t = cos(2. * params->theta);
636  REAL8 k[PLUS_MINUS_DIM] = { ckp[1] - 1., ckp[1] + 1. };
637  REAL8 s2k1 = sin(2. * params->kappa1);
638  REAL8 c2k1 = cos(2. * params->kappa1);
639 
640  REAL8 vP[OMEGA_POWER_DIM];
641  vP[0] = 1.;
642  vP[1] = cbrt(params->totalmass * expr[OMEGA] * G_CP2 / LAL_C_SI);
643  for (int i = 2; i < OMEGA_POWER_DIM; ++i) {
644  vP[i] = vP[1] * vP[i - 1];
645  }
646  params->eps = vP[2];
647  REAL8 eta = params->nu / (1. + params->nu) / (1. + params->nu);
648  REAL8 eps_corr[PN_ORDER][PN_ORDER];
649  eps_corr[PN00][PN10] = params->eps
650  * (1. + vP[2] * (1. - eta / 3.) );
651  eps_corr[PN05][PN10] = eps_corr[PN00][PN10];
652  eps_corr[PN00][PN15] = params->eps
653  * (1. + vP[2] * (1. - eta / 3.) + vP[3] * (eta + (2. / 3.) / (eta / params->nu)));
654  REAL8 epssqrt = sqrt(params->eps);
655  params->xi = params->nu / epssqrt; // second small parameter
656 
657  REAL8 phin[PHI_PSI_DIM], psin[PHI_PSI_DIM];
658  phin[0] = psin[0] = 1.;
659  for (int i = 1; i < PHI_PSI_DIM; ++i) {
660  phin[i] = i * expr[PHI];
661  psin[i] = i * expr[PSI];
662 
663  }
664  REAL8 *waveampcoeffs = XLALDmatrix(AMPCOEFF_DIM, PLUS_MINUS_DIM);
665  switch (params->pnamp) {
666  case (PNDEF): //default value -1 contains all corrections
667  case (PN15): {
668  waveampcoeffs[PN15_C + AMPCOEFF_DIM * PLUS_] = LAL_PI_2 * ( //
669  6. * skp[2] * stp[2] * cos(psin[2]) + ( //
670  (stp[2] - 2.) * ( //
671  /* */cos(phin[2] + psin[2]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
672  /**/+ cos(phin[2] - psin[2]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
673  ) - 2. * skp[1] * s2t * ( //
674  /* */sin(phin[1] + psin[2]) * k[PLUS_] //
675  /**/+ sin(phin[1] - psin[2]) * k[MINUS] //
676  )//
677  )//
678  ) + 3. * log(eps_corr[PN00][PN15] / omega0) * ( //
679  4. * ctp[1] * skp[1] * stp[1] * ( //
680  -k[MINUS] * cos(phin[1] - psin[2])
681  + k[PLUS_] * cos(phin[1] + psin[2]) //
682  )//
683  + sin(phin[2] - psin[2]) * (2. * ckp[1] - skp[2] + 2.) * (-stp[2] + 2.) //
684  + sin(phin[2] + psin[2]) * (2. * ckp[1] + skp[2] - 2.) //
685  + 6. * skp[2] * sin(psin[2]) * stp[2] //
686  );
687  waveampcoeffs[PN15_C + AMPCOEFF_DIM * MINUS] = LAL_PI * ( //
688  ctp[1] * ( //
689  /* */sin(phin[2] + psin[2]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
690  /**/+ sin(phin[2] - psin[2]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
691  ) - 2. * skp[1] * stp[1] * ( //
692  /* */cos(phin[1] + psin[2]) * k[PLUS_]
693  /**/+ cos(phin[1] - psin[2]) * k[MINUS] //
694  )//
695  ) + 6. * log(eps_corr[PN00][PN15] / omega0) * ( //
696  (2. * ckp[1] - skp[2] + 2.) * ctp[1] * cos(phin[2] - psin[2])
697  + (2. * ckp[1] + skp[2] - 2.) * ctp[1] * cos(phin[2] + psin[2])
698  + 2. * k[MINUS] * skp[1] * stp[1] * sin(phin[1] - psin[2])
699  - 2. * k[PLUS_] * skp[1] * stp[1] * sin(phin[1] + psin[2]) //
700  );
701 
702  waveampcoeffs[PN15_B + AMPCOEFF_DIM * PLUS_] = chi1_1 / 2. * (4. * skp[1] * ( //
703  /**/ckp[1] * skp[1] * cos(phin[2]) - c2k1 * ctp[1] * stp[1] * sin(phin[1])
704  /**/+ ckp[1] * skp[1] * stp[2] * ( //
705  /* */6. * sin(psin[1]) * sin(psin[1]) - 2. + sin(phin[1]) * sin(phin[1]) //
706  )//
707  ) + ( //
708  2. * ctp[1] * skp[1] * stp[1] * ( //
709  /* */(-3. * k[PLUS_] - 4. * skp[2]) * sin(phin[1] + psin[2])
710  /**/+ (+3. * k[MINUS] - 4. * skp[2]) * sin(phin[1] - psin[2]) //
711  ) + (stp[2] - 2.) * ( //
712  /* */(-2. * k[PLUS_] + (2. * ckp[1] - 3.) * skp[2]) * cos(phin[2] + psin[2])
713  /**/+ (-2. * k[MINUS] + (2. * ckp[1] + 3.) * skp[2]) * cos(phin[2] - psin[2]) //
714  )//
715  ));
716  waveampcoeffs[PN15_B + AMPCOEFF_DIM * MINUS] = chi1_1 * ( //
717  -2. * cos(phin[1]) * skp[1] * (stp[1] * c2k1 + ctp[1] * s2k1 * sin(phin[1])) + ( //
718  ctp[1] * ( //
719  /* */(-2. * k[PLUS_] + (2. * ckp[1] - 3.) * skp[2]) * sin(phin[2] + psin[2])
720  /**/+ (-2. * k[MINUS] + (2. * ckp[1] + 3.) * skp[2]) * sin(phin[2] - psin[2]) //
721  ) + skp[1] * stp[1] * ( //
722  /* */(-3. * k[PLUS_] - 4. * skp[2]) * cos(phin[1] + psin[2])
723  /**/+ (+3. * k[MINUS] - 4. * skp[2]) * cos(phin[1] - psin[2]) //
724  )//
725  )//
726  );
727  waveampcoeffs[PN15_A + AMPCOEFF_DIM * PLUS_] = 1. / 12288. * (12. * ctp[1] * skp[1] * stp[2] * ( //
728  cos(psin[3]) * ( //
729  /**/1701. * (2. - 3. * stp[2]) * skp[4] + 72. * skp[2] * (63. * stp[2] + 178.) //
730  ) + cos(psin[1]) * ( //
731  /**/-14. * (2. - 3. * stp[2]) * skp[4] //
732  /**/- 8. * skp[2] * (7. * stp[2] + 162.) + 16. * (stp[2] + 66.) //
733  ) - 4375. * (2. - 3. * stp[2]) * skp[4] * cos(psin[5]) //
734  ) + ( //
735  2. * (stp[2] - 2.) * skp[4] * stp[3] * ( //
736  /* */sin(phin[5] + psin[1]) * k[PLUS_]
737  /**/+ sin(phin[5] - psin[1]) * k[MINUS] //
738  ) + 4. * ctp[1] * skp[3] * stp[2] * ( //
739  /* */cos(phin[4] + psin[1]) * params->ccoeff15pn[0 + PN15DIM * PLUS_]
740  /**/+ cos(phin[4] - psin[1]) * params->ccoeff15pn[0 + PN15DIM * MINUS] //
741  ) + 16. * ctp[1] * skp[1] * ( //
742  /* */cos(phin[2] + psin[1]) * params->ccoeff15pn[1 + PN15DIM * PLUS_]
743  /**/+ cos(phin[2] - psin[1]) * params->ccoeff15pn[1 + PN15DIM * MINUS] //
744  ) + 1250. * skp[4] * stp[1] * (105. * stp[4] - 126. * stp[2] + 28.) * ( //
745  /* */sin(phin[1] + psin[5]) * k[PLUS_]
746  /**/+ sin(phin[1] - psin[5]) * k[MINUS] //
747  ) + 625. * (stp[2] - 2.) * stp[3] * ( //
748  /* */sin(phin[5] + psin[5]) * params->ccoeff15pn[2 + PN15DIM * PLUS_]
749  /**/+ sin(phin[5] - psin[5]) * params->ccoeff15pn[2 + PN15DIM * MINUS] //
750  ) + 6. * skp[2] * stp[1] * ( //
751  /* */sin(phin[3] + psin[1]) * params->ccoeff15pn[3 + PN15DIM * PLUS_]
752  /**/+ sin(phin[3] - psin[1]) * params->ccoeff15pn[3 + PN15DIM * MINUS] //
753  ) + 243. * (stp[2] - 2.) * skp[2] * stp[3] * ( //
754  /* */sin(phin[5] + psin[3]) * params->ccoeff05pn[1 + PN05DIM * PLUS_]
755  /**/+ sin(phin[5] - psin[3]) * params->ccoeff05pn[1 + PN05DIM * MINUS] //
756  ) + 4. * stp[1] * ( //
757  /* */sin(phin[1] + psin[1]) * params->ccoeff15pn[4 + PN15DIM * PLUS_]
758  /**/+ sin(phin[1] - psin[1]) * params->ccoeff15pn[4 + PN15DIM * MINUS] //
759  ) + 5000. * ctp[1] * skp[3] * (15. * stp[4] - 12. * stp[2] + 2.) * ( //
760  /* */cos(phin[2] + psin[5]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
761  /**/+ cos(phin[2] - psin[5]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
762  ) - 1250. * ctp[1] * skp[1] * stp[2] * (5. * stp[2] - 6.) * ( //
763  /* */cos(phin[4] + psin[5]) * params->ccoeff10pn[0 + PN10DIM * PLUS_]
764  /**/+ cos(phin[4] - psin[5]) * params->ccoeff10pn[0 + PN10DIM * MINUS] //
765  ) + 1875. * skp[2] * stp[1] * (8. - 22. * stp[2] + 15. * stp[4]) * ( //
766  /* */sin(phin[3] + psin[5]) * params->ccoeff05pn[1 + PN05DIM * PLUS_]
767  /**/+ sin(phin[3] - psin[5]) * params->ccoeff05pn[1 + PN05DIM * MINUS] //
768  ) + 216. * ctp[1] * skp[1] * ( //
769  /* */cos(phin[2] + psin[3]) * params->ccoeff15pn[5 + PN15DIM * PLUS_]
770  /**/+ cos(phin[2] - psin[3]) * params->ccoeff15pn[5 + PN15DIM * MINUS] //
771  ) + 27. * stp[1] * ( //
772  /* */sin(phin[3] + psin[3]) * params->ccoeff15pn[6 + PN15DIM * PLUS_]
773  /**/+ sin(phin[3] - psin[3]) * params->ccoeff15pn[6 + PN15DIM * MINUS] //
774  ) + 54. * ctp[1] * skp[1] * stp[2] * ( //
775  /* */cos(phin[4] + psin[3]) * params->ccoeff15pn[7 + PN15DIM * PLUS_]
776  /**/+ cos(phin[4] - psin[3]) * params->ccoeff15pn[7 + PN15DIM * MINUS] //
777  ) + 54. * skp[2] * stp[1] * ( //
778  /* */sin(phin[1] + psin[3]) * params->ccoeff15pn[8 + PN15DIM * PLUS_]
779  /**/+ sin(phin[1] - psin[3]) * params->ccoeff15pn[8 + PN15DIM * MINUS] //
780  )//
781  ));
782  waveampcoeffs[PN15_A + AMPCOEFF_DIM * MINUS] = 1. / 6144. * (192. * ckp[1] * skp[1] * stp[2] * ( //
783  /* */sin(psin[1]) * (64. - skp[2] * (7. * stp[2] - 6.) + 4. * stp[2])
784  /**/+ 27. * sin(psin[3]) * skp[2] * (7. * stp[2] - 6.) //
785  ) + ( //
786  4. * skp[3] * stp[2] * ( //
787  /* */sin(phin[4] + psin[1]) * params->ccoeff15pn[9 + PN15DIM * PLUS_]
788  /**/+ sin(phin[4] - psin[1]) * params->ccoeff15pn[9 + PN15DIM * MINUS] //
789  ) - 2. * ctp[1] * skp[4] * stp[3] * ( //
790  /* */cos(phin[5] + psin[1]) * k[PLUS_]
791  /**/+ cos(phin[5] - psin[1]) * k[MINUS] //
792  ) - 243. * ctp[1] * skp[2] * stp[3] * ( //
793  /* */cos(phin[5] + psin[3]) * params->ccoeff05pn[1 + PN05DIM * PLUS_]
794  /**/+ cos(phin[5] - psin[3]) * params->ccoeff05pn[1 + PN05DIM * MINUS] //
795  ) - 625. * ctp[1] * stp[3] * ( //
796  /* */cos(phin[5] + psin[5]) * params->ccoeff15pn[2 + PN15DIM * PLUS_]
797  /**/+ cos(phin[5] - psin[5]) * params->ccoeff15pn[2 + PN15DIM * MINUS] //
798  ) + 3. * skp[2] * s2t * ( //
799  /* */cos(phin[3] + psin[1]) * params->ccoeff15pn[10 + PN15DIM * PLUS_]
800  /**/+ cos(phin[3] - psin[1]) * params->ccoeff15pn[10 + PN15DIM * MINUS] //
801  ) + 27. * ctp[1] * stp[1] * ( //
802  /* */cos(phin[3] + psin[3]) * params->ccoeff15pn[11 + PN15DIM * PLUS_]
803  /**/+ cos(phin[3] - psin[3]) * params->ccoeff15pn[11 + PN15DIM * MINUS] //
804  ) + 1875. * ctp[1] * skp[2] * stp[1] * (4. - 9. * stp[2]) * ( //
805  /* */cos(phin[3] + psin[5]) * params->ccoeff05pn[1 + PN05DIM * PLUS_]
806  /**/+ cos(phin[3] - psin[5]) * params->ccoeff05pn[1 + PN05DIM * MINUS] //
807  ) + 2. * s2t * ( //
808  /* */cos(phin[1] + psin[1]) * params->ccoeff15pn[12 + PN15DIM * PLUS_]
809  /**/+ cos(phin[1] - psin[1]) * params->ccoeff15pn[12 + PN15DIM * MINUS] //
810  ) + 27. * skp[2] * s2t * ( //
811  /* */cos(phin[1] + psin[3]) * params->ccoeff15pn[13 + PN15DIM * PLUS_]
812  /**/+ cos(phin[1] - psin[3]) * params->ccoeff15pn[13 + PN15DIM * MINUS] //
813  ) + 8. * skp[1] * ( //
814  /* */sin(phin[2] + psin[1]) * params->ccoeff15pn[14 + PN15DIM * PLUS_]
815  /**/+ sin(phin[2] - psin[1]) * params->ccoeff15pn[14 + PN15DIM * MINUS] //
816  ) + 4375. * (2. - 3. * stp[2]) * skp[4] * s2t * ( //
817  /* */cos(phin[1] + psin[5]) * k[PLUS_]
818  /**/+ cos(phin[1] - psin[5]) * k[MINUS] //
819  ) + 108. * skp[1] * ( //
820  /* */sin(phin[2] + psin[3]) * params->ccoeff15pn[15 + PN15DIM * PLUS_]
821  /**/+ sin(phin[2] - psin[3]) * params->ccoeff15pn[15 + PN15DIM * MINUS] //
822  ) + 162. * skp[1] * stp[2] * ( //
823  /* */sin(phin[4] + psin[3]) * params->ccoeff15pn[16 + PN15DIM * PLUS_]
824  /**/+ sin(phin[4] - psin[3]) * params->ccoeff15pn[16 + PN15DIM * MINUS] //
825  ) - 2500. * skp[3] * (2. - 13. * stp[2] + 12. * stp[4]) * ( //
826  /* */sin(phin[2] + psin[5]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
827  /**/+ sin(phin[2] - psin[5]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
828  ) - 1250. * skp[1] * stp[2] * (3. - 4. * stp[2]) * ( //
829  /* */sin(phin[4] + psin[5]) * params->ccoeff10pn[0 + PN10DIM * PLUS_]
830  /**/+ sin(phin[4] - psin[5]) * params->ccoeff10pn[0 + PN10DIM * MINUS] //
831  )//
832  ));
833  }
834 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
835  __attribute__ ((fallthrough)); /* no break */
836 #endif
837  case (PN10): {
838  REAL8 phi1_1 = LAL_PI_2;
839  waveampcoeffs[PN10_D + AMPCOEFF_DIM * PLUS_] = chi1_1 / 2. * ( //
840  ( //
841  /* */-k[PLUS_] * ctp[1] * sin(phin[2] + psin[1] - phi1_1)
842  /**/+ k[MINUS] * ctp[1] * sin(phin[2] - psin[1] - phi1_1) //
843  ) + skp[1] * stp[1] * ( //
844  /* */sin(phin[1] + psin[1]) - sin(phin[1] - psin[1])
845  /**/+ cos(phin[1] + psin[1] - phi1_1) - cos(phin[1] - psin[1] - phi1_1) //
846  )//
847  );
848  waveampcoeffs[PN10_D + AMPCOEFF_DIM * MINUS] = chi1_1 / 4. * (-stp[2] * ( //
849  2. * (sin(phi1_1) + 2.) * ckp[1] * sin(psin[1]) + 2. * cos(phi1_1) * cos(psin[1]) //
850  ) + ( //
851  2. * ctp[1] * skp[1] * stp[1] * ( //
852  /**/cos(phin[1] + psin[1]) - cos(phin[1] - psin[1])
853  /**/- sin(phin[1] + psin[1] - phi1_1) + sin(phin[1] - psin[1] - phi1_1) //
854  ) + (stp[2] - 2.) * ( //
855  /* */k[PLUS_] * cos(phin[2] + psin[1] - phi1_1)
856  /**/- k[MINUS] * cos(phin[2] - psin[1] - phi1_1) //
857  )//
858  ));
859  waveampcoeffs[PN10_C + AMPCOEFF_DIM * PLUS_] = chi1_1 / 2. * stp[1]
860  * (k[PLUS_] * sin(phin[1] + psin[1]) - k[MINUS] * sin(phin[1] - psin[1]));
861  waveampcoeffs[PN10_C + AMPCOEFF_DIM * MINUS] = chi1_1 / 2. * stp[1]
862  * (2. * skp[1] * sin(psin[1]) * stp[1]
863  + ctp[1] * (k[PLUS_] * cos(phin[1] + psin[1]) - k[MINUS] * cos(phin[1] - psin[1])));
864  waveampcoeffs[PN10_B + AMPCOEFF_DIM * PLUS_] = 1. / 24. * ( //
865  ckp[1] * skp[1] * stp[2] * ( //
866  4. * (15. * stp[2] + 51.) * cos(psin[2])
867  + 20. * skp[2] * (7. * stp[2] - 6.) * (4. * cos(psin[4]) - cos(psin[2])) //
868  ) + ( //
869  s2t * ( //
870  /* */sin(phin[3] + psin[2]) * params->ccoeff10pn[7 + PN10DIM * PLUS_]
871  /**/+ sin(phin[3] - psin[2]) * params->ccoeff10pn[7 + PN10DIM * MINUS] //
872  ) + s2t * ( //
873  /* */sin(phin[1] + psin[2]) * params->ccoeff10pn[8 + PN10DIM * PLUS_]
874  /**/+ sin(phin[1] - psin[2]) * params->ccoeff10pn[8 + PN10DIM * MINUS] //
875  ) + ctp[2] * s2t * 8. * ( //
876  /* */sin(phin[3] + psin[4]) * params->ccoeff10pn[9 + PN10DIM * PLUS_]
877  /**/+ sin(phin[3] - psin[4]) * params->ccoeff10pn[9 + PN10DIM * MINUS] //
878  ) + 2. * skp[1] * ( //
879  /* */cos(phin[2] + psin[2]) * params->ccoeff10pn[10 + PN10DIM * PLUS_]
880  /**/+ cos(phin[2] - psin[2]) * params->ccoeff10pn[10 + PN10DIM * MINUS] //
881  ) + 8. * skp[2] * s2t * (7. * stp[2] - 3.) * ( //
882  /* */sin(phin[1] + psin[4]) * (+3. * k[PLUS_] + 4. * skp[2])
883  /**/+ sin(phin[1] - psin[4]) * (-3. * k[MINUS] + 4. * skp[2]) //
884  ) + 16. * skp[1] * (2. - 8. * stp[2] + 7. * stp[4]) * ( //
885  /* */cos(phin[2] + psin[4]) * params->ccoeff10pn[11 + PN10DIM * PLUS_]
886  /**/+ cos(phin[2] - psin[4]) * params->ccoeff10pn[11 + PN10DIM * MINUS] //
887  ) + skp[1] * stp[2] * (stp[2] - 2.) * ( //
888  /* */cos(phin[4] + psin[2]) * params->ccoeff10pn[11 + PN10DIM * PLUS_]
889  /**/+ cos(phin[4] - psin[2]) * params->ccoeff10pn[11 + PN10DIM * MINUS] - 8. * ( //
890  /* */cos(phin[4] + psin[4]) * params->ccoeff05pn[1 + PN05DIM * PLUS_]
891  /**/+ cos(phin[4] - psin[4]) * params->ccoeff05pn[1 + PN05DIM * MINUS] //
892  )//
893  )//
894  )//
895  );
896  waveampcoeffs[PN10_B + AMPCOEFF_DIM * MINUS] = 1. / 12. * ( //
897  30. * ctp[1] * skp[1] * sin(psin[2]) * stp[2] * (3. * skp[2] - 2.) + (stp[1] * ( //
898  /* */cos(phin[1] + psin[2]) * params->ccoeff10pn[12 + PN10DIM * PLUS_]
899  /**/+ cos(phin[1] - psin[2]) * params->ccoeff10pn[12 + PN10DIM * MINUS] //
900  ) + 2. * ctp[1] * skp[1] * ( //
901  ( //
902  /* */sin(phin[2] + psin[2]) * params->ccoeff10pn[13 + PN10DIM * PLUS_]
903  /**/+ sin(phin[2] - psin[2]) * params->ccoeff10pn[13 + PN10DIM * MINUS] //
904  ) + 4. * (7. * stp[2] - 2.) * ( //
905  /* */sin(phin[2] + psin[4]) * params->ccoeff10pn[11 + PN10DIM * PLUS_]
906  /**/+ sin(phin[2] - psin[4]) * params->ccoeff10pn[11 + PN10DIM * MINUS] //
907  )//
908  ) + ctp[1] * skp[1] * stp[2] * ( //
909  ( //
910  /* */sin(phin[4] + psin[2]) * params->ccoeff10pn[11 + PN10DIM * PLUS_]
911  /**/+ sin(phin[4] - psin[2]) * params->ccoeff10pn[11 + PN10DIM * MINUS] //
912  ) - 8. * ( //
913  /* */sin(phin[4] + psin[4]) * params->ccoeff05pn[1 + PN05DIM * PLUS_]
914  /**/+ sin(phin[4] - psin[4]) * params->ccoeff05pn[1 + PN05DIM * MINUS] //
915  )//
916  ) + stp[1] * ( //
917  ( //
918  /* */cos(phin[3] + psin[2]) * params->ccoeff10pn[14 + PN10DIM * PLUS_]
919  /**/+ cos(phin[3] - psin[2]) * params->ccoeff10pn[14 + PN10DIM * MINUS] //
920  ) + 4. * (2. - 3. * stp[2]) * ( //
921  /* */cos(phin[3] + psin[4]) * params->ccoeff10pn[9 + PN10DIM * PLUS_]
922  /**/+ cos(phin[3] - psin[4]) * params->ccoeff10pn[9 + PN10DIM * MINUS] //
923  )//
924  ) + 4. * skp[2] * stp[1] * (6. - 7. * stp[2]) * ( //
925  /* */(-3. * k[PLUS_] - 4. * skp[2]) * cos(phin[1] + psin[4])
926  /**/+ (+3. * k[MINUS] - 4. * skp[2]) * cos(phin[1] - psin[4]) //
927  ))//
928  );
929  waveampcoeffs[PN10_A + AMPCOEFF_DIM * PLUS_] = 1. / 48. * (2. * skp[2] * stp[2] * ( //
930  /* */5. * skp[2] * (7. * stp[2] - 6.) * (cos(psin[2]) - 4. * cos(psin[4]))
931  /**/- 2. * (15. * stp[2] + 51.) * cos(psin[2]) //
932  ) + ( //
933  /* */16. * skp[3] * s2t * (7. * stp[2] - 3.) * ( //
934  /* */sin(phin[1] + psin[4]) * k[PLUS_]
935  /**/+ sin(phin[1] - psin[4]) * k[MINUS] //
936  ) - (stp[2] - 2.) * skp[2] * stp[2] * ( //
937  /* */cos(phin[4] + psin[2]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
938  /**/+ cos(phin[4] - psin[2]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
939  ) + 4. * (stp[2] - 2.) * stp[2] * ( //
940  /* */cos(phin[4] + psin[4]) * params->ccoeff10pn[0 + PN10DIM * PLUS_]
941  /**/+ cos(phin[4] - psin[4]) * params->ccoeff10pn[0 + PN10DIM * MINUS] //
942  ) + 2. * skp[1] * s2t * ( //
943  ( //
944  /* */sin(phin[1] + psin[2]) * params->ccoeff10pn[1 + PN10DIM * PLUS_]
945  /**/+ sin(phin[1] - psin[2]) * params->ccoeff10pn[1 + PN10DIM * MINUS] //
946  ) + (
947  /* */sin(phin[3] + psin[2]) * params->ccoeff10pn[2 + PN10DIM * PLUS_]
948  /**/+ sin(phin[3] - psin[2]) * params->ccoeff10pn[2 + PN10DIM * MINUS] //
949  ) - 8. * ctp[2] * (
950  /* */sin(phin[3] + psin[4]) * params->ccoeff05pn[1 + PN05DIM * PLUS_]
951  /**/+ sin(phin[3] - psin[4]) * params->ccoeff05pn[1 + PN05DIM * MINUS] //
952  )//
953  ) + 2. * ( //
954  /* */cos(phin[2] + psin[2]) * params->ccoeff10pn[3 + PN10DIM * PLUS_]
955  /**/+ cos(phin[2] - psin[2]) * params->ccoeff10pn[3 + PN10DIM * MINUS] //
956  ) - 16. * skp[2] * (7. * stp[4] - 2. * (4. * stp[2] - 1.)) * ( //
957  /* */cos(phin[2] + psin[4]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
958  /**/+ cos(phin[2] - psin[4]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
959  )//
960  ));
961  waveampcoeffs[PN10_A + AMPCOEFF_DIM * MINUS] = 1. / 24. * ( //
962  60. * ckp[1] * ctp[1] * skp[2] * sin(psin[2]) * stp[2] + ( //
963  8. * skp[3] * stp[1] * (7. * stp[2] - 6.) * ( //
964  /* */cos(phin[1] + psin[4]) * k[PLUS_]
965  /**/+ cos(phin[1] - psin[4]) * k[MINUS] //
966  ) - ctp[1] * skp[2] * stp[2] * ( //
967  /* */sin(phin[4] + psin[2]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
968  /**/+ sin(phin[4] - psin[2]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
969  ) + 2. * skp[1] * stp[1] * ( //
970  /* */cos(phin[1] + psin[2]) * params->ccoeff10pn[4 + PN10DIM * PLUS_]
971  /**/+ cos(phin[1] - psin[2]) * params->ccoeff10pn[4 + PN10DIM * MINUS] //
972  ) + 4. * ctp[1] * stp[2] * ( //
973  /* */sin(phin[4] + psin[4]) * params->ccoeff10pn[0 + PN10DIM * PLUS_]
974  /**/+ sin(phin[4] - psin[4]) * params->ccoeff10pn[0 + PN10DIM * MINUS] //
975  ) + 2. * ctp[1] * ( //
976  /* */sin(phin[2] + psin[2]) * params->ccoeff10pn[5 + PN10DIM * PLUS_]
977  /**/+ sin(phin[2] - psin[2]) * params->ccoeff10pn[5 + PN10DIM * MINUS] //
978  ) - 8. * ctp[1] * skp[2] * (7. * stp[2] - 2.) * ( //
979  /* */sin(phin[2] + psin[4]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
980  /**/+ sin(phin[2] - psin[4]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
981  ) - 8. * (2. - 3. * stp[2]) * skp[1] * stp[1] * ( //
982  /* */cos(phin[3] + psin[4]) * params->ccoeff05pn[1 + PN05DIM * PLUS_]
983  /**/+ cos(phin[3] - psin[4]) * params->ccoeff05pn[1 + PN05DIM * MINUS] //
984  ) + 2. * skp[1] * stp[1] * ( //
985  /* */cos(phin[3] + psin[2]) * params->ccoeff10pn[6 + PN10DIM * PLUS_]
986  /**/+ cos(phin[3] - psin[2]) * params->ccoeff10pn[6 + PN10DIM * MINUS] //
987  )//
988  )//
989  );
990  }
991 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
992  __attribute__ ((fallthrough)); /* no break */
993 #endif
994  case (PN05): {
995  waveampcoeffs[PN05_B + AMPCOEFF_DIM * PLUS_] = 1. / 64. * ( //
996  4. * ckp[1] * ctp[1] * stp[2] * (135. * cos(psin[3]) * skp[2] + cos(psin[1]) * (4. - 15. * skp[2]) //
997  ) + ( //
998  2. * ctp[1] * ( //
999  9. * (2. - 3. * stp[2]) * ( //
1000  /* */cos(phin[2] + psin[3]) * params->ccoeff05pn[5 + PN05DIM * PLUS_]
1001  /**/+ cos(phin[2] - psin[3]) * params->ccoeff05pn[5 + PN05DIM * MINUS] //
1002  )//
1003  /* */+ cos(phin[2] + psin[1]) * params->ccoeff05pn[6 + PN05DIM * PLUS_]
1004  + cos(phin[2] - psin[1]) * params->ccoeff05pn[6 + PN05DIM * MINUS] //
1005  ) - skp[1] * stp[1] * (stp[2] - 2.) * ( //
1006  27. * ( //
1007  /* */sin(phin[3] + psin[3]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
1008  /**/+ sin(phin[3] - psin[3]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
1009  )//
1010  /**/+ sin(phin[3] + psin[1]) * params->ccoeff05pn[7 + PN05DIM * PLUS_]
1011  /**/+ sin(phin[3] - psin[1]) * params->ccoeff05pn[7 + PN05DIM * MINUS] //
1012  ) + skp[1] * stp[1] * ( //
1013  45. * (2. - 3. * stp[2]) * ( //
1014  /* */sin(phin[1] + psin[3]) * params->ccoeff05pn[7 + PN05DIM * PLUS_]
1015  /**/+ sin(phin[1] - psin[3]) * params->ccoeff05pn[7 + PN05DIM * MINUS] //
1016  )//
1017  /**/+ sin(phin[1] + psin[1]) * params->ccoeff05pn[8 + PN05DIM * PLUS_]
1018  /**/+ sin(phin[1] - psin[1]) * params->ccoeff05pn[8 + PN05DIM * MINUS] //
1019  ))//
1020  );
1021  waveampcoeffs[PN05_B + AMPCOEFF_DIM * MINUS] = 1. / 32. * (32. * sin(psin[1]) * stp[2] * c2k1 + ( //
1022  ctp[1] * skp[1] * stp[1] * (27. * ( //
1023  /* */cos(phin[3] + psin[3]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
1024  /**/+ cos(phin[3] - psin[3]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
1025  ) + ( //
1026  /**/cos(phin[3] + psin[1]) + 45. * cos(phin[1] + psin[3]) //
1027  ) * params->ccoeff05pn[7 + PN05DIM * PLUS_] + ( //
1028  /**/cos(phin[3] - psin[1]) + 45. * cos(phin[1] - psin[3]) //
1029  ) * params->ccoeff05pn[7 + PN05DIM * MINUS] + ( //
1030  /* */params->ccoeff05pn[9 + PN05DIM * PLUS_] * cos(phin[1] + psin[1])
1031  /**/+ params->ccoeff05pn[9 + PN05DIM * MINUS] * cos(phin[1] - psin[1]) //
1032  )) - (18. * c2t) * ( //
1033  /* */sin(phin[2] + psin[3]) * params->ccoeff05pn[5 + PN05DIM * PLUS_]
1034  /**/+ sin(phin[2] - psin[3]) * params->ccoeff05pn[5 + PN05DIM * MINUS] //
1035  ) - 2. * ( //
1036  /* */sin(phin[2] + psin[1]) * params->ccoeff05pn[10 + PN05DIM * PLUS_]
1037  /**/+ sin(phin[2] - psin[1]) * params->ccoeff05pn[10 + PN05DIM * MINUS] //
1038  )//
1039  ));
1040  waveampcoeffs[PN05_A + AMPCOEFF_DIM * PLUS_] = 1. / 64. * (4. * ctp[1] * skp[1] * stp[2] * ( //
1041  /**/-45. * skp[2] * cos(psin[3]) + cos(psin[1]) * (5. * skp[2] - 4.) //
1042  ) - skp[2] * stp[1] * ((stp[2] - 2.) * ( //
1043  /* */sin(phin[3] + psin[1]) * k[PLUS_]
1044  /**/+ sin(phin[3] - psin[1]) * k[MINUS] //
1045  ) - 45. * (2. - 3. * stp[2]) * ( //
1046  /* */sin(phin[1] + psin[3]) * k[PLUS_]
1047  /**/+ sin(phin[1] - psin[3]) * k[MINUS] //
1048  )) + stp[1] * ( //
1049  /* */sin(phin[1] + psin[1]) * params->ccoeff05pn[0 + PN05DIM * PLUS_]
1050  /**/+ sin(phin[1] - psin[1]) * params->ccoeff05pn[0 + PN05DIM * MINUS] //
1051  - 9. * (stp[2] - 2.) * ( //
1052  /* */sin(phin[3] + psin[3]) * params->ccoeff05pn[1 + PN05DIM * PLUS_]
1053  /**/+ sin(phin[3] - psin[3]) * params->ccoeff05pn[1 + PN05DIM * MINUS] //
1054  )//
1055  ) + 2. * ctp[1] * skp[1] * ( //
1056  /* */cos(phin[2] + psin[1]) * params->ccoeff05pn[2 + PN05DIM * PLUS_]
1057  /**/+ cos(phin[2] - psin[1]) * params->ccoeff05pn[2 + PN05DIM * MINUS] //
1058  + 9. * (2. - 3. * stp[2]) * ( //
1059  /* */cos(phin[2] + psin[3]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
1060  /**/+ cos(phin[2] - psin[3]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
1061  )//
1062  ));
1063  waveampcoeffs[PN05_A + AMPCOEFF_DIM * MINUS] = 1. / 32. * ( //
1064  -16. * s2k1 * stp[2] * sin(psin[1]) + ctp[1] * stp[1] * skp[2] * ( //
1065  /* */cos(phin[3] + psin[1]) * k[PLUS_]
1066  /**/+ cos(phin[3] - psin[1]) * k[MINUS] //
1067  + 45. * ( //
1068  /* */cos(phin[1] + psin[3]) * k[PLUS_]
1069  /**/+ cos(phin[1] - psin[3]) * k[MINUS] //
1070  )//
1071  ) + 0.5 * s2t * ( //
1072  /* */cos(phin[1] + psin[1]) * params->ccoeff05pn[3 + PN05DIM * PLUS_]
1073  /**/+ cos(phin[1] - psin[1]) * params->ccoeff05pn[3 + PN05DIM * MINUS] //
1074  + 9. * ( //
1075  /* */cos(phin[3] + psin[3]) * params->ccoeff05pn[1 + PN05DIM * PLUS_]
1076  /**/+ cos(phin[3] - psin[3]) * params->ccoeff05pn[1 + PN05DIM * MINUS] //
1077  )//
1078  ) + 2. * skp[1] * ( //
1079  /* */sin(phin[2] + psin[1]) * params->ccoeff05pn[4 + PN05DIM * PLUS_]
1080  /**/+ sin(phin[2] - psin[1]) * params->ccoeff05pn[4 + PN05DIM * MINUS] //
1081  - 9. * c2t * ( //
1082  /* */sin(phin[2] + psin[3]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
1083  /**/+ sin(phin[2] - psin[3]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
1084  )//
1085  )//
1086  );
1087  }
1088 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
1089  __attribute__ ((fallthrough)); /* no break */
1090 #endif
1091  case (PN00): {
1092  waveampcoeffs[PN00_B + AMPCOEFF_DIM * PLUS_] = 0.5 * (s2t * ( //
1093  /* */sin(phin[1] + psin[2]) * params->ccoeff00pn[1 + PN00DIM * PLUS_]
1094  /**/+ sin(phin[1] - psin[2]) * params->ccoeff00pn[1 + PN00DIM * MINUS] //
1095  ) + skp[1] * (stp[2] - 2.) * ( //
1096  /* */cos(phin[2] + psin[2]) * k[PLUS_]
1097  /**/+ cos(phin[2] - psin[2]) * k[MINUS] //
1098  ) - 3. * s2k1 * stp[2] * cos(psin[2]));
1099  waveampcoeffs[PN00_B + AMPCOEFF_DIM * MINUS] = ctp[1] * skp[1] * ( //
1100  /* */sin(phin[2] + psin[2]) * k[PLUS_]
1101  /**/+ sin(phin[2] - psin[2]) * k[MINUS] //
1102  ) + stp[1] * ( //
1103  /* */cos(phin[1] + psin[2]) * params->ccoeff00pn[1 + PN00DIM * PLUS_]
1104  /**/+ cos(phin[1] - psin[2]) * params->ccoeff00pn[1 + PN00DIM * MINUS]);
1105 
1106  waveampcoeffs[PN00_A + AMPCOEFF_DIM * PLUS_] = 0.25 * ((stp[2] - 2.) * ( //
1107  /* */cos(phin[2] + psin[2]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
1108  /**/+ cos(phin[2] - psin[2]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
1109  ) - 2. * skp[1] * s2t * ( //
1110  /* */sin(phin[1] + psin[2]) * k[PLUS_]
1111  /**/+ sin(phin[1] - psin[2]) * k[MINUS] //
1112  ) + 6. * skp[2] * stp[2] * cos(psin[2]));
1113  waveampcoeffs[PN00_A + AMPCOEFF_DIM * MINUS] = 0.5 * (ctp[1] * ( //
1114  /* */sin(phin[2] + psin[2]) * params->ccoeff00pn[0 + PN00DIM * PLUS_]
1115  /**/+ sin(phin[2] - psin[2]) * params->ccoeff00pn[0 + PN00DIM * MINUS] //
1116  ) - 2. * stp[1] * skp[1] * ( //
1117  /* */cos(phin[1] + psin[2]) * k[PLUS_]
1118  /**/+ cos(phin[1] - psin[2]) * k[MINUS]));
1119 
1120  }
1121  }
1122 
1123 
1124  REAL8 h[PLUS_CROSS_DIM] = { 0., 0. };
1125  switch (params->pnamp) {
1126  case (PNDEF): //default value -1 contains all corrections
1127  case (PN15):
1128  // Highest order, only leading order of eps(omega) is needed.
1129  for (int i = PLUS_; i < PLUS_CROSS_DIM; ++i) {
1130  h[i] += params->eps * epssqrt// eps_corr[PN00][PN15] * sqrt(eps_corr[PN00][PN15])
1131  * (waveampcoeffs[PN15_A + AMPCOEFF_DIM * i] + waveampcoeffs[PN15_B + AMPCOEFF_DIM * i]
1132  + waveampcoeffs[PN15_C + AMPCOEFF_DIM * i]);
1133  }
1134 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
1135  __attribute__ ((fallthrough)); /* no break */
1136 #endif
1137  case (PN10):
1138  // Since highest order is 1.5 PN and there is no 0.5 PN order correction to eps(omega), leading order eps is enough.
1139  for (int i = PLUS_; i < PLUS_CROSS_DIM; ++i) {
1140  h[i] += params->eps // eps_corr[PN00][PN10]
1141  * (4. * params->xi * waveampcoeffs[PN05_A + AMPCOEFF_DIM * i]
1142  + waveampcoeffs[PN10_A + AMPCOEFF_DIM * i] + waveampcoeffs[PN10_C + AMPCOEFF_DIM * i]
1143  + params->beta1 * waveampcoeffs[PN10_B + AMPCOEFF_DIM * i]
1144  + params->beta1 * waveampcoeffs[PN10_D + AMPCOEFF_DIM * i]);
1145  }
1146 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
1147  __attribute__ ((fallthrough)); /* no break */
1148 #endif
1149  case (PN05):
1150  //The 0.5 PN correction needs to include 1 PN correction of eps(omega) the amplitude is taken to 1.5 PN order
1151  for (int i = PLUS_; i < PLUS_CROSS_DIM; ++i) {
1152  REAL8 temp = waveampcoeffs[PN05_A + AMPCOEFF_DIM * i]
1153  + params->beta1 * waveampcoeffs[PN05_B + AMPCOEFF_DIM * i]
1154  - 2. * params->xi * waveampcoeffs[PN00_A + AMPCOEFF_DIM * i];
1155  if (params->pnamp == PN15) {
1156  h[i] += 1. / (params->eps * epssqrt) * (eps_corr[PN05][PN10] * sqrt(eps_corr[PN05][PN10])) * epssqrt * temp;
1157  } else {
1158  h[i] +=epssqrt * temp;
1159  }
1160  }
1161 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
1162  __attribute__ ((fallthrough)); /* no break */
1163 #endif
1164  case (PN00):
1165  // If the amplitude is taken to 1 PN order, the eps(omega) needs to include the 1 PN correction, if amplitude is taken to 1.5 PN order, that eps(omega) needs to include corrections up to 1.5 PN order
1166  for (int i = PLUS_; i < PLUS_CROSS_DIM; ++i) {
1167  REAL8 temp = waveampcoeffs[PN00_A + AMPCOEFF_DIM * i]
1168  + params->beta1 * waveampcoeffs[PN00_B + AMPCOEFF_DIM * i];
1169  if (params->pnamp == PN15) {
1170  h[i] += 1. / (params->eps * epssqrt) * eps_corr[PN00][PN15] * sqrt(eps_corr[PN00][PN15]) * temp;
1171  } else if (params->pnamp == PN10) {
1172  h[i] += 1. / (params->eps * epssqrt) * eps_corr[PN00][PN10] * sqrt(eps_corr[PN00][PN10]) * temp;
1173  } else {
1174  h[i] += /* */temp;
1175  }
1176  }
1177  break;
1178  }
1179  REAL8 ampcoeff = 2. * G_CP2 * params->totalmass * params->eps * epssqrt * params->xi / params->dist;
1180  // rotating h+ and hx into the LALSimulation precessing radiation frame
1181  (*hplus)->data->data[idx] = ampcoeff * (h[PLUS_]*cos(params->polarizationangle)-h[MINUS]*sin(params->polarizationangle));
1182  (*hcross)->data->data[idx] = ampcoeff * (h[MINUS]*cos(params->polarizationangle)+h[PLUS_]*sin(params->polarizationangle));
1183  XLALFreeDmatrix(waveampcoeffs);
1184  return XLAL_SUCCESS;
1185 }
1186 
1187 /**
1188  * Interface routine, calculating the prefered variables for the Spin-dominated waveforms
1189  */
1190 int XLALSimInspiralSpinDominatedWaveformInterfaceTD(REAL8TimeSeries **hplus, /**< +-polarization waveform */
1191 REAL8TimeSeries **hcross, /**< x-polarization waveform */
1192 REAL8 deltaT, /**< sampling interval (s) */
1193 REAL8 m1, /**< mass of companion 1 (kg) */
1194 REAL8 m2, /**< mass of companion 2 (kg) */
1195 REAL8 fStart, /**< start GW frequency (Hz) */
1196 REAL8 fRef, /**< end GW frequency (Hz) */
1197 REAL8 D, /**< distance of source (m) */
1198 REAL8 s1x, /**< initial value of S1x */
1199 REAL8 s1y, /**< initial value of S1y */
1200 REAL8 s1z, /**< initial value of S1z */
1201 REAL8 lnhatx, /**< initial value of LNhatx */
1202 REAL8 lnhaty, /**< initial value of LNhaty */
1203 REAL8 lnhatz, /**< initial value of LNhatz */
1204 REAL8 incl, /**< inclination, angle between L and line of sight N */
1205 int phaseO, /**< twice PN phase order */
1206 int amplitudeO, /**< twice PN amplitude order */
1207 REAL8 phiRef /**< Reference phase at the Reference Frequency */
1208 ) {
1209  enum {
1210  X, Y, Z, DIM_XYZ,
1211  };
1212  REAL8 totalmass, nu, chi1, beta1, kappa1, totalJ, S1, omega, eta, romega, v, LN, theta, polarizationangle;
1213  REAL8 phin0;
1214  totalmass = m1 + m2;
1215  if (m1 > m2) {
1216  nu = m2 / m1;
1217  } else {
1218  nu = m1 / m2;
1219  }
1220  if (LAL_SDW_MAX_PN_PARAM < 100. * nu * nu) {
1221  XLALPrintError(
1222  "XLAL Error: Spin-dominated waveforms error: Please make sure that the total mass is higher than 45 solar mass, and mass ratio is lower than 0.03125. Also above 130 solar mass be aware that high starting frequency may result in termination right after start, due to high value of the pn parameter. \n");
1223  XLAL_ERROR(XLAL_EDOM);
1224  } //too high mass ratio for the waveform, abort
1225  omega = fStart * LAL_PI;
1226  eta = nu / (1. + nu) / (1. + nu);
1227  chi1 = sqrt(s1x * s1x + s1y * s1y + s1z * s1z);
1228  if (chi1 < 0.5) {
1229  XLALPrintError(
1230  "XLAL Error: Spin-dominated waveforms error: Please make sure that the dimensionless spin parameter is higher than 0.5 \n");
1231  XLAL_ERROR(XLAL_EDOM);
1232  }
1233  REAL8 Lnh[DIM_XYZ] = { lnhatx, lnhaty, lnhatz };
1234  REAL8 LNHdotS1 = (Lnh[X] * s1x + Lnh[Y] * s1y + Lnh[Z] * s1z) / chi1;
1235  if (LNHdotS1 - 1.0 > 0 && LNHdotS1 - 1.0 < 1.0e-12)
1236  {kappa1 = acos(1.0);
1237  } else {
1238  kappa1 = acos((Lnh[X] * s1x + Lnh[Y] * s1y + Lnh[Z] * s1z) / chi1);
1239  }
1240 
1241  // Calculate the orbital angular momentum, up to 1.5 PN, with SO corrections
1242  v = cbrt(G_CP2 * totalmass * omega / LAL_C_SI);
1243  romega = G_CP2 * totalmass / v / v * (1.);
1244  LN = eta * totalmass * romega * romega * omega;
1245  // Calculate Spin magnitude, and the total angular momentum J
1246  S1 = chi1 * LAL_G_SI / LAL_C_SI * totalmass * totalmass * eta / nu;
1247 
1248  REAL8 JxN[DIM_XYZ], JxL[DIM_XYZ], Nvec[DIM_XYZ];
1249  REAL8 J[DIM_XYZ] = { LN * Lnh[X] + S1 * s1x / chi1, LN * Lnh[Y] + S1 * s1y / chi1, LN * Lnh[Z] + S1 * s1z / chi1 };
1250  totalJ = sqrt(J[X] * J[X] + J[Y] * J[Y] + J[Z] * J[Z]);
1251  // calculate the remaining angles
1252  if (kappa1 < 1e-7) {
1253  kappa1 = 0.;
1254  phin0 = LAL_PI_2;
1255  beta1 = 0.;
1256  polarizationangle=LAL_PI;
1257  theta = incl;
1258  } else if (kappa1 - LAL_PI > 0 && kappa1 - LAL_PI < 1.0e-12) {
1259  kappa1 = LAL_PI;
1260  phin0 = LAL_PI_2;
1261  beta1 = LAL_PI;
1262  polarizationangle=LAL_PI;
1263  theta = incl;
1264  } else {
1265  REAL8 JdotS = (J[X] * s1x + J[Y] * s1y + J[Z] * s1z) / totalJ / chi1 ;
1266  if (JdotS - 1.0 > 0 && JdotS - 1.0 < 1.0e-12){
1267  beta1 = acos(1.0);
1268  } else{
1269  beta1 = acos((J[X] * s1x + J[Y] * s1y + J[Z] * s1z) / totalJ / chi1);
1270  }
1271  // Line of sight vectorProd
1272  Nvec[X]=0.;
1273  Nvec[Y]=cos(incl);
1274  Nvec[Z]=sin(incl);
1275  theta = acos((J[X] * Nvec[X] + J[Y] * Nvec[Y] + J[Z] * Nvec[Z]) / totalJ );
1276 
1277 
1278 
1279  vectorProd(J, Lnh, totalJ, JxL);
1280  vectorProd(J, Nvec, totalJ, JxN);
1281  REAL8 JxNxN[DIM_XYZ], NxLxN[DIM_XYZ], LxN[DIM_XYZ], NLNxJNN[DIM_XYZ];
1282 
1283  REAL8 JxNxNamp, NxLxNamp, polarizationanglesign, JxLamp;
1284  vectorProd(Lnh, Nvec, 1.0, LxN);
1285  vectorProd(JxN, Nvec, 1.0, JxNxN);
1286  vectorProd(Nvec, LxN, 1.0, NxLxN);
1287 
1288  JxLamp = sqrt(JxL[X]*JxL[X]+JxL[Y]*JxL[Y]+JxL[Z]*JxL[Z]);
1289  JxNxNamp=sqrt(JxNxN[X]*JxNxN[X]+JxNxN[Y]*JxNxN[Y]+JxNxN[Z]*JxNxN[Z]);
1290  NxLxNamp=sqrt(NxLxN[X]*NxLxN[X]+NxLxN[Y]*NxLxN[Y]+NxLxN[Z]*NxLxN[Z]);
1291 
1292  REAL8 JxNxNdotNxLxN = (JxNxN[X]*NxLxN[X]+JxNxN[Y]*NxLxN[Y]+JxNxN[Z]*NxLxN[Z])/JxNxNamp/NxLxNamp ;
1293 
1294  if (JxNxNdotNxLxN - 1.0 > 0 && JxNxNdotNxLxN - 1.0 < 1.0e-12){
1295  polarizationangle = acos(1.0);
1296  } else if (JxNxNdotNxLxN + 1.0 < 0 && fabs(JxNxNdotNxLxN + 1.0) < 1.0e-12){
1297  polarizationangle = acos(-1.0);
1298  }
1299  else {
1300  polarizationangle = acos((JxNxN[X]*NxLxN[X]+JxNxN[Y]*NxLxN[Y]+JxNxN[Z]*NxLxN[Z])/JxNxNamp/NxLxNamp);
1301  }
1302 
1303  vectorProd(NxLxN, JxNxN, JxNxNamp*NxLxNamp, NLNxJNN);
1304 
1305  if (beta1 < 1e-7){
1306  polarizationangle =LAL_PI;
1307  } else{
1308  polarizationanglesign = NLNxJNN[X]*Nvec[X]+NLNxJNN[Y]*Nvec[Y]+NLNxJNN[Z]*Nvec[Z];
1309 
1310  if (polarizationanglesign < 0.) {
1311  polarizationangle *= -1.0;
1312  }
1313  }
1314 
1315  REAL8 JxNxJ[DIM_XYZ];
1316  REAL8 alpha0;
1317  vectorProd(JxN, J, totalJ, JxNxJ);
1318 
1319  REAL8 JxNxJamp;
1320  JxNxJamp = sqrt(JxNxJ[X]*JxNxJ[X]+JxNxJ[Y]*JxNxJ[Y]+JxNxJ[Z]*JxNxJ[Z]);
1321 
1322  alpha0 = acos((JxNxJ[X] * JxL[X] + JxNxJ[Y] * JxL[Y] + JxNxJ[Z] * JxL[Z])/JxNxJamp/JxLamp);
1323 
1324  phin0 = alpha0;
1325 
1326  REAL8 JNJxJL[DIM_XYZ];
1327  vectorProd(JxNxJ, JxL, JxNxJamp*JxLamp, JNJxJL);
1328  REAL8 JNJxJLsign = JNJxJL[X]*J[X]+JNJxJL[Y]*J[Y]+JNJxJL[Z]*J[Z];
1329 
1330  if (JNJxJLsign > 0.) {
1331  phin0 *= -1.0;
1332  }
1333  phiRef += LAL_PI_2;
1334  }
1335  // calling the SDW driver with the prefered variables
1336  int n = XLALSimInspiralSpinDominatedWaveformDriver(hplus, hcross, totalmass, nu, chi1, D, kappa1, beta1, theta,
1337  fStart, fRef, phaseO, amplitudeO, deltaT, phiRef, phin0, polarizationangle);
1338  return n;
1339 }
1340 
1341 /**
1342  * Function calculating the Spin-Dominated waveforms
1343  * This waveform is an inspiral only, 1 spin, precessing waveform.
1344  * For the formulae see the appendix of Arxiv:1209.1722
1345  */
1346 int XLALSimInspiralSpinDominatedWaveformDriver(REAL8TimeSeries **hplus, /**< +-polarization waveform */
1347 REAL8TimeSeries **hcross, /**< x-polarization waveform */
1348 REAL8 totalmass, /**< total mass of the binary */
1349 REAL8 nu, /**< mass ratio */
1350 REAL8 chi1, /**< dimensionless spin paramter */
1351 REAL8 D, /**< Distance to the source */
1352 REAL8 kappa1, /**< Angle span by S_1 and L */
1353 REAL8 beta1, /**< Angle span by J and S_1 */
1354 REAL8 theta, /**< Angle span by the line of sight and J */
1355 REAL8 fStart, /**< Starting gravitational wave frequency*/
1356 REAL8 fRef, /**< Ending gravitational wave frequency*/
1357 int phaseO, /**< twice PN phase order */
1358 int amplitudeO, /**< twice PN amplitude order */
1359 REAL8 deltaT, /**< Sampling time interval */
1360 REAL8 phiRef, /**< Reference phase at the Reference Frequency */
1361 REAL8 phin0, /**< Starting value of the phi_n parameter */
1362 REAL8 polarizationangle /**< Angle to rotate the radiaton frame to the default LALSimulation radiation frame */
1363 ) {
1364  int idx;
1365  int n;
1366  unsigned int i;
1367  REAL8 phiShift;
1368  LIGOTimeGPS tStart = LIGOTIMEGPSZERO;
1369  /* check inputs for sanity */
1370  if (*hplus) {
1371  XLAL_ERROR(XLAL_EFAULT);
1372  }
1373  if (*hcross) {
1374  XLAL_ERROR(XLAL_EFAULT);
1375  }
1376  if (deltaT <= 0) {
1377  XLAL_ERROR(XLAL_EDOM);
1378  }
1379  if (totalmass < 0) {
1380  XLAL_ERROR(XLAL_EDOM);
1381  }
1382  if (fStart <= 0) {
1383  XLAL_ERROR(XLAL_EDOM);
1384  }
1385  /* set up the integrator*/
1386  LALAdaptiveRungeKuttaIntegrator *integrator = XLALAdaptiveRungeKutta4Init(LAL_SDW_NUM_VARIABLES,
1387  XLALSpinDominatedWaveformDerivatives, XLALSpinDominatedWaveformStoppingTest, LAL_SDW_ABSOLUTE_TOLERANCE,
1388  LAL_SDW_RELATIVE_TOLERANCE);
1389  if (!integrator) {
1390  XLALPrintError("XLAL Error - %s: Cannot allocate integrator\n", __func__);
1391  XLAL_ERROR(XLAL_EFUNC);
1392  }
1393  /* stop the integration only when the test is true */
1394  integrator->stopontestonly = 1;
1395  LALSDWaveformParams params;
1396  params.totalmass = totalmass;
1397  params.nu = nu;
1398  params.chi1 = chi1;
1399  params.dist = D;
1400  params.kappa1 = kappa1;
1401  params.beta1 = beta1;
1402  params.theta = theta;
1403  params.eps = 0.;
1404  params.xi = 0.;
1405  params.pnamp = amplitudeO;
1406  params.pnphase = phaseO;
1407  params.prevdomega = 0.;
1408  params.polarizationangle = polarizationangle;
1409  n = XLALSpinDominatedWaveformConstantCoefficients(&params);
1410  if (n < 0) {
1411  XLAL_ERROR(XLAL_EFUNC);
1412  }
1413  REAL8 yin[LAL_SDW_NUM_VARIABLES];
1414  yin[PHI] = phin0;
1415  yin[OMEGA] = fStart * LAL_PI;
1416  yin[PSI] = 0.;
1417  REAL8Array *yout;
1418  // estimating the length of the waveform
1419  REAL8 length = 5. / 256. * pow(fStart * LAL_PI, -8. / 3.) * (1. + params.nu) * (1. + params.nu) / params.nu
1420  * pow(G_CP2 * params.totalmass / LAL_C_SI, -5. / 3.);
1421  INT4 intLen = XLALAdaptiveRungeKutta4Hermite(integrator, (void *) &params, yin, 0.0, length, deltaT, &yout);
1422  UNUSED INT4 intReturn = integrator->returncode;
1423  XLALAdaptiveRungeKuttaFree(integrator);
1424  REAL8TimeSeries *phin = XLALCreateREAL8TimeSeries("PHI_N", &tStart, 0., deltaT, &lalDimensionlessUnit, intLen);
1425  REAL8TimeSeries *omega = XLALCreateREAL8TimeSeries("OMEGA", &tStart, 0., deltaT, &lalDimensionlessUnit, intLen);
1426  REAL8TimeSeries *psi = XLALCreateREAL8TimeSeries("ORBITAL_PHASE", &tStart, 0., deltaT, &lalDimensionlessUnit,
1427  intLen);
1428  for (idx = 0; idx < intLen; idx++) {
1429  phin->data->data[idx] = yout->data[intLen + idx];
1430  omega->data->data[idx] = yout->data[2 * intLen + idx];
1431  psi->data->data[idx] = yout->data[3 * intLen + idx];
1432  }
1433 
1434  if (fRef == 0) {
1435  phiShift = phiRef - psi->data->data[0];
1436  for (i = 0; i < psi->data->length; i++) {
1437  psi->data->data[i] += phiShift;
1438  }
1439  } else if (fRef == fStart) {
1440  phiShift = phiRef - psi->data->data[0];
1441  for (i = 0; i < psi->data->length; i++) {
1442  psi->data->data[i] += phiShift;
1443  }
1444  } else {
1445  XLALPrintError(
1446  "XLAL Error: Spin-dominated waveforms error: Please set the reference frequency as the starting frequency, Setting 0 will default to the starting frequency. \n");
1447  XLAL_ERROR(XLAL_EDOM);
1448  }
1449  if ((*hplus) && (*hcross)) {
1450  if ((*hplus)->data->length != (*hcross)->data->length) {
1451  XLALPrintError("*** h+ and hx differ in length\n");
1452  XLAL_ERROR(XLAL_EFAILED);
1453  } else {
1454  if ((int) (*hplus)->data->length < intLen) {
1455  XLALPrintError("*** ERROR: h+ and hx too short\n");
1456  XLAL_ERROR(XLAL_EFAILED);
1457  } else {
1458  XLALGPSAdd(&((*hplus)->epoch), -intLen * deltaT);
1459  XLALGPSAdd(&((*hcross)->epoch), -intLen * deltaT);
1460  }
1461  }
1462  } else {
1463  XLALGPSAdd(&tStart, -intLen * deltaT);
1464  *hplus = XLALCreateREAL8TimeSeries("H+", &tStart, 0.0, deltaT, &lalStrainUnit, intLen);
1465  *hcross = XLALCreateREAL8TimeSeries("Hx", &tStart, 0.0, deltaT, &lalStrainUnit, intLen);
1466  if (*hplus == NULL || *hcross == NULL) {
1467  XLAL_ERROR(XLAL_ENOMEM);
1468  }
1469  }
1470 
1471 
1472  REAL8 expr[3];
1473  for (idx = 0; idx < intLen; idx++) {
1474  expr[PHI] = phin->data->data[idx];
1475  expr[OMEGA] = omega->data->data[idx];
1476  expr[PSI] = psi->data->data[idx];
1477 
1478  n = XLALSpinDominatedWaveformBuild(&params, expr, hplus, hcross, idx);
1479  if (n < 0) {
1480  XLAL_ERROR(XLAL_EFUNC);
1481  }
1482 
1483  }
1484 
1485  XLALDestroyREAL8Array(yout);
1486  XLALDestroyREAL8TimeSeries(phin);
1487  XLALDestroyREAL8TimeSeries(omega);
1488  XLALDestroyREAL8TimeSeries(psi);
1489  return intLen;
1490 }
1491 
1492 /**
1493  * Function calculating the derivatives of the three time dependent variables of the Spin-Dominated waveforms (SDW)
1494  * The first paramter is phi_n, Eq 27 of Arxiv:1005.5330, taken for 1 spin case, and integrated over an orbital period.
1495  * The second parameter is omega, the derivative is taken from Arxiv: astro-ph/0504538, up to 2 PN orders with 1 spin. (In order to stay consistent with SDW)
1496  * The thirs parameter is the phase.
1497  */
1498 static INT4 XLALSpinDominatedWaveformDerivatives(UNUSED REAL8 t, const REAL8 values[], REAL8 dvalues[], void *mparams) {
1499  LALSDWaveformParams *params = (LALSDWaveformParams *) mparams;
1500  // parameters required for the time derivatives
1501  REAL8 vP[OMEGA_POWER_DIM];
1502  vP[0] = 1.;
1503  vP[1] = cbrt(LAL_G_SI * params->totalmass * values[OMEGA] / LAL_C_SI / LAL_C_SI / LAL_C_SI);
1504  for (int i = 2; i < OMEGA_POWER_DIM; ++i) {
1505  vP[i] = vP[1] * vP[i - 1];
1506  }
1507 
1508  params->eps = vP[2];
1509  REAL8 epsP3 = params->eps * params->eps * params->eps;
1510  params->xi = params->nu / sqrt(params->eps);
1511  REAL8 eta = params->nu / (1. + params->nu) / (1. + params->nu);
1512  REAL8 phasecoeff = 96. / 5. * eta * vP[5] * values[OMEGA] * values[OMEGA];
1513  REAL8 sinKappa1 = sin(params->kappa1);
1514  REAL8 cosKappa1 = cos(params->kappa1);
1515  // Calculating the derivatives
1516 
1517  dvalues[PHI] = 0;
1518  if (params->kappa1 != 0 ) {
1519  switch (params->pnphase) {
1520  case (PNDEF): //default value -1 contains all corrections
1521  case (PN20):
1522  dvalues[PHI] += +3. / 2. / LAL_G_SI / params->totalmass / sinKappa1 * params->chi1 * params->chi1 * LAL_C_SI
1523  * LAL_C_SI * LAL_C_SI * epsP3 * sqrt(params->eps)
1524  * (1. - 2. * params->xi * sqrt(params->eps)) * (sinKappa1 * cosKappa1
1525  + params->beta1 * cosKappa1 * cosKappa1);
1526 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
1527  __attribute__ ((fallthrough)); /* no break */
1528 #endif
1529  case (PN15):
1530  dvalues[PHI] += params->chi1 * LAL_C_SI * LAL_C_SI * LAL_C_SI * epsP3 / 2. / LAL_G_SI / params->totalmass
1531  * ( 5. * sqrt(params->eps) * params->xi - 4.)*(1. + cosKappa1 * params->beta1 / sinKappa1 );
1532 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
1533  __attribute__ ((fallthrough)); /* no break */
1534 #endif
1535  case (PN10):
1536  /* no break */
1537  case (PN05):
1538  /* no break */
1539  case (PN00):
1540  break;
1541  }
1542  }
1543 
1544  dvalues[OMEGA] = 0;
1545  switch (params->pnphase) {
1546  case (PNDEF): //default value -1 contains all corrections
1547  case (PN20):
1548  dvalues[OMEGA] += phasecoeff * vP[PN20] * (34103. / 18144. + 13661. / 2016. * eta + 59. / 18. * eta * eta);
1549  // QM and Self-Spin components taken together
1550  dvalues[OMEGA] += phasecoeff * vP[PN20]
1551  * (5. / 2. * params->chi1 * params->chi1 * eta / params->nu * (3. * cosKappa1 * cosKappa1 - 1.)
1552  + 1. / 96. * params->chi1 * params->chi1 * eta / params->nu * (6. + sinKappa1 * sinKappa1));
1553 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
1554  __attribute__ ((fallthrough)); /* no break */
1555 #endif
1556  case (PN15):
1557  dvalues[OMEGA] += phasecoeff * vP[PN15] * (4. * LAL_PI);
1558  // SO component
1559  dvalues[OMEGA] += phasecoeff * vP[PN15]
1560  * (-1 / 12. * cosKappa1 * params->chi1 * (113. * eta / params->nu + 75. * eta));
1561 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
1562  __attribute__ ((fallthrough)); /* no break */
1563 #endif
1564  case (PN10):
1565  dvalues[OMEGA] += -phasecoeff * vP[PN10] * (743. / 336. + 11. / 4. * eta);
1566 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
1567  __attribute__ ((fallthrough)); /* no break */
1568 #endif
1569  case (PN05):
1570  /* no break */
1571  case (PN00):
1572  dvalues[OMEGA] += phasecoeff * 1.;
1573  break;
1574  }
1575  dvalues[PSI] = values[OMEGA]-dvalues[PHI]*cos(params->kappa1-params->beta1);
1576 
1577  return GSL_SUCCESS;
1578 } /* end of XLALSpinDominatedWaveformDerivatives */
1579 
1580 /**
1581  * Stopping test for the Spin-Dominated waveforms. Using MECO, or the desired ending frequency. The maximum value of the PN parameter is set to 0.8.
1582  */
1583 static INT4 XLALSpinDominatedWaveformStoppingTest(UNUSED REAL8 t, const REAL8 values[], REAL8 dvalues[],
1584 UNUSED void *mparams) {
1585  LALSDWaveformParams *params = (LALSDWaveformParams *) mparams;
1586  REAL8 vP[OMEGA_POWER_DIM];
1587  vP[0] = 1.;
1588  vP[1] = cbrt(G_CP2 * params->totalmass * values[OMEGA] / LAL_C_SI);
1589  for (int i = 2; i < OMEGA_POWER_DIM; ++i) {
1590  vP[i] = vP[1] * vP[i - 1];
1591  }
1592  REAL8 cosKappa1 = cos(params->kappa1);
1593  REAL8 eta = params->nu / (1. + params->nu) / (1. + params->nu);
1594  REAL8 omega = values[OMEGA];
1595  REAL8 domega = dvalues[OMEGA];
1596  REAL8 d2omega = dvalues[OMEGA] - params->prevdomega;
1597  params->prevdomega = dvalues[OMEGA];
1598  REAL8 mecotest;
1599  mecotest = 0.;
1600 
1601  switch (params->pnphase) {
1602  case (PNDEF): //default value -1 contains all corrections
1603  case (PN20):
1604  mecotest += +6. * vP[PN20]
1605  * (1. / 8. * (-27. + 19. * eta - eta * eta / 3.)
1606  - (3. * cosKappa1 * cosKappa1 - 1.) / 2. * params->chi1 * params->chi1 * eta / params->nu);
1607 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
1608  __attribute__ ((fallthrough)); /* no break */
1609 #endif
1610  case (PN15):
1611  mecotest += +5. * vP[PN15] * (8. / 3. * eta / params->nu + 2. * eta) * cosKappa1 * params->chi1;
1612 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
1613  __attribute__ ((fallthrough)); /* no break */
1614 #endif
1615  case (PN10):
1616  mecotest += -4. * vP[PN10] * (3. + eta / 3.) / 4.;
1617 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
1618  __attribute__ ((fallthrough)); /* no break */
1619 #endif
1620  case (PN05):
1621  /* no break */
1622  case (PN00):
1623  mecotest += +2.;
1624  break;
1625  }
1626 
1627  /* Check value of SEOBNRv2 stopping frequency
1628 
1629  REAL8 mcheck1 = params->totalmass / (1. + params->nu);
1630  REAL8 mcheck2 = params->totalmass - mcheck1;
1631  REAL8 spincheckz1 = params->chi1 * cos(params->kappa1);
1632  REAL8 spincheckz2 = 0.0;
1633  REAL8 seobnr_stop_freq = XLALSimIMRSpinAlignedEOBPeakFrequency(mcheck1, mcheck2, spincheckz1, spincheckz2, 2);
1634 */
1635 
1636  if (mecotest < 0) {
1637  XLALPrintWarning("** LALSimInspiralSDW WARNING **: MECO reached\n");
1638  return -1;
1639 
1640  } /*else if (omega > seobnr_stop_freq*LAL_PI) {
1641  XLALPrintWarning("** LALSimInspiralSDW WARNING **: SEOBNR stopping frequency reached\n");
1642  return -1;
1643 
1644  } */else if (isnan(omega)) {
1645  XLALPrintWarning("** LALSimInspiralSDW WARNING **: omega is NAN\n");
1646  return -1;
1647  } else if (vP[1] >= 1.) {
1648  XLALPrintWarning("** LALSimInspiralSDW WARNING **: PN parameter is too large\n");
1649  return -1;
1650  } else if (domega < 0.0) {
1651  XLALPrintWarning("** LALSimInspiralSDW WARNING **: domega < 0\n");
1652  return -1;
1653  } else if (d2omega <= 0.) {
1654  XLALPrintWarning("** LALSimInspiralSDW WARNING **: d2omega < 0\n");
1655  return -1;
1656  }
1657  return GSL_SUCCESS;
1658 } /* End of XLALSpinDominatedWaveformStoppingTest */
LALMalloc
#define LALMalloc(n)
LALFree
#define LALFree(p)
LAL_SDW_ABSOLUTE_TOLERANCE
#define LAL_SDW_ABSOLUTE_TOLERANCE
Definition: LALSimInspiralSpinDominatedWaveform.c:34
LAL_SDW_RELATIVE_TOLERANCE
#define LAL_SDW_RELATIVE_TOLERANCE
Definition: LALSimInspiralSpinDominatedWaveform.c:35
XLALSpinDominatedWaveformConstantCoefficients
int XLALSpinDominatedWaveformConstantCoefficients(LALSDWaveformParams *params)
Function for calculating the constant coefficients of Spin-Dominated waveforms See tables 1 to 5 in t...
Definition: LALSimInspiralSpinDominatedWaveform.c:132
XLALFreeDmatrix
static void XLALFreeDmatrix(REAL8 *m)
Function for freeing memory for a matrix.
Definition: LALSimInspiralSpinDominatedWaveform.c:124
LAL_SDW_NUM_VARIABLES
#define LAL_SDW_NUM_VARIABLES
Definition: LALSimInspiralSpinDominatedWaveform.c:36
G_CP2
static const REAL8 G_CP2
Definition: LALSimInspiralSpinDominatedWaveform.c:47
XLALSpinDominatedWaveformBuild
int XLALSpinDominatedWaveformBuild(LALSDWaveformParams *params, REAL8 expr[], REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, int idx)
Function building the wavefrom from the calculated parameters at a given time For the formulae see th...
Definition: LALSimInspiralSpinDominatedWaveform.c:612
XLALSpinDominatedWaveformDerivatives
static INT4 XLALSpinDominatedWaveformDerivatives(UNUSED REAL8 t, const REAL8 values[], REAL8 dvalues[], void *mparams)
Function calculating the derivatives of the three time dependent variables of the Spin-Dominated wave...
Definition: LALSimInspiralSpinDominatedWaveform.c:1498
XLALDmatrix
static REAL8 * XLALDmatrix(INT8 nrh, INT8 nch)
Function for allocating memory for a matrix.
Definition: LALSimInspiralSpinDominatedWaveform.c:111
XLALSpinDominatedWaveformStoppingTest
static INT4 XLALSpinDominatedWaveformStoppingTest(UNUSED REAL8 t, const REAL8 values[], REAL8 dvalues[], UNUSED void *mparams)
Stopping test for the Spin-Dominated waveforms.
Definition: LALSimInspiralSpinDominatedWaveform.c:1583
XLALSimInspiralSpinDominatedWaveformDriver
int XLALSimInspiralSpinDominatedWaveformDriver(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, REAL8 totalmass, REAL8 nu, REAL8 chi1, REAL8 D, REAL8 kappa1, REAL8 beta1, REAL8 theta, REAL8 fStart, REAL8 fRef, int phaseO, int amplitudeO, REAL8 deltaT, REAL8 phiRef, REAL8 phin0, REAL8 polarizationangle)
Function calculating the Spin-Dominated waveforms This waveform is an inspiral only,...
Definition: LALSimInspiralSpinDominatedWaveform.c:1346
XLALSimInspiralSpinDominatedWaveformInterfaceTD
int XLALSimInspiralSpinDominatedWaveformInterfaceTD(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 fStart, REAL8 fRef, REAL8 D, REAL8 s1x, REAL8 s1y, REAL8 s1z, REAL8 lnhatx, REAL8 lnhaty, REAL8 lnhatz, REAL8 incl, int phaseO, int amplitudeO, REAL8 phiRef)
Interface routine, calculating the prefered variables for the Spin-dominated waveforms.
Definition: LALSimInspiralSpinDominatedWaveform.c:1190
LAL_SDW_MAX_PN_PARAM
#define LAL_SDW_MAX_PN_PARAM
Definition: LALSimInspiralSpinDominatedWaveform.c:37
vectorProd
#define vectorProd(lhs, rhs, denominator, result)
Definition: LALSimInspiralSpinDominatedWaveform.c:69
i
double i
Definition: bh_ringdown.c:118
e
double e
Definition: bh_ringdown.c:117
theta
double theta
Definition: bh_sphwf.c:118
__attribute__
#define __attribute__(x)
XLALDestroyREAL8Array
void XLALDestroyREAL8Array(REAL8Array *)
XLALAdaptiveRungeKuttaFree
void XLALAdaptiveRungeKuttaFree(LALAdaptiveRungeKuttaIntegrator *integrator)
XLALAdaptiveRungeKutta4Hermite
int XLALAdaptiveRungeKutta4Hermite(LALAdaptiveRungeKuttaIntegrator *integrator, void *params, REAL8 *yinit, REAL8 tinit, REAL8 tend_in, REAL8 deltat, REAL8Array **yout)
XLALAdaptiveRungeKutta4Init
LALAdaptiveRungeKuttaIntegrator * XLALAdaptiveRungeKutta4Init(int dim, int(*dydt)(double t, const double y[], double dydt[], void *params), int(*stop)(double t, const double y[], double dydt[], void *params), double eps_abs, double eps_rel)
LAL_PI_2
#define LAL_PI_2
LAL_C_SI
#define LAL_C_SI
LAL_PI
#define LAL_PI
LAL_G_SI
#define LAL_G_SI
LIGOTIMEGPSZERO
#define LIGOTIMEGPSZERO
REAL8
double REAL8
INT8
int64_t INT8
INT4
int32_t INT4
m
static const INT4 m
XLALCreateREAL8TimeSeries
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
XLALDestroyREAL8TimeSeries
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
lalStrainUnit
const LALUnit lalStrainUnit
lalDimensionlessUnit
const LALUnit lalDimensionlessUnit
XLAL_ERROR
#define XLAL_ERROR(...)
XLALPrintError
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLALPrintWarning
int int XLALPrintWarning(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_ENOMEM
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFAULT
XLAL_EFUNC
XLAL_EFUNC
XLAL_EDOM
XLAL_EDOM
XLAL_EFAILED
XLAL_EFAILED
XLALGPSAdd
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
LALSDWaveformParams
Structure containing the prefered variabloes for Spin-Dominated waveforms.
Definition: LALSimInspiralSpinDominatedWaveform.c:77
REAL8Array::data
REAL8 * data
REAL8TimeSeries::data
REAL8Sequence * data
REAL8Vector::data
REAL8 * data
params
Definition: burst.c:245
deltaT
double deltaT
Definition: unicorn.c:24