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
LALSimInspiralSpinTaylorT5duplicate.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2012 P. Ajith
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 <lal/Units.h>
22#include <lal/LALConstants.h>
23#include <lal/LALAdaptiveRungeKuttaIntegrator.h>
24#include <lal/TimeSeries.h>
25#include <lal/LALSimInspiral.h>
26#include <lal/VectorOps.h>
27#include "check_series_macros.h"
28#define LAL_PISQR 9.869604401089358
29
30#define UNUSED(expr) do { (void)(expr); } while (0)
31
32/* use error codes above 1024 to avoid conflicts with GSL */
33#define LALSIMINSPIRAL_ST5_TEST_ENERGY 1025
34#define LALSIMINSPIRAL_ST5_TEST_VDOT 1026
35#define LALSIMINSPIRAL_ST5_TEST_COORDINATE 1027
36#define LALSIMINSPIRAL_ST5_TEST_VNAN 1028
37#define LALSIMINSPIRAL_ST5_TEST_FREQBOUND 1029
38#define LALSIMINSPIRAL_ST5_DERIVATIVE_OMEGANONPOS 1030
39
40/* (2x) Highest available PN order - UPDATE IF NEW ORDERS ADDED!!*/
41#define LAL_MAX_PN_ORDER 8
42/* Number of variables used for precessing waveforms */
43#define LAL_NUM_ST5_VARIABLES 10
44/* absolute and relative tolerance for adaptive Runge-Kutta ODE integrator */
45#define LAL_ST5_ABSOLUTE_TOLERANCE 1.e-12
46#define LAL_ST5_RELATIVE_TOLERANCE 1.e-12
47
48/* terminate the ODE integration at or below this value of the PN parameter */
49#define PN_V_MAX 0.6
50
51/**
52 * Structure containing the non-dynamical coefficients needed
53 * to evolve a spinning, precessing binary and produce a waveform.
54 * This struct is passed to the 2 static functions below
55 */
56
57typedef struct SpinTaylorT5StructParams {
58 double mass1;
59 double mass2;
60 double totalMass;
61 double eta;
62 double delta;
63 int phaseO;
64 double massRatio;
65 double m1Sqr;
66 double m2Sqr;
67 double mSqr;
68 double etaSqr;
69 double deltaEta;
70 double dEbF0NonSpin;
71 double dEbF1NonSpin;
72 double dEbF2NonSpin;
73 double dEbF3NonSpin;
74 double dEbF4NonSpin;
75 double dEbF5NonSpin;
76 double dEbF6NonSpin;
77 double dEbF6LogNonSpin;
78 double dEbF7NonSpin;
79 double phiOrb0NonSpin;
80 double phiOrb1NonSpin;
81 double phiOrb2NonSpin;
82 double phiOrb3NonSpin;
83 double phiOrb4NonSpin;
84 double phiOrb5LogNonSpin;
85 double phiOrb6NonSpin;
86 double phiOrb6LogNonSpin;
87 double phiOrb7NonSpin;
88 double v0;
89 double vMax;
90 double phiRef;
91} SpinTaylorT5Params;
92
93/* declaration of static functions */
94static int XLALSimInspiralSpinTaylorT5Derivatives(
95 double t,
96 const double y[],
97 double dy[],
98 void *mParams
99 );
100
101static REAL8 XLALdEnergyByFluxSpinPrec(
102 double v,
103 double *chi1,
104 double *chi2,
105 double *LNh,
106 SpinTaylorT5Params *params
107 );
108
109static int polarizationsInRadiationFrame(
110 REAL8TimeSeries **hplus, /**< +-polarization waveform [returned] */
111 REAL8TimeSeries **hcross, /**< x-polarization waveform [returned] */
112 REAL8TimeSeries *V, /**< post-Newtonian parameter */
113 REAL8TimeSeries *Phi, /**< orbital phase */
114 REAL8TimeSeries *LNhxVec, /**< unit orbital ang. mom. x comp. */
115 REAL8TimeSeries *LNhyVec, /**< unit orbital ang. mom. y comp. */
116 REAL8TimeSeries *LNhzVec, /**< unit orbital ang. mom. z comp. */
117 REAL8 m1, /**< mass1 (KG) */
118 REAL8 m2, /**< mass2 (KG) */
119 REAL8 r, /**< distance (meter) */
120 REAL8 Theta /**< angle between the initial total ang momentum and line of sight. */
121 );
122
123static int computeOrbitalPhase(
124 SpinTaylorT5Params *params,
125 REAL8TimeSeries *V,
126 REAL8TimeSeries *LNhxVec,
127 REAL8TimeSeries *LNhyVec,
128 REAL8TimeSeries *LNhzVec,
129 REAL8TimeSeries *S1xVec,
130 REAL8TimeSeries *S1yVec,
131 REAL8TimeSeries *S1zVec,
132 REAL8TimeSeries *S2xVec,
133 REAL8TimeSeries *S2yVec,
134 REAL8TimeSeries *S2zVec,
135 REAL8TimeSeries *orbPhase);
136
137static int spinTaylorT5Init(
138 REAL8 m1,
139 REAL8 m2,
140 REAL8 fStart,
141 REAL8 deltaT,
142 REAL8 phiRef,
143 UINT4 phaseO,
144 SpinTaylorT5Params *mParams
145 );
146
147static int XLALSimInspiralSpinTaylorT5StoppingTest(
148 double t,
149 const double y[],
150 double dy[],
151 void *mParams
152 );
153
154static double dotProduct(double a[3], double b[3]);
155static void crossProduct(double a[3], double b[3], double c[3]);
156static void rotateVector(double a[3], double b[3], double phi, double theta, double psi);
157
158/**
159 * @addtogroup LALSimInspiralSpinTaylor_c
160 * @{
161 */
162
163/**
164 * Generate time-domain generic spinning PN waveforms in the SpinTaylorT5 approximaton.
165 */
166int XLALSimInspiralSpinTaylorT5duplicate(
167 REAL8TimeSeries **hplus, /**< +-polarization waveform */
168 REAL8TimeSeries **hcross, /**< x-polarization waveform */
169 REAL8 phiRef, /**< orbital phase at reference pt. */
170 REAL8 deltaT, /**< sampling interval (s) */
171 REAL8 m1, /**< mass of companion 1 (kg) */
172 REAL8 m2, /**< mass of companion 2 (kg) */
173 REAL8 fStart, /**< start GW frequency (Hz) */
174 REAL8 r, /**< distance of source (m) */
175 REAL8 s1x, /**< initial value of S1x */
176 REAL8 s1y, /**< initial value of S1y */
177 REAL8 s1z, /**< initial value of S1z */
178 REAL8 s2x, /**< initial value of S2x */
179 REAL8 s2y, /**< initial value of S2y */
180 REAL8 s2z, /**< initial value of S2z */
181 REAL8 incAngle, /**< inclination angle with J_ini */
182 int phaseO, /**< twice PN phase order */
183 int amplitudeO /**< twice PN amplitude order */
184 ) {
185
186 XLAL_PRINT_DEPRECATION_WARNING("XLALSimInspiralSpinTaylorT5");
187 UINT4 i, intStatus, lenReturn;
188 LIGOTimeGPS tStart = LIGOTIMEGPSZERO;
189 REAL8TimeSeries *orbPhase=NULL, *V=NULL, *LNhxVec=NULL, *LNhyVec=NULL, *LNhzVec=NULL;
190 REAL8TimeSeries *S1xVec=NULL, *S1yVec=NULL, *S1zVec=NULL, *S2xVec=NULL, *S2yVec=NULL, *S2zVec=NULL;
191 LALAdaptiveRungeKuttaIntegrator *integrator = NULL; /* GSL integrator object */
192 SpinTaylorT5Params *mParams;
193 SpinTaylorT5Params SpinTaylorParameters;
194 REAL8 LNh[3], S1[3], S2[3], J[3], Jh[3], Lmag, Jmag, v;
195
196 mParams = &SpinTaylorParameters; /* structure containing various parameters and coefficients */
197 REAL8 yinit[LAL_NUM_ST5_VARIABLES]; /* initial values of parameters */
198 REAL8Array *yout=NULL; /* time series of variables returned from integrator */
199
200 /* initialize the mParams containing the the PN expansion parameters */
201 if (spinTaylorT5Init(m1, m2, fStart, deltaT, phiRef, phaseO, mParams) != XLAL_SUCCESS) {
202 XLALPrintError("XLAL Error - %s: Unable to initialize the mParams structure.\n", __func__);
203 return XLAL_FAILURE;
204 }
205
206 /* Estimate length of waveform using Newtonian chirp time formula. give 10% extra time */
207 UINT4 dataLength = pow(2, ceil(log2(1.1*5*mParams->totalMass/(256.*pow(mParams->v0,8.)*mParams->eta)/deltaT)));
208
209 /* Adjust tStart so last sample is at time=0 */
210 XLALGPSAdd(&tStart, -1.0*(dataLength-1)*deltaT);
211
212 /* binary parameters in a coordinate system whose z axis is along the initial
213 orbital angular momentum */
214 LNh[0] = 0.; /* unit vector along Newt. ang. momtm */
215 LNh[1] = 0.; /* unit vector along Newt. ang. momtm */
216 LNh[2] = 1.; /* unit vector along Newt. ang. momtm */
217 S1[0] = s1x*mParams->m1Sqr; /* spin vector of the first object */
218 S1[1] = s1y*mParams->m1Sqr; /* spin vector of the first object */
219 S1[2] = s1z*mParams->m1Sqr; /* spin vector of the first object */
220 S2[0] = s2x*mParams->m2Sqr; /* spin vector of the second object */
221 S2[1] = s2y*mParams->m2Sqr; /* spin vector of the second object */
222 S2[2] = s2z*mParams->m2Sqr; /* spin vector of the second object */
223
224 /* compute the initial total angular momentum */
225 Lmag = (mParams->eta*mParams->mSqr/mParams->v0)*(1. + (1.5+mParams->eta/6.)*mParams->v0*mParams->v0); // magnitde of orb ang momentum (1PN)
226 //Lmag = (mParams->eta*mParams->mSqr/mParams->v0); /* magnitude of orb ang momentum (Newtonian) */
227 for (i=0; i<3; i++) J[i] = Lmag*LNh[i] + S1[i] + S2[i]; // total angular momentum vec
228 Jmag = sqrt(J[0]*J[0] + J[1]*J[1] + J[2]*J[2]); // magnitude of tot ang momentum
229 for (i=0; i<3; i++) Jh[i] = J[i]/Jmag; // unit vector along J
230
231 /* transform to the coordinate system defined by the initial total angular momentum */
232 REAL8 JIniTheta = -acos(Jh[2]);
233 REAL8 JIniPhi = -atan2(Jh[1], Jh[0]);
234
235 /* if both objects are non-spinning, the coordinate system defined by the total angular
236 momentum (= orbital angular momentum) will cause the dynamical variables to hit
237 coordinate singularities. Thus, for the case of non-spinning binaries, we avoid this
238 situlation by choosing a coordinate system which is misaligned by 45 degs with the
239 total angular momentum direction. We also shift the inclination angle by this amount
240 so that, at the end, the waveforms are not affected by this choice of coordinate
241 system */
242 if (JIniTheta == 0) {
243 JIniTheta = LAL_PI/4.;
244 incAngle = incAngle+JIniTheta;
245 }
246
247 rotateVector(LNh, LNh, JIniPhi, JIniTheta, 0.);
248 rotateVector(S1, S1, JIniPhi, JIniTheta, 0.);
249 rotateVector(S2, S2, JIniPhi, JIniTheta, 0.);
250 rotateVector(Jh, Jh, JIniPhi, JIniTheta, 0.);
251
252 // initialise the evolution variables.
253 yinit[0] = v = mParams->v0; // PN expansion parameter
254 yinit[1] = LNh[0]; // unit vector along Newt. ang. momtm
255 yinit[2] = LNh[1]; // unit vector along Newt. ang. momtm
256 yinit[3] = LNh[2]; // unit vector along Newt. ang. momtm
257 yinit[4] = S1[0]; // spin vector of the first object
258 yinit[5] = S1[1]; // spin vector of the first object
259 yinit[6] = S1[2]; // spin vector of the first object
260 yinit[7] = S2[0]; // spin vector of the second object
261 yinit[8] = S2[1]; // spin vector of the second object
262 yinit[9] = S2[2]; // spin vector of the second object
263
264 /* initialize the integrator */
265 integrator = XLALAdaptiveRungeKutta4Init(LAL_NUM_ST5_VARIABLES,
266 XLALSimInspiralSpinTaylorT5Derivatives,
267 XLALSimInspiralSpinTaylorT5StoppingTest,
268 LAL_ST5_ABSOLUTE_TOLERANCE, LAL_ST5_RELATIVE_TOLERANCE);
269
270 if( !integrator ) {
271 XLALPrintError("XLAL Error - %s: Cannot allocate integrator\n", __func__);
272 XLAL_ERROR(XLAL_EFUNC);
273 }
274
275 /* stop the integration only when the test is true */
276 integrator->stopontestonly = 0;
277
278 /* run the ntegration; note: time is measured in \hat{t} = t / M */
279 lenReturn = XLALAdaptiveRungeKutta4Hermite(integrator, (void *) mParams, yinit,
280 0.0, dataLength*deltaT, deltaT, &yout);
281
282 intStatus = integrator->returncode;
283 XLALAdaptiveRungeKuttaFree(integrator);
284
285 if (!lenReturn) {
286 XLALPrintError("XLAL Error - %s: integration failed with errorcode %d.\n", __func__, intStatus);
287 XLAL_ERROR(XLAL_EFUNC);
288 }
289
290 /* Print warning about abnormal termination */
291 if (intStatus != 0 && intStatus != LALSIMINSPIRAL_ST5_TEST_ENERGY
292 && intStatus != LALSIMINSPIRAL_ST5_TEST_FREQBOUND) {
293 XLALPrintWarning("XLAL Warning - %s: integration terminated with code %d.\n \
294 Waveform parameters were m1 = %e, m2 = %e, S1 = [%e,%e,%e], S2 = [%e,%e,%e], LNh = [%e,%e,%e]\n",
295 __func__, intStatus, m1/LAL_MSUN_SI, m2/LAL_MSUN_SI, S1[0], S1[1], S1[2], S2[0],
296 S2[1], S2[2], LNh[0], LNh[1], LNh[2]);
297 }
298
299 /* allocate memory for vectors storing the PN parameter, orbital phase
300 and the unit vector along the Newtonian orbital angular momentum */
301 V = XLALCreateREAL8TimeSeries( "PN Parameter", &tStart, 0.0, deltaT, &lalStrainUnit, lenReturn);
302 orbPhase = XLALCreateREAL8TimeSeries( "Orbital Phase", &tStart, 0.0, deltaT, &lalStrainUnit, lenReturn);
303 LNhxVec = XLALCreateREAL8TimeSeries( "Uvec along Newt ang mom (X comp)", &tStart, 0., deltaT, &lalStrainUnit, lenReturn);
304 LNhyVec = XLALCreateREAL8TimeSeries( "Uvec along Newt ang mom (Y comp)", &tStart, 0., deltaT, &lalStrainUnit, lenReturn);
305 LNhzVec = XLALCreateREAL8TimeSeries( "Uvec along Newt ang mom (Z comp)", &tStart, 0., deltaT, &lalStrainUnit, lenReturn);
306 S1xVec = XLALCreateREAL8TimeSeries( "Spin ang mom of body 1 (X comp)", &tStart, 0., deltaT, &lalStrainUnit, lenReturn);
307 S1yVec = XLALCreateREAL8TimeSeries( "Spin ang mom of body 1 (Y comp)", &tStart, 0., deltaT, &lalStrainUnit, lenReturn);
308 S1zVec = XLALCreateREAL8TimeSeries( "Spin ang mom of body 1 (Z comp)", &tStart, 0., deltaT, &lalStrainUnit, lenReturn);
309 S2xVec = XLALCreateREAL8TimeSeries( "Spin ang mom of body 2 (X comp)", &tStart, 0., deltaT, &lalStrainUnit, lenReturn);
310 S2yVec = XLALCreateREAL8TimeSeries( "Spin ang mom of body 2 (Y comp)", &tStart, 0., deltaT, &lalStrainUnit, lenReturn);
311 S2zVec = XLALCreateREAL8TimeSeries( "Spin ang mom of body 2 (Z comp)", &tStart, 0., deltaT, &lalStrainUnit, lenReturn);
312
313 memset(V->data->data, 0, lenReturn*sizeof(REAL8));
314 memset(orbPhase->data->data, 0, lenReturn*sizeof(REAL8));
315 memset(LNhxVec->data->data, 0, lenReturn*sizeof(REAL8));
316 memset(LNhyVec->data->data, 0, lenReturn*sizeof(REAL8));
317 memset(LNhzVec->data->data, 0, lenReturn*sizeof(REAL8));
318 memset(S1xVec->data->data, 0, lenReturn*sizeof(REAL8));
319 memset(S1yVec->data->data, 0, lenReturn*sizeof(REAL8));
320 memset(S1zVec->data->data, 0, lenReturn*sizeof(REAL8));
321 memset(S2xVec->data->data, 0, lenReturn*sizeof(REAL8));
322 memset(S2yVec->data->data, 0, lenReturn*sizeof(REAL8));
323 memset(S2zVec->data->data, 0, lenReturn*sizeof(REAL8));
324
325 /* Copy dynamical variables from yout array to output time series.
326 note that yout[0:dataLength=1] is time */
327 for (i = 0; i<lenReturn; i++) {
328 V->data->data[i] = yout->data[lenReturn+ i];
329 LNhxVec->data->data[i] = yout->data[2*lenReturn+ i];
330 LNhyVec->data->data[i] = yout->data[3*lenReturn+ i];
331 LNhzVec->data->data[i] = yout->data[4*lenReturn+ i];
332 S1xVec->data->data[i] = yout->data[5*lenReturn+ i];
333 S1yVec->data->data[i] = yout->data[6*lenReturn+ i];
334 S1zVec->data->data[i] = yout->data[7*lenReturn+ i];
335 S2xVec->data->data[i] = yout->data[8*lenReturn+ i];
336 S2yVec->data->data[i] = yout->data[9*lenReturn+ i];
337 S2zVec->data->data[i] = yout->data[10*lenReturn+ i];
338 }
339
340 XLALDestroyREAL8Array(yout);
341
342 /* compute the orbital phase */
343 if (computeOrbitalPhase(mParams, V, LNhxVec, LNhyVec, LNhzVec, S1xVec, S1yVec, S1zVec,
344 S2xVec, S2yVec, S2zVec, orbPhase) != XLAL_SUCCESS) {
345 XLALPrintError("XLAL Error - %s: Unable to compute the orbital phase .\n", __func__);
346 return XLAL_FAILURE;
347 }
348
349 /* project the polarizations into the radiation frame */
350 switch (amplitudeO) {
351 case 0:
352 if (polarizationsInRadiationFrame(hplus, hcross, V, orbPhase, LNhxVec,
353 LNhyVec, LNhzVec, m1, m2, r, incAngle) != XLAL_SUCCESS) {
354 XLALPrintError("XLAL Error - %s: Unable to project the waveforms into radiation frame.\n", __func__);
355 return XLAL_FAILURE;
356 }
357 break;
358 default:
359 XLALPrintError("XLAL Error - %s: Polarizations are currently computed only in amplitudeO = 0\n", __func__);
360 return XLAL_FAILURE;
361 }
362
363 /* clear memory */
364 XLALDestroyREAL8TimeSeries(V);
365 XLALDestroyREAL8TimeSeries(orbPhase);
366 XLALDestroyREAL8TimeSeries(LNhxVec);
367 XLALDestroyREAL8TimeSeries(LNhyVec);
368 XLALDestroyREAL8TimeSeries(LNhzVec);
369 XLALDestroyREAL8TimeSeries(S1xVec);
370 XLALDestroyREAL8TimeSeries(S1yVec);
371 XLALDestroyREAL8TimeSeries(S1zVec);
372 XLALDestroyREAL8TimeSeries(S2xVec);
373 XLALDestroyREAL8TimeSeries(S2yVec);
374 XLALDestroyREAL8TimeSeries(S2zVec);
375
376 return XLAL_SUCCESS;
377}
378
379/** @} */
380
381/*
382 * Evolution of dynamical variables in the SpinTaylorT5
383 */
384static int XLALSimInspiralSpinTaylorT5Derivatives(
385 double t,
386 const double y[],
387 double dydt[],
388 void *mParams
389 ) {
390
391 REAL8 LNh[3], S1[3], S2[3], chi1[3], chi2[3], Omega1[3], Omega2[3];
392 REAL8 dS1byDt[3], dS2byDt[3], dLNhByDt[3];
393 REAL8 m, eta, q, delta, v, dEbyF, om0, v2, S1dotLNh, S2dotLNh, Lmag;
394 UINT4 phaseO, i, onePFpnFlag = 1, onePNFlag = 1;
395 SpinTaylorT5Params *params = (SpinTaylorT5Params*) mParams;
396 UNUSED(t);
397
398 /* binary parameters */
399 m = params->totalMass; /* in seconds */
400 eta = params->eta;
401 delta = params->delta;
402 phaseO = params->phaseO;
403 q = params->massRatio;
404
405 /* evolution variables. */
406 v = y[0]; /* PN expansion parameter: (m \omega_orb)^(1/3) */
407 LNh[0] = y[1]; /* unit vector in the direction of Newt. ang. momentum (LNh_x) */
408 LNh[1] = y[2]; /* unit vector in the direction of Newt. ang. momentum (LNh_y) */
409 LNh[2] = y[3]; /* unit vector in the direction of Newt. ang. momentum (LNh_z) */
410 S1[0] = y[4]; /* spin vector of the first object */
411 S1[1] = y[5]; /* spin vector of the first object */
412 S1[2] = y[6]; /* spin vector of the first object */
413 S2[0] = y[7]; /* spin vector of the second object */
414 S2[1] = y[8]; /* spin vector of the second object */
415 S2[2] = y[9]; /* spin vector of the second object */
416
417 /* comptute the dimensionless spins */
418 chi1[0] = S1[0]/params->m1Sqr;
419 chi1[1] = S1[1]/params->m1Sqr;
420 chi1[2] = S1[2]/params->m1Sqr;
421 chi2[0] = S2[0]/params->m2Sqr;
422 chi2[1] = S2[1]/params->m2Sqr;
423 chi2[2] = S2[2]/params->m2Sqr;
424
425 /* compute energy and flux function */
426 dEbyF = XLALdEnergyByFluxSpinPrec(v, chi1, chi2, LNh, params);
427
428 /* spin evolution equations - BBF Eq.(7.5 - 7.8) and Kesden et al 2010 Eq (2.2) */
429 v2 = v*v;
430 om0 = v2*v2*v/m;
431 S1dotLNh = dotProduct(S1, LNh);
432 S2dotLNh = dotProduct(S2, LNh);
433
434 /* if low PN orders are used, set the higher order PN terms to zero */
435 if (phaseO < 6) onePFpnFlag = 0; /* 1.5 PN correction to leading order spin (3PN) */
436 if (phaseO < 5) onePNFlag = 0; /* next to leading order in spin (2.5PN) */
437 if (phaseO < 3) om0 = 0; /* leading order spin (1.5PN) */
438
439 for (i=0; i<3; i++) {
440 Omega1[i] = om0*((0.75 + eta/2. - 0.75*delta) * LNh[i]
441 + onePNFlag * v *(-3.*(S2dotLNh + S1dotLNh*q) * LNh[i] + S2[i])/(2.*params->mSqr)
442 + onePFpnFlag * v2 *(0.5625 + 1.25*eta - params->etaSqr/24. - 0.5625*delta
443 + 0.625*params->deltaEta)*LNh[i]);
444
445 Omega2[i] = om0*((0.75 + eta/2. + 0.75*delta) * LNh[i]
446 + onePNFlag * v *(-3.*(S1dotLNh + S2dotLNh/q) * LNh[i] + S1[i])/(2.*params->mSqr)
447 + onePFpnFlag * v2 *(0.5625 + 1.25*eta - params->etaSqr/24. + 0.5625*delta
448 - 0.625*params->deltaEta)*LNh[i]);
449 }
450
451 crossProduct(Omega1, S1, dS1byDt);
452 crossProduct(Omega2, S2, dS2byDt);
453
454 /* angular momentum evolution eqn (BCV2, Eq (9) */
455 Lmag = (eta*params->mSqr/v)*(1. + (1.5+eta/6.)*v2);
456 dLNhByDt[0] = -(dS1byDt[0] + dS2byDt[0])/Lmag;
457 dLNhByDt[1] = -(dS1byDt[1] + dS2byDt[1])/Lmag;
458 dLNhByDt[2] = -(dS1byDt[2] + dS2byDt[2])/Lmag;
459
460 /* evolve the variables */
461 dydt[0] = -1./(dEbyF*m);
462 dydt[1] = dLNhByDt[0];
463 dydt[2] = dLNhByDt[1];
464 dydt[3] = dLNhByDt[2];
465 dydt[4] = dS1byDt[0];
466 dydt[5] = dS1byDt[1];
467 dydt[6] = dS1byDt[2];
468 dydt[7] = dS2byDt[0];
469 dydt[8] = dS2byDt[1];
470 dydt[9] = dS2byDt[2];
471
472 if (dEbyF > 0.0) /* energy test fails! (dE/F > 0) */
473 return LALSIMINSPIRAL_ST5_TEST_ENERGY;
474 else
475 return GSL_SUCCESS;
476
477}
478
479
480/*
481 * Compute the re-expanded (dEnergy/dv)/Flux for generic spinning binaries
482 */
483static REAL8 XLALdEnergyByFluxSpinPrec(
484 double v,
485 double *chi1,
486 double *chi2,
487 double *LNh,
488 SpinTaylorT5Params *params
489 ) {
490
491 REAL8 chi_s[3], chi_a[3];
492 REAL8 chisDotLNh, chiaDotLNh, chisDotChia, chisSqr, chiaSqr;
493 REAL8 dEbF0 = 0., dEbF1 = 0., dEbF2 = 0., dEbF3 = 0., dEbF4 = 0., dEbF5 = 0., dEbF6 = 0.;
494 REAL8 dEbF6L = 0., dEbF7 = 0.;
495 REAL8 v2, v3, v4, v5, v6, v7, v9;
496
497 REAL8 eta = params->eta;
498 REAL8 delta = params->delta;
499 UINT4 phaseO = params->phaseO;
500 REAL8 eta_p2 = eta*eta;
501
502 /* symmetric and anti-symmetric combinations of spin vectors */
503 chi_s[0] = 0.5*(chi1[0]+chi2[0]);
504 chi_s[1] = 0.5*(chi1[1]+chi2[1]);
505 chi_s[2] = 0.5*(chi1[2]+chi2[2]);
506 chi_a[0] = 0.5*(chi1[0]-chi2[0]);
507 chi_a[1] = 0.5*(chi1[1]-chi2[1]);
508 chi_a[2] = 0.5*(chi1[2]-chi2[2]);
509
510 /* dot products */
511 chisDotLNh = dotProduct(chi_s, LNh);
512 chiaDotLNh = dotProduct(chi_a, LNh);
513 chisDotChia = dotProduct(chi_s, chi_a);
514 chisSqr = dotProduct(chi_s, chi_s);
515 chiaSqr = dotProduct(chi_a, chi_a);
516
517 switch (phaseO) {
518 case -1: /* highest available PN order */
519 case 8: /* pseudo 4PN */
520 case 7: /* 3.5 PN */
521 dEbF7 = params->dEbF7NonSpin;
522#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
523 __attribute__ ((fallthrough));
524#endif
525 case 6: /* 3 PN */
526 dEbF6 = params->dEbF6NonSpin;
527 dEbF6L = params->dEbF6LogNonSpin;
528#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
529 __attribute__ ((fallthrough));
530#endif
531 case 5: /* 2.5 PN */
532 dEbF5 = params->dEbF5NonSpin
533 + chiaDotLNh*delta*(72.71676587301587 + (7*eta)/2.)
534 + chisDotLNh*(72.71676587301587 - (1213*eta)/18. - (17*eta_p2)/2.);
535#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
536 __attribute__ ((fallthrough));
537#endif
538 case 4: /* 2 PN */
539 dEbF4 = params->dEbF4NonSpin
540 + (233*chisDotChia*delta)/48. - (719*chiaDotLNh*chisDotLNh*delta)/48.
541 + chiaSqr*(2.4270833333333335 - 10*eta) + chisDotLNh*chisDotLNh*(-7.489583333333333 - eta/24.)
542 + chisSqr*(2.4270833333333335 + (7*eta)/24.) + chiaDotLNh*chiaDotLNh*(-7.489583333333333 + 30*eta);
543#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
544 __attribute__ ((fallthrough));
545#endif
546 case 3: /* 1.5 PN */
547 dEbF3 = params->dEbF3NonSpin
548 + (113*chiaDotLNh*delta)/12. + chisDotLNh*(9.416666666666666 - (19*eta)/3.);
549#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
550 __attribute__ ((fallthrough));
551#endif
552 case 2: /* 1 PN */
553 dEbF2 = params->dEbF2NonSpin;
554#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
555 __attribute__ ((fallthrough));
556#endif
557 case 1: /* 0.5 PN */
558 dEbF1 = params->dEbF1NonSpin;
559#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
560 __attribute__ ((fallthrough));
561#endif
562 case 0: /* Newtonian */
563 dEbF0 = params->dEbF0NonSpin;
564 break;
565 default:
566 XLALPrintError("XLAL Error - %s: Invalid phase. PN order %d\n", __func__, phaseO);
567 XLAL_ERROR(XLAL_EINVAL);
568 break;
569 }
570
571 v2 = v*v; v3 = v2*v; v4 = v2*v2; v5 = v4*v; v6 = v3*v3; v7 = v6*v; v9 = v7*v2;
572
573 /* compute the re-expanded dEnerby/Flux function */
574 return (dEbF0/v9)*(1. + dEbF1*v + dEbF2*v2 + dEbF3*v3 + dEbF4*v4 + dEbF5*v5
575 + (dEbF6 + dEbF6L*log(v))*v6 + dEbF7*v7);
576
577}
578
579/*
580 * Compute the GW polarizations in the radiation frame
581 */
582static int polarizationsInRadiationFrame(
583 REAL8TimeSeries **hplus, /**< +-polarization waveform [returned] */
584 REAL8TimeSeries **hcross, /**< x-polarization waveform [returned] */
585 REAL8TimeSeries *V, /**< post-Newtonian parameter */
586 REAL8TimeSeries *Phi, /**< orbital phase */
587 REAL8TimeSeries *LNhxVec, /**< unit orbital ang. mom. x comp. */
588 REAL8TimeSeries *LNhyVec, /**< unit orbital ang. mom. y comp. */
589 REAL8TimeSeries *LNhzVec, /**< unit orbital ang. mom. z comp. */
590 REAL8 m1, /**< mass1 (KG) */
591 REAL8 m2, /**< mass2 (KG) */
592 REAL8 r, /**< distance (meter) */
593 REAL8 Theta /**< angle between the initial total ang momentum and line of sight. */
594 ) {
595
596 REAL8 deltaPhi, twoPhiS, alphaDot;
597 REAL8 LNx, LNy, LNz, v2, LNx_p2, LNy_p2, LNz_p2, LN_xz;
598 UINT4 i;
599 REAL8Vector *alphaVec=NULL, *iVec=NULL;
600 REAL8 cos_TwoPhiS, sin_TwoPhiS;
601
602 UINT4 dataLength = V->data->length;
603 REAL8 deltaT = V->deltaT;
604 REAL8 m = (m1+m2)*LAL_MTSUN_SI/LAL_MSUN_SI;
605 REAL8 eta = m1*m2/(m1+m2)/(m1+m2);
606 REAL8 amp = 4.*eta*m/(r/LAL_C_SI);
607
608 /* allocate memory for the arrays */
609 alphaVec = XLALCreateREAL8Vector(dataLength);
610 iVec = XLALCreateREAL8Vector(dataLength);
611 memset(alphaVec->data, 0, alphaVec->length * sizeof( REAL8 ));
612 memset(iVec->data, 0, iVec->length * sizeof( REAL8 ));
613
614 /* angles describing the time-evolution of the orientation of the orbital*/
615 /* angular momentum in the source frame */
616 for (i=0; i<dataLength; i++) {
617 if (V->data->data[i]) {
618 alphaVec->data[i] = atan2(LNhyVec->data->data[i], LNhxVec->data->data[i]); // called alpha in BCV2
619 iVec->data[i] = acos(LNhzVec->data->data[i]); // called i in BCV2
620 }
621 }
622
623 /* unwrap the angles */
624 XLALREAL8VectorUnwrapAngle(alphaVec, alphaVec);
625 XLALREAL8VectorUnwrapAngle(iVec, iVec);
626
627 deltaPhi = 0.;
628
629 /* Allocate polarization vectors and set to 0 */
630 *hplus = XLALCreateREAL8TimeSeries( "H_PLUS", &V->epoch,
631 0.0, V->deltaT, &lalStrainUnit, V->data->length );
632 *hcross = XLALCreateREAL8TimeSeries( "H_CROSS", &V->epoch,
633 0.0, V->deltaT, &lalStrainUnit, V->data->length );
634 if ( ! hplus || ! hcross ) XLAL_ERROR(XLAL_EFUNC);
635 memset((*hplus)->data->data, 0, (*hplus)->data->length*sizeof(*(*hplus)->data->data));
636 memset((*hcross)->data->data, 0, (*hcross)->data->length*sizeof(*(*hcross)->data->data));
637
638 REAL8 cos_Theta = cos(Theta);
639 REAL8 cos_TwoTheta = cos(2.*Theta);
640 REAL8 sin_Theta = sin(Theta);
641 REAL8 sin_TwoTheta = sin(2.*Theta);
642
643 for (i=0; i<dataLength; i++) {
644
645 if (V->data->data[i]) {
646
647 LNx = LNhxVec->data->data[i];
648 LNy = LNhyVec->data->data[i];
649 LNz = LNhzVec->data->data[i];
650 v2 = V->data->data[i]*V->data->data[i];
651
652 LNx_p2 = LNx*LNx;
653 LNy_p2 = LNy*LNy;
654 LNz_p2 = LNz*LNz;
655 LN_xz = LNx*LNz;
656
657 /* compute d alpha/dt*/
658 if (i == dataLength-1) alphaDot = (alphaVec->data[i] - alphaVec->data[i-1])/deltaT;
659 else alphaDot = (alphaVec->data[i+1] - alphaVec->data[i])/deltaT;
660
661 /* \int \cos(i) \alpha_dot deltaT -- integrating Eq.(18) of BCV2 */
662 deltaPhi += LNz*alphaDot*deltaT;
663
664 /* the phase measured in the detector */
665 twoPhiS = 2.*(Phi->data->data[i] - deltaPhi);
666
667 cos_TwoPhiS = cos(twoPhiS);
668 sin_TwoPhiS = sin(twoPhiS);
669
670 /* compute the polarizations. For derivation, see the Mathematica notebook Precession_ProjectiontoRadiationFrame.nb*/
671 /* (based on Sec IIC of BCV2)*/
672 (*hplus)->data->data[i] = -(0.5*amp*v2*(cos_TwoPhiS*(-LNx_p2 + LNy_p2*LNz_p2 + ((LNy
673 + LN_xz)*cos_Theta + sin_Theta
674 - LNz_p2*sin_Theta)*((LNy - LN_xz)*cos_Theta + (-1. + LNz_p2)*sin_Theta))
675 + LNy*sin_TwoPhiS*(LN_xz*(3. + cos_TwoTheta) - (-1. + LNz_p2)*sin_TwoTheta)))/(-1. + LNz_p2);
676
677 (*hcross)->data->data[i] = -(amp*v2*(cos_Theta*(-(LNx*LNy*(1 + LNz_p2)*cos_TwoPhiS)
678 + (-LNx_p2 + LNy_p2)*LNz*sin_TwoPhiS)
679 + (-1. + LNz_p2)*(LNy*LNz*cos_TwoPhiS + LNx*sin_TwoPhiS)*sin_Theta))/(-1. + LNz_p2);
680
681 }
682
683 }
684
685 XLALDestroyREAL8Vector(alphaVec);
686 XLALDestroyREAL8Vector(iVec);
687
688 return XLAL_SUCCESS;
689}
690
691/*
692 * Compute the orbital phase as an explicit function of v (TaylorT2 approximant)
693 * Ref. Eq.(3.4) of http://arxiv.org/abs/1107.1267
694 */
695static int computeOrbitalPhase(
696 SpinTaylorT5Params *params,
697 REAL8TimeSeries *V,
698 REAL8TimeSeries *LNhxVec,
699 REAL8TimeSeries *LNhyVec,
700 REAL8TimeSeries *LNhzVec,
701 REAL8TimeSeries *S1xVec,
702 REAL8TimeSeries *S1yVec,
703 REAL8TimeSeries *S1zVec,
704 REAL8TimeSeries *S2xVec,
705 REAL8TimeSeries *S2yVec,
706 REAL8TimeSeries *S2zVec,
707 REAL8TimeSeries *orbPhase
708 ) {
709
710 REAL8 chi_s[3], chi_a[3], LNh[3];
711 REAL8 chisDotLNh, chiaDotLNh, chisDotChia, chisSqr, chiaSqr;
712 REAL8 v, v2, v3, v4, v5, v6, v7;
713 REAL8 phi0 = 0., phi1 = 0., phi2 = 0., phi3 = 0., phi4 = 0., phi5 = 0.;
714 REAL8 phi6 = 0., phi6L = 0., phi7 = 0.;
715 UINT4 i;
716
717 REAL8 eta = params->eta;
718 REAL8 delta = params->delta;
719 UINT4 phaseO = params->phaseO;
720 REAL8 eta_p2 = eta*eta;
721
722 for (i = 0; i<V->data->length; i++) {
723
724 if (V->data->data[i]) {
725
726 /* symmetric and anti-symmetric combinations of spin vectors */
727 chi_s[0] = 0.5*(S1xVec->data->data[i]/params->m1Sqr + S2xVec->data->data[i]/params->m2Sqr);
728 chi_s[1] = 0.5*(S1yVec->data->data[i]/params->m1Sqr + S2yVec->data->data[i]/params->m2Sqr);
729 chi_s[2] = 0.5*(S1zVec->data->data[i]/params->m1Sqr + S2zVec->data->data[i]/params->m2Sqr);
730 chi_a[0] = 0.5*(S1xVec->data->data[i]/params->m1Sqr - S2xVec->data->data[i]/params->m2Sqr);
731 chi_a[1] = 0.5*(S1yVec->data->data[i]/params->m1Sqr - S2yVec->data->data[i]/params->m2Sqr);
732 chi_a[2] = 0.5*(S1zVec->data->data[i]/params->m1Sqr - S2zVec->data->data[i]/params->m2Sqr);
733
734 LNh[0] = LNhxVec->data->data[i];
735 LNh[1] = LNhyVec->data->data[i];
736 LNh[2] = LNhzVec->data->data[i];
737
738 /* dot products */
739 chisDotLNh = dotProduct(chi_s, LNh);
740 chiaDotLNh = dotProduct(chi_a, LNh);
741 chisDotChia = dotProduct(chi_s, chi_a);
742 chisSqr = dotProduct(chi_s, chi_s);
743 chiaSqr = dotProduct(chi_a, chi_a);
744
745 v = V->data->data[i];
746
747 switch (phaseO) {
748 case -1: /* highest available PN order */
749 case 8: /* pseudo 4PN */
750 case 7: /* 3.5 PN */
751 phi7 = params->phiOrb7NonSpin;
752#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
753 __attribute__ ((fallthrough));
754#endif
755 case 6: /* 3 PN */
756 phi6 = params->phiOrb6NonSpin;
757 phi6L = params->phiOrb6LogNonSpin;
758#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
759 __attribute__ ((fallthrough));
760#endif
761 case 5: /* 2.5 PN */
762 phi5 = params->phiOrb5LogNonSpin
763 + chiaDotLNh*((-732985*delta)/2016. - (35*delta*eta)/2.)
764 + chisDotLNh*(-363.58382936507934 + (6065*eta)/18. + (85*eta_p2)/2.);
765#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
766 __attribute__ ((fallthrough));
767#endif
768 case 4: /* 2 PN */
769 phi4 = params->phiOrb4NonSpin
770 + (1165*chisDotChia*delta)/48. - (3595*chiaDotLNh*chisDotLNh*delta)/48.
771 + chiaSqr*(12.135416666666666 - 50*eta) + chisDotLNh*chisDotLNh*(-37.447916666666664 - (5*eta)/24.)
772 + chisSqr*(12.135416666666666 + (35*eta)/24.) + chiaDotLNh*chiaDotLNh*(-37.447916666666664 + 150*eta);
773#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
774 __attribute__ ((fallthrough));
775#endif
776 case 3: /* 1.5 PN */
777 phi3 = params->phiOrb3NonSpin
778 + (565*chiaDotLNh*delta)/24. + chisDotLNh*(23.541666666666668 - (95*eta)/6.);
779#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
780 __attribute__ ((fallthrough));
781#endif
782 case 2: /* 1 PN */
783 phi2 = params->phiOrb2NonSpin;
784#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
785 __attribute__ ((fallthrough));
786#endif
787 case 1: /* 0.5 PN */
788 phi1 = params->phiOrb1NonSpin;
789#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
790 __attribute__ ((fallthrough));
791#endif
792 case 0: /* Newtonian */
793 phi0 = params->phiOrb0NonSpin;
794 break;
795 default:
796 XLALPrintError("XLAL Error - %s: Invalid phase. PN order %d\n", __func__, phaseO);
797 XLAL_ERROR(XLAL_EINVAL);
798 break;
799 }
800
801 v2 = v*v; v3 = v2*v; v4 = v2*v2; v5 = v4*v; v6 = v3*v3; v7 = v6*v;
802
803 /* compute the orbital phase */
804 orbPhase->data->data[i] = (phi0/v5)*(1. + phi1*v + phi2*v2 + phi3*v3 + phi4*v4 + phi5*log(v)*v5
805 + (phi6 + phi6L*log(v))*v6 + phi7*v7);
806 }
807 }
808
809 /* apply a constant phase-shift such that phi_orb(v) = phiRef/2 at v = v_0 */
810 REAL8 phiOrb0 = orbPhase->data->data[0]-params->phiRef/2.;
811 for (i = 0; i<orbPhase->data->length; i++) orbPhase->data->data[i] -= phiOrb0;
812
813 return XLAL_SUCCESS;
814}
815
816/**
817 * Initialize the non-dynamical variables required for the SpinTaylorT5 evolution
818 */
819static int spinTaylorT5Init(
820 REAL8 m1,
821 REAL8 m2,
822 REAL8 fStart,
823 REAL8 deltaT,
824 REAL8 phiRef,
825 UINT4 phaseO,
826 SpinTaylorT5Params *mParams
827 ) {
828
829 /* various parameters describing the binary */
830 mParams->mass1 = m1*LAL_MTSUN_SI/LAL_MSUN_SI; /* m1 in seconds */
831 mParams->mass2 = m2*LAL_MTSUN_SI/LAL_MSUN_SI; /* m2 in seconds */
832 mParams->totalMass = mParams->mass1+mParams->mass2; /* m1+m2 in seconds */
833 mParams->eta = mParams->mass1*mParams->mass2/(mParams->totalMass*mParams->totalMass);
834 mParams->delta = (mParams->mass1-mParams->mass2)/mParams->totalMass; /* (m1-m2)/(m1+m2) */
835 mParams->phaseO = phaseO; /* twice the PN phase order */
836 mParams->massRatio = mParams->mass2/mParams->mass1; /* q = m2/m1, where m2 > m1 */
837 mParams->m1Sqr = mParams->mass1*mParams->mass1; /* m1^2 in seconds^2 */
838 mParams->m2Sqr = mParams->mass2*mParams->mass2; /* m2^2 in seconds^2 */
839 mParams->mSqr = mParams->totalMass*mParams->totalMass; /* m^2 in seconds^2 */
840 mParams->etaSqr = mParams->eta*mParams->eta; /* eta^2 */
841 mParams->deltaEta = mParams->delta*mParams->eta; /* delta * eta */
842 mParams->v0 = cbrt(LAL_PI*mParams->totalMass*fStart); /* starting value for the PN parameter */
843 mParams->vMax = cbrt(0.4*LAL_PI*mParams->totalMass/deltaT); /* set an emperical maximum on the PN parameter (0.9 f_Nyquist) */
844 mParams->phiRef = phiRef;
845
846 /* coefficients of the reexpanded dEnergy/flux function (non-spinning) */
847 mParams->dEbF0NonSpin = -5./(32.*mParams->eta);
848 mParams->dEbF1NonSpin = 0.;
849 mParams->dEbF2NonSpin = 2.2113095238095237 + (11*mParams->eta)/4.;
850 mParams->dEbF3NonSpin = -4*LAL_PI;
851 mParams->dEbF4NonSpin = 3.010315295099521 + (5429*mParams->eta)/1008. + (617*mParams->etaSqr)/144.;
852 mParams->dEbF5NonSpin = (-7729*LAL_PI)/672. + (13*mParams->eta*LAL_PI)/8.;
853 mParams->dEbF6NonSpin = -115.2253249962622 + (3147553127*mParams->eta)/1.2192768e7 - (15211*mParams->etaSqr)/6912.
854 + (25565*mParams->etaSqr*mParams->eta)/5184. + (1712*LAL_GAMMA)/105. + (32*LAL_PISQR)/3.
855 - (451*mParams->eta*LAL_PISQR)/48. + (1712*log(4))/105.;
856 mParams->dEbF6LogNonSpin = 16.304761904761904;
857 mParams->dEbF7NonSpin = (-15419335*LAL_PI)/1.016064e6 - (75703*mParams->eta*LAL_PI)/6048.
858 + (14809*mParams->etaSqr*LAL_PI)/3024.;
859
860 /* coefficients of the PN expansion of the orbital phase (non-spinning) */
861 mParams->phiOrb0NonSpin = -1/(32.*mParams->eta);
862 mParams->phiOrb1NonSpin = 0.;
863 mParams->phiOrb2NonSpin = 3.685515873015873 + (55*mParams->eta)/12.;
864 mParams->phiOrb3NonSpin = -10*LAL_PI;
865 mParams->phiOrb4NonSpin = 15.051576475497606 + (27145*mParams->eta)/1008. + (3085*mParams->etaSqr)/144.;
866 mParams->phiOrb5LogNonSpin = (38645*LAL_PI)/672. - (65*mParams->eta*LAL_PI)/8.;
867 mParams->phiOrb6NonSpin = 657.6504345051205 - (15737765635*mParams->eta)/1.2192768e7
868 + (76055*mParams->etaSqr)/6912. - (127825*mParams->etaSqr*mParams->eta)/5184.
869 - (1712*LAL_GAMMA)/21. - (160*LAL_PISQR)/3. + (2255*mParams->eta*LAL_PISQR)/48.
870 - (1712*log(4))/21.;
871 mParams->phiOrb6LogNonSpin = -81.52380952380952;
872 mParams->phiOrb7NonSpin = (77096675*LAL_PI)/2.032128e6 + (378515*mParams->eta*LAL_PI)/12096.
873 - (74045*mParams->etaSqr*LAL_PI)/6048.;
874
875 return XLAL_SUCCESS;
876}
877
878/*
879 * Internal function called by the integration routine.
880 * Stops the integration if
881 * 1) The energy decreases with increasing orbital frequency
882 * 3) The orbital frequency becomes infinite
883 */
884static int XLALSimInspiralSpinTaylorT5StoppingTest(
885 double t,
886 const double y[],
887 double dy[],
888 void *mParams
889 ) {
890
891 SpinTaylorT5Params *params = (SpinTaylorT5Params*) mParams;
892 UNUSED(t);
893 REAL8 LNh[3], S1[3], S2[3], chi1[3], chi2[3], v, dEbyF;
894
895 /* evolution variables. */
896 v = y[0]; /* PN expansion parameter: (m \omega_orb)^(1/3) */
897 LNh[0] = y[1]; /* unit vector in the direction of Newt. ang. momentum (LNh_x) */
898 LNh[1] = y[2]; /* unit vector in the direction of Newt. ang. momentum (LNh_y) */
899 LNh[2] = y[3]; /* unit vector in the direction of Newt. ang. momentum (LNh_z) */
900 S1[0] = y[4]; /* spin vector of the first object */
901 S1[1] = y[5]; /* spin vector of the first object */
902 S1[2] = y[6]; /* spin vector of the first object */
903 S2[0] = y[7]; /* spin vector of the second object */
904 S2[1] = y[8]; /* spin vector of the second object */
905 S2[2] = y[9]; /* spin vector of the second object */
906
907 /* comptute the dimensionless spins */
908 chi1[0] = S1[0]/params->m1Sqr;
909 chi1[1] = S1[1]/params->m1Sqr;
910 chi1[2] = S1[2]/params->m1Sqr;
911 chi2[0] = S2[0]/params->m2Sqr;
912 chi2[1] = S2[1]/params->m2Sqr;
913 chi2[2] = S2[2]/params->m2Sqr;
914
915 /* compute energy and flux function */
916 dEbyF = XLALdEnergyByFluxSpinPrec(v, chi1, chi2, LNh, params);
917
918 if (dEbyF > -1.e-10) /* energy test fails! (dE/F > 0) */
919 return LALSIMINSPIRAL_ST5_TEST_ENERGY;
920 else if (dy[0] < 0.0) /* dv/dt < 0! */
921 return LALSIMINSPIRAL_ST5_TEST_VDOT;
922 else if (isnan(v) || isinf(v)) /* v is nan! */
923 return LALSIMINSPIRAL_ST5_TEST_VNAN;
924 else if ((v < params->v0) || (v > params->vMax) || (v < 0.) || (v >= PN_V_MAX))
925 return LALSIMINSPIRAL_ST5_TEST_FREQBOUND;
926 else /* Step successful, continue integrating */
927 return GSL_SUCCESS;
928}
929
930static double dotProduct(double a[3], double b[3]) {
931
932 return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
933}
934
935static void crossProduct(double a[3], double b[3], double c[3]) {
936
937 c[0] = a[1]*b[2] - a[2]*b[1];
938 c[1] = a[2]*b[0] - a[0]*b[2];
939 c[2] = a[0]*b[1] - a[1]*b[0];
940}
941
942static void rotateVector(double a[3], double b[3], double phi, double theta, double psi) {
943/* rotate vector "a" by three Euler angles phi, theta, psi around axes
944 * z, y and x, respectively. b = R_x(psi) R_y(theta) R_z(phi) a */
945
946 double x = a[0];
947 double y = a[1];
948 double z = a[2];
949
950 b[0] = x*cos(phi)*cos(theta) - y*cos(theta)*sin(phi)
951 + z*sin(theta);
952
953 b[1] = -(z*cos(theta)*sin(psi)) + x*(cos(psi)*sin(phi)
954 + cos(phi)*sin(psi)*sin(theta))
955 + y*(cos(phi)*cos(psi) - sin(phi)*sin(psi)*sin(theta));
956
957 b[2] = z*cos(psi)*cos(theta) + x*(sin(phi)*sin(psi)
958 - cos(phi)*cos(psi)*sin(theta))
959 + y*(cos(phi)*sin(psi) + cos(psi)*sin(phi)*sin(theta));
960
961}
delta
static double double delta
Definition: LALSimIMRPhenomTHM_fits.h:105
c
#define c
Definition: LALSimInspiralHGimri.c:35
LALSIMINSPIRAL_ST5_TEST_VNAN
#define LALSIMINSPIRAL_ST5_TEST_VNAN
Definition: LALSimInspiralSpinTaylorT5duplicate.c:36
polarizationsInRadiationFrame
static int polarizationsInRadiationFrame(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, REAL8TimeSeries *V, REAL8TimeSeries *Phi, REAL8TimeSeries *LNhxVec, REAL8TimeSeries *LNhyVec, REAL8TimeSeries *LNhzVec, REAL8 m1, REAL8 m2, REAL8 r, REAL8 Theta)
Definition: LALSimInspiralSpinTaylorT5duplicate.c:582
LALSIMINSPIRAL_ST5_TEST_VDOT
#define LALSIMINSPIRAL_ST5_TEST_VDOT
Definition: LALSimInspiralSpinTaylorT5duplicate.c:34
spinTaylorT5Init
static int spinTaylorT5Init(REAL8 m1, REAL8 m2, REAL8 fStart, REAL8 deltaT, REAL8 phiRef, UINT4 phaseO, SpinTaylorT5Params *mParams)
Initialize the non-dynamical variables required for the SpinTaylorT5 evolution.
Definition: LALSimInspiralSpinTaylorT5duplicate.c:819
XLALdEnergyByFluxSpinPrec
static REAL8 XLALdEnergyByFluxSpinPrec(double v, double *chi1, double *chi2, double *LNh, SpinTaylorT5Params *params)
Definition: LALSimInspiralSpinTaylorT5duplicate.c:483
LAL_ST5_ABSOLUTE_TOLERANCE
#define LAL_ST5_ABSOLUTE_TOLERANCE
Definition: LALSimInspiralSpinTaylorT5duplicate.c:45
LALSIMINSPIRAL_ST5_TEST_ENERGY
#define LALSIMINSPIRAL_ST5_TEST_ENERGY
Definition: LALSimInspiralSpinTaylorT5duplicate.c:33
LAL_NUM_ST5_VARIABLES
#define LAL_NUM_ST5_VARIABLES
Definition: LALSimInspiralSpinTaylorT5duplicate.c:43
LALSIMINSPIRAL_ST5_TEST_FREQBOUND
#define LALSIMINSPIRAL_ST5_TEST_FREQBOUND
Definition: LALSimInspiralSpinTaylorT5duplicate.c:37
computeOrbitalPhase
static int computeOrbitalPhase(SpinTaylorT5Params *params, REAL8TimeSeries *V, REAL8TimeSeries *LNhxVec, REAL8TimeSeries *LNhyVec, REAL8TimeSeries *LNhzVec, REAL8TimeSeries *S1xVec, REAL8TimeSeries *S1yVec, REAL8TimeSeries *S1zVec, REAL8TimeSeries *S2xVec, REAL8TimeSeries *S2yVec, REAL8TimeSeries *S2zVec, REAL8TimeSeries *orbPhase)
Definition: LALSimInspiralSpinTaylorT5duplicate.c:695
crossProduct
static void crossProduct(double a[3], double b[3], double c[3])
Definition: LALSimInspiralSpinTaylorT5duplicate.c:935
dotProduct
static double dotProduct(double a[3], double b[3])
Definition: LALSimInspiralSpinTaylorT5duplicate.c:930
LAL_ST5_RELATIVE_TOLERANCE
#define LAL_ST5_RELATIVE_TOLERANCE
Definition: LALSimInspiralSpinTaylorT5duplicate.c:46
XLALSimInspiralSpinTaylorT5StoppingTest
static int XLALSimInspiralSpinTaylorT5StoppingTest(double t, const double y[], double dy[], void *mParams)
Definition: LALSimInspiralSpinTaylorT5duplicate.c:884
XLALSimInspiralSpinTaylorT5Derivatives
static int XLALSimInspiralSpinTaylorT5Derivatives(double t, const double y[], double dy[], void *mParams)
Definition: LALSimInspiralSpinTaylorT5duplicate.c:384
rotateVector
static void rotateVector(double a[3], double b[3], double phi, double theta, double psi)
Definition: LALSimInspiralSpinTaylorT5duplicate.c:942
i
double i
Definition: bh_ringdown.c:118
theta
double theta
Definition: bh_sphwf.c:118
check_series_macros.h
__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_C_SI
#define LAL_C_SI
LAL_MSUN_SI
#define LAL_MSUN_SI
LAL_PI
#define LAL_PI
LAL_MTSUN_SI
#define LAL_MTSUN_SI
LAL_GAMMA
#define LAL_GAMMA
LIGOTIMEGPSZERO
#define LIGOTIMEGPSZERO
REAL8
double REAL8
UINT4
uint32_t UINT4
XLALSimInspiralSpinTaylorT5duplicate
int XLALSimInspiralSpinTaylorT5duplicate(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, REAL8 phiRef, REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 fStart, REAL8 r, REAL8 s1x, REAL8 s1y, REAL8 s1z, REAL8 s2x, REAL8 s2y, REAL8 s2z, REAL8 incAngle, int phaseO, int amplitudeO)
Generate time-domain generic spinning PN waveforms in the SpinTaylorT5 approximaton.
Definition: LALSimInspiralSpinTaylorT5duplicate.c:166
r
static const INT4 r
m
static const INT4 m
q
static const INT4 q
a
static const INT4 a
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
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLALREAL8VectorUnwrapAngle
int XLALREAL8VectorUnwrapAngle(REAL8Vector *out, const REAL8Vector *in)
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_PRINT_DEPRECATION_WARNING
#define XLAL_PRINT_DEPRECATION_WARNING(replacement)
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EFUNC
XLAL_EINVAL
XLAL_EINVAL
XLAL_FAILURE
XLAL_FAILURE
XLALGPSAdd
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
y
list y
x
x
REAL8Array::data
REAL8 * data
REAL8TimeSeries::data
REAL8Sequence * data
REAL8Vector::data
REAL8 * data
SpinTaylorT5Params
Structure containing the non-dynamical coefficients needed to evolve a spinning, precessing binary an...
Definition: LALSimInspiralSpinTaylorT5duplicate.c:57
params
Definition: burst.c:245
params::phiRef
double phiRef
Definition: inspiral.c:231
V
double V
Definition: unicorn.c:25
deltaT
double deltaT
Definition: unicorn.c:24