Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-5e288d3
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
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
47static const REAL8 G_CP2 = LAL_G_SI / LAL_C_SI / LAL_C_SI;
48
49enum {
50 PNDEF = -1 ,PN00 = 0, PN05 = 1, PN10 = 2, PN15 =3, PN20 =4, PN25 =5, PN3 =6, PN_ORDER =7,
51};
52
53typedef enum {
54 PHI, OMEGA, PSI, PHASE = PSI,
55} TIME_DEPENDANT;
56
57typedef enum {
58 PLUS_ = 0, MINUS = 1, CROSS_ = 1, PLUS_MINUS_DIM = 2, PLUS_CROSS_DIM = 2,
59} COMPONENTS;
60
61typedef enum {
62 PN00DIM = 2, PN05DIM = 11, PN10DIM = 15, PN15DIM = 17,
63} COEFFICIENT_DIMENSIONS;
64
65typedef 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 */
77typedef 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
98static INT4 XLALSpinDominatedWaveformStoppingTest(UNUSED REAL8 t, const REAL8 values[], REAL8 dvalues[],
99UNUSED void *mparams);
100
101static INT4 XLALSpinDominatedWaveformDerivatives(UNUSED REAL8 t, const REAL8 values[], REAL8 dvalues[], void *mparams);
102
103int XLALSpinDominatedWaveformBuild(LALSDWaveformParams *params, REAL8 expr[], REAL8TimeSeries **hplus,
104 REAL8TimeSeries **hcross, int idx);
105
106int XLALSpinDominatedWaveformConstantCoefficients(LALSDWaveformParams * params);
107
108/**
109 * Function for allocating memory for a matrix
110 */
111static 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 */
124static 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 */
132int 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 */
612int XLALSpinDominatedWaveformBuild(LALSDWaveformParams *params, /**< The SDW parameters */
613REAL8 expr[], /**< The 3 time dependent variables of the waveform at the time indexed by idx */
614REAL8TimeSeries **hplus, /**< +-polarization waveform */
615REAL8TimeSeries **hcross, /**< x-polarization waveform */
616int 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 */
1190int XLALSimInspiralSpinDominatedWaveformInterfaceTD(REAL8TimeSeries **hplus, /**< +-polarization waveform */
1191REAL8TimeSeries **hcross, /**< x-polarization waveform */
1192REAL8 deltaT, /**< sampling interval (s) */
1193REAL8 m1, /**< mass of companion 1 (kg) */
1194REAL8 m2, /**< mass of companion 2 (kg) */
1195REAL8 fStart, /**< start GW frequency (Hz) */
1196REAL8 fRef, /**< end GW frequency (Hz) */
1197REAL8 D, /**< distance of source (m) */
1198REAL8 s1x, /**< initial value of S1x */
1199REAL8 s1y, /**< initial value of S1y */
1200REAL8 s1z, /**< initial value of S1z */
1201REAL8 lnhatx, /**< initial value of LNhatx */
1202REAL8 lnhaty, /**< initial value of LNhaty */
1203REAL8 lnhatz, /**< initial value of LNhatz */
1204REAL8 incl, /**< inclination, angle between L and line of sight N */
1205int phaseO, /**< twice PN phase order */
1206int amplitudeO, /**< twice PN amplitude order */
1207REAL8 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 */
1346int XLALSimInspiralSpinDominatedWaveformDriver(REAL8TimeSeries **hplus, /**< +-polarization waveform */
1347REAL8TimeSeries **hcross, /**< x-polarization waveform */
1348REAL8 totalmass, /**< total mass of the binary */
1349REAL8 nu, /**< mass ratio */
1350REAL8 chi1, /**< dimensionless spin paramter */
1351REAL8 D, /**< Distance to the source */
1352REAL8 kappa1, /**< Angle span by S_1 and L */
1353REAL8 beta1, /**< Angle span by J and S_1 */
1354REAL8 theta, /**< Angle span by the line of sight and J */
1355REAL8 fStart, /**< Starting gravitational wave frequency*/
1356REAL8 fRef, /**< Ending gravitational wave frequency*/
1357int phaseO, /**< twice PN phase order */
1358int amplitudeO, /**< twice PN amplitude order */
1359REAL8 deltaT, /**< Sampling time interval */
1360REAL8 phiRef, /**< Reference phase at the Reference Frequency */
1361REAL8 phin0, /**< Starting value of the phi_n parameter */
1362REAL8 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 */
1498static 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 */
1583static INT4 XLALSpinDominatedWaveformStoppingTest(UNUSED REAL8 t, const REAL8 values[], REAL8 dvalues[],
1584UNUSED 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
XLALDmatrix
static REAL8 * XLALDmatrix(INT8 nrh, INT8 nch)
Function for allocating memory for a matrix.
Definition: LALSimInspiralSpinDominatedWaveform.c:111
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
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)
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)
XLALAdaptiveRungeKutta4Hermite
int XLALAdaptiveRungeKutta4Hermite(LALAdaptiveRungeKuttaIntegrator *integrator, void *params, REAL8 *yinit, REAL8 tinit, REAL8 tend_in, REAL8 deltat, REAL8Array **yout)
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
XLALDestroyREAL8TimeSeries
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
XLALCreateREAL8TimeSeries
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
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