LALSimulation  5.4.0.1-fe68b98
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 
57 typedef 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 */
94 static int XLALSimInspiralSpinTaylorT5Derivatives(
95  double t,
96  const double y[],
97  double dy[],
98  void *mParams
99  );
100 
101 static REAL8 XLALdEnergyByFluxSpinPrec(
102  double v,
103  double *chi1,
104  double *chi2,
105  double *LNh,
106  SpinTaylorT5Params *params
107  );
108 
109 static 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 
123 static 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 
137 static 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 
147 static int XLALSimInspiralSpinTaylorT5StoppingTest(
148  double t,
149  const double y[],
150  double dy[],
151  void *mParams
152  );
153 
154 static double dotProduct(double a[3], double b[3]);
155 static void crossProduct(double a[3], double b[3], double c[3]);
156 static 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  */
166 int 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  */
384 static 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  */
483 static 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  */
582 static 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  */
695 static 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  */
819 static 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  */
884 static 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 
930 static 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 
935 static 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 
942 static 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)
XLALAdaptiveRungeKutta4Hermite
int XLALAdaptiveRungeKutta4Hermite(LALAdaptiveRungeKuttaIntegrator *integrator, void *params, REAL8 *yinit, REAL8 tinit, REAL8 tend_in, REAL8 deltat, REAL8Array **yout)
XLALAdaptiveRungeKutta4Init
LALAdaptiveRungeKuttaIntegrator * XLALAdaptiveRungeKutta4Init(int dim, int(*dydt)(double t, const double y[], double dydt[], void *params), int(*stop)(double t, const double y[], double dydt[], void *params), double eps_abs, double eps_rel)
LAL_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
XLALCreateREAL8TimeSeries
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
XLALDestroyREAL8TimeSeries
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
lalStrainUnit
const LALUnit lalStrainUnit
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)
x
list x
y
list y
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