LALSimulation  5.4.0.1-fe68b98
LALSimIMRPSpinInspiralRD.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2011 Riccardo Sturani, John Veitch, Drew Keppel
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 <stdlib.h>
21 #include <math.h>
22 #include <gsl/gsl_linalg.h>
23 #include <gsl/gsl_interp.h>
24 #include <gsl/gsl_spline.h>
25 #include <lal/LALStdlib.h>
26 #include <lal/AVFactories.h>
27 #include <lal/SeqFactories.h>
28 #include <lal/Units.h>
29 #include <lal/TimeSeries.h>
30 #include <lal/LALConstants.h>
31 #include <lal/SeqFactories.h>
32 #include <lal/RealFFT.h>
33 #include <lal/SphericalHarmonics.h>
34 #include <lal/LALAdaptiveRungeKuttaIntegrator.h>
35 #include <lal/LALSimInspiral.h>
36 #include <lal/LALSimIMR.h>
37 #include <lal/LALSimSphHarmMode.h>
38 #include <lal/Date.h>
39 
40 //#include "LALSimIMRPhenSpin.h"
41 #include "LALSimPhenSpinRingDown.c"
42 #include "LALSimInspiralPNCoefficients.c"
43 
44 #define LALSIMINSPIRAL_PST4_TEST_ENERGY 1025
45 #define LALSIMINSPIRAL_PST4_TEST_OMEGADOT 1026
46 #define LALSIMINSPIRAL_PST4_TEST_COORDINATE 1027
47 #define LALSIMINSPIRAL_PST4_TEST_OMEGANAN 1028
48 #define LALSIMINSPIRAL_PST4_TEST_FREQBOUND 1029
49 #define LALSIMINSPIRAL_PST4_DERIVATIVE_OMEGANONPOS 1030
50 #define LALSIMINSPIRAL_PST4_TEST_OMEGAMATCH 1031
51 
52 #define LAL_NUM_PST4_VARIABLES 12
53 
54 #define LAL_PST4_ABSOLUTE_TOLERANCE 1.e-12
55 #define LAL_PST4_RELATIVE_TOLERANCE 1.e-12
56 
57 /* Minimum integration length */
58 #define minIntLen 16
59 /* For turning on debugging messages*/
60 #define DEBUG_RD 0
61 
62 #define nModes 8
63 #define RD_EFOLDS 10
64 
65 static REAL8 dsign(INT4 i){
66  if (i>=0) return 1.;
67  else return -1.;
68 }
69 
70 typedef struct tagLALSimInspiralPhenSpinTaylorT4Coeffs {
71  REAL8 M; ///< total mass in seconds
72  REAL8 eta; ///< symmetric mass ratio
73  REAL8 m1ByM; ///< m1 / M
74  REAL8 m2ByM; ///< m2 / M
75  REAL8 dt; ///< UNDOCUMENTED
76  REAL8 wdotnewt; ///< leading order coefficient of wdot = \f$\dot{\omega}\f$
77  REAL8 wdotcoeff[LAL_MAX_PN_ORDER]; ///< coeffs. of PN corrections to wdot
78  REAL8 wdotlogcoeff; ///< coefficient of log term in wdot
79  REAL8 Enewt; ///< coeffs. of PN corrections to energy
80  REAL8 Ecoeff[LAL_MAX_PN_ORDER]; ///< coeffs. of PN corrections to energy
81  REAL8 wdot3S1O, wdot3S2O; ///< non-dynamical 1.5PN SO corrections
82  REAL8 wdot4S1S2; ///< non-dynamical 2PN SS correction
83  REAL8 wdot4S1OS2O; ///< non-dynamical 2PN SS correction
84  REAL8 wdot4S1S1,wdot4S2S2; ///< non-dynamical 2PN SS correction
85  REAL8 wdot4S1OS1O,wdot4S2OS2O; ///< non-dynamical 2PN SS correction
86  REAL8 wdot4QMS1; ///< non-dynamical S1^2 2PN quadrupole-monopole correction
87  REAL8 wdot4QMS1O; ///< non-dynamical (S1.L)^2 2PN quadrupole-monopole correction
88  REAL8 wdot4QMS2; ///< non-dynamical S2^2 2PN quadrupole-monopole correction
89  REAL8 wdot4QMS2O; ///< non-dynamical (S2.L)^2 2PN quadrupole-monopole correction
90  REAL8 wdot5S1O, wdot5S2O; ///< non-dynamical 2.5PN SO corrections
91  REAL8 wdot6S1O, wdot6S2O; ///< non-dynamical 3PN SO corrections
92  REAL8 E3S1O, E3S2O; ///< non-dynamical 1.5PN SO corrections
93  REAL8 E4S1S2,E4S1OS2O; ///< non-dynamical 2PN SS correction
94  REAL8 E4QMS1; ///< non-dynamical S1^2 2PN quadrupole-monopole correction
95  REAL8 E4QMS1O;///< non-dynamical (S1.L)^2 2PN quadrupole-monopole correction
96  REAL8 E4QMS2; ///< non-dynamical S2^2 2PN quadrupole-monopole correction
97  REAL8 E4QMS2O;///< non-dynamical (S2.L)^2 2PN quadrupole-monopole correction
98  REAL8 E5S1O, E5S2O; ///< non-dynamical 2.5PN SO corrections
99  REAL8 E6S1O, E6S2O; ///< non-dynamical 2.5PN SO corrections
100  REAL8 LNhat3S1O, LNhat3S2O; ///< non-dynamical 1.5PN SO corrections
101  REAL8 LNhat4SS; ///< non-dynamical 2PN SS correction
102  REAL8 wdottidal5pn; ///< leading order tidal correction
103  REAL8 wdottidal6pn; ///< next to leading order tidal correction
104  REAL8 Etidal5pn; ///< leading order tidal correction to energy
105  REAL8 Etidal6pn; ///< next to leading order tidal correction to energy
106  REAL8 S1dot3, S2dot3; ///< UNDOCUMENTED
107  REAL8 Sdot4S2,Sdot4S2O; ///< UNDOCUMENTED
108  REAL8 S1dot4QMS1O,S2dot4QMS2O; ///< UNDOCUMENTED
109  REAL8 S1dot5, S2dot5; ///< UNDOCUMENTED
110  REAL8 S1dot6, S2dot6; ///< UNDOCUMENTED
111  REAL8 fStart; ///< starting GW frequency of integration
112  REAL8 fEnd; ///< ending GW frequency of integration
113 } LALSimInspiralPhenSpinTaylorT4Coeffs;
114 
115 static REAL8 OmMatch(REAL8 LNhS1, REAL8 LNhS2, REAL8 S1S1, REAL8 S1S2, REAL8 S2S2) {
116 
117  const REAL8 omM = 0.05;
118  const REAL8 omMsz12 = 9.97e-4;
119  const REAL8 omMs1d2 = -2.032e-3;
120  const REAL8 omMssq = 5.629e-3;
121  const REAL8 omMsz1d2 = 8.646e-3;
122  const REAL8 omMszsq = -5.909e-3;
123  const REAL8 omM3s1d2 = 1.801e-3;
124  const REAL8 omM3ssq = -1.4059e-2;
125  const REAL8 omM3sz1d2 = 1.5483e-2;
126  const REAL8 omM3szsq = 8.922e-3;
127 
128  return omM + /*6.05e-3 * sqrtOneMinus4Eta +*/
129  omMsz12 * (LNhS1 + LNhS2) +
130  omMs1d2 * (S1S2) +
131  omMssq * (S1S1 + S2S2) +
132  omMsz1d2 * (LNhS1 * LNhS2) +
133  omMszsq * (LNhS1 * LNhS1 + LNhS2 * LNhS2) +
134  omM3s1d2 * (LNhS1 + LNhS2) * (S1S2) +
135  omM3ssq * (LNhS1 + LNhS2) * (S1S1+S2S2) +
136  omM3sz1d2 * (LNhS1 + LNhS2) * (LNhS1*LNhS2) +
137  omM3szsq * (LNhS1 + LNhS2) * (LNhS1*LNhS1+LNhS2*LNhS2);
138 } /* End of OmMatch */
139 
140 static REAL8 fracRD(REAL8 LNhS1, REAL8 LNhS2, REAL8 S1S1, REAL8 S1S2, REAL8 S2S2) {
141 
142  const double frac0 = 0.840;
143  const double fracsz12 = -2.145e-2;
144  const double fracs1d2 = -4.421e-2;
145  const double fracssq = -2.643e-2;
146  const double fracsz1d2 = -5.876e-2;
147  const double fracszsq = -2.215e-2;
148 
149  return frac0 +
150  fracsz12 * (LNhS1 + LNhS2) +
151  fracs1d2 * (S1S2) +
152  fracssq * (S1S1 + S2S2) +
153  fracsz1d2 * (LNhS1 * LNhS2) +
154  fracszsq * (LNhS1 * LNhS1 + LNhS2 * LNhS2);
155 } /* End of fracRD */
156 
157 /**
158  * Convenience function to set up XLALSimInspiralSpinTaylotT4Coeffs struct
159  */
160 static INT4 XLALSimIMRPhenSpinParamsSetup(LALSimInspiralPhenSpinTaylorT4Coeffs *params, /**< PN params [returned] */
161  REAL8 dt, /**< Sampling in secs */
162  REAL8 fStart, /**< Starting frequency of integration*/
163  REAL8 fEnd, /**< Ending frequency of integration*/
164  REAL8 mass1, /**< Mass 1 in solar mass units */
165  REAL8 mass2, /**< Mass 2 in solar mass units */
166  REAL8 lambda1, /**< Tidal par1*/
167  REAL8 lambda2, /**< Tidal par2*/
168  REAL8 quadparam1, /**< Quad-monop par1*/
169  REAL8 quadparam2, /**< Quad-monop par2*/
170  int spinO, /**< Spin interaction */
171  int tideO, /**< Tidal iteraction interaction */
172  LALDict *LALparams, /**< Extra parameters */
173  UINT4 order /**< twice PN Order in Phase */
174 )
175 {
176  /* Zero the coefficients */
177  memset(params, 0, sizeof(LALSimInspiralPhenSpinTaylorT4Coeffs));
178  params->eta = mass1*mass2/(mass1+mass2)/(mass1+mass2);
179  params->m1ByM = mass1 / (mass1+mass2);
180  params->m2ByM = mass2 / (mass1+mass2);
181  params->M = (mass1+mass2)*LAL_MTSUN_SI;
182  REAL8 unitHz = params->M *((REAL8) LAL_PI);
183 
184  params->fEnd = fEnd*unitHz; /*On the left side there is actually an omega*/
185  params->fStart = fStart*unitHz; /*On the left side there is actually an omega*/
186  if (fEnd>0.)
187  params->dt = dt*(fEnd-fStart)/fabs(fEnd-fStart);
188  else
189  params->dt = dt;
190 
191  REAL8 phi1 = XLALSimInspiralWaveformParamsLookupNonGRPhi1(LALparams);
192 
193  params->wdotnewt = XLALSimInspiralTaylorT4wdot_0PNCoeff(params->eta);
194  params->Enewt = XLALSimInspiralPNEnergy_0PNCoeff(params->eta);
195 
196  switch (order) {
197 
198  case -1: // Use the highest PN order available.
199  case 7:
200  params->wdotcoeff[7] = XLALSimInspiralTaylorT4wdot_7PNCoeff(params->eta);
201 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
202  __attribute__ ((fallthrough));
203 #endif
204 
205  case 6:
206  params->Ecoeff[6] = XLALSimInspiralPNEnergy_6PNCoeff(params->eta);
207  params->wdotcoeff[6] = XLALSimInspiralTaylorT4wdot_6PNCoeff(params->eta);
208  params->wdotlogcoeff = XLALSimInspiralTaylorT4wdot_6PNLogCoeff(params->eta);
209  params->wdot6S1O = XLALSimInspiralTaylorT4wdot_6PNSOCoeff(params->m1ByM);
210  params->wdot6S2O = XLALSimInspiralTaylorT4wdot_6PNSOCoeff(params->m2ByM);
211 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
212  __attribute__ ((fallthrough));
213 #endif
214 
215  case 5:
216  params->Ecoeff[5] = 0.;
217  params->wdotcoeff[5] = XLALSimInspiralTaylorT4wdot_5PNCoeff(params->eta);
218  params->E5S1O = XLALSimInspiralPNEnergy_5PNSOCoeff(params->m1ByM);
219  params->E5S2O = XLALSimInspiralPNEnergy_5PNSOCoeff(params->m1ByM);
220  params->wdot5S1O = XLALSimInspiralTaylorT4wdot_5PNSOCoeff(params->m1ByM);
221  params->wdot5S2O = XLALSimInspiralTaylorT4wdot_5PNSOCoeff(params->m2ByM);
222  params->S1dot5 = XLALSimInspiralSpinDot_5PNCoeff(params->m1ByM);
223  params->S2dot5 = XLALSimInspiralSpinDot_5PNCoeff(params->m2ByM);
224 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
225  __attribute__ ((fallthrough));
226 #endif
227 
228  case 4:
229  params->wdotcoeff[4] = XLALSimInspiralTaylorT4wdot_4PNCoeff(params->eta);
230  params->Ecoeff[4] = XLALSimInspiralPNEnergy_4PNCoeff(params->eta);
231  params->wdot4S1S2 = XLALSimInspiralTaylorT4wdot_4PNS1S2CoeffAvg(params->eta);
232  params->wdot4S1OS2O = XLALSimInspiralTaylorT4wdot_4PNS1OS2OCoeffAvg(params->eta);
233  params->E4S1S2 = XLALSimInspiralPNEnergy_4PNS1S2CoeffAvg(params->eta);
234  params->E4S1OS2O = XLALSimInspiralPNEnergy_4PNS1OS2OCoeffAvg(params->eta);
235  params->wdot4S1S1 = XLALSimInspiralTaylorT4wdot_4PNS1S1CoeffAvg(params->m1ByM);
236  params->wdot4S1OS1O = XLALSimInspiralTaylorT4wdot_4PNS1OS1OCoeffAvg(params->m1ByM);
237  params->wdot4S2S2 = XLALSimInspiralTaylorT4wdot_4PNS1S1CoeffAvg(params->m2ByM);
238  params->wdot4S2OS2O = XLALSimInspiralTaylorT4wdot_4PNS1OS1OCoeffAvg(params->m2ByM);
239  params->wdot4QMS1 = quadparam1*XLALSimInspiralTaylorT4wdot_4PNQMS1S1CoeffAvg(params->m1ByM);
240  params->wdot4QMS1O = quadparam1*XLALSimInspiralTaylorT4wdot_4PNQMS1OS1OCoeffAvg(params->m1ByM);
241  params->wdot4QMS2 = quadparam2*XLALSimInspiralTaylorT4wdot_4PNQMS1S1CoeffAvg(params->m2ByM);
242  params->wdot4QMS2O = quadparam2*XLALSimInspiralTaylorT4wdot_4PNQMS1OS1OCoeffAvg(params->m2ByM);
243  params->E4QMS1 = quadparam1*XLALSimInspiralPNEnergy_4PNQMS1S1CoeffAvg(params->m1ByM);
244  params->E4QMS2 = quadparam2*XLALSimInspiralPNEnergy_4PNQMS1S1CoeffAvg(params->m2ByM);
245  params->E4QMS1O = quadparam1*XLALSimInspiralPNEnergy_4PNQMS1OS1OCoeffAvg(params->m1ByM);
246  params->E4QMS2O = quadparam2*XLALSimInspiralPNEnergy_4PNQMS1OS1OCoeffAvg(params->m2ByM);
247  params->Sdot4S2 = XLALSimInspiralSpinDot_4PNS2CoeffAvg;
248  params->Sdot4S2O = XLALSimInspiralSpinDot_4PNS2OCoeffAvg;
249  params->S1dot4QMS1O = quadparam1*XLALSimInspiralSpinDot_4PNQMSOCoeffAvg(params->m1ByM);
250  params->S2dot4QMS2O = quadparam2*XLALSimInspiralSpinDot_4PNQMSOCoeffAvg(params->m2ByM);
251 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
252  __attribute__ ((fallthrough));
253 #endif
254 
255  case 3:
256  params->Ecoeff[3] = 0.;
257  params->wdotcoeff[3] = XLALSimInspiralTaylorT4wdot_3PNCoeff(params->eta);
258  params->wdot3S1O = XLALSimInspiralTaylorT4wdot_3PNSOCoeff(params->m1ByM);
259  params->wdot3S2O = XLALSimInspiralTaylorT4wdot_3PNSOCoeff(params->m1ByM);
260  params->E3S1O = XLALSimInspiralPNEnergy_3PNSOCoeff(params->m1ByM);
261  params->E3S2O = XLALSimInspiralPNEnergy_3PNSOCoeff(params->m1ByM);
262  params->S1dot3 = XLALSimInspiralSpinDot_3PNCoeff(params->m1ByM);
263  params->S2dot3 = XLALSimInspiralSpinDot_3PNCoeff(params->m2ByM);
264 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
265  __attribute__ ((fallthrough));
266 #endif
267 
268  case 2:
269  params->Ecoeff[2] = XLALSimInspiralPNEnergy_2PNCoeff(params->eta);
270  params->wdotcoeff[2] = XLALSimInspiralTaylorT4wdot_2PNCoeff(params->eta);
271 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
272  __attribute__ ((fallthrough));
273 #endif
274 
275  case 1:
276  params->Ecoeff[1] = 0.;
277  params->wdotcoeff[1] = phi1;
278 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
279  __attribute__ ((fallthrough));
280 #endif
281 
282  case 0:
283  params->Ecoeff[0] = 1.;
284  params->wdotcoeff[0] = 1.;
285  break;
286 
287  case 8:
288  XLALPrintError("*** LALSimIMRPhenSpinInspiralRD ERROR: PhenSpin approximant not available at pseudo4PN order\n");
289  XLAL_ERROR(XLAL_EDOM);
290  break;
291 
292  default:
293  XLALPrintError("*** LALSimIMRPhenSpinInspiralRD ERROR: Impossible to create waveform with %d order\n",order);
294  XLAL_ERROR(XLAL_EFAILED);
295  break;
296  }
297 
298  switch (spinO) {
299 
300  case LAL_SIM_INSPIRAL_SPIN_ORDER_0PN:
301  case LAL_SIM_INSPIRAL_SPIN_ORDER_05PN:
302  case LAL_SIM_INSPIRAL_SPIN_ORDER_1PN:
303  /*This kills all spin effects in the phase. Still there are spin effects
304  in the waveform due to orbital plane precession*/
305  params->E6S1O = 0.;
306  params->E6S2O = 0.;
307  params->wdot6S1O = 0.;
308  params->wdot6S1O = 0.;
309  params->S1dot6 = 0.;
310  params->S2dot6 = 0.;
311 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
312  __attribute__ ((fallthrough));
313 #endif
314 
315  case LAL_SIM_INSPIRAL_SPIN_ORDER_15PN:
316  /* This keeps only the leading spin-orbit interactions*/
317  params->wdot4S1S2 = 0.;
318  params->wdot4S1OS2O = 0.;
319  params->E4QMS1 = 0.;
320  params->E4QMS2 = 0.;
321  params->E4QMS1O = 0.;
322  params->E4QMS2O = 0.;
323 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
324  __attribute__ ((fallthrough));
325 #endif
326 
327  case LAL_SIM_INSPIRAL_SPIN_ORDER_2PN:
328  /* This kills all spin interaction intervening at 2.5PN order or higher*/
329  params->E5S1O = 0.;
330  params->E5S2O = 0.;
331  params->wdot5S1O = 0.;
332  params->wdot5S2O = 0.;
333  params->S1dot5 = 0.;
334  params->S2dot5 = 0.;
335 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
336  __attribute__ ((fallthrough));
337 #endif
338 
339  case LAL_SIM_INSPIRAL_SPIN_ORDER_25PN:
340  /* This kills all spin interaction intervening at 3PN order or higher*/
341  params->wdot6S1O = 0.;
342  params->wdot6S2O = 0.;
343 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
344  __attribute__ ((fallthrough));
345 #endif
346 
347  case LAL_SIM_INSPIRAL_SPIN_ORDER_3PN:
348  case LAL_SIM_INSPIRAL_SPIN_ORDER_ALL:
349  default:
350  break;
351  }
352 
353  switch( tideO ) {
354  case LAL_SIM_INSPIRAL_TIDAL_ORDER_ALL:
355  case LAL_SIM_INSPIRAL_TIDAL_ORDER_6PN:
356  params->wdottidal6pn = lambda1 * XLALSimInspiralTaylorT4wdot_12PNTidalCoeff(params->m1ByM) + lambda2 * XLALSimInspiralTaylorT4wdot_12PNTidalCoeff(params->m2ByM);
357  params->Etidal6pn = lambda1*XLALSimInspiralPNEnergy_12PNTidalCoeff(params->m1ByM) + lambda2*XLALSimInspiralPNEnergy_12PNTidalCoeff(params->m2ByM);
358 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
359  __attribute__ ((fallthrough));
360 #endif
361  case LAL_SIM_INSPIRAL_TIDAL_ORDER_5PN:
362  params->wdottidal5pn = lambda1 * XLALSimInspiralTaylorT4wdot_10PNTidalCoeff(params->m1ByM) + lambda2 * XLALSimInspiralTaylorT4wdot_10PNTidalCoeff(params->m2ByM);
363  params->Etidal5pn = lambda1*XLALSimInspiralPNEnergy_10PNTidalCoeff(params->m1ByM) + lambda2*XLALSimInspiralPNEnergy_10PNTidalCoeff(params->m2ByM);
364 #if __GNUC__ >= 7 && !defined __INTEL_COMPILER
365  __attribute__ ((fallthrough));
366 #endif
367  case LAL_SIM_INSPIRAL_TIDAL_ORDER_0PN:
368  break;
369  default:
370  XLALPrintError("XLAL Error - %s: Invalid tidal PN order %d\n",
371  __func__, tideO );
372  XLAL_ERROR(XLAL_EINVAL);
373  break;
374  }
375 
376  return XLAL_SUCCESS;
377 } /* End of XLALSimIMRPhenSpinParamsSetup */
378 
379 static INT4 XLALSpinInspiralDerivatives(UNUSED double t,
380  const double values[],
381  double dvalues[],
382  void *mparams)
383 {
384  REAL8 omega; // time-derivative of the orbital phase
385  REAL8 LNhx, LNhy, LNhz; // orbital angular momentum unit vector
386  REAL8 S1x, S1y, S1z; // dimension-less spin variable S/M^2
387  REAL8 S2x, S2y, S2z;
388  REAL8 LNhS1, LNhS2; // scalar products
389  REAL8 domega; // derivative of omega
390  REAL8 dLNhx, dLNhy, dLNhz; // derivatives of \f$\hat L_N\f$ components
391  REAL8 dS1x, dS1y, dS1z; // derivative of \f$S_i\f$
392  REAL8 dS2x, dS2y, dS2z;
393  REAL8 energy,energyold;
394 
395  /* auxiliary variables*/
396  REAL8 S1S2, S1S1, S2S2; // Scalar products
397  REAL8 alphadotcosi; // alpha is the right ascension of L, i(iota) the angle between L and J
398  REAL8 v, v2, v4, v5, v6, v7;
399  REAL8 tmpx, tmpy, tmpz, cross1x, cross1y, cross1z, cross2x, cross2y, cross2z, LNhxy;
400 
401  LALSimInspiralPhenSpinTaylorT4Coeffs *params = (LALSimInspiralPhenSpinTaylorT4Coeffs *) mparams;
402 
403  /* --- computation start here --- */
404  omega = values[1];
405 
406  LNhx = values[2];
407  LNhy = values[3];
408  LNhz = values[4];
409 
410  S1x = values[5];
411  S1y = values[6];
412  S1z = values[7];
413 
414  S2x = values[8];
415  S2y = values[9];
416  S2z = values[10];
417 
418  energyold = values[11];
419 
420  v = cbrt(omega);
421  v2 = v * v;
422  v4 = omega * v;
423  v5 = omega * v2;
424  v6 = omega * omega;
425  v7 = omega * omega * v;
426 
427  // Omega derivative without spin effects up to 3.5 PN
428  // this formula does not include the 1.5PN shift mentioned in arXiv:0810.5336, five lines below (3.11)
429  domega = params->wdotcoeff[0]
430  + v * (params->wdotcoeff[1]
431  + v * (params->wdotcoeff[2]
432  + v * (params->wdotcoeff[3]
433  + v * (params->wdotcoeff[4]
434  + v * (params->wdotcoeff[5]
435  + v * (params->wdotcoeff[6] + params->wdotlogcoeff * log(v)
436  + v * params->wdotcoeff[7]))))));
437 
438  energy = (params->Ecoeff[0] + v2 * (params->Ecoeff[2] +
439  v2 * (params->Ecoeff[4] +
440  v2 * params->Ecoeff[6])));
441 
442  domega+= omega*v7* ( params->wdottidal5pn + v2 * ( params->wdottidal6pn ) );
443  energy+= omega*v7* ( params->Etidal5pn + v2 * params->Etidal6pn);
444 
445  // Adding spin effects
446  // L dot S1,2
447  LNhS1 = (LNhx * S1x + LNhy * S1y + LNhz * S1z);
448  LNhS2 = (LNhx * S2x + LNhy * S2y + LNhz * S2z);
449 
450  // wdotSO15si = -1/12 (...)
451  domega += omega * (params->wdot3S1O * LNhS1 + params->wdot3S2O * LNhS2); // see e.g. Buonanno et al. gr-qc/0211087
452 
453  energy += omega * (params->E3S1O * LNhS1 + params->E3S2O * LNhS2); // see e.g. Blanchet et al. gr-qc/0605140
454 
455  // wdotSS2 = -1/48 eta ...
456  S1S1 = (S1x * S1x + S1y * S1y + S1z * S1z);
457  S2S2 = (S2x * S2x + S2y * S2y + S2z * S2z);
458  S1S2 = (S1x * S2x + S1y * S2y + S1z * S2z);
459  domega += v4 * ( params->wdot4S1S2 * S1S2 + params->wdot4S1OS2O * LNhS1 * LNhS2); // see e.g. Buonanno et al. arXiv:0810.5336
460  domega += v4 * ( params->wdot4QMS1O * LNhS1*LNhS1 + params->wdot4QMS2O * LNhS2*LNhS2 + params->wdot4QMS1 * S1S1 + params->wdot4QMS2 * S2S2 );
461  // see Racine et al. arXiv:0812.4413
462 
463  energy += v4 * (params->E4S1S2 * S1S2 + params->E4S1OS2O * LNhS1 * LNhS2); // see e.g. Buonanno et al. as above
464  energy += v4 * (params->E4QMS1 * S1S1 + params->E4QMS2 * S2S2 + params->E4QMS1O * LNhS1 * LNhS1 + params->E4QMS2O * LNhS2 * LNhS2); // see Racine et al. as above
465 
466  // wdotspin25SiLNh = see below
467  domega += v5 * (params->wdot5S1O * LNhS1 + params->wdot5S2O * LNhS2); //see (8.3) of Blanchet et al.
468  energy += v5 * (params->E5S1O * LNhS1 + params->E5S2O * LNhS2); //see (7.9) of Blanchet et al.
469 
470  domega += omega*omega * (params->wdot6S1O * LNhS1 + params->wdot6S2O * LNhS2); // see (6.5) of arXiv:1104.5659
471 
472  // Setting the right pre-factor
473  domega *= params->wdotnewt * v5 * v6;
474  energy *= params->Enewt * v2;
475 
476  /*Derivative of the angular momentum and spins */
477 
478  /* dS1, 1.5PN */
479  cross1x = (LNhy * S1z - LNhz * S1y);
480  cross1y = (LNhz * S1x - LNhx * S1z);
481  cross1z = (LNhx * S1y - LNhy * S1x);
482 
483  dS1x = params->S1dot3 * v5 * cross1x;
484  dS1y = params->S1dot3 * v5 * cross1y;
485  dS1z = params->S1dot3 * v5 * cross1z;
486 
487  /* dS1, 2PN */
488  tmpx = S1z * S2y - S1y * S2z;
489  tmpy = S1x * S2z - S1z * S2x;
490  tmpz = S1y * S2x - S1x * S2y;
491 
492  // S1S2 contribution see. eq. 2.23 of arXiv:0812.4413
493  dS1x += v6 * (params->Sdot4S2*tmpx + params->Sdot4S2O * LNhS2 * cross1x);
494  dS1y += v6 * (params->Sdot4S2*tmpy + params->Sdot4S2O * LNhS2 * cross1y);
495  dS1z += v6 * (params->Sdot4S2*tmpz + params->Sdot4S2O * LNhS2 * cross1z);
496  // S1S1 contribution
497  dS1x += v6 * LNhS1 * cross1x * params->S1dot4QMS1O;
498  dS1y += v6 * LNhS1 * cross1y * params->S1dot4QMS1O;
499  dS1z += v6 * LNhS1 * cross1z * params->S1dot4QMS1O;
500 
501  // dS1, 2.5PN, eq. 7.8 of Blanchet et al. gr-qc/0605140
502  dS1x += params->S1dot5 * v7 * cross1x;
503  dS1y += params->S1dot5 * v7 * cross1y;
504  dS1z += params->S1dot5 * v7 * cross1z;
505 
506  /* dS2, 1.5PN */
507  cross2x = (LNhy * S2z - LNhz * S2y);
508  cross2y = (LNhz * S2x - LNhx * S2z);
509  cross2z = (LNhx * S2y - LNhy * S2x);
510 
511  dS2x = params->S2dot3 * v5 * cross2x;
512  dS2y = params->S2dot3 * v5 * cross2y;
513  dS2z = params->S2dot3 * v5 * cross2z;
514 
515  /* dS2, 2PN */
516  dS2x += v6 * (-params->Sdot4S2*tmpx + params->Sdot4S2O * LNhS1 * cross2x);
517  dS2y += v6 * (-params->Sdot4S2*tmpy + params->Sdot4S2O * LNhS1 * cross2y);
518  dS2z += v6 * (-params->Sdot4S2*tmpz + params->Sdot4S2O * LNhS1 * cross2z);
519  // S2S2 contribution
520  dS2x += v6 * LNhS2 * cross2x * params->S2dot4QMS2O;
521  dS2y += v6 * LNhS2 * cross2y * params->S2dot4QMS2O;
522  dS2z += v6 * LNhS2 * cross2z * params->S2dot4QMS2O;
523 
524  // dS2, 2.5PN, eq. 7.8 of Blanchet et al. gr-qc/0605140
525  dS2x += params->S2dot5 * v7 * cross2x;
526  dS2y += params->S2dot5 * v7 * cross2y;
527  dS2z += params->S2dot5 * v7 * cross2z;
528 
529  dLNhx = -(dS1x + dS2x) * v / params->eta;
530  dLNhy = -(dS1y + dS2y) * v / params->eta;
531  dLNhz = -(dS1z + dS2z) * v / params->eta;
532 
533  /* dphi */
534  LNhxy = LNhx * LNhx + LNhy * LNhy;
535 
536  if (LNhxy > 0.0)
537  alphadotcosi = LNhz * (LNhx * dLNhy - LNhy * dLNhx) / LNhxy;
538  else
539  {
540  //XLALPrintWarning("*** LALSimIMRPSpinInspiralRD WARNING ***: alphadot set to 0, LNh:(%12.4e %12.4e %12.4e)\n",LNhx,LNhy,LNhz);
541  alphadotcosi = 0.;
542  }
543 
544  /* dvalues->data[0] is the phase derivative */
545  /* omega is the derivative of the orbital phase omega \neq dvalues->data[0] */
546  dvalues[0] = omega - alphadotcosi;
547  dvalues[1] = domega;
548 
549  dvalues[2] = dLNhx;
550  dvalues[3] = dLNhy;
551  dvalues[4] = dLNhz;
552 
553  dvalues[5] = dS1x;
554  dvalues[6] = dS1y;
555  dvalues[7] = dS1z;
556 
557  dvalues[8] = dS2x;
558  dvalues[9] = dS2y;
559  dvalues[10] = dS2z;
560 
561  dvalues[11] = (energy-energyold)/params->dt*params->M;
562 
563  return GSL_SUCCESS;
564 } /* end of XLALSpinInspiralDerivatives */
565 
566 static INT4 XLALGenerateWaveDerivative (REAL8Vector *dwave,
567  REAL8Vector *wave,
568  REAL8 dt
569 )
570 {
571  /* XLAL error handling */
572  INT4 errcode = XLAL_SUCCESS;
573 
574  /* For checking GSL return codes */
575  INT4 gslStatus;
576 
577  UINT4 j;
578  double *x, *y;
579  double dy;
580  gsl_interp_accel *acc;
581  gsl_spline *spline;
582 
583  if (wave->length!=dwave->length)
584  XLAL_ERROR( XLAL_EFUNC );
585 
586  /* Getting interpolation and derivatives of the waveform using gsl spline routine */
587  /* Initialize arrays and supporting variables for gsl */
588 
589  x = (double *) LALMalloc(wave->length * sizeof(double));
590  y = (double *) LALMalloc(wave->length * sizeof(double));
591 
592  if ( !x || !y )
593  {
594  if ( x ) LALFree (x);
595  if ( y ) LALFree (y);
596  XLAL_ERROR( XLAL_ENOMEM );
597  }
598 
599  for (j = 0; j < wave->length; ++j)
600  {
601  x[j] = j;
602  y[j] = wave->data[j];
603  }
604 
605  XLAL_CALLGSL( acc = (gsl_interp_accel*) gsl_interp_accel_alloc() );
606  XLAL_CALLGSL( spline = (gsl_spline*) gsl_spline_alloc(gsl_interp_cspline, wave->length) );
607  if ( !acc || !spline )
608  {
609  if ( acc ) gsl_interp_accel_free(acc);
610  if ( spline ) gsl_spline_free(spline);
611  LALFree( x );
612  LALFree( y );
613  XLAL_ERROR( XLAL_ENOMEM );
614  }
615 
616  /* Gall gsl spline interpolation */
617  XLAL_CALLGSL( gslStatus = gsl_spline_init(spline, x, y, wave->length) );
618  if ( gslStatus != GSL_SUCCESS )
619  {
620  gsl_spline_free(spline);
621  gsl_interp_accel_free(acc);
622  LALFree( x );
623  LALFree( y );
624  XLAL_ERROR( XLAL_EFUNC );
625  }
626 
627  /* Getting first and second order time derivatives from gsl interpolations */
628  for (j = 0; j < wave->length; ++j)
629  {
630  XLAL_CALLGSL(gslStatus = gsl_spline_eval_deriv_e( spline, j, acc, &dy ) );
631  if (gslStatus != GSL_SUCCESS )
632  {
633  gsl_spline_free(spline);
634  gsl_interp_accel_free(acc);
635  LALFree( x );
636  LALFree( y );
637  XLAL_ERROR( XLAL_EFUNC );
638  }
639  dwave->data[j] = (REAL8)(dy / dt);
640  }
641 
642  /* Free gsl variables */
643  gsl_spline_free(spline);
644  gsl_interp_accel_free(acc);
645  LALFree(x);
646  LALFree(y);
647 
648  return errcode;
649 } /* End of XLALGenerateWaveDerivative */
650 
651 static INT4 XLALSimSpinInspiralTest(UNUSED double t, const double values[], double dvalues[], void *mparams) {
652 
653  LALSimInspiralPhenSpinTaylorT4Coeffs *params = (LALSimInspiralPhenSpinTaylorT4Coeffs *) mparams;
654 
655  REAL8 omega = values[1];
656  REAL8 energy = values[11];
657  REAL8 denergy = dvalues[11];
658 
659  if ( (energy > 0.0) || (( denergy*params->dt/params->M > - 0.001*energy )&&(energy<0.) ) ) {
660  if (energy>0.) XLALPrintWarning("*** Test: LALSimIMRPSpinInspiral WARNING **: Bounding energy >ve!\n");
661  else
662  XLALPrintWarning("*** Test: LALSimIMRPSpinInspiral WARNING **: Energy increases dE %12.6e dE*dt %12.6e 1pMEn %12.4e M: %12.4e, eta: %12.4e om %12.6e \n", denergy, denergy*params->dt/params->M, - 0.001*energy, params->M/LAL_MTSUN_SI, params->eta, omega);
663  return LALSIMINSPIRAL_PST4_TEST_ENERGY;
664  }
665  else if (omega < 0.0) {
666  XLALPrintWarning("** LALSimIMRPSpinInspiral WARNING **: omega < 0 M: %12.4e, eta: %12.4e om %12.6e\n",params->M, params->eta, omega);
667  return LALSIMINSPIRAL_PST4_DERIVATIVE_OMEGANONPOS;
668  }
669  else if (dvalues[1] < 0.0) {
670  /* omegadot < 0 */
671  return LALSIMINSPIRAL_PST4_TEST_OMEGADOT;
672  }
673  else if (isnan(omega)) {
674  /* omega is nan */
675  return LALSIMINSPIRAL_PST4_TEST_OMEGANAN;
676  }
677  else if ( params->fEnd > 0. && params->fStart > params->fEnd && omega < params->fEnd) {
678  /* freq. below bound in backward integration */
679  return LALSIMINSPIRAL_PST4_TEST_FREQBOUND;
680  }
681  else if ( params->fEnd > params->fStart && omega > params->fEnd) {
682  /* freq. above bound in forward integration */
683  return LALSIMINSPIRAL_PST4_TEST_FREQBOUND;
684  }
685  else
686  return GSL_SUCCESS;
687 } /* End of XLALSimSpinInspiralTest */
688 
689 
690 static INT4 XLALSimIMRPhenSpinTest(UNUSED double t, const double values[], double dvalues[], void *mparams) {
691 
692  LALSimInspiralPhenSpinTaylorT4Coeffs *params = (LALSimInspiralPhenSpinTaylorT4Coeffs *) mparams;
693 
694  REAL8 omega = values[1];
695  REAL8 energy = values[11];
696  REAL8 denergy = dvalues[11];
697 
698  REAL8 LNhS1=(values[2]*values[5]+values[3]*values[6]+values[4]*values[7])/params->m1ByM/params->m1ByM;
699  REAL8 LNhS2=(values[2]*values[8]+values[3]*values[9]+values[4]*values[10])/params->m2ByM/params->m2ByM;
700  REAL8 S1sq =(values[5]*values[5]+values[6]*values[6]+values[7]*values[7])/pow(params->m1ByM,4);
701  REAL8 S2sq =(values[8]*values[8]+values[9]*values[9]+values[10]*values[10])/pow(params->m2ByM,4);
702  REAL8 S1S2 =(values[5]*values[8]+values[6]*values[9]+values[7]*values[10])/pow(params->m1ByM*params->m2ByM,2);
703 
704  REAL8 omegaMatch=OmMatch(LNhS1,LNhS2,S1sq,S1S2,S2sq)+0.0005;
705 
706  if ( (energy > 0.0) || (( denergy*params->dt/params->M > - 0.001*energy )&&(energy<0.) ) ) {
707  if (energy>0.) XLALPrintWarning("*** Test: LALSimIMRPSpinInspiralRD WARNING **: Bounding energy >ve!\n");
708  else
709  XLALPrintWarning("*** Test: LALSimIMRPSpinInspiralRD WARNING **: Energy increases dE %12.6e dE*dt %12.6e 1pMEn %12.4e M: %12.4e, eta: %12.4e om %12.6e \n", denergy, denergy*params->dt/params->M, - 0.001*energy, params->M/LAL_MTSUN_SI, params->eta, omega);
710  return LALSIMINSPIRAL_PST4_TEST_ENERGY;
711  }
712  else if (omega < 0.0) {
713  XLALPrintWarning("** LALSimIMRPSpinInspiralRD WARNING **: omega < 0 M: %12.4e, eta: %12.4e om %12.6e\n",params->M, params->eta, omega);
714  return LALSIMINSPIRAL_PST4_DERIVATIVE_OMEGANONPOS;
715  }
716  else if (dvalues[1] < 0.0) {
717  /* omegadot < 0 */
718  return LALSIMINSPIRAL_PST4_TEST_OMEGADOT;
719  }
720  else if (isnan(omega)) {
721  /* omega is nan */
722  return LALSIMINSPIRAL_PST4_TEST_OMEGANAN;
723  }
724  else if ( params->fEnd > 0. && params->fStart > params->fEnd && omega < params->fEnd) {
725  /* freq. below bound in backward integration */
726  return LALSIMINSPIRAL_PST4_TEST_FREQBOUND;
727  }
728  else if ( params->fEnd > params->fStart && omega > params->fEnd) {
729  /* freq. above bound in forward integration */
730  return LALSIMINSPIRAL_PST4_TEST_FREQBOUND;
731  }
732  else if (omega>omegaMatch) {
733  return LALSIMINSPIRAL_PST4_TEST_OMEGAMATCH;
734  }
735  else
736  return GSL_SUCCESS;
737 } /* End of XLALSimIMRPhenSpinTest */
738 
739 typedef struct tagLALSimInspiralInclAngle {
740  REAL8 cHi;
741  REAL8 sHi;
742  REAL8 ci;
743  REAL8 si;
744  REAL8 ci2;
745  REAL8 si2;
746  REAL8 cHi2;
747  REAL8 sHi2;
748  REAL8 cHi3;
749  REAL8 sHi3;
750  REAL8 cHi4;
751  REAL8 sHi4;
752  REAL8 cHi5;
753  REAL8 sHi5;
754  REAL8 cHi6;
755  REAL8 sHi6;
756  REAL8 cHi8;
757  REAL8 sHi8;
758  REAL8 cDi;
759  REAL8 sDi;
760 } LALSimInspiralInclAngle;
761 
762 static INT4 XLALSimSpinInspiralFillL2Modes(COMPLEX16Vector *hL2,
763  REAL8 v,
764  REAL8 eta,
765  REAL8 dm,
766  REAL8 Psi,
767  REAL8 alpha,
768  LALSimInspiralInclAngle *an
769  )
770 {
771  const INT4 os=2;
772  REAL8 amp20 = sqrt(1.5);
773  REAL8 v2 = v*v;
774  REAL8 damp = 1.;
775 
776  hL2->data[2+os] = ( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
777  ( cos(2.*(Psi+alpha)) * an->cHi4 + cos(2.*(Psi-alpha)) * an->sHi4 ) +
778  v * dm/3.*an->si * ( cos(Psi-2.*alpha) * an->sHi2 + cos(Psi + 2.*alpha) * an->cHi2 ) );
779 
780  hL2->data[2+os]+=I*( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
781  (-sin(2.*(Psi+alpha)) * an->cHi4 + sin(2.*(Psi-alpha)) * an->sHi4 ) +
782  v * dm/3.*an->si * ( sin(Psi-2.*alpha) * an->sHi2 - sin(Psi + 2. * alpha) * an->cHi2 ) );
783 
784  hL2->data[-2+os] = ( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
785  ( cos(2. * (Psi + alpha)) * an->cHi4 + cos(2. * (Psi - alpha)) * an->sHi4 ) -
786  v * dm / 3. * an->si * ( cos(Psi - 2. * alpha) * an->sHi2 + cos(Psi + 2. * alpha) * an->cHi2 ) );
787 
788  hL2->data[-2+os]+=I*( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
789  ( sin(2.*(Psi + alpha))*an->cHi4 - sin(2.*(Psi-alpha)) * an->sHi4 ) +
790  v*dm/3.*an->si * ( sin(Psi-2.*alpha) * an->sHi2 - sin(Psi+2.*alpha) * an->cHi2 ) );
791 
792  hL2->data[1+os] = an->si * ( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
793  ( -cos(2. * Psi - alpha) * an->sHi2 + cos(2. * Psi + alpha) * an->cHi2 ) +
794  v * dm / 3. * ( -cos(Psi + alpha) * (an->ci + an->cDi)/2. - cos(Psi - alpha) * an->sHi2 * (1. + 2. * an->ci) ) );
795 
796  hL2->data[1+os]+= an->si *I*( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
797  ( -sin(2.*Psi-alpha ) * an->sHi2 - sin(2.*Psi + alpha) * an->cHi2 ) +
798  v * dm / 3. * (sin(Psi + alpha) * (an->ci + an->cDi)/2. - sin(Psi - alpha) * an->sHi2 * (1.+2.*an->ci) ) );
799 
800  hL2->data[-1+os] = an->si * ( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
801  ( cos(2.*Psi-alpha) * an->sHi2 - cos(2.*Psi+alpha)*an->cHi2) +
802  v * dm / 3. * ( -cos(Psi + alpha) * (an->ci + an->cDi)/2. - cos(Psi - alpha) * an->sHi2 * (1. + 2. * an->ci) ) );
803 
804  hL2->data[-1+os]+= an->si *I*( 1./( 1. + damp * v2 / 42. * (107. - 55. * eta) ) *
805  ( -sin(2. * Psi - alpha) * an->sHi2 - sin(2. * Psi + alpha) * an->cHi2 ) -
806  v * dm / 3. * ( sin(Psi + alpha) * (an->ci + an->cDi)/2. - sin(Psi - alpha) * an->sHi2 * (1. + 2. * an->ci) ) );
807 
808  hL2->data[os] = amp20 * ( an->si2/( 1. + damp *v2/42.*(107.-55.*eta) )*cos(2.*Psi) + I*v*dm/3.*an->sDi*sin(Psi) );
809 
810  return XLAL_SUCCESS;
811 } /* End of XLALSimSpinInspiralFillL2Modes*/
812 
813 static INT4 XLALSimSpinInspiralFillL3Modes(COMPLEX16Vector *hL3,
814  REAL8 v,
815  REAL8 eta,
816  REAL8 dm,
817  REAL8 Psi,
818  REAL8 alpha,
819  LALSimInspiralInclAngle *an)
820 {
821  const INT4 os=3;
822  REAL8 amp32 = sqrt(1.5);
823  REAL8 amp31 = sqrt(0.15);
824  REAL8 amp30 = 1. / sqrt(5)/2.;
825  REAL8 v2 = v*v;
826 
827  hL3->data[3+os] = (v * dm * (-9.*cos(3.*(Psi-alpha))*an->sHi6 - cos(Psi-3.*alpha)*an->sHi4*an->cHi2 + cos(Psi+3.*alpha)*an->sHi2*an->cHi4 + 9.*cos(3.*(Psi+alpha))*an->cHi6) +
828  v2 * 4. * an->si *(1.-3.*eta)* ( -cos(2.*Psi-3.*alpha)*an->sHi4 + cos(2.*Psi+3.*alpha)*an->cHi4) );
829 
830  hL3->data[3+os]+= I*(v * dm * (-9.*sin(3.*(Psi-alpha))*an->sHi6 - sin(Psi-3.*alpha)*an->sHi4*an->cHi2 - sin(Psi+3.*alpha)*an->sHi2*an->cHi4 - 9.*sin(3.*(Psi+alpha))* an->cHi6) +
831  v2 * 4. * an->si *(1.-3.*eta)* ( -sin(2.*Psi-3.*alpha)*an->sHi4 -sin(2.*Psi+3.*alpha)*an->cHi4) );
832 
833  hL3->data[-3+os] = (-v * dm * (-9.*cos(3.*(Psi-alpha))*an->sHi6 - cos(Psi-3.*alpha)*an->sHi4*an->cHi2 + cos(Psi+3.*alpha)*an->sHi2*an->cHi4 + 9.*cos(3.*(Psi+alpha))*an->cHi6) +
834  v2 * 4. * an->si *(1.-3.*eta)*( -cos(2.*Psi-3.*alpha)*an->sHi4 + cos(2.*Psi+3.*alpha)*an->cHi4) );
835 
836  hL3->data[-3+os]+=I*(v * dm *(-9.*sin(3.*(Psi-alpha))*an->sHi6 - sin(Psi-3.*alpha)*an->sHi4*an->cHi2 - sin(Psi+3.*alpha)*an->sHi2*an->cHi4 - 9.*sin(3.*(Psi+alpha))* an->cHi6) -
837  v2 * 4. * an->si * (1.-3.*eta)*( -sin(2.*Psi-3.*alpha)*an->sHi4 - sin(2.*Psi+3.*alpha)*an->cHi4 ) );
838 
839  hL3->data[2+os] = amp32 * ( v * dm/3. * (27.*cos(3.*Psi-2.*alpha)*an->si*an->sHi4 + 27.*cos(3.*Psi+2.*alpha)*an->si*an->cHi4 + cos(Psi+2.*alpha)*an->cHi3*(5.*an->sHi-3.*an->si*an->cHi-3.*an->ci*an->sHi) /2. + cos(Psi-2.*alpha)*an->sHi3*(5.*an->cHi+3.*an->ci*an->cHi-3.*an->si*an->sHi) /2. ) +
840  v2*(1./3.-eta) * (-8.*an->cHi4*(3.*an->ci-2.)*cos(2.*(Psi+alpha)) + 8.*an->sHi4*(3.*an->ci+2.)*cos(2.*(Psi-alpha)) ) );
841 
842  hL3->data[2+os]+= amp32*I*( v * dm/3. * ( 27.*sin(3.*Psi-2.*alpha)*an->si*an->sHi4 - 27.*cos(3.*Psi+2.*alpha)*an->si*an->cHi4 - sin(Psi+2.*alpha)*an->cHi3*(5.*an->sHi-3.*an->si*an->cHi-3.*an->ci*an->sHi) /2. + sin(Psi-2.*alpha)*an->sHi3*(5.*an->cHi+3.*an->ci*an->cHi-3.*an->si*an->sHi)/2. ) +
843  v2*(1./3.-eta) * ( 8.*an->cHi4*(3.*an->ci-2.)*sin(2.*(Psi+alpha)) + 8.*an->sHi4*(3.*an->ci+2.)*sin(2.*(Psi-alpha)) ) );
844 
845  hL3->data[-2+os] = amp32 * ( v * dm/3. * (27.*cos(3.*Psi-2.*alpha)*an->si*an->sHi4 + 27.*cos(3.*Psi+2.*alpha)*an->si*an->cHi4 + cos(Psi+2.*alpha)*an->cHi3*(5.*an->sHi-3.*an->si*an->cHi-3.*an->ci*an->sHi) /2. + cos(Psi-2.*alpha)*an->sHi3*(5.*an->cHi+3.*an->ci*an->cHi-3.*an->si*an->sHi) /2. ) -
846  v2*(1./3.-eta) * ( 8.*an->cHi4*(3.*an->ci-2.)*cos(2.*(Psi+alpha)) - 8.*an->sHi4*(3.*an->ci+2.)*cos(2.*(Psi-alpha)) ) );
847 
848  hL3->data[-2+os]+= amp32*I*(-v * dm/3. * (27.*sin(3.*Psi-2.*alpha)*an->si*an->sHi4 - 27.*cos(3.*Psi+2.*alpha)*an->si*an->cHi4 - sin(Psi+2.*alpha)*an->cHi3*(5.*an->sHi-3.*an->si*an->cHi-3.*an->ci*an->sHi) /2.+ sin(Psi-2.*alpha)*an->sHi3*(5.*an->cHi+3.*an->ci*an->cHi-3.*an->si*an->sHi) /2.) +
849  v2*(1./3.-eta) * (8.*an->cHi4*(3.*an->ci-2.)*sin(2.*(Psi+alpha)) + 8.*an->sHi4*(3.*an->ci+2.)*sin(2.*(Psi-alpha)) ) );
850 
851  hL3->data[1+os] = amp31 * ( v * dm/6. * ( -135.*cos(3.*Psi-alpha)*an->sHi*an->sHi2 + 135.*cos(3.*Psi+alpha)*an->sHi*an->cHi2 + cos(Psi+alpha)*an->cHi2*(15.*an->cDi-20.*an->ci+13.)/2. - cos(Psi-alpha)*an->sHi2*(15.*an->cDi+20.*an->ci+13.)/2.)
852  + v2*(1./3.-eta) * ( 20.*an->cHi3*cos(2.*Psi+alpha)*(3.*(an->sHi*an->ci+an->cHi*an->si)-5.*an->sHi) + 20.*an->sHi3*cos(2.*Psi-alpha)*(3.*(an->cHi2*an->ci-an->sHi*an->si)+5.*an->cHi) ) );
853 
854  hL3->data[1+os]+= amp31*I*(-v * dm/6. * ( -135.*cos(3.*Psi-alpha)*an->si2*an->sHi2 + 135.*cos(3.*Psi+alpha)*an->si2*an->cHi2 + cos(Psi+alpha)*an->cHi2*(15.*an->cDi-20.*an->ci+13.)/2. - cos(Psi-alpha)*an->sHi2*(15.*an->cDi+20.*an->ci+13.)/2. )
855  - v2*(1./3.-eta) * ( 20.*an->cHi3*cos(2.*Psi+alpha)*(3.*(an->sHi*an->ci+an->cHi*an->si)-5.*an->sHi) + 20.*an->sHi3*cos(2.*Psi-alpha)*(3.*(an->cHi*an->ci-an->sHi*an->si)+5.*an->cHi) ) );
856 
857  hL3->data[-1+os] = amp31 * (-v * dm/6. * ( -135.*cos(3.*Psi-alpha)*an->si2*an->sHi2 + 135.*cos(3.*Psi+alpha)*an->si2*an->cHi2 + cos(Psi+alpha)*an->cHi2*(15.*an->cDi-20.*an->ci+13.)/2.- cos(Psi-alpha) * an->sHi2*(15.*an->cDi+20.*an->ci+13.)/2. ) -
858  v2 * (1./3.-eta)* ( 20.*an->cHi3*cos(2.*Psi+alpha)*(3.*(an->sHi*an->ci+an->cHi*an->si)-5.*an->sHi) + 20.*an->sHi3*cos(2.*Psi-alpha)*(3.*(an->cHi*an->ci-an->sHi*an->si)+5.*an->cHi) ) );
859 
860  hL3->data[-1+os]+= amp31*I*(v * dm/6. * ( -135.*sin(3.*Psi-alpha)*an->si2*an->sHi2 - 135.*sin(3.*Psi+alpha)*an->si2*an->cHi2 - sin(Psi+alpha)*an->cHi2*(15.*an->cDi-20.*an->ci+13.)/2. - sin(Psi-alpha)*an->sHi2*(15.*an->cDi+20.*an->ci+13.)/2.)
861  -v2 * (1./3.-eta)* ( 20.*an->cHi3*sin(2.*Psi+alpha)*(3.*(an->sHi*an->ci+an->ci2*an->si)-5.*an->si2) - 20.*an->sHi3*sin(2.*Psi-alpha)*(3.*(an->ci2*an->ci-an->si2*an->si)+5.*an->ci2) ) );
862 
863  hL3->data[os] = amp30 * I * ( v * dm * ( cos(Psi)*an->si*(cos(2.*Psi)*(45.*an->si2)-(25.*an->cDi-21.) ) ) +
864  v2*(1.-3.*eta) * (80.*an->si2*an->cHi*sin(2.*Psi) ) );
865 
866  return XLAL_SUCCESS;
867 
868 } /*End of XLALSimSpinInspiralFillL3Modes*/
869 
870 static INT4 XLALSimSpinInspiralFillL4Modes(COMPLEX16Vector *hL4,
871  UNUSED REAL8 v,
872  REAL8 eta,
873  UNUSED REAL8 dm,
874  REAL8 Psi,
875  REAL8 alpha,
876  LALSimInspiralInclAngle *an
877  )
878 {
879  const INT4 os=4;
880  REAL8 amp43 = - sqrt(2.);
881  REAL8 amp42 = sqrt(7.)/2.;
882  REAL8 amp41 = sqrt(3.5)/4.;
883  REAL8 amp40 = sqrt(17.5)/16.;
884 
885  hL4->data[4+os] = (1. - 3.*eta) * ( 4.*an->sHi8*cos(4.*(Psi-alpha)) + cos(2.*Psi-4.*alpha)*an->sHi6*an->cHi2 + an->sHi2*an->cHi6*cos(2.*Psi+4.*alpha) + 4.*an->cHi8*cos(4.*(Psi+alpha)) );
886 
887  hL4->data[4+os]+= (1. - 3.*eta)*I*( 4.*an->sHi8*sin(4.*(Psi-alpha)) + sin(2.*Psi-4.*alpha)*an->sHi6*an->cHi2 - an->sHi2*an->cHi6*sin(2.*Psi+4.*alpha) - 4.*an->cHi8*sin(4.*(Psi+alpha)) );
888 
889  hL4->data[-4+os] = (1. - 3.*eta) * (4.*an->sHi8*cos(4.*(Psi-alpha)) + cos(2.*Psi-4.*alpha)*an->sHi6*an->cHi2 + an->sHi2*an->cHi6*cos(2.*Psi+4.*alpha) + 4.*an->cHi8*cos(4.*(Psi+alpha) ) );
890 
891  hL4->data[-4+os]+=-(1. - 3.*eta) *I*(4.*an->sHi8*sin(4.*(Psi-alpha)) + sin(2*Psi-4.*alpha)*an->sHi6*an->cHi2 - an->sHi2*an->cHi6*sin(2.*Psi+4.*alpha) - 4.*an->cHi8*sin(4.*(Psi+alpha)) );
892 
893  hL4->data[3+os] = amp43 * (1. - 3.*eta) * an->si * ( 4.*an->sHi6*cos(4.*Psi-3.*alpha) - 4.*an->cHi6*cos(4.*Psi+3.*alpha) - an->sHi4*(an->ci+0.5)/2.*cos(2.*Psi-3.*alpha) + an->cHi4*(an->ci-0.5)*cos(2.*Psi+3.*alpha) ); /****/
894 
895  hL4->data[3+os]+= amp43*I*(1. - 3.*eta) * an->si * ( 4.*an->sHi6*sin(4.*Psi-3.*alpha) + 4.*an->cHi6*sin(4.*Psi+3.*alpha) - an->sHi4*(an->ci+0.5)/2.*sin(2.*Psi-3.*alpha) + an->cHi4*(an->ci-0.5)*sin(2.*Psi+3.*alpha) ); /****/
896 
897  hL4->data[-3+os] = -amp43 * (1. - 3.*eta) * an->si * ( 4.*an->sHi6*cos(4.*Psi-3.*alpha) - 4.*an->cHi6*cos(4.*Psi+3.*alpha) - an->sHi4*(an->ci+0.5)/2.*cos(2.*Psi-3.*alpha) + an->cHi4*(an->ci-0.5)*cos(2.*Psi+3.*alpha) ); /****/
898 
899  hL4->data[-3+os]+= amp43*I*(1. - 3.*eta) * an->si * ( 4.*an->sHi6*sin(4.*Psi-3.*alpha) + 4.*an->cHi6*sin(4.*Psi+3.*alpha) - an->sHi4*(an->ci+0.5)/2.*sin(2.*Psi-3.*alpha) + an->cHi4*(an->ci-0.5)*sin(2.*Psi+3.*alpha) ); /****/
900 
901  hL4->data[2+os] = amp42 * (1. - 3.*eta) * ( 16.*an->sHi6*an->cHi2*cos(4.*Psi-2.*alpha) + 16.*an->cHi6*an->sHi2*cos(4.*Psi+2.*alpha) - an->cHi4*cos(2.*(Psi+alpha))*(an->cDi-2.*an->ci+9./7.)/2. - an->sHi4*cos(2.*(Psi-alpha))*(an->cDi+2.*an->ci+9./7.)/2. );
902 
903  hL4->data[2+os]+= amp42 *I*(1. - 3.*eta) * ( 16.*an->sHi6*an->cHi2 * sin(4.*Psi-2.*alpha) - 16.*an->cHi6*an->sHi2*sin(4.*Psi+2.*alpha) + an->cHi4*sin(2.*(Psi+alpha))*(an->cDi-2.*an->ci+9./7.)/2. - an->sHi4*sin(2.*(Psi-alpha))*(an->cDi+2.*an->ci+9./7.)/2. );
904 
905  hL4->data[-2+os] = amp42 * (1. - 3.*eta) * ( 16.*an->sHi6*an->cHi2*cos(4.*Psi-2.*alpha) + 16.*an->cHi6*an->sHi2*cos(4.*Psi+2.*alpha) - an->cHi4*cos(2.*(Psi+alpha))*(an->cDi-2.*an->ci+9./7.)/2. - an->sHi4*cos(2.*(Psi-alpha))*(an->cDi+2.*an->ci+9./7.)/2. );
906 
907  hL4->data[-2+os]+=-amp42 *I*(1. - 3.*eta) * ( 16.*an->sHi6*an->cHi2*sin(4.*Psi-2.*alpha) - 16.*an->cHi6*an->sHi2*sin(4.*Psi+2.*alpha) + an->cHi4*sin(2.*(Psi+alpha))*(an->cDi-2.*an->ci+9./7.)/2. - an->sHi4*sin(2.*(Psi-alpha))*(an->cDi+2.*an->ci+9./7.)/2. );
908 
909  hL4->data[1+os] = amp41 * (1. - 3.*eta) * ( -64.*an->sHi5*an->cHi3*cos(4.*Psi-alpha) + 64.*an->sHi3*an->cHi5*cos(4.*Psi+alpha) - an->sHi3*cos(2.*Psi-alpha)*((an->cDi*an->cHi-an->sDi*an->sHi) + 2.*(an->cHi*an->ci-an->sHi*an->si) + 19./7.*an->cHi) + an->cHi3*cos(2.*Psi+alpha)*((an->cDi*an->sHi+an->sDi*an->cHi) - 2.*(an->si*an->cHi+an->ci*an->si2) +19./7.*an->cHi) );
910 
911  hL4->data[1+os]+= amp41*I*(1. - 3.*eta) * ( -64.*an->sHi5*an->cHi3 * sin(4.*Psi-alpha) - 64.*an->sHi3*an->cHi5 * sin(4.*Psi+alpha) - an->sHi3*sin(2.*Psi-alpha)*((an->cDi*an->cHi-an->sDi*an->sHi) + 2.*(an->cHi*an->ci-an->sHi*an->si) + 19./7.*an->cHi) - an->cHi3*sin(2.*Psi+alpha)*((an->cDi*an->sHi+an->sDi*an->cHi) - 2.*(an->si*an->cHi+an->ci*an->sHi) + 19./7.*an->cHi) );
912 
913  hL4->data[-1+os] = -amp41 * (1. - 3.*eta) * ( -64*an->sHi5*an->cHi3 * cos(4.*Psi-alpha) + 64.*an->sHi3*an->cHi5*cos(4.*Psi+alpha) - an->sHi3*cos(2.*Psi-alpha)*((an->cDi*an->cHi-an->sDi*an->sHi) + 2.*(an->cHi*an->ci-an->sHi*an->si) + 19./7.*an->ci2) + an->cHi3*cos(2.*Psi+alpha)*((an->cDi*an->sHi+an->sDi*an->cHi) - 2.*(an->si*an->cHi+an->ci*an->sHi) + 19./7.*an->ci2) );
914 
915  hL4->data[-1+os]+= amp41 *I*(1. - 3.*eta) *I*( -64.*an->sHi5*an->cHi3 * sin(4.*Psi-alpha) - 64.*an->sHi3*an->cHi5 * sin(4.*Psi+alpha) - an->sHi3*sin(2.*Psi-alpha)*((an->cDi*an->cHi-an->sDi*an->sHi) + 2.*(an->cHi*an->ci-an->sHi*an->si) + 19./7.*an->ci2) - an->cHi3*sin(2.*Psi+alpha)*((an->cDi*an->sHi+an->sDi*an->cHi) - 2.*(an->si*an->cHi+an->ci*an->sHi) + 19./7.*an->cHi) );
916 
917  hL4->data[os] = amp40 * (1.-3.*eta) * an->si2 * (8.*an->si2*cos(4.*Psi) + cos(2.*Psi)*(an->cDi+5./7.) );
918 
919  return XLAL_SUCCESS;
920 } /* End of XLALSimSpinInspiralFillL4Modes*/
921 
922 static INT4 XLALSimInspiralSpinTaylorT4Engine(REAL8TimeSeries **omega, /**< post-Newtonian parameter [returned]*/
923  REAL8TimeSeries **Phi, /**< orbital phase [returned]*/
924  REAL8TimeSeries **LNhatx, /**< LNhat vector x component [returned]*/
925  REAL8TimeSeries **LNhaty, /**< " " " y component [returned]*/
926  REAL8TimeSeries **LNhatz, /**< " " " z component [returned]*/
927  REAL8TimeSeries **S1x, /**< Spin1 vector x component [returned]*/
928  REAL8TimeSeries **S1y, /**< " " " y component [returned]*/
929  REAL8TimeSeries **S1z, /**< " " " z component [returned]*/
930  REAL8TimeSeries **S2x, /**< Spin2 vector x component [returned]*/
931  REAL8TimeSeries **S2y, /**< " " " y component [returned]*/
932  REAL8TimeSeries **S2z, /**< " " " z component [returned]*/
933  REAL8TimeSeries **Energy, /**< Energy [returned]*/
934  const REAL8 yinit[], /**< UNDOCUMENTED */
935  const INT4 lengthH, /**< UNDOCUMENTED */
936  const Approximant approx, /**< Allow to choose w/o ringdown */
937  LALSimInspiralPhenSpinTaylorT4Coeffs *params /**< UNDOCUMENTED */
938  )
939 {
940  UINT4 idx;
941  INT4 jdx;
942  UINT4 intLen;
943  INT4 intReturn;
944 
945  REAL8 S1x0,S1y0,S1z0,S2x0,S2y0,S2z0; /** Used to store initial spin values */
946  REAL8Array *yout; /** Used to store integration output */
947 
948  LALAdaptiveRungeKuttaIntegrator *integrator;
949 
950  /* allocate the integrator */
951  if (approx == PhenSpinTaylor)
952  integrator = XLALAdaptiveRungeKutta4Init(LAL_NUM_PST4_VARIABLES,XLALSpinInspiralDerivatives,XLALSimSpinInspiralTest,LAL_PST4_ABSOLUTE_TOLERANCE,LAL_PST4_RELATIVE_TOLERANCE);
953  else
954  integrator = XLALAdaptiveRungeKutta4Init(LAL_NUM_PST4_VARIABLES,XLALSpinInspiralDerivatives,XLALSimIMRPhenSpinTest,LAL_PST4_ABSOLUTE_TOLERANCE,LAL_PST4_RELATIVE_TOLERANCE);
955 
956  if (!integrator) {
957  XLALPrintError("XLAL Error - %s: Cannot allocate integrator\n", __func__);
958  XLAL_ERROR(XLAL_EFUNC);
959  }
960 
961  /* stop the integration only when the test is true */
962  integrator->stopontestonly = 1;
963 
964  REAL8 *yin = (REAL8 *) LALMalloc(sizeof(REAL8) * LAL_NUM_PST4_VARIABLES);
965  for (idx=0; idx<LAL_NUM_PST4_VARIABLES; idx++) yin[idx]=yinit[idx];
966  S1x0=yinit[5];
967  S1y0=yinit[6];
968  S1z0=yinit[7];
969  S2x0=yinit[8];
970  S2y0=yinit[9];
971  S2z0=yinit[10];
972 
973  //REAL8 dtInt=1./OmMatch(0,0,0,0,0)/50.*fabs(params->dt)/params->dt;
974  REAL8 length=((REAL8)lengthH)*fabs(params->dt)/params->M;
975  intLen = XLALAdaptiveRungeKutta4Hermite(integrator,(void *)params,yin,0.0,length,params->dt/params->M,&yout);
976 
977  intReturn = integrator->returncode;
978  XLALAdaptiveRungeKuttaFree(integrator);
979 
980  if (intReturn == XLAL_FAILURE) {
981  XLALPrintError("** LALSimIMRPSpinInspiralRD Error **: Adaptive Integrator\n");
982  XLALPrintError(" m: %12.4e %12.4e Mom %12.4e\n",params->m1ByM*params->M,params->m2ByM*params->M,params->fStart);
983  XLALPrintError(" S1: %12.4e %12.4e %12.4e\n",S1x0,S1y0,S1z0);
984  XLALPrintError(" S2: %12.4e %12.4e %12.4e\n",S2x0,S2y0,S2z0);
985  XLAL_ERROR(XLAL_EFUNC);
986  }
987  /* End integration*/
988 
989  /* Start of the integration checks*/
990  if (intLen<minIntLen) {
991  XLALPrintError("** LALSimIMRPSpinInspiralRD ERROR **: integration too short! intReturnCode %d, integration length %d, at least %d required\n",intReturn,intLen,minIntLen);
992  if (XLALClearErrno() == XLAL_ENOMEM) {
993  XLAL_ERROR( XLAL_ENOMEM);
994  } else {
995  XLAL_ERROR( XLAL_EFAILED);
996  }
997  }
998 
999  const LIGOTimeGPS tStart=LIGOTIMEGPSZERO;
1000  *omega = XLALCreateREAL8TimeSeries( "OMEGA", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
1001  *Phi = XLALCreateREAL8TimeSeries( "ORBITAL_PHASE", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
1002  *LNhatx = XLALCreateREAL8TimeSeries( "LNHAT_X_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
1003  *LNhaty = XLALCreateREAL8TimeSeries( "LNHAT_Y_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
1004  *LNhatz = XLALCreateREAL8TimeSeries( "LNHAT_Z_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
1005  *S1x = XLALCreateREAL8TimeSeries( "SPIN1_X_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
1006  *S1y = XLALCreateREAL8TimeSeries( "SPIN1_Y_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
1007  *S1z = XLALCreateREAL8TimeSeries( "SPIN1_Z_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
1008  *S2x = XLALCreateREAL8TimeSeries( "SPIN2_X_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
1009  *S2y = XLALCreateREAL8TimeSeries( "SPIN2_Y_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
1010  *S2z = XLALCreateREAL8TimeSeries( "SPIN2_Z_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
1011  *Energy = XLALCreateREAL8TimeSeries( "LNHAT_Z_COMPONENT", &tStart, 0., params->dt, &lalDimensionlessUnit, intLen);
1012  if ( !omega || !Phi || !S1x || !S1y || !S1z || !S2x || !S2y || !S2z || !LNhatx || !LNhaty || !LNhatz || !Energy ) {
1013  XLALDestroyREAL8Array(yout);
1014  XLAL_ERROR(XLAL_EFUNC);
1015  }
1016 
1017  /* Copy dynamical variables from yout array to output time series.
1018  * Note the first 'len' members of yout are the time steps.
1019  */
1020  INT4 sign=params->dt > 0. ? 1 : -1;
1021  jdx= (intLen-1)*(-sign+1)/2;
1022 
1023  for (idx=0;idx<intLen;idx++) {
1024  (*Phi)->data->data[idx] = yout->data[intLen+jdx];
1025  (*omega)->data->data[idx] = yout->data[2*intLen+jdx];
1026  (*LNhatx)->data->data[idx] = yout->data[3*intLen+jdx];
1027  (*LNhaty)->data->data[idx] = yout->data[4*intLen+jdx];
1028  (*LNhatz)->data->data[idx] = yout->data[5*intLen+jdx];
1029  (*S1x)->data->data[idx] = yout->data[6*intLen+jdx];
1030  (*S1y)->data->data[idx] = yout->data[7*intLen+jdx];
1031  (*S1z)->data->data[idx] = yout->data[8*intLen+jdx];
1032  (*S2x)->data->data[idx] = yout->data[9*intLen+jdx];
1033  (*S2y)->data->data[idx] = yout->data[10*intLen+jdx];
1034  (*S2z)->data->data[idx] = yout->data[11*intLen+jdx];
1035  (*Energy)->data->data[idx] = yout->data[12*intLen+jdx];
1036  jdx+=sign;
1037  }
1038 
1039  XLALDestroyREAL8Array(yout);
1040  return intReturn;
1041 } /* End of XLALSimInspiralSpinTaylorT4Engine */
1042 
1043 static INT4 XLALSimInspiralComputeInclAngle(REAL8 ciota, LALSimInspiralInclAngle *angle){
1044  angle->ci=ciota;
1045  angle->si=sqrt(1.-ciota*ciota);
1046  angle->ci2=angle->ci*angle->ci;
1047  angle->si2=angle->si*angle->si;
1048  angle->cDi=angle->ci*angle->ci-angle->si*angle->si;
1049  angle->sDi=2.*angle->ci*angle->si;
1050  angle->cHi=sqrt((1.+angle->ci)/2.);
1051  angle->sHi=sqrt((1.-angle->ci)/2.);
1052  angle->cHi2=(1.+angle->ci)/2.;
1053  angle->sHi2=(1.-angle->ci)/2.;
1054  angle->cHi3=angle->cHi*angle->cHi2;
1055  angle->sHi3=angle->sHi*angle->sHi2;
1056  angle->cHi4=angle->cHi2*angle->cHi2;
1057  angle->sHi4=angle->sHi2*angle->sHi2;
1058  angle->cHi6=angle->cHi2*angle->cHi4;
1059  angle->sHi6=angle->sHi2*angle->sHi4;
1060  angle->cHi8=angle->cHi4*angle->cHi4;
1061  angle->sHi8=angle->sHi4*angle->sHi4;
1062 
1063  return XLAL_SUCCESS;
1064 
1065 } /* End of XLALSimInspiralComputeInclAngle*/
1066 
1067 /**
1068  * The following lines are necessary in the case L is initially parallel to
1069  * N so that alpha is undefined at the beginning but different from zero at the first
1070  * step (this happens if the spins are not aligned with L).
1071  * Such a discontinuity of alpha would induce
1072  * a discontinuity of the waveform between its initial value and its value after the
1073  * first integration step. This does not happen during the integration as in that
1074  * case alpha can be safely set to the previous value, just before L becomes parallel
1075  * to N. In the case L stays all the time parallel to N than alpha can be
1076  * safely set to zero, as it is.
1077  */
1078 
1079 static INT4 XLALSimInspiralComputeAlpha(LALSimInspiralPhenSpinTaylorT4Coeffs params, REAL8 LNhx, REAL8 LNhy, REAL8 S1x, REAL8 S1y, REAL8 S2x, REAL8 S2y,REAL8 *alpha){
1080  if ((LNhy*LNhy+LNhx*LNhx)==0.) {
1081  REAL8 S1xy=S1x*S1x+S1y*S1y;
1082  REAL8 S2xy=S2x*S2x+S2y*S2y;
1083  if ((S1xy+S2xy)==0.) {
1084  *alpha=0.;
1085  }
1086  else {
1087  REAL8 c1=0.75+params.eta/2-0.75*(params.m1ByM-params.m2ByM);
1088  REAL8 c2=0.75+params.eta/2+0.75*(params.m1ByM-params.m2ByM);
1089  *alpha=atan2(-c1*S1x-c2*S2x,c1*S1y+c2*S2y);
1090  }
1091  }
1092  else {
1093  *alpha=atan2(LNhy,LNhx);
1094  }
1095  return XLAL_SUCCESS;
1096 } /*End of XLALSimInspiralComputeAlpha*/
1097 
1098 /**
1099  * Here we use the following convention:
1100  * the coordinates of the spin vectors spin1,2 and the inclination
1101  * variable refers to different physical parameters according to the value of
1102  * axisChoice:
1103  *
1104  * * LAL_SIM_INSPIRAL_FRAME_AXIS_TOTAL_J: inclination denotes the angle
1105  * between reference J and the view direction N
1106  * (for N=(0,0,1), Jhat=(sin(inc),0,cos(inc)) )
1107  * and the spins are given with respect to initial J.
1108  * Note that not all values of the spin will correspond
1109  * to physically well defined configurations.
1110  *
1111  * * LAL_SIM_INSPIRAL_FRAME_AXIS_ORBITAL_L: inclination denotes the angle
1112  * between the reference L and the view direction N
1113  * (for N=(0,0,1), LNhat=(sin(inc),0,cos(inc)) )
1114  * and the spins are given wrt initial L.
1115  *
1116  * * LAL_SIM_INSPIRAL_FRAME_AXIS_VIEW: inclination denotes the angle
1117  * between the reference L and the view direction N
1118  * (for N=(0,0,1), LNhat=(sin(inc),0,cos(inc)) )
1119  * and the spins are given wrt to N.
1120  * (default)
1121  *
1122  * The spin magnitudes are normalized to the individual mass^2, i.e.
1123  * they are dimension-less.
1124  * The modulus of the initial angular momentum is fixed by m1,m2 and
1125  * initial frequency.
1126  * The polarization angle is not used here, it enters the pattern
1127  * functions along with the angles marking the sky position of the
1128  * source.
1129  */
1130 
1131 /*static void rotateX(REAL8 phi,REAL8 *vx, REAL8 *vy, REAL8 *vz){
1132  REAL8 tmp[3]={*vx,*vy,*vz};
1133  *vx=*vy=*vz=0.;
1134  REAL8 rotX[3][3]={{1.,0.,0.},{0,cos(phi),-sin(phi)},{0,sin(phi),cos(phi)}};
1135  INT4 idx;
1136  for (idx=0;idx<3;idx++) {
1137  *vx+=rotX[0][idx]*tmp[idx];
1138  *vy+=rotX[1][idx]*tmp[idx];
1139  *vz+=rotX[2][idx]*tmp[idx];
1140  }
1141  }*/
1142 static void rotateY(REAL8 phi,REAL8 *vx, REAL8 *vy, REAL8 *vz){
1143  REAL8 rotY[3][3]={{cos(phi),0.,sin(phi)},{0.,1.,0.},{-sin(phi),0.,cos(phi)}};
1144  REAL8 tmp[3]={*vx,*vy,*vz};
1145  *vx=*vy=*vz=0.;
1146  INT4 idx;
1147  for (idx=0;idx<3;idx++) {
1148  *vx+=rotY[0][idx]*tmp[idx];
1149  *vy+=rotY[1][idx]*tmp[idx];
1150  *vz+=rotY[2][idx]*tmp[idx];
1151  }
1152 }
1153 static void rotateZ(REAL8 phi,REAL8 *vx, REAL8 *vy, REAL8 *vz){
1154  REAL8 tmp[3]={*vx,*vy,*vz};
1155  REAL8 rotZ[3][3]={{cos(phi),-sin(phi),0.},{sin(phi),cos(phi),0.},{0.,0.,1.}};
1156  *vx=*vy=*vz=0.;
1157  INT4 idx;
1158  for (idx=0;idx<3;idx++) {
1159  *vx+=rotZ[0][idx]*tmp[idx];
1160  *vy+=rotZ[1][idx]*tmp[idx];
1161  *vz+=rotZ[2][idx]*tmp[idx];
1162  }
1163 }
1164 
1165 static INT4 XLALSimIMRPhenSpinInspiralSetAxis(REAL8 **yinitOut, /* [returned] */
1166  REAL8 *iota, /* [returned] */
1167  REAL8 *phiN, /* azimuthal angle of the view direction [returned]*/
1168  REAL8 *yinitIn, /* initial values of dynamical variables */
1169  const REAL8 inclination, /* input*/
1170  const REAL8 mass1, /* in MSun units */
1171  const REAL8 mass2, /* in MSun units */
1172  LALSimInspiralFrameAxis axisChoice)
1173 {
1174  // Magnitude of the Newtonian orbital angular momentum
1175  REAL8 omega=yinitIn[1];
1176  REAL8 LNmag = mass1*mass2 / cbrt(omega);
1177  REAL8 Mass = mass1 + mass2;
1178  REAL8 S1[3],S2[3];
1179  REAL8 phiJ=0.;
1180  REAL8 thetaJ=0.;
1181  REAL8 LNh[3],N[3],J[3];
1182  REAL8 Jmag=0.;
1183  INT4 idx;
1184 
1185  // Physical values of the spins
1186  S1[0] = yinitIn[5] * mass1 * mass1;
1187  S1[1] = yinitIn[6] * mass1 * mass1;
1188  S1[2] = yinitIn[7] * mass1 * mass1;
1189  S2[0] = yinitIn[8] * mass2 * mass2;
1190  S2[1] = yinitIn[9] * mass2 * mass2;
1191  S2[2] = yinitIn[10] * mass2 * mass2;
1192 
1193  switch (axisChoice) {
1194 
1195  case LAL_SIM_INSPIRAL_FRAME_AXIS_TOTAL_J:
1196  LNh[0] = -(S1[0]+S2[0])/LNmag;
1197  LNh[1] = -(S1[1]+S2[1])/LNmag;
1198  REAL8 LNhxy2 = LNh[0]*LNh[0]+LNh[1]*LNh[1];
1199  LNh[2]=0.;
1200  if (LNhxy2<=1.) {
1201  LNh[2]=sqrt(1.-sqrt(LNhxy2));
1202  if ( (LNh[2]*LNmag)<(S1[2]+S2[2]) ) {
1203  XLALPrintError("** LALSimIMRPSpinInspiralRD error *** for s1 (%12.4e %12.4e %12.4e)\n",S1[0],S1[1],S1[2]);
1204  XLALPrintError(" s2 (%12.4e %12.4e %12.4e)\n",S2[0],S2[1],S2[2]);
1205  XLALPrintError(" wrt to J for m: (%12.4e %12.4e) and v= %12.4e\n",mass1,mass2,cbrt(omega));
1206  XLALPrintError(" it is impossible to determine the sign of LNhz\n");
1207  XLAL_ERROR(XLAL_EDOM);
1208  }
1209  }
1210  else {
1211  XLALPrintError("** LALSimIMRPSpinInspiralRD error *** unphysical values of s1 (%12.4e %12.4e %12.4e)\n",S1[0],S1[1],S1[2]);
1212  XLALPrintError(" s2 (%12.4e %12.4e %12.4e)\n",S2[0],S2[1],S2[2]);
1213  XLALPrintError(" wrt to J for m: (%12.4e %12.4e) and v= %12.4e\n",mass1,mass2,cbrt(omega));
1214  XLAL_ERROR(XLAL_EDOM);
1215  }
1216  *iota=inclination;
1217  *phiN=LAL_PI;
1218  break;
1219 
1220  case LAL_SIM_INSPIRAL_FRAME_AXIS_ORBITAL_L:
1221  J[0]=S1[0]+S2[0];
1222  J[1]=S1[1]+S2[1];
1223  J[2]=S1[2]+S2[2]+LNmag;
1224  N[0]=-sin(inclination);
1225  N[1]=0.;
1226  N[2]=cos(inclination);
1227  LNh[0]=0.;
1228  LNh[1]=0.;
1229  LNh[2]=1.;
1230  Jmag=sqrt(J[0]*J[0]+J[1]*J[1]+J[2]*J[2]);
1231  if (Jmag>0.) {
1232  phiJ=atan2(J[1],J[0]);
1233  thetaJ=acos(J[2]/Jmag);
1234  }
1235  rotateZ(-phiJ,&N[0],&N[1],&N[2]);
1236  rotateY(-thetaJ,&N[0],&N[1],&N[2]);
1237  *iota=acos(N[2]);
1238  *phiN=atan2(N[1],N[0]);
1239 
1240  rotateZ(-phiJ,&J[0],&J[1],&J[2]);
1241  rotateY(-thetaJ,&J[0],&J[1],&J[2]);
1242  printf("Check J: %12.4e %12.4e %12.4e\n",J[0],J[1],J[2]);
1243  printf("Check Jmag: %12.4e\n",Jmag);
1244 
1245  break;
1246 
1247  case LAL_SIM_INSPIRAL_FRAME_AXIS_VIEW:
1248  default:
1249  LNh[0] = sin(inclination);
1250  LNh[1] = 0.;
1251  LNh[2] = cos(inclination);
1252  J[0]=S1[0]+S2[0]+LNh[0]*LNmag;
1253  J[1]=S1[1]+S2[1];
1254  J[2]=S1[2]+S2[2]+LNh[2]*LNmag;
1255  N[0]=0.;
1256  N[1]=0.;
1257  N[2]=1.;
1258  Jmag=sqrt(J[0]*J[0]+J[1]*J[1]+J[2]*J[2]);
1259  if (Jmag>0.) {
1260  phiJ=atan2(J[1],J[0]);
1261  thetaJ=acos(J[2]/Jmag);
1262  }
1263  rotateZ(-phiJ,&N[0],&N[1],&N[2]);
1264  rotateY(-thetaJ,&N[0],&N[1],&N[2]);
1265  *iota=acos(N[2]);
1266  *phiN=atan2(N[1],N[0]);
1267 
1268  rotateZ(-phiJ,&J[0],&J[1],&J[2]);
1269  rotateY(-thetaJ,&J[0],&J[1],&J[2]);
1270  printf("Check J: %12.4e %12.4e %12.4e\n",J[0],J[1],J[2]);
1271  printf("Check Jmag: %12.4e\n",Jmag);
1272 
1273  break;
1274  }
1275 
1276  rotateZ(-phiJ, &S1[0], &S1[1], &S1[2]);
1277  rotateZ(-phiJ, &S2[0], &S2[1], &S2[2]);
1278  rotateZ(-phiJ, &LNh[0],&LNh[1],&LNh[2]);
1279  rotateY(-thetaJ,&S1[0], &S1[1], &S1[2]);
1280  rotateY(-thetaJ,&S2[0], &S2[1], &S2[2]);
1281  rotateY(-thetaJ,&LNh[0],&LNh[1],&LNh[2]);
1282 
1283  *yinitOut = (REAL8 *) LALMalloc(sizeof(REAL8) * LAL_NUM_PST4_VARIABLES);
1284  *yinitOut[0]=yinitIn[0];
1285  *yinitOut[1]=yinitIn[1];
1286  *yinitOut[2]=LNh[0];
1287  *yinitOut[3]=LNh[1];
1288  *yinitOut[4]=LNh[2];
1289  *yinitOut[5]=S1[0]/Mass/Mass;
1290  *yinitOut[6]=S1[1]/Mass/Mass;
1291  *yinitOut[7]=S1[2]/Mass/Mass;
1292  *yinitOut[8]=S2[0]/Mass/Mass;
1293  *yinitOut[9]=S2[1]/Mass/Mass;
1294  *yinitOut[10]=S2[2]/Mass/Mass;
1295  for (idx=11;idx<LAL_NUM_PST4_VARIABLES;idx++)
1296  *yinitOut[idx]=yinitIn[idx];
1297 
1298  return XLAL_SUCCESS;
1299 
1300 } /* End of XLALSimIMRPhenSpinInspiralSetAxis*/
1301 
1302 /**
1303  * PhenSpin Initialization
1304  */
1305 
1306 static INT4 XLALSimIMRPhenSpinInitialize(REAL8 mass1, /**< in Msun units */
1307  REAL8 mass2, /**< in Msun units */
1308  REAL8 lambda1, /**< Tidal par1*/
1309  REAL8 lambda2, /**< Tidal par2*/
1310  REAL8 quadparam1, /**< Quad-monopole par1*/
1311  REAL8 quadparam2, /**< Quad-monopole par2*/
1312  REAL8 *yinit, /**< Initial values*/
1313  REAL8 fStart, /**< in Hz*/
1314  REAL8 fEnd, /**< in Hz*/
1315  REAL8 deltaT, /**< sampling time (sec)*/
1316  INT4 phaseO, /**< (twice) phase order*/
1317  LALSimInspiralPhenSpinTaylorT4Coeffs *params, /**< Coefficient holder*/
1318  LALDict *LALparams, /**< Extra parameters */
1319  Approximant approx /**< Approximant*/)
1320 {
1321  if (fStart<=0.) {
1322  XLALPrintError("** LALSimIMRPSpinInspiralRD error *** non >ve value of fMin %12.4e\n",fStart);
1323  XLAL_ERROR(XLAL_EDOM);
1324  }
1325 
1326  REAL8 S1x=yinit[5];
1327  REAL8 S1y=yinit[6];
1328  REAL8 S1z=yinit[7];
1329  REAL8 S2x=yinit[8];
1330  REAL8 S2y=yinit[9];
1331  REAL8 S2z=yinit[10];
1332 
1333  REAL8 LNhS1 = S1z;
1334  REAL8 LNhS2 = S2z;
1335  REAL8 S1S1 = S1x*S1x + S1y*S1y + S1z*S1z;
1336  REAL8 S1S2 = S1x*S2x + S1y*S2y + S1z*S2z;
1337  REAL8 S2S2 = S2x*S2x + S2y*S2y + S2z*S2z;
1338  REAL8 unitHz = (mass1+mass2)*LAL_MTSUN_SI; /* convert m from msun to seconds */
1339  REAL8 initOmega = fStart*unitHz * (REAL8) LAL_PI;
1340  REAL8 omegaMatch = OmMatch(LNhS1,LNhS2,S1S1,S1S2,S2S2);
1341  yinit[1]=initOmega;
1342 
1343  if (approx==PhenSpinTaylorRD) {
1344  if ( initOmega > omegaMatch ) {
1345  if ((S1x==S1y)&&(S1x==0)&&(S2x==S2y)&&(S2y==0.)) {
1346  initOmega = 0.95*omegaMatch;
1347  yinit[1]=initOmega;
1348  XLALPrintWarning("*** LALSimIMRPSpinInspiralRD WARNING ***: Initial frequency reset from %12.6e to %12.6e Hz, m:(%12.4e,%12.4e)\n",fStart,initOmega/unitHz/LAL_PI,mass1,mass2);
1349  }
1350  else {
1351  XLALPrintError("*** LALSimIMRPSpinInspiralRD ERROR ***: Initial frequency %12.6e Hz too high, as fMatch estimated %12.6e Hz, m:(%12.4e,%12.4e)\n",fStart,omegaMatch/unitHz/LAL_PI,mass1,mass2);
1352  XLAL_ERROR(XLAL_EFAILED);
1353  }
1354  }
1355  }
1356 
1357  /* setup coefficients for PN equations */
1358  if(XLALSimIMRPhenSpinParamsSetup(params,deltaT,fStart,fEnd,mass1,mass2,lambda1,lambda2,quadparam1,quadparam2,XLALSimInspiralWaveformParamsLookupPNSpinOrder(LALparams),XLALSimInspiralWaveformParamsLookupPNTidalOrder(LALparams),LALparams,phaseO)) {
1359  XLAL_ERROR(XLAL_ENOMEM);
1360  }
1361 
1362  return XLAL_SUCCESS;
1363 
1364 } /* End of XLALSimIMRPhenSpinInitialize*/
1365 
1366 /* Appends the start and end time series together, skipping the redundant last
1367  * sample of begin. Frees end before returning a pointer to the result, which is
1368  * the resized start series. */
1369 static REAL8TimeSeries *appendTSandFree(REAL8TimeSeries *start, REAL8TimeSeries *end) {
1370  UINT4 origlen = start->data->length;
1371  start = XLALResizeREAL8TimeSeries(start, 0,
1372  start->data->length + end->data->length - 1);
1373 
1374  memcpy(start->data->data + origlen -2, end->data->data,
1375  (end->data->length)*sizeof(REAL8));
1376 
1377  XLALGPSAdd(&(start->epoch), -end->deltaT*(end->data->length - 1));
1378 
1379  XLALDestroyREAL8TimeSeries(end);
1380 
1381  return start;
1382 }
1383 
1384 
1385 static INT4 XLALSimIMRHybridRingdownWave(
1386  REAL8Vector *rdwave1, /**< Real part of ringdown */
1387  REAL8Vector *rdwave2, /**< Imaginary part of ringdown */
1388  const REAL8 dt, /**< Sampling interval */
1389  const REAL8 mass1, /**< First component mass (in Solar masses) */
1390  const REAL8 mass2, /**< Second component mass (in Solar masses) */
1391  REAL8VectorSequence *inspwave1, /**< Values and derivatives of real part of inspiral waveform */
1392  REAL8VectorSequence *inspwave2, /**< Values and derivatives of Imaginary part of inspiral waveform */
1393  COMPLEX16Vector *modefreqs, /**< Complex frequencies of ringdown (scaled by total mass) */
1394  REAL8Vector *matchrange /**< Times which determine the comb size for ringdown attachment */
1395 )
1396 {
1397 
1398  /* XLAL error handling */
1399  INT4 errcode = XLAL_SUCCESS;
1400 
1401  /* For checking GSL return codes */
1402  INT4 gslStatus;
1403 
1404  UINT4 i, j, k, nmodes = 8;
1405 
1406  /* Sampling rate from input */
1407  REAL8 t1, t2, t3, t4, t5, rt;
1408  gsl_matrix *coef;
1409  gsl_vector *hderivs;
1410  gsl_vector *x;
1411  gsl_permutation *p;
1412  REAL8Vector *modeamps;
1413  INT4 s;
1414  REAL8 tj=0.;
1415  REAL8 m;
1416 
1417  /* mass in geometric units */
1418  m = (mass1 + mass2) * LAL_MTSUN_SI;
1419  t5 = (matchrange->data[0] - matchrange->data[1]) * m;
1420  rt = -t5 / 5.;
1421 
1422  t4 = t5 + rt;
1423  t3 = t4 + rt;
1424  t2 = t3 + rt;
1425  t1 = t2 + rt;
1426 
1427  // printf(" ** t1 %12.4e t5 %12.4e\n",t1,t5);
1428 
1429  if ( inspwave1->length != 2 || inspwave2->length != 2 ||
1430  modefreqs->length != nmodes )
1431  {
1432  XLAL_ERROR( XLAL_EBADLEN );
1433  }
1434 
1435  /* Solving the linear system for QNMs amplitude coefficients using gsl routine */
1436  /* Initiate matrices and supporting variables */
1437  XLAL_CALLGSL( coef = (gsl_matrix *) gsl_matrix_alloc(2 * nmodes, 2 * nmodes) );
1438  XLAL_CALLGSL( hderivs = (gsl_vector *) gsl_vector_alloc(2 * nmodes) );
1439  XLAL_CALLGSL( x = (gsl_vector *) gsl_vector_alloc(2 * nmodes) );
1440  XLAL_CALLGSL( p = (gsl_permutation *) gsl_permutation_alloc(2 * nmodes) );
1441 
1442  /* Check all matrices and variables were allocated */
1443  if ( !coef || !hderivs || !x || !p )
1444  {
1445  if (coef) gsl_matrix_free(coef);
1446  if (hderivs) gsl_vector_free(hderivs);
1447  if (x) gsl_vector_free(x);
1448  if (p) gsl_permutation_free(p);
1449 
1450  XLAL_ERROR( XLAL_ENOMEM );
1451  }
1452 
1453  /* Define the linear system Ax=y */
1454  /* Matrix A (2*n by 2*n) has block symmetry. Define half of A here as "coef" */
1455  /* Define y here as "hderivs" */
1456  for (i = 0; i < nmodes; ++i)
1457  {
1458  gsl_matrix_set(coef, 0, i, 1);
1459  gsl_matrix_set(coef, 1, i, - cimag(modefreqs->data[i]));
1460  gsl_matrix_set(coef, 2, i, exp(-cimag(modefreqs->data[i])*t1) * cos(creal(modefreqs->data[i])*t1));
1461  gsl_matrix_set(coef, 3, i, exp(-cimag(modefreqs->data[i])*t2) * cos(creal(modefreqs->data[i])*t2));
1462  gsl_matrix_set(coef, 4, i, exp(-cimag(modefreqs->data[i])*t3) * cos(creal(modefreqs->data[i])*t3));
1463  gsl_matrix_set(coef, 5, i, exp(-cimag(modefreqs->data[i])*t4) * cos(creal(modefreqs->data[i])*t4));
1464  gsl_matrix_set(coef, 6, i, exp(-cimag(modefreqs->data[i])*t5) * cos(creal(modefreqs->data[i])*t5));
1465  gsl_matrix_set(coef, 7, i, exp(-cimag(modefreqs->data[i])*t5) *
1466  (-cimag(modefreqs->data[i]) * cos(creal(modefreqs->data[i])*t5)
1467  -creal(modefreqs->data[i]) * sin(creal(modefreqs->data[i])*t5) ));
1468  gsl_matrix_set(coef, 8, i, 0);
1469  gsl_matrix_set(coef, 9, i, -creal(modefreqs->data[i]));
1470  gsl_matrix_set(coef, 10, i, -exp(-cimag(modefreqs->data[i])*t1) * sin(creal(modefreqs->data[i])*t1));
1471  gsl_matrix_set(coef, 11, i, -exp(-cimag(modefreqs->data[i])*t2) * sin(creal(modefreqs->data[i])*t2));
1472  gsl_matrix_set(coef, 12, i, -exp(-cimag(modefreqs->data[i])*t3) * sin(creal(modefreqs->data[i])*t3));
1473  gsl_matrix_set(coef, 13, i, -exp(-cimag(modefreqs->data[i])*t4) * sin(creal(modefreqs->data[i])*t4));
1474  gsl_matrix_set(coef, 14, i, -exp(-cimag(modefreqs->data[i])*t5) * sin(creal(modefreqs->data[i])*t5));
1475  gsl_matrix_set(coef, 15, i, -exp(-cimag(modefreqs->data[i])*t5) *
1476  ( cimag(modefreqs->data[i]) * sin(creal(modefreqs->data[i])*t5)
1477  -creal(modefreqs->data[i]) * cos(creal(modefreqs->data[i])*t5)));
1478  }
1479 
1480  gsl_vector_set(hderivs, 0, inspwave1->data[5]);
1481  gsl_vector_set(hderivs, 0 + nmodes, inspwave2->data[5]);
1482  gsl_vector_set(hderivs, 1, inspwave1->data[11]);
1483  gsl_vector_set(hderivs, 1 + nmodes, inspwave2->data[11]);
1484  gsl_vector_set(hderivs, 2, inspwave1->data[4]);
1485  gsl_vector_set(hderivs, 2 + nmodes, inspwave2->data[4]);
1486  gsl_vector_set(hderivs, 3, inspwave1->data[3]);
1487  gsl_vector_set(hderivs, 3 + nmodes, inspwave2->data[3]);
1488  gsl_vector_set(hderivs, 4, inspwave1->data[2]);
1489  gsl_vector_set(hderivs, 4 + nmodes, inspwave2->data[2]);
1490  gsl_vector_set(hderivs, 5, inspwave1->data[1]);
1491  gsl_vector_set(hderivs, 5 + nmodes, inspwave2->data[1]);
1492  gsl_vector_set(hderivs, 6, inspwave1->data[0]);
1493  gsl_vector_set(hderivs, 6 + nmodes, inspwave2->data[0]);
1494  gsl_vector_set(hderivs, 7, inspwave1->data[6]);
1495  gsl_vector_set(hderivs, 7 + nmodes, inspwave2->data[6]);
1496 
1497  /* Complete the definition for the rest half of A */
1498  for (i = 0; i < nmodes; ++i)
1499  {
1500  for (k = 0; k < nmodes; ++k)
1501  {
1502  gsl_matrix_set(coef, i, k + nmodes, - gsl_matrix_get(coef, i + nmodes, k));
1503  gsl_matrix_set(coef, i + nmodes, k + nmodes, gsl_matrix_get(coef, i, k));
1504  }
1505  }
1506 
1507 #if DEBUG_RD
1508  printf("\nRingdown matching matrix:\n");
1509  for (i = 0; i < 16; ++i) {
1510  for (j = 0; j < 16; ++j) {
1511  printf("%8.1e ",gsl_matrix_get(coef,i,j));
1512  }
1513  printf(" | %8.1e\n",gsl_vector_get(hderivs,i));
1514  }
1515 #endif
1516 
1517  /* Call gsl LU decomposition to solve the linear system */
1518  XLAL_CALLGSL( gslStatus = gsl_linalg_LU_decomp(coef, p, &s) );
1519  if ( gslStatus == GSL_SUCCESS )
1520  {
1521  XLAL_CALLGSL( gslStatus = gsl_linalg_LU_solve(coef, p, hderivs, x) );
1522  }
1523  if ( gslStatus != GSL_SUCCESS )
1524  {
1525  gsl_matrix_free(coef);
1526  gsl_vector_free(hderivs);
1527  gsl_vector_free(x);
1528  gsl_permutation_free(p);
1529  XLAL_ERROR( XLAL_EFUNC );
1530  }
1531 
1532  /* Putting solution to an XLAL vector */
1533  modeamps = XLALCreateREAL8Vector(2 * nmodes);
1534 
1535  if ( !modeamps )
1536  {
1537  gsl_matrix_free(coef);
1538  gsl_vector_free(hderivs);
1539  gsl_vector_free(x);
1540  gsl_permutation_free(p);
1541  XLAL_ERROR( XLAL_ENOMEM );
1542  }
1543 
1544  for (i = 0; i < nmodes; ++i)
1545  {
1546  modeamps->data[i] = gsl_vector_get(x, i);
1547  modeamps->data[i + nmodes] = gsl_vector_get(x, i + nmodes);
1548  }
1549 
1550 #if DEBUG_RD
1551  for (i = 0; i < nmodes; ++i)
1552  {
1553  printf("%d: om %12.4e 1/tau %12.4e A %12.4e B %12.4e \n",i,creal(modefreqs->data[i]),cimag(modefreqs->data[i]),modeamps->data[i],modeamps->data[i + nmodes]);
1554  }
1555 #endif
1556 
1557  /* Free all gsl linear algebra objects */
1558  gsl_matrix_free(coef);
1559  gsl_vector_free(hderivs);
1560  gsl_vector_free(x);
1561  gsl_permutation_free(p);
1562 
1563  //double tOffset=(matchrange->data[2]-matchrange->data[1])*m;
1564 
1565  /* Build ring-down waveforms */
1566 
1567  FILE *frd=fopen("checkrdPS.dat","w");
1568  double a1=0.;
1569  double a2=0.;
1570  INT4 jdx;
1571  for (jdx = -5; jdx < 0; ++jdx) {
1572  tj = jdx * dt;
1573  a1=0.;
1574  a2=0.;
1575  for (i = 0; i < nmodes; ++i) {
1576  a1 += exp(- tj * cimag(modefreqs->data[i]))
1577  * ( modeamps->data[i] * cos(tj * creal(modefreqs->data[i]))
1578  + modeamps->data[i + nmodes] * sin(tj * creal(modefreqs->data[i])) );
1579  a2 += exp(- tj * cimag(modefreqs->data[i]))
1580  * (- modeamps->data[i] * sin(tj * creal(modefreqs->data[i]))
1581  + modeamps->data[i + nmodes] * cos(tj * creal(modefreqs->data[i])) );
1582  }
1583  fprintf(frd," %d %12.4e %12.4e %12.4e\n",jdx,matchrange->data[1]*m+tj,.631*a1,.631*a2);
1584  }
1585  for (j = 0; j < rdwave1->length; ++j) {
1586  tj = j * dt;
1587  rdwave1->data[j] = 0;
1588  rdwave2->data[j] = 0;
1589  for (i = 0; i < nmodes; ++i) {
1590  rdwave1->data[j] += exp(- tj * cimag(modefreqs->data[i]))
1591  * ( modeamps->data[i] * cos(tj * creal(modefreqs->data[i]))
1592  + modeamps->data[i + nmodes] * sin(tj * creal(modefreqs->data[i])) );
1593  rdwave2->data[j] += exp(- tj * cimag(modefreqs->data[i]))
1594  * (- modeamps->data[i] * sin(tj * creal(modefreqs->data[i]))
1595  + modeamps->data[i + nmodes] * cos(tj * creal(modefreqs->data[i])) );
1596  }
1597  if (j<20) fprintf(frd," %d %12.4e %12.4e %12.4e\n",j,matchrange->data[1]*m+tj,.631*rdwave1->data[j],.631*rdwave2->data[j]);
1598  }
1599  fclose(frd);
1600 
1601  XLALDestroyREAL8Vector(modeamps);
1602  return errcode;
1603 }
1604 
1605 static INT4 XLALUpSampling(REAL8Vector* vHi, REAL8 dtHi, REAL8Vector* v, REAL8 dt)
1606 {
1607  UINT4 idx;
1608  gsl_interp_accel *acc;
1609  gsl_spline *spline;
1610  INT4 gslStatus;
1611 
1612  double *x = (double *) LALMalloc(v->length * sizeof(double));
1613  double *xHi = (double *) LALMalloc(vHi->length * sizeof(double));
1614 
1615  if ( !x || !xHi ) {
1616  XLALPrintError("** LALSimIMRPSpinInspiralRD ERROR **: allocation failed in interpolation routine\n");
1617  XLAL_ERROR( XLAL_ENOMEM );
1618  }
1619 
1620  for (idx = 0; idx < v->length; idx++) x[idx] = idx*dt;
1621  for (idx = 0; idx < vHi->length; idx++) xHi[idx] = idx*dtHi;
1622  /* printf("First point %12.4e %12.4e\n",x[0],xHi[0]);
1623  printf("Last point %12.4e %12.4e\n",x[v->length-1],xHi[vHi->length-1]);
1624  printf(" dt %12.4e %d, dt %12.4e %d\n",dt,v->length,dtHi,vHi->length);*/
1625 
1626  XLAL_CALLGSL( acc = (gsl_interp_accel*) gsl_interp_accel_alloc() );
1627  XLAL_CALLGSL( spline = (gsl_spline*) gsl_spline_alloc(gsl_interp_cspline, v->length) );
1628  if ( !acc || !spline )
1629  {
1630  if ( acc ) gsl_interp_accel_free(acc);
1631  if ( spline ) gsl_spline_free(spline);
1632  LALFree( x );
1633  XLAL_ERROR( XLAL_ENOMEM );
1634  }
1635 
1636  /* Gall gsl spline interpolation */
1637  XLAL_CALLGSL( gslStatus = gsl_spline_init(spline, x, v->data, v->length) );
1638  if ( gslStatus != GSL_SUCCESS )
1639  {
1640  gsl_spline_free(spline);
1641  gsl_interp_accel_free(acc);
1642  LALFree( x );
1643  XLAL_ERROR( XLAL_EFUNC );
1644  }
1645 
1646  /* Getting first and second order time derivatives from gsl interpolations */
1647  for (idx = 0; idx < vHi->length; idx++) {
1648  vHi->data[idx]=gsl_spline_eval(spline, xHi[idx], acc);
1649  // printf(" vHi[%d] %12.4e\n",idx,vHi->data[idx]);
1650  }
1651 
1652  /* Free gsl variables */
1653  gsl_spline_free(spline);
1654  gsl_interp_accel_free(acc);
1655  LALFree(x);
1656  LALFree(xHi);
1657 
1658  return 0;
1659 }
1660 
1661 
1662 /**
1663  * @addtogroup LALSimIMRPSpinInspiralRD_c
1664  * @{
1665  */
1666 
1667 /**
1668  * Driver routine to compute the PhenSpin Inspiral waveform
1669  * without ring-down
1670  *
1671  * All units are SI units.
1672  */
1673 INT4 XLALSimSpinInspiralGenerator(REAL8TimeSeries **hPlus, /**< +-polarization waveform [returned] */
1674  REAL8TimeSeries **hCross, /**< x-polarization waveform [returned] */
1675  REAL8 phi_start, /**< start phase */
1676  REAL8 deltaT, /**< sampling interval */
1677  REAL8 m1, /**< mass of companion 1 */
1678  REAL8 m2, /**< mass of companion 2 */
1679  REAL8 f_start, /**< start frequency */
1680  REAL8 f_ref, /**< reference frequency */
1681  REAL8 r, /**< distance of source */
1682  REAL8 iota, /**< incination of source (rad) */
1683  REAL8 s1x, /**< x-component of dimensionless spin for object 1 */
1684  REAL8 s1y, /**< y-component of dimensionless spin for object 1 */
1685  REAL8 s1z, /**< z-component of dimensionless spin for object 1 */
1686  REAL8 s2x, /**< x-component of dimensionless spin for object 2 */
1687  REAL8 s2y, /**< y-component of dimensionless spin for object 2 */
1688  REAL8 s2z, /**< z-component of dimensionless spin for object 2 */
1689  INT4 phaseO, /**< twice post-Newtonian phase order */
1690  INT4 UNUSED ampO, /**< twice post-Newtonian amplitude order */
1691  REAL8 lambda1, /**< Tidal par1*/
1692  REAL8 lambda2, /**< Tidal par2*/
1693  REAL8 quadparam1, /**< Quad-monopole par1*/
1694  REAL8 quadparam2, /**< Quad-monopole par1*/
1695  LALDict *LALparams /**< Extra parameters */ )
1696 {
1697 
1698  INT4 errcode=0;
1699  INT4 errcodeInt=0;
1700  INT4 intLen; /* Length of arrays after integration*/
1701  INT4 lengthH;
1702  INT4 idx,kdx;
1703  LALSimInspiralPhenSpinTaylorT4Coeffs params;
1704  REAL8 mass1=m1/LAL_MSUN_SI;
1705  REAL8 mass2=m2/LAL_MSUN_SI;
1706  REAL8 phiN=0.;
1707 
1708  REAL8 *yinitOut=NULL;
1709  REAL8 yinit[LAL_NUM_PST4_VARIABLES];
1710  yinit[0] = phi_start;
1711  yinit[1] = 0.;
1712  yinit[2] = 0.;
1713  yinit[3] = 0.;
1714  yinit[4] = cos(iota);
1715  yinit[5] = s1x;
1716  yinit[6] = s1y;
1717  yinit[7] = s1z;
1718  yinit[8] = s2x;
1719  yinit[9] = s2y;
1720  yinit[10]= s2z;
1721  yinit[11]= 0.;
1722 
1723  REAL8 tn = XLALSimInspiralTaylorLength(deltaT, m1, m2, f_start, phaseO);
1724  REAL8 x = 1.1 * (tn + 1. ) / deltaT;
1725  INT4 length = ceil(log10(x)/log10(2.));
1726  lengthH = pow(2, length);
1727  REAL8TimeSeries *omega=NULL;
1728  REAL8TimeSeries *Phi=NULL;
1729  REAL8TimeSeries *LNhatx=NULL;
1730  REAL8TimeSeries *LNhaty=NULL;
1731  REAL8TimeSeries *LNhatz=NULL;
1732  REAL8TimeSeries *S1x=NULL;
1733  REAL8TimeSeries *S1y=NULL;
1734  REAL8TimeSeries *S1z=NULL;
1735  REAL8TimeSeries *S2x=NULL;
1736  REAL8TimeSeries *S2y=NULL;
1737  REAL8TimeSeries *S2z=NULL;
1738  REAL8TimeSeries *Energy=NULL;
1739 
1740  REAL8 iotaTmp=iota;
1741  if (f_ref<=f_start) {
1742  if (XLALSimIMRPhenSpinInitialize(mass1,mass2,lambda1,lambda2,quadparam1,quadparam2,yinit,f_start,-1.,deltaT,phaseO,&params,LALparams,XLALGetApproximantFromString("PhenSpinTaylor")))
1743  XLAL_ERROR(XLAL_EFUNC);
1744  if(XLALSimIMRPhenSpinInspiralSetAxis(&yinitOut,&iota,&phiN,yinit,iotaTmp,mass1,mass2,XLALSimInspiralWaveformParamsLookupFrameAxis(LALparams)))
1745  XLAL_ERROR(XLAL_EFUNC);
1746  for (idx=0;idx<LAL_NUM_PST4_VARIABLES;idx++)
1747  yinit[idx] = yinitOut[idx];
1748  errcodeInt=XLALSimInspiralSpinTaylorT4Engine(&omega,&Phi,&LNhatx,&LNhaty,&LNhatz,&S1x,&S1y,&S1z,&S2x,&S2y,&S2z,&Energy,yinit,lengthH,PhenSpinTaylor,&params);
1749  intLen=Phi->data->length;
1750  }
1751  else {
1752  REAL8TimeSeries *Phi1,*omega1,*LNhatx1,*LNhaty1,*LNhatz1,*S1x1,*S1y1,*S1z1,*S2x1,*S2y1,*S2z1,*Energy1;
1753  if (XLALSimIMRPhenSpinInitialize(mass1,mass2,lambda1,lambda2,quadparam1,quadparam2,yinit,f_ref,f_start,deltaT,phaseO,&params,LALparams,XLALGetApproximantFromString("PhenSpinTaylor")))
1754  XLAL_ERROR(XLAL_EFUNC);
1755  if(XLALSimIMRPhenSpinInspiralSetAxis(&yinitOut,&iota,&phiN,yinit,iotaTmp,mass1,mass2,XLALSimInspiralWaveformParamsLookupFrameAxis(LALparams)))
1756  XLAL_ERROR(XLAL_EFUNC);
1757  for (idx=0;idx<LAL_NUM_PST4_VARIABLES;idx++)
1758  yinit[idx] = yinitOut[idx];
1759  REAL8 dyTmp[LAL_NUM_PST4_VARIABLES];
1760  REAL8 energy;
1761  XLALSpinInspiralDerivatives(0., yinit,dyTmp,&params);
1762  energy=dyTmp[11]*params.dt/params.M+yinit[11];
1763  yinit[11]=energy;
1764 
1765  errcodeInt=XLALSimInspiralSpinTaylorT4Engine(&omega1,&Phi1,&LNhatx1,&LNhaty1,&LNhatz1,&S1x1,&S1y1,&S1z1,&S2x1,&S2y1,&S2z1,&Energy1,yinit,lengthH,PhenSpinTaylor,&params);
1766 
1767  INT4 intLen1=Phi1->data->length;
1768  /* report on abnormal termination*/
1769  if ( (errcodeInt != LALSIMINSPIRAL_PST4_TEST_FREQBOUND ) ) {
1770  XLALPrintError("** LALSimIMRPSpinInspiralRD WARNING **: integration terminated with code %d.\n",errcode);
1771  XLALPrintError(" 1025: Energy increases\n 1026: Omegadot -ve\n 1028: Omega NAN\n 1029: Freqbound\n 1030: Omega -ve\n");
1772  XLALPrintError(" Waveform parameters were m1 = %14.6e, m2 = %14.6e, inc = %10.6f, fref %10.4f Hz\n", m1, m2, iota, f_ref);
1773  XLALPrintError(" S1 = (%10.6f,%10.6f,%10.6f)\n", s1x, s1y, s1z);
1774  XLALPrintError(" S2 = (%10.6f,%10.6f,%10.6f)\n", s2x, s2y, s2z);
1775  }
1776 
1777  yinit[0] = Phi1->data->data[intLen1-1];
1778  yinit[1] = omega1->data->data[intLen1-1];
1779  yinit[2] = LNhatx1->data->data[intLen1-1];
1780  yinit[3] = LNhaty1->data->data[intLen1-1];
1781  yinit[4] = LNhatz1->data->data[intLen1-1];
1782  yinit[5] = S1x1->data->data[intLen1-1];
1783  yinit[6] = S1y1->data->data[intLen1-1];
1784  yinit[7] = S1z1->data->data[intLen1-1];
1785  yinit[8] = S2x1->data->data[intLen1-1];
1786  yinit[9] = S2y1->data->data[intLen1-1];
1787  yinit[10]= S2z1->data->data[intLen1-1];
1788  yinit[11]= Energy1->data->data[intLen1-1];
1789 
1790  REAL8TimeSeries *omega2,*Phi2,*LNhatx2,*LNhaty2,*LNhatz2,*S1x2,*S1y2,*S1z2,*S2x2,*S2y2,*S2z2,*Energy2;
1791 
1792  params.fEnd=-1.;
1793  params.dt*=-1.;
1794  errcodeInt=XLALSimInspiralSpinTaylorT4Engine(&omega2,&Phi2,&LNhatx2,&LNhaty2,&LNhatz2,&S1x2,&S1y2,&S1z2,&S2x2,&S2y2,&S2z2,&Energy2,yinit,lengthH,PhenSpinTaylor,&params);
1795 
1796  REAL8 phiRef=Phi1->data->data[Phi1->data->length-1];
1797 
1798  omega =appendTSandFree(omega1,omega2);
1799  Phi =appendTSandFree(Phi1,Phi2);
1800  LNhatx=appendTSandFree(LNhatx1,LNhatx2);
1801  LNhaty=appendTSandFree(LNhaty1,LNhaty2);
1802  LNhatz=appendTSandFree(LNhatz1,LNhatz2);
1803  S1x =appendTSandFree(S1x1,S1x2);
1804  S1y =appendTSandFree(S1y1,S1y2);
1805  S1z =appendTSandFree(S1z1,S1z2);
1806  S2x =appendTSandFree(S2x1,S2x2);
1807  S2y =appendTSandFree(S2y1,S2y2);
1808  S2z =appendTSandFree(S2z1,S2z2);
1809  Energy=appendTSandFree(Energy1,Energy2);
1810  intLen=Phi->data->length;
1811  for (idx=0;idx<intLen;idx++) Phi->data->data[idx]-=phiRef;
1812 
1813  }
1814 
1815  /* report on abnormal termination*/
1816  if ( (errcodeInt != LALSIMINSPIRAL_PST4_TEST_ENERGY) ) {
1817  XLALPrintWarning("** LALSimIMRPSpinInspiralRD WARNING **: integration terminated with code %d.\n",errcode);
1818  XLALPrintWarning(" Waveform parameters were m1 = %14.6e, m2 = %14.6e, inc = %10.6f,\n", m1, m2, iota);
1819  XLALPrintWarning(" S1 = (%10.6f,%10.6f,%10.6f)\n", s1x, s1y, s1z);
1820  XLALPrintWarning(" S2 = (%10.6f,%10.6f,%10.6f)\n", s2x, s2y, s2z);
1821  }
1822 
1823  LIGOTimeGPS tStart=LIGOTIMEGPSZERO;
1824  COMPLEX16Vector* hL2tmp=XLALCreateCOMPLEX16Vector(5);
1825  COMPLEX16Vector* hL3tmp=XLALCreateCOMPLEX16Vector(7);
1826  COMPLEX16Vector* hL4tmp=XLALCreateCOMPLEX16Vector(9);
1827  COMPLEX16TimeSeries* hL2=XLALCreateCOMPLEX16TimeSeries( "hL2", &tStart, 0., deltaT, &lalDimensionlessUnit, 5*intLen);
1828  COMPLEX16TimeSeries* hL3=XLALCreateCOMPLEX16TimeSeries( "hL3", &tStart, 0., deltaT, &lalDimensionlessUnit, 7*intLen);
1829  COMPLEX16TimeSeries* hL4=XLALCreateCOMPLEX16TimeSeries( "hL4", &tStart, 0., deltaT, &lalDimensionlessUnit, 9*intLen);
1830  for (idx=0;idx<(int)hL2->data->length;idx++) hL2->data->data[idx]=0.;
1831  for (idx=0;idx<(int)hL3->data->length;idx++) hL3->data->data[idx]=0.;
1832  for (idx=0;idx<(int)hL4->data->length;idx++) hL4->data->data[idx]=0.;
1833 
1834  REAL8TimeSeries *hPtmp=XLALCreateREAL8TimeSeries( "hPtmp", &tStart, 0., deltaT, &lalDimensionlessUnit, intLen);
1835  REAL8TimeSeries *hCtmp=XLALCreateREAL8TimeSeries( "hCtmp", &tStart, 0., deltaT, &lalDimensionlessUnit, intLen);
1836  COMPLEX16TimeSeries *hLMtmp=XLALCreateCOMPLEX16TimeSeries( "hLMtmp", &tStart, 0., deltaT, &lalDimensionlessUnit, intLen);
1837  for (idx=0;idx<(int)hPtmp->data->length;idx++) {
1838  hPtmp->data->data[idx]=0.;
1839  hCtmp->data->data[idx]=0.;
1840  hLMtmp->data->data[idx]=0.;
1841  }
1842 
1843  LALSimInspiralInclAngle trigAngle;
1844 
1845  REAL8 amp22ini = -2.0 * m1*m2/(m1+m2) * LAL_G_SI/pow(LAL_C_SI,2.) / r * sqrt(16. * LAL_PI / 5.);
1846  REAL8 amp33ini = -amp22ini * sqrt(5./42.)/4.;
1847  REAL8 amp44ini = amp22ini * sqrt(5./7.) * 2./9.;
1848  REAL8 alpha,v,v2,Psi,om;
1849  REAL8 eta=mass1*mass2/(mass1+mass2)/(mass1+mass2);
1850  REAL8 dm=(mass1-mass2)/(mass1+mass2);
1851 
1852  for (idx=0;idx<intLen;idx++) {
1853  om=omega->data->data[idx];
1854  v=cbrt(om);
1855  v2=v*v;
1856  Psi=Phi->data->data[idx] -2.*om*(1.-eta*v2)*log(om);
1857  errcode =XLALSimInspiralComputeAlpha(params,LNhatx->data->data[idx],LNhaty->data->data[idx],S1x->data->data[idx],S1y->data->data[idx],S2x->data->data[idx],S2y->data->data[idx],&alpha);
1858 
1859  errcode+=XLALSimInspiralComputeInclAngle(LNhatz->data->data[idx],&trigAngle);
1860  errcode+=XLALSimSpinInspiralFillL2Modes(hL2tmp,v,eta,dm,Psi,alpha,&trigAngle);
1861  for (kdx=0;kdx<5;kdx++) hL2->data->data[5*idx+kdx]=hL2tmp->data[kdx]*amp22ini*v2;
1862  errcode+=XLALSimSpinInspiralFillL3Modes(hL3tmp,v,eta,dm,Psi,alpha,&trigAngle);
1863  for (kdx=0;kdx<7;kdx++) hL3->data->data[7*idx+kdx]=hL3tmp->data[kdx]*amp33ini*v2;
1864  errcode+=XLALSimSpinInspiralFillL4Modes(hL4tmp,v,eta,dm,Psi,alpha,&trigAngle);
1865  for (kdx=0;kdx<9;kdx++) hL4->data->data[9*idx+kdx]=hL4tmp->data[kdx]*amp44ini*v2*v2;
1866  }
1867  XLALDestroyCOMPLEX16Vector(hL2tmp);
1868  XLALDestroyCOMPLEX16Vector(hL3tmp);
1869  XLALDestroyCOMPLEX16Vector(hL4tmp);
1870 
1871  XLALDestroyREAL8TimeSeries(omega);
1872  XLALDestroyREAL8TimeSeries(Phi);
1873  XLALDestroyREAL8TimeSeries(LNhatx);
1874  XLALDestroyREAL8TimeSeries(LNhaty);
1875  XLALDestroyREAL8TimeSeries(LNhatz);
1876  XLALDestroyREAL8TimeSeries(S1x);
1877  XLALDestroyREAL8TimeSeries(S1y);
1878  XLALDestroyREAL8TimeSeries(S1z);
1879  XLALDestroyREAL8TimeSeries(S2x);
1880  XLALDestroyREAL8TimeSeries(S2y);
1881  XLALDestroyREAL8TimeSeries(S2z);
1882  XLALDestroyREAL8TimeSeries(Energy);
1883 
1884  INT4 m,l;
1885  int modesChoice=XLALSimInspiralWaveformParamsLookupModesChoice(LALparams);
1886  if ( ( modesChoice & LAL_SIM_INSPIRAL_MODES_CHOICE_RESTRICTED) == LAL_SIM_INSPIRAL_MODES_CHOICE_RESTRICTED ) {
1887  l=2;
1888  for (m=-l;m<=l;m++) {
1889  for (idx=0;idx<intLen;idx++) hLMtmp->data->data[idx]=hL2->data->data[(m+l)+idx*(2*l+1)];
1890  XLALSimAddMode(hPtmp,hCtmp,hLMtmp,iota,phiN,l,m,0);
1891  }
1892  }
1893  XLALDestroyCOMPLEX16TimeSeries(hL2);
1894  if ( ( modesChoice & LAL_SIM_INSPIRAL_MODES_CHOICE_3L) == LAL_SIM_INSPIRAL_MODES_CHOICE_3L ) {
1895  l=3;
1896  for (m=-l;m<=l;m++) {
1897  for (idx=0;idx<intLen;idx++)
1898  hLMtmp->data->data[idx]=hL3->data->data[(m+l)+idx*(2*l+1)];
1899  XLALSimAddMode(hPtmp,hCtmp,hLMtmp,iota,phiN,l,m,0);
1900  }
1901  }
1902  XLALDestroyCOMPLEX16TimeSeries(hL3);
1903  if ( ( modesChoice & LAL_SIM_INSPIRAL_MODES_CHOICE_3L) == LAL_SIM_INSPIRAL_MODES_CHOICE_3L ) {
1904  l=4;
1905  for (m=-l;m<=l;m++) {
1906  for (idx=0;idx<intLen;idx++)
1907  hLMtmp->data->data[idx]=hL4->data->data[(m+l)+idx*(2*l+1)];
1908  XLALSimAddMode(hPtmp,hCtmp,hLMtmp,iota,phiN,l,m,0);
1909  }
1910  }
1911  XLALDestroyCOMPLEX16TimeSeries(hL4);
1912  XLALDestroyCOMPLEX16TimeSeries(hLMtmp);
1913 
1914  REAL8 tPeak=intLen*deltaT;
1915  if ((*hPlus) && (*hCross)) {
1916  if ((*hPlus)->data->length!=(*hCross)->data->length) {
1917  XLALPrintError("*** LALSimIMRPSpinInspiralRD ERROR: h+ and hx differ in length: %d vs. %d\n",(*hPlus)->data->length,(*hCross)->data->length);
1918  XLAL_ERROR(XLAL_EFAILED);
1919  }
1920  else {
1921  if ((int)(*hPlus)->data->length<intLen) {
1922  XLALPrintError("*** LALSimIMRPSpinInspiralRD ERROR: h+ and hx too short: %d vs. %d\n",(*hPlus)->data->length,intLen);
1923  XLAL_ERROR(XLAL_EFAILED);
1924  }
1925  else {
1926  XLALGPSAdd(&((*hPlus)->epoch),-tPeak);
1927  XLALGPSAdd(&((*hCross)->epoch),-tPeak);
1928  }
1929  }
1930  }
1931  else {
1932  XLALGPSAdd(&tStart,-tPeak);
1933  *hPlus = XLALCreateREAL8TimeSeries("H+", &tStart, 0.0, deltaT, &lalDimensionlessUnit, intLen);
1934  *hCross = XLALCreateREAL8TimeSeries("Hx", &tStart, 0.0, deltaT, &lalDimensionlessUnit, intLen);
1935  if(*hPlus == NULL || *hCross == NULL)
1936  XLAL_ERROR(XLAL_ENOMEM);
1937  }
1938 
1939  INT4 minLen=hPtmp->data->length < (*hPlus)->data->length ? hPtmp->data->length : (*hPlus)->data->length;
1940  for (idx=0;idx<minLen;idx++) {
1941  (*hPlus)->data->data[idx] =hPtmp->data->data[idx];
1942  (*hCross)->data->data[idx]=hCtmp->data->data[idx];
1943  }
1944  for (idx=minLen;idx<(int)(*hPlus)->data->length;idx++) {
1945  (*hPlus)->data->data[idx] =0.;
1946  (*hCross)->data->data[idx]=0.;
1947  }
1948 
1949  XLALDestroyREAL8TimeSeries(hPtmp);
1950  XLALDestroyREAL8TimeSeries(hCtmp);
1951 
1952  return errcode;
1953 
1954 } /* End of XLALSimSpinInspiralGenerator*/
1955 
1956 INT4 XLALSimIMRPhenSpinFinalMassSpin(REAL8 *finalMass,
1957  REAL8 *finalSpin,
1958  REAL8 mass1,
1959  REAL8 mass2,
1960  REAL8 s1s1,
1961  REAL8 s2s2,
1962  REAL8 s1L,
1963  REAL8 s2L,
1964  REAL8 s1s2,
1965  REAL8 energy)
1966 {
1967  /* XLAL error handling */
1968  INT4 errcode = XLAL_SUCCESS;
1969  REAL8 qq,ll,eta;
1970 
1971  /* See eq.(6) in arXiv:0904.2577 */
1972  REAL8 ma1,ma2,a12,a12l;
1973  REAL8 cosa1=0.;
1974  REAL8 cosa2=0.;
1975  REAL8 cosa12=0.;
1976 
1977  REAL8 t0=-2.9;
1978  REAL8 t3=2.6;
1979  REAL8 s4=-0.123;
1980  REAL8 s5=0.45;
1981  REAL8 t2=16.*(0.6865-t3/64.-sqrt(3.)/2.);
1982 
1983  /* get a local copy of the intrinstic parameters */
1984  qq = mass2/mass1;
1985  eta = mass1*mass2/((mass1+mass2)*(mass1+mass2));
1986  /* done */
1987  ma1 = sqrt( s1s1 );
1988  ma2 = sqrt( s2s2 );
1989 
1990  if (ma1>0.) cosa1 = s1L/ma1;
1991  else cosa1=0.;
1992  if (ma2>0.) cosa2 = s2L/ma2;
1993  else cosa2=0.;
1994  if ((ma1>0.)&&(ma2>0.)) {
1995  cosa12 = s1s2/ma1/ma2;
1996  }
1997  else cosa12=0.;
1998 
1999  a12 = ma1*ma1 + ma2*ma2*qq*qq*qq*qq + 2.*ma1*ma2*qq*qq*cosa12 ;
2000  a12l = ma1*cosa1 + ma2*cosa2*qq*qq ;
2001  ll = 2.*sqrt(3.)+ t2*eta + t3*eta*eta + s4*a12/(1.+qq*qq)/(1.+qq*qq) + (s5*eta+t0+2.)/(1.+qq*qq)*a12l;
2002 
2003  /* Estimate final mass by adding the negative binding energy to the rest mass*/
2004  *finalMass = 1. + energy;
2005 
2006  /* Estimate final spin */
2007  *finalSpin = sqrt( a12 + 2.*ll*qq*a12l + ll*ll*qq*qq)/(1.+qq)/(1.+qq);
2008 
2009  /* Check value of finalMass */
2010  if (*finalMass < 0.) {
2011  XLALPrintWarning("*** LALSimIMRPSpinInspiralRD ERROR: Estimated final mass <0 : %12.6f\n ",*finalMass);
2012  XLAL_ERROR( XLAL_ERANGE);
2013  }
2014 
2015  /* Check value of finalSpin */
2016  if ((*finalSpin > 1.)||(*finalSpin < 0.)) {
2017  if ((*finalSpin>=1.)&&(*finalSpin<1.01)) {
2018  XLALPrintWarning("*** LALSimIMRPSpinInspiralRD WARNING: Estimated final Spin slightly >1 : %11.3e\n ",*finalSpin);
2019  XLALPrintWarning(" (m1=%8.3f m2=%8.3f s1sq=%8.3f s2sq=%8.3f s1L=%8.3f s2L=%8.3f s1s2=%8.3f ) final spin set to 1 and code goes on\n",mass1,mass2,s1s1,s2s2,s1L,s2L,s1s2);
2020  *finalSpin = .99999;
2021  }
2022  else {
2023  XLALPrintError("*** LALSimIMRPSpinInspiralRD ERROR: Unphysical estimation of final Spin : %11.3e\n ",*finalSpin);
2024  XLALPrintWarning(" (m1=%8.3f m2=%8.3f s1sq=%8.3f s2sq=%8.3f s1L=%8.3f s2L=%8.3f s1s2=%8.3f )\n",mass1,mass2,s1s1,s2s2,s1L,s2L,s1s2);
2025  XLALPrintError("*** Code aborts\n");
2026  XLAL_ERROR( XLAL_ERANGE);
2027  }
2028  }
2029 
2030  return errcode;
2031 } /* End of XLALSimIMRPhenSpinFinalMassSpin*/
2032 
2033 /**
2034  * Driver routine for generating PhenSpinRD waveforms
2035  */
2036 INT4 XLALSimIMRPhenSpinInspiralRDGenerator(REAL8TimeSeries **hPlus, /**< +-polarization waveform [returned] */
2037  REAL8TimeSeries **hCross, /**< x-polarization waveform [returned] */
2038  REAL8 phi_start, /**< start phase */
2039  REAL8 deltaT, /**< sampling interval */
2040  REAL8 m1, /**< mass of companion 1 in SI units */
2041  REAL8 m2, /**< mass of companion 2 in SI units */
2042  REAL8 f_start, /**< start frequency */
2043  REAL8 f_ref, /**< reference frequency */
2044  REAL8 r, /**< distance of source */
2045  REAL8 iota, /**< inclination of source (rad) */
2046  REAL8 s1x, /**< x-component of dimensionless spin for object 1 */
2047  REAL8 s1y, /**< y-component of dimensionless spin for object 1 */
2048  REAL8 s1z, /**< z-component of dimensionless spin for object 1 */
2049  REAL8 s2x, /**< x-component of dimensionless spin for object 2 */
2050  REAL8 s2y, /**< y-component of dimensionless spin for object 2 */
2051  REAL8 s2z, /**< z-component of dimensionless spin for object 2 */
2052  INT4 phaseO, /**< twice post-Newtonian phase order */
2053  INT4 UNUSED ampO, /**< twice post-Newtonian amplitude order */
2054  REAL8 lambda1, /**< Tidal par1*/
2055  REAL8 lambda2, /**< Tidal par2*/
2056  REAL8 quadparam1, /**< Quad-monopole par1*/
2057  REAL8 quadparam2, /**< Quad-monopole par2*/
2058  LALDict *LALparams /**< Extra parameters */ )
2059 {
2060 
2061  const double rateHi=16384;
2062  if (1./deltaT>rateHi) {
2063  XLALPrintError("** LALSimIMRPSpinInspiralRD ERROR **: rate must be smaller than %8.0f Hz, value %8.0f Hz given\n",rateHi,1./deltaT);
2064  XLAL_ERROR(XLAL_EFUNC);
2065  }
2066  INT4 errcode=0;
2067  INT4 errcodeInt=0;
2068  UINT4 lengthH=0; /* Length of hPlus and hCross passed, 0 if NULL*/
2069  UINT4 intLen; /* Length of arrays after integration*/
2070  INT4 idx,jdx,kdx;
2071  LALSimInspiralPhenSpinTaylorT4Coeffs params;
2072  REAL8 S1S1=s1x*s1x+s1y*s1y+s1z*s1z;
2073  REAL8 S2S2=s1x*s1x+s1y*s1y+s1z*s1z;
2074  REAL8 mass1=m1/LAL_MSUN_SI;
2075  REAL8 mass2=m2/LAL_MSUN_SI;
2076  REAL8 Mass=mass1+mass2;
2077  REAL8 Mtime=Mass*LAL_MTSUN_SI;
2078  REAL8 phiN=0.;
2079 
2080  REAL8 *yinitOut=NULL;
2081  REAL8 yinit[LAL_NUM_PST4_VARIABLES];
2082  yinit[0] = phi_start;
2083  yinit[1] = 0.;
2084  yinit[2] = 0.;
2085  yinit[3] = 0.;
2086  yinit[4] = cos(iota);
2087  yinit[5] = s1x;
2088  yinit[6] = s1y;
2089  yinit[7] = s1z;
2090  yinit[8] = s2x;
2091  yinit[9] = s2y;
2092  yinit[10]= s2z;
2093  yinit[11]= 0.;
2094 
2095  REAL8TimeSeries *omega, *Phi, *LNhatx, *LNhaty, *LNhatz;
2096  REAL8TimeSeries *S1x, *S1y, *S1z, *S2x, *S2y, *S2z, *Energy;
2097 
2098  REAL8 tn = XLALSimInspiralTaylorLength(deltaT, m1, m2, f_start, phaseO);
2099  REAL8 x = 1.1 * (tn + 1. ) / deltaT;
2100  INT4 length = ceil(log10(x)/log10(2.));
2101  lengthH = pow(2, length);
2102 
2103  REAL8 iotaTmp=iota;
2104  if (f_ref<=f_start) {
2105  if (XLALSimIMRPhenSpinInitialize(mass1,mass2,lambda1,lambda2,quadparam1,quadparam2,yinit,f_start,-1.,deltaT,phaseO,&params,LALparams,XLALGetApproximantFromString("PhenSpinTaylorRD")))
2106  XLAL_ERROR(XLAL_EFUNC);
2107  if(XLALSimIMRPhenSpinInspiralSetAxis(&yinitOut,&iota,&phiN,yinit,iotaTmp,mass1,mass2,XLALSimInspiralWaveformParamsLookupFrameAxis(LALparams)))
2108  XLAL_ERROR(XLAL_EFUNC);
2109  for (idx=0;idx<LAL_NUM_PST4_VARIABLES;idx++)
2110  yinit[idx]=yinitOut[idx];
2111  errcodeInt=XLALSimInspiralSpinTaylorT4Engine(&omega,&Phi,&LNhatx,&LNhaty,&LNhatz,&S1x,&S1y,&S1z,&S2x,&S2y,&S2z,&Energy,yinit,lengthH,PhenSpinTaylorRD,&params);
2112  intLen=Phi->data->length;
2113  }
2114  else /* do both forward and backward integration*/ {
2115  REAL8TimeSeries *Phi1, *omega1, *LNhatx1, *LNhaty1, *LNhatz1;
2116  REAL8TimeSeries *S1x1, *S1y1, *S1z1, *S2x1, *S2y1, *S2z1, *Energy1;
2117  if (XLALSimIMRPhenSpinInitialize(mass1,mass2,lambda1,lambda2,quadparam1,quadparam2,yinit,f_ref,f_start,deltaT,phaseO,&params,LALparams,XLALGetApproximantFromString("PhenSpinTaylorRD")))
2118  XLAL_ERROR(XLAL_EFUNC);
2119  if(XLALSimIMRPhenSpinInspiralSetAxis(&yinitOut,&iota,&phiN,yinit,iotaTmp,mass1,mass2,XLALSimInspiralWaveformParamsLookupFrameAxis(LALparams)))
2120  XLAL_ERROR(XLAL_EFUNC);
2121  for (idx=0;idx<LAL_NUM_PST4_VARIABLES;idx++)
2122  yinit[idx]=yinitOut[idx];
2123  REAL8 dyTmp[LAL_NUM_PST4_VARIABLES];
2124  REAL8 energy;
2125  XLALSpinInspiralDerivatives(0., yinit,dyTmp,&params);
2126  energy=dyTmp[11]*params.dt/params.M+yinit[11];
2127  yinit[11]=energy;
2128 
2129  errcode=XLALSimInspiralSpinTaylorT4Engine(&omega1,&Phi1,&LNhatx1,&LNhaty1,&LNhatz1,&S1x1,&S1y1,&S1z1,&S2x1,&S2y1,&S2z1,&Energy1,yinit,lengthH,PhenSpinTaylorRD,&params);
2130  /* report on abnormal termination*/
2131  if ( errcode != LALSIMINSPIRAL_PST4_TEST_FREQBOUND ) {
2132  XLALPrintError("** LALSimIMRPSpinInspiralRD WARNING **: integration terminated with code %d.\n",errcode);
2133  XLALPrintError(" 1025: Energy increases\n 1026: Omegadot -ve\n 1027: Freqbound\n 1028: Omega NAN\n 1030: Omega -ve\n 1031: Omega > OmegaMatch\n");
2134  XLALPrintError(" Waveform parameters were m1 = %14.6e, m2 = %14.6e, inc = %10.6f, fref %10.4f Hz\n", mass1, mass2, iota, f_ref);
2135  XLALPrintError(" S1 = (%10.6f,%10.6f,%10.6f)\n", s1x, s1y, s1z);
2136  XLALPrintError(" S2 = (%10.6f,%10.6f,%10.6f)\n", s2x, s2y, s2z);
2137  }
2138 
2139  INT4 intLen1=Phi1->data->length;
2140 
2141  yinit[0] = Phi1->data->data[intLen1-1];
2142  yinit[1] = omega1->data->data[intLen1-1];
2143  yinit[2] = LNhatx1->data->data[intLen1-1];
2144  yinit[3] = LNhaty1->data->data[intLen1-1];
2145  yinit[4] = LNhatz1->data->data[intLen1-1];
2146  yinit[5] = S1x1->data->data[intLen1-1];
2147  yinit[6] = S1y1->data->data[intLen1-1];
2148  yinit[7] = S1z1->data->data[intLen1-1];
2149  yinit[8] = S2x1->data->data[intLen1-1];
2150  yinit[9] = S2y1->data->data[intLen1-1];
2151  yinit[10]= S2z1->data->data[intLen1-1];
2152  yinit[11]= Energy1->data->data[intLen1-1];
2153 
2154  REAL8TimeSeries *omega2, *Phi2, *LNhatx2, *LNhaty2, *LNhatz2;
2155  REAL8TimeSeries *S1x2, *S1y2, *S1z2, *S2x2, *S2y2, *S2z2, *Energy2;
2156 
2157  params.fEnd=-1.;
2158  params.dt*=-1.;
2159  errcodeInt=XLALSimInspiralSpinTaylorT4Engine(&omega2,&Phi2,&LNhatx2,&LNhaty2,&LNhatz2,&S1x2,&S1y2,&S1z2,&S2x2,&S2y2,&S2z2,&Energy2,yinit,lengthH,PhenSpinTaylorRD,&params);
2160 
2161  REAL8 phiRef=Phi1->data->data[omega1->data->length-1];
2162  omega =appendTSandFree(omega1,omega2);
2163  Phi =appendTSandFree(Phi1,Phi2);
2164  LNhatx=appendTSandFree(LNhatx1,LNhatx2);
2165  LNhaty=appendTSandFree(LNhaty1,LNhaty2);
2166  LNhatz=appendTSandFree(LNhatz1,LNhatz2);
2167  S1x =appendTSandFree(S1x1,S1x2);
2168  S1y =appendTSandFree(S1y1,S1y2);
2169  S1z =appendTSandFree(S1z1,S1z2);
2170  S2x =appendTSandFree(S2x1,S2x2);
2171  S2y =appendTSandFree(S2y1,S2y2);
2172  S2z =appendTSandFree(S2z1,S2z2);
2173  Energy=appendTSandFree(Energy1,Energy2);
2174  intLen=Phi->data->length;
2175  for (idx=0;idx<(int)intLen;idx++) Phi->data->data[idx]-=phiRef;
2176 
2177  } /* End of int forward + integration backward*/
2178 
2179  /* report on abnormal termination*/
2180  if ( (errcodeInt != LALSIMINSPIRAL_PST4_TEST_OMEGAMATCH) ) {
2181  XLALPrintError("** LALSimIMRPSpinInspiralRD WARNING **: integration terminated with code %d.\n",errcode);
2182  XLALPrintError(" 1025: Energy increases\n 1026: Omegadot -ve\n 1027: Freqbound\n 1028: Omega NAN\n 1030: Omega -ve\n");
2183  XLALPrintError(" Waveform parameters were m1 = %14.6e, m2 = %14.6e, inc = %10.6f,\n", m1, m2, iota);
2184  XLALPrintError(" S1 = (%10.6f,%10.6f,%10.6f)\n", s1x, s1y, s1z);
2185  XLALPrintError(" S2 = (%10.6f,%10.6f,%10.6f)\n", s2x, s2y, s2z);
2186  }
2187 
2188  INT4 count=intLen;
2189  double tPeak=0.,tm=0.;
2190  LIGOTimeGPS tStart=LIGOTIMEGPSZERO;
2191  COMPLEX16Vector* hL2tmp=XLALCreateCOMPLEX16Vector(5);
2192  COMPLEX16Vector* hL3tmp=XLALCreateCOMPLEX16Vector(7);
2193  COMPLEX16Vector* hL4tmp=XLALCreateCOMPLEX16Vector(9);
2194  COMPLEX16TimeSeries* hL2=XLALCreateCOMPLEX16TimeSeries( "hL2", &tStart, 0., deltaT, &lalDimensionlessUnit, 5*lengthH);
2195  COMPLEX16TimeSeries* hL3=XLALCreateCOMPLEX16TimeSeries( "hL3", &tStart, 0., deltaT, &lalDimensionlessUnit, 7*lengthH);
2196  COMPLEX16TimeSeries* hL4=XLALCreateCOMPLEX16TimeSeries( "hL4", &tStart, 0., deltaT, &lalDimensionlessUnit, 9*lengthH);
2197  for (idx=0;idx<(int)hL2->data->length;idx++) hL2->data->data[idx]=0.;
2198  for (idx=0;idx<(int)hL3->data->length;idx++) hL3->data->data[idx]=0.;
2199  for (idx=0;idx<(int)hL4->data->length;idx++) hL4->data->data[idx]=0.;
2200 
2201  REAL8TimeSeries *hPtmp=XLALCreateREAL8TimeSeries( "hPtmp", &tStart, 0., deltaT, &lalDimensionlessUnit, lengthH);
2202  REAL8TimeSeries *hCtmp=XLALCreateREAL8TimeSeries( "hCtmp", &tStart, 0., deltaT, &lalDimensionlessUnit, lengthH);
2203  COMPLEX16TimeSeries *hLMtmp=XLALCreateCOMPLEX16TimeSeries( "hLMtmp", &tStart, 0., deltaT, &lalDimensionlessUnit, lengthH);
2204  for (idx=0;idx<(int)hPtmp->data->length;idx++) {
2205  hPtmp->data->data[idx]=0.;
2206  hCtmp->data->data[idx]=0.;
2207  hLMtmp->data->data[idx]=0.;
2208  }
2209 
2210  LALSimInspiralInclAngle trigAngle;
2211 
2212  REAL8 amp22ini = -2.0 * m1*m2/(m1+m2) * LAL_G_SI/pow(LAL_C_SI,2.) / r * sqrt(16. * LAL_PI / 5.);
2213  REAL8 amp33ini = -amp22ini * sqrt(5./42.)/4.;
2214  REAL8 amp44ini = amp22ini * sqrt(5./7.) * 2./9.;
2215  REAL8 alpha=0.,v=0.,v2=0.,Psi,om;
2216  REAL8 eta=m1*m2/(m1+m2)/(m1+m2);
2217  REAL8 dm=(m1-m2)/(m1+m2);
2218  int modesChoice=XLALSimInspiralWaveformParamsLookupModesChoice(LALparams);
2219 
2220  for (idx=0;idx<(int)intLen;idx++) {
2221  om=omega->data->data[idx];
2222  v=cbrt(om);
2223  v2=v*v;
2224  Psi=Phi->data->data[idx] -2.*om*(1.-eta*v2)*log(om);
2225  errcode =XLALSimInspiralComputeAlpha(params,LNhatx->data->data[idx],LNhaty->data->data[idx],S1x->data->data[idx],S1y->data->data[idx],S2x->data->data[idx],S2y->data->data[idx],&alpha);
2226  errcode+=XLALSimInspiralComputeInclAngle(LNhatz->data->data[idx],&trigAngle);
2227  errcode+=XLALSimSpinInspiralFillL2Modes(hL2tmp,v,eta,dm,Psi,alpha,&trigAngle);
2228  for (kdx=0;kdx<5;kdx++) hL2->data->data[5*idx+kdx]=hL2tmp->data[kdx]*amp22ini*v2;
2229  errcode+=XLALSimSpinInspiralFillL3Modes(hL3tmp,v,eta,dm,Psi,alpha,&trigAngle);
2230  for (kdx=0;kdx<7;kdx++) hL3->data->data[7*idx+kdx]=hL3tmp->data[kdx]*amp33ini*v2;
2231  errcode+=XLALSimSpinInspiralFillL4Modes(hL4tmp,v,eta,dm,Psi,alpha,&trigAngle);
2232  for (kdx=0;kdx<9;kdx++) hL4->data->data[9*idx+kdx]=hL4tmp->data[kdx]*amp44ini*v2*v2;
2233  }
2234 
2235  if (errcodeInt==LALSIMINSPIRAL_PST4_TEST_OMEGAMATCH) {
2236  REAL8 LNhS1,LNhS2,S1S2,omegaMatch;
2237  REAL8 m1Msq=pow(params.m1ByM,2.);
2238  REAL8 m2Msq=pow(params.m2ByM,2.);
2239 
2240  INT4 iMatchUp=0;
2241  INT4 iMatch=0;
2242  const double dtHi=1./rateHi;
2243  INT4 stkLength,stkLenHi;
2244 
2245  REAL8Vector *Phi_s = NULL;
2246  REAL8Vector *omega_s = NULL;
2247  REAL8Vector *LNhx_s = NULL;
2248  REAL8Vector *LNhy_s = NULL;
2249  REAL8Vector *LNhz_s = NULL;
2250  REAL8Vector *S1x_s = NULL;
2251  REAL8Vector *S1y_s = NULL;
2252  REAL8Vector *S1z_s = NULL;
2253  REAL8Vector *S2x_s = NULL;
2254  REAL8Vector *S2y_s = NULL;
2255  REAL8Vector *S2z_s = NULL;
2256  REAL8Vector *Energy_s = NULL;
2257 
2258  REAL8Vector *PhiHi = NULL;
2259  REAL8Vector *omegaHi = NULL;
2260  REAL8Vector *LNhxHi = NULL;
2261  REAL8Vector *LNhyHi = NULL;
2262  REAL8Vector *LNhzHi = NULL;
2263  REAL8Vector *S1xHi = NULL;
2264  REAL8Vector *S1yHi = NULL;
2265  REAL8Vector *S1zHi = NULL;
2266  REAL8Vector *S2xHi = NULL;
2267  REAL8Vector *S2yHi = NULL;
2268  REAL8Vector *S2zHi = NULL;
2269  REAL8Vector *EnergyHi = NULL;
2270 
2271  idx=omega->data->length-2-( (int) ( ((double)minIntLen)*dtHi/deltaT ) );
2272  do {
2273  idx--;
2274  LNhS1=(LNhatx->data->data[idx]*S1x->data->data[idx]+LNhaty->data->data[idx]*S1y->data->data[idx]+LNhatz->data->data[idx]*S1z->data->data[idx])/m1Msq;
2275  LNhS2=(LNhatx->data->data[idx]*S2x->data->data[idx]+LNhaty->data->data[idx]*S2y->data->data[idx]+LNhatz->data->data[idx]*S2z->data->data[idx])/m2Msq;
2276  S1S2=(S1x->data->data[idx]*S2x->data->data[idx]+S1y->data->data[idx]*S2y->data->data[idx]+S1z->data->data[idx]*S2z->data->data[idx])/m1Msq/m2Msq;
2277  omegaMatch=OmMatch(LNhS1,LNhS2,S1S1,S1S2,S2S2);
2278  } while ((idx>0)&&(omega->data->data[abs(idx)]>omegaMatch));
2279  if (idx<0) {
2280  XLALPrintError(" *** XLALSimIMRPSpinInspiralRD ERROR ***: impossible to attach phen part\n");
2281  XLAL_ERROR(XLAL_EFAILED);
2282  }
2283  iMatch=idx;
2284  stkLength=omega->data->length-iMatch;
2285  Phi_s = XLALCreateREAL8Vector(stkLength);
2286  omega_s = XLALCreateREAL8Vector(stkLength);
2287  LNhx_s = XLALCreateREAL8Vector(stkLength);
2288  LNhy_s = XLALCreateREAL8Vector(stkLength);
2289  LNhz_s = XLALCreateREAL8Vector(stkLength);
2290  S1x_s = XLALCreateREAL8Vector(stkLength);
2291  S1y_s = XLALCreateREAL8Vector(stkLength);
2292  S1z_s = XLALCreateREAL8Vector(stkLength);
2293  S2x_s = XLALCreateREAL8Vector(stkLength);
2294  S2y_s = XLALCreateREAL8Vector(stkLength);
2295  S2z_s = XLALCreateREAL8Vector(stkLength);
2296  Energy_s = XLALCreateREAL8Vector(stkLength);
2297  for (jdx=0; jdx < ( ( (int)omega->data->length ) - iMatch); jdx++) {
2298  kdx=jdx+iMatch;
2299  Phi_s->data[jdx] = Phi->data->data[kdx];
2300  omega_s->data[jdx] = omega->data->data[kdx];
2301  LNhx_s->data[jdx] = LNhatx->data->data[kdx];
2302  LNhy_s->data[jdx] = LNhaty->data->data[kdx];
2303  LNhz_s->data[jdx] = LNhatz->data->data[kdx];
2304  S1x_s->data[jdx] = S1x->data->data[kdx];
2305  S1y_s->data[jdx] = S1y->data->data[kdx];
2306  S1z_s->data[jdx] = S1z->data->data[kdx];
2307  S2x_s->data[jdx] = S2x->data->data[kdx];
2308  S2y_s->data[jdx] = S2y->data->data[kdx];
2309  S2z_s->data[jdx] = S2z->data->data[kdx];
2310  Energy_s->data[jdx] = Energy->data->data[kdx];
2311  LNhS1=(LNhatx->data->data[kdx]*S1x->data->data[kdx]+LNhaty->data->data[kdx]*S1y->data->data[kdx]+LNhatz->data->data[kdx]*S1z->data->data[kdx])/m1Msq;
2312  LNhS2=(LNhatx->data->data[kdx]*S2x->data->data[kdx]+LNhaty->data->data[kdx]*S2y->data->data[kdx]+LNhatz->data->data[kdx]*S2z->data->data[kdx])/m2Msq;
2313  S1S2=(S1x->data->data[kdx]*S2x->data->data[kdx]+S1y->data->data[kdx]*S2y->data->data[kdx]+S1z->data->data[kdx]*S2z->data->data[kdx])/m1Msq/m2Msq;
2314  }
2315  XLALDestroyREAL8TimeSeries(Phi);
2316  XLALDestroyREAL8TimeSeries(omega);
2317  XLALDestroyREAL8TimeSeries(LNhatx);
2318  XLALDestroyREAL8TimeSeries(LNhaty);
2319  XLALDestroyREAL8TimeSeries(LNhatz);
2320  XLALDestroyREAL8TimeSeries(S1x);
2321  XLALDestroyREAL8TimeSeries(S1y);
2322  XLALDestroyREAL8TimeSeries(S1z);
2323  XLALDestroyREAL8TimeSeries(S2x);
2324  XLALDestroyREAL8TimeSeries(S2y);
2325  XLALDestroyREAL8TimeSeries(S2z);
2326  XLALDestroyREAL8TimeSeries(Energy);
2327 
2328  stkLenHi=((int) (deltaT/dtHi))*(stkLength-1)+1;
2329  PhiHi = XLALCreateREAL8Vector(stkLenHi);
2330  omegaHi = XLALCreateREAL8Vector(stkLenHi);
2331  REAL8Vector* omMHi = XLALCreateREAL8Vector(stkLenHi);
2332  LNhxHi = XLALCreateREAL8Vector(stkLenHi);
2333  LNhyHi = XLALCreateREAL8Vector(stkLenHi);
2334  LNhzHi = XLALCreateREAL8Vector(stkLenHi);
2335  S1xHi = XLALCreateREAL8Vector(stkLenHi);
2336  S1yHi = XLALCreateREAL8Vector(stkLenHi);
2337  S1zHi = XLALCreateREAL8Vector(stkLenHi);
2338  S2xHi = XLALCreateREAL8Vector(stkLenHi);
2339  S2yHi = XLALCreateREAL8Vector(stkLenHi);
2340  S2zHi = XLALCreateREAL8Vector(stkLenHi);
2341  EnergyHi = XLALCreateREAL8Vector(stkLenHi);
2342 
2343  XLALUpSampling(PhiHi, dtHi, Phi_s, deltaT);
2344  XLALUpSampling(omegaHi, dtHi, omega_s, deltaT);
2345  XLALUpSampling(LNhxHi, dtHi, LNhx_s, deltaT);
2346  XLALUpSampling(LNhyHi, dtHi, LNhy_s, deltaT);
2347  XLALUpSampling(LNhzHi, dtHi, LNhz_s, deltaT);
2348  XLALUpSampling(S1xHi, dtHi, S1x_s, deltaT);
2349  XLALUpSampling(S1yHi, dtHi, S1y_s, deltaT);
2350  XLALUpSampling(S1zHi, dtHi, S1z_s, deltaT);
2351  XLALUpSampling(S2xHi, dtHi, S2x_s, deltaT);
2352  XLALUpSampling(S2yHi, dtHi, S2y_s, deltaT);
2353  XLALUpSampling(S2zHi, dtHi, S2z_s, deltaT);
2354  XLALUpSampling(EnergyHi, dtHi, Energy_s, deltaT);
2355 
2356  XLALDestroyREAL8Vector(Phi_s);
2357  XLALDestroyREAL8Vector(omega_s);
2358  XLALDestroyREAL8Vector(LNhx_s);
2359  XLALDestroyREAL8Vector(LNhy_s);
2360  XLALDestroyREAL8Vector(LNhz_s);
2361  XLALDestroyREAL8Vector(S1x_s);
2362  XLALDestroyREAL8Vector(S1y_s);
2363  XLALDestroyREAL8Vector(S1z_s);
2364  XLALDestroyREAL8Vector(S2x_s);
2365  XLALDestroyREAL8Vector(S2y_s);
2366  XLALDestroyREAL8Vector(S2z_s);
2367  XLALDestroyREAL8Vector(Energy_s);
2368  XLALDestroyREAL8Vector(omMHi);
2369 
2370  idx=omegaHi->length;
2371  do {
2372  idx--;
2373  LNhS1=(LNhxHi->data[idx]*S1xHi->data[idx]+LNhyHi->data[idx]*S1yHi->data[idx]+LNhzHi->data[idx]*S1zHi->data[idx])/m1Msq;
2374  LNhS2=(LNhxHi->data[idx]*S2xHi->data[idx]+LNhyHi->data[idx]*S2yHi->data[idx]+LNhzHi->data[idx]*S2zHi->data[idx])/m2Msq;
2375  S1S2=(S1xHi->data[idx]*S2xHi->data[idx]+S1yHi->data[idx]*S2yHi->data[idx]+S1zHi->data[idx]*S2zHi->data[idx])/m1Msq/m2Msq;
2376  omegaMatch=OmMatch(LNhS1,LNhS2,S1S1,S1S2,S2S2);
2377  if ((omegaMatch>omegaHi->data[idx])&&(omegaHi->data[idx]<0.1)) {
2378  if (omegaHi->data[idx-1]<omegaHi->data[idx]) iMatchUp=idx;
2379  // The numerical integrator sometimes stops and stores twice the last
2380  // omega value, this 'if' instruction avoids keeping two identical
2381  // values of omega at the end of the integration.
2382  }
2383  } while ((idx>0)&&(iMatchUp==0));
2384 
2385  REAL8Vector *domegaHi = XLALCreateREAL8Vector(stkLenHi);
2386  REAL8Vector *dLNhxHi = XLALCreateREAL8Vector(stkLenHi);
2387  REAL8Vector *dLNhyHi = XLALCreateREAL8Vector(stkLenHi);
2388  REAL8Vector *dLNhzHi = XLALCreateREAL8Vector(stkLenHi);
2389  REAL8Vector *dalphaHi = XLALCreateREAL8Vector(stkLenHi);
2390 
2391  errcode = XLALGenerateWaveDerivative(domegaHi,omegaHi,dtHi);
2392  errcode += XLALGenerateWaveDerivative(dLNhxHi,LNhxHi,dtHi);
2393  errcode += XLALGenerateWaveDerivative(dLNhyHi,LNhyHi,dtHi);
2394  errcode += XLALGenerateWaveDerivative(dLNhzHi,LNhzHi,dtHi);
2395 
2396  if ( (errcode != 0) || (domegaHi->data[iMatchUp]<0.) ) {
2397  XLALPrintError("**** LALSimIMRPhenSpinInspiralRD ERROR ****: error generating derivatives");
2398  XLALPrintError(" m: : %12.5f %12.5f\n",mass1,mass2);
2399  XLALPrintError(" S1: : %12.5f %12.5f %12.5f\n",s1x,s1y,s1z);
2400  XLALPrintError(" S2: : %12.5f %12.5f %12.5f\n",s2x,s2y,s2z);
2401  XLALPrintError(" omM %12.5f om[%d] %12.5f\n",omegaMatch,iMatch,omega->data->data[iMatch]);
2402  XLAL_ERROR(XLAL_EFAILED);
2403  }
2404  else {
2405  REAL8 LNhxy;
2406  for (idx=0;idx<stkLenHi;idx++) {
2407  LNhxy = LNhxHi->data[idx] * LNhxHi->data[idx] + LNhyHi->data[idx] * LNhyHi->data[idx];
2408  if (LNhxy > 0.) {
2409  dalphaHi->data[idx] = (LNhxHi->data[idx] * dLNhyHi->data[idx] - LNhyHi->data[idx] * dLNhxHi->data[idx]) / LNhxy;
2410  } else {
2411  dalphaHi->data[idx] = 0.;
2412  }
2413  }
2414  }
2415 
2416  count = iMatch;
2417  tm=((REAL8) iMatch)*deltaT;
2418  const INT4 fac=2;
2419  REAL8 dt=((REAL8)fac)*dtHi;
2420  REAL8 t0 = tm + ((REAL8) iMatchUp )*dtHi;
2421  REAL8 tm1 = t0 - dt;
2422  REAL8 dom = omegaHi->data[iMatchUp] - omegaHi->data[iMatchUp-fac];
2423  REAL8 tAs = (t0 * domegaHi->data[iMatchUp] - tm1 * dom/dt)/ (domegaHi->data[iMatchUp] - dom/dt);
2424  REAL8 om1 = dom * (tAs -t0)*(tAs-tm1)/dt/tAs;
2425  REAL8 om0 = omegaHi->data[iMatchUp] - om1 / (1. - t0 / tAs);
2426 
2427  REAL8 dalpha1 = (dalphaHi->data[iMatchUp]-dalphaHi->data[iMatchUp-fac]) * (tAs - t0) * (tAs - tm1)/dt/tAs;
2428  REAL8 dalpha0 = dalphaHi->data[iMatchUp] - dalpha1 / (1. - t0 / tAs);
2429 
2430  while ((tm+deltaT)<=t0) {
2431  count++;
2432  tm+=deltaT;
2433  }
2434 
2435  if ((tAs < t0) || (om1 < 0.)) {
2436  XLALPrintError("**** LALSimIMRPhenSpinInspiralRD ERROR ****: Could not attach phen part for:\n");
2437  XLALPrintError(" tAs %12.6e t0 %12.6e om1 %12.6e\n",tAs,t0,om1);
2438  XLALPrintError(" m1 = %14.6e, m2 = %14.6e, inc = %10.6f,\n", mass1, mass2, iota);
2439  XLALPrintError(" S1 = (%10.6f,%10.6f,%10.6f)\n", s1x, s1y, s1z);
2440  XLALPrintError(" S2 = (%10.6f,%10.6f,%10.6f)\n", s2x, s2y, s2z);
2441  XLAL_ERROR(XLAL_EFAILED);
2442  }
2443  else /*if Phen part is sane go for this*/ {
2444 
2445  REAL8 Psi0;
2446  REAL8 alpha0,energy;
2447 
2448  XLALSimInspiralComputeInclAngle(LNhzHi->data[iMatchUp],&trigAngle);
2449  om = omegaHi->data[iMatchUp];
2450  Psi = PhiHi->data[iMatchUp];
2451  Psi0 = Psi + tAs * (om1/Mtime -dalpha1*trigAngle.ci) * log(1. - t0 / tAs);
2452  errcode =XLALSimInspiralComputeAlpha(params,LNhxHi->data[iMatchUp],LNhyHi->data[iMatchUp],S1xHi->data[iMatchUp],S1yHi->data[iMatchUp],S2xHi->data[iMatchUp],S2yHi->data[iMatchUp],&alpha);
2453  alpha0 = alpha + tAs * dalpha1 * log(1. - t0 / tAs);
2454  energy = EnergyHi->data[iMatchUp];
2455 
2456  XLALDestroyREAL8Vector(PhiHi);
2457  XLALDestroyREAL8Vector(omegaHi);
2458  XLALDestroyREAL8Vector(LNhxHi);
2459  XLALDestroyREAL8Vector(LNhyHi);
2460  XLALDestroyREAL8Vector(LNhzHi);
2461  XLALDestroyREAL8Vector(S1xHi);
2462  XLALDestroyREAL8Vector(S1yHi);
2463  XLALDestroyREAL8Vector(S1zHi);
2464  XLALDestroyREAL8Vector(S2xHi);
2465  XLALDestroyREAL8Vector(S2yHi);
2466  XLALDestroyREAL8Vector(S2zHi);
2467  XLALDestroyREAL8Vector(EnergyHi);
2468  XLALDestroyREAL8Vector(domegaHi);
2469  XLALDestroyREAL8Vector(dLNhxHi);
2470  XLALDestroyREAL8Vector(dLNhyHi);
2471  XLALDestroyREAL8Vector(dLNhzHi);
2472  XLALDestroyREAL8Vector(dalphaHi);
2473 
2474  /* Estimate final mass and spin*/
2475  REAL8 finalMass,finalSpin;
2476  errcode=XLALSimIMRPhenSpinFinalMassSpin(&finalMass,&finalSpin,m1,m2,S1S1,S2S2,LNhS1,LNhS2,S1S2,energy);
2477 
2478  /* Get QNM frequencies */
2479  COMPLEX16Vector *modefreqs=XLALCreateCOMPLEX16Vector(1);
2480  errcode+=XLALSimIMRPhenSpinGenerateQNMFreq(modefreqs, 2, 2, finalMass, finalSpin, Mass);
2481  if (errcode) {
2482  XLALPrintError("**** LALSimIMRPhenSpinInspiralRD ERROR ****: impossible to generate RingDown frequency\n");
2483  XLALPrintError( " m (%11.4e %11.4e) f0 %11.4e\n",mass1, mass2, f_start);
2484  XLALPrintError( " S1 (%8.4f %8.4f %8.4f)\n", s1x,s1y,s1z);
2485  XLALPrintError( " S2 (%8.4f %8.4f %8.4f)\n", s2x,s2y,s2z);
2486  XLALDestroyCOMPLEX16Vector(modefreqs);
2487  XLAL_ERROR(XLAL_EFAILED);
2488  }
2489 
2490  REAL8 omegaRD = creal(modefreqs->data[0])*Mass*LAL_MTSUN_SI/2.;
2491  REAL8 frOmRD = fracRD(LNhS1,LNhS2,S1S1,S1S2,S2S2)*omegaRD;
2492 
2493  INT4 up=((int)(deltaT/dtHi));
2494  INT4 upcntP=0,upcnt=0;
2495  INT4 cntI=count;
2496  do {
2497  count++;
2498  for (idx=0;idx<up;idx++) {
2499  tm += dtHi;
2500  om = om1 / (1. - tm / tAs) + om0;
2501  if ((om>=frOmRD)&&(upcntP==0)) {
2502  upcntP=upcnt;
2503  }
2504  Psi = Psi0 + (- tAs * (om1/Mtime-dalpha1*trigAngle.ci) * log(1. - tm / tAs) + (om0/Mtime-dalpha0*trigAngle.ci) * (tm - t0) ) - 2.*om*(1.-eta*pow(om,2./3.))*log(om);
2505  alpha = alpha0 + ( dalpha0 * (tm - t0) - dalpha1 * tAs * log(1. - tm / tAs) );
2506  v = cbrt(om);
2507  v2 = v*v;
2508  upcnt++;
2509  }
2510  errcode =XLALSimSpinInspiralFillL2Modes(hL2tmp,v,eta,dm,Psi,alpha,&trigAngle);
2511  for (kdx=0;kdx<5;kdx++) hL2->data->data[5*count+kdx]=hL2tmp->data[kdx]*amp22ini*v2;
2512  errcode+=XLALSimSpinInspiralFillL3Modes(hL3tmp,v,eta,dm,Psi,alpha,&trigAngle);
2513  for (kdx=0;kdx<7;kdx++) hL3->data->data[7*count+kdx]=hL3tmp->data[kdx]*amp33ini*v2;
2514  errcode+=XLALSimSpinInspiralFillL4Modes(hL4tmp,v,eta,dm,Psi,alpha,&trigAngle);
2515  for (kdx=0;kdx<9;kdx++) hL4->data->data[9*count+kdx]=hL4tmp->data[kdx]*amp44ini*v2*v2;
2516  } while ( (om < frOmRD) && (tm < tAs) );
2517  tPeak=cntI*deltaT+upcntP*dtHi;
2518 
2519  /*--------------------------------------------------------------
2520  * Attach the ringdown waveform to the end of inspiral/merger
2521  -------------------------------------------------------------*/
2522 
2523  static const INT4 nPtsComb=6;
2524  if (upcntP<nPtsComb) {
2525  XLALPrintError("*** LALSimIMRPSpinInspiralRD ERROR: impossible to attach RD");
2526  XLAL_ERROR( XLAL_ENOMEM );
2527  }
2528 
2529  REAL8Vector *PsiMat = XLALCreateREAL8Vector( nPtsComb+2 );
2530  REAL8Vector *alpMat = XLALCreateREAL8Vector( nPtsComb+2 );
2531  REAL8Vector *velMat = XLALCreateREAL8Vector( nPtsComb+2 );
2532 
2533  REAL8Vector *waveR = XLALCreateREAL8Vector( nPtsComb+2 );
2534  REAL8Vector *dwaveR = XLALCreateREAL8Vector( nPtsComb+2 );
2535  REAL8Vector *waveI = XLALCreateREAL8Vector( nPtsComb+2 );
2536  REAL8Vector *dwaveI = XLALCreateREAL8Vector( nPtsComb+2 );
2537  REAL8VectorSequence *inspWaveR = XLALCreateREAL8VectorSequence( 2, nPtsComb );
2538  REAL8VectorSequence *inspWaveI = XLALCreateREAL8VectorSequence( 2, nPtsComb );
2539 
2540  INT4 nRDWave = (INT4) (RD_EFOLDS / fabs(cimag(modefreqs->data[0])) / deltaT);
2541 
2542  REAL8Vector *matchrange=XLALCreateREAL8Vector(3);
2543  matchrange->data[2]=count*deltaT/Mtime;
2544  matchrange->data[0]=tPeak/Mtime-12.;
2545  matchrange->data[1]=tPeak/Mtime;
2546  double dtMat=(matchrange->data[1]-matchrange->data[0])*Mtime/(nPtsComb-1);
2547  tm=tPeak+dtMat;
2548  REAL8Vector *tmArray=XLALCreateREAL8Vector(nPtsComb+2);
2549  for (idx=nPtsComb+1;idx>=0;idx--) {
2550  tmArray->data[idx]=tm;
2551  om = om1 / (1. - tm / tAs) + om0;
2552  PsiMat->data[idx] = Psi0 + (- tAs * (om1/Mtime-dalpha1*trigAngle.ci) * log(1. - tm / tAs) + (om0/Mtime-dalpha0*trigAngle.ci) * (tm - t0) ) - 2.*om*(1.-eta*pow(om,2./3.))*log(om);
2553  alpMat->data[idx] = alpha0 + ( dalpha0 * (tm - t0) - dalpha1 * tAs * log(1. - tm / tAs) );
2554  velMat->data[idx] = cbrt(om);
2555  tm -= dtMat;
2556  }
2557 
2558  /* Check memory was allocated */
2559  if ( !waveR || !dwaveR || !waveI || !dwaveI || !inspWaveR || !inspWaveI ) {
2560  XLALDestroyCOMPLEX16Vector( modefreqs );
2561  if (waveR) XLALDestroyREAL8Vector( waveR );
2562  if (dwaveR) XLALDestroyREAL8Vector( dwaveR );
2563  if (waveI) XLALDestroyREAL8Vector( waveI );
2564  if (dwaveI) XLALDestroyREAL8Vector( dwaveI );
2565  if (inspWaveR) XLALDestroyREAL8VectorSequence( inspWaveR );
2566  if (inspWaveI) XLALDestroyREAL8VectorSequence( inspWaveI );
2567  XLAL_ERROR( XLAL_ENOMEM );
2568  }
2569 
2570  INT4 l,m;
2571  if ( ( modesChoice & LAL_SIM_INSPIRAL_MODES_CHOICE_RESTRICTED) == LAL_SIM_INSPIRAL_MODES_CHOICE_RESTRICTED ) {
2572  REAL8Vector *rdwave1l2 = XLALCreateREAL8Vector( nRDWave );
2573  REAL8Vector *rdwave2l2 = XLALCreateREAL8Vector( nRDWave );
2574  memset( rdwave1l2->data, 0, rdwave1l2->length * sizeof( REAL8 ) );
2575  memset( rdwave2l2->data, 0, rdwave2l2->length * sizeof( REAL8 ) );
2576  l=2;
2577  for (m=-l;m<=l;m++) {
2578  for (idx=0;idx<nPtsComb+2;idx++) {
2579  errcode =XLALSimSpinInspiralFillL2Modes(hL2tmp,velMat->data[idx],eta,dm,PsiMat->data[idx],alpMat->data[idx],&trigAngle);
2580  waveR->data[idx]=creal(hL2tmp->data[m+l])*amp22ini*pow(velMat->data[idx],2);
2581  waveI->data[idx]=cimag(hL2tmp->data[m+l])*amp22ini*pow(velMat->data[idx],2);
2582  }
2583 
2584  errcode+=XLALGenerateWaveDerivative(dwaveR,waveR,dtMat);
2585  errcode+=XLALGenerateWaveDerivative(dwaveI,waveI,dtMat);
2586  for (idx=0;idx<nPtsComb;idx++) {
2587  inspWaveR->data[idx] =waveR->data[idx+1];
2588  inspWaveR->data[idx+ nPtsComb] =dwaveR->data[idx+1];
2589  inspWaveI->data[idx] =waveI->data[idx+1];
2590  inspWaveI->data[idx+ nPtsComb] =dwaveI->data[idx+1];
2591  }
2592 
2593  FILE *out=fopen("checkiwPS.dat","w");
2594  for (idx=0;idx<(int)tmArray->length;idx++) {
2595  fprintf(out," %d %12.4e %12.4e %12.4e\n",idx,tmArray->data[idx],.631*waveR->data[idx],.631*waveI->data[idx]);
2596  }
2597  fclose(out);
2598 
2599  if (modefreqs) XLALDestroyCOMPLEX16Vector(modefreqs);
2600  modefreqs=XLALCreateCOMPLEX16Vector(nModes);
2601  errcode+=XLALSimIMRPhenSpinGenerateQNMFreq(modefreqs, l, abs(m), finalMass, dsign(m)*finalSpin, Mass);
2602  errcode+=XLALSimIMRHybridRingdownWave(rdwave1l2,rdwave2l2,deltaT,mass1,mass2,inspWaveR,inspWaveI,modefreqs,matchrange);
2603  for (idx=0;idx<=count;idx++) {
2604  hLMtmp->data->data[idx]=hL2->data->data[5*idx+(l+m)];
2605  }
2606  for (idx=count; idx<count+nRDWave-1;idx++) {
2607  hLMtmp->data->data[idx]=rdwave1l2->data[idx-count+1]+I*rdwave2l2->data[idx-count+1];
2608  }
2609  XLALSimAddMode(hPtmp,hCtmp,hLMtmp,iota,phiN,l,m,0);
2610  }
2611  XLALDestroyREAL8Vector(rdwave1l2);
2612  XLALDestroyREAL8Vector(rdwave2l2);
2613  }
2614  XLALDestroyCOMPLEX16Vector(hL2tmp);
2615  XLALDestroyCOMPLEX16TimeSeries(hL2);
2616  XLALDestroyREAL8Vector(tmArray);
2617 
2618  if ( ( modesChoice & LAL_SIM_INSPIRAL_MODES_CHOICE_3L) == LAL_SIM_INSPIRAL_MODES_CHOICE_3L ) {
2619  printf("Aggiunto modo l=3\n");
2620  REAL8Vector *rdwave1l3 = XLALCreateREAL8Vector( nRDWave );
2621  REAL8Vector *rdwave2l3 = XLALCreateREAL8Vector( nRDWave );
2622  l=3;
2623  for (m=-l;m<=l;m++) {
2624  for (idx=0;idx<nPtsComb+2;idx++) {
2625  errcode =XLALSimSpinInspiralFillL3Modes(hL3tmp,velMat->data[idx],eta,dm,PsiMat->data[idx],alpMat->data[idx],&trigAngle);
2626  waveR->data[idx]=creal(hL3tmp->data[m+l])*amp33ini*velMat->data[idx]*velMat->data[idx];
2627  waveI->data[idx]=cimag(hL3tmp->data[m+l])*amp33ini*velMat->data[idx]*velMat->data[idx];
2628  }
2629  errcode+=XLALGenerateWaveDerivative(waveR,dwaveR,dtMat);
2630  errcode+=XLALGenerateWaveDerivative(waveI,dwaveI,dtMat);
2631  for (idx=0;idx<nPtsComb;idx++) {
2632  inspWaveR->data[idx] = waveR->data[idx+1];
2633  inspWaveR->data[idx+ nPtsComb] = dwaveR->data[idx+1];
2634  inspWaveI->data[idx] = waveI->data[idx+1];
2635  inspWaveI->data[idx+ nPtsComb] = dwaveI->data[idx+1];
2636  }
2637  errcode+=XLALSimIMRPhenSpinGenerateQNMFreq(modefreqs, l, abs(m), finalMass, dsign(m)*finalSpin, Mass);
2638  errcode+=XLALSimIMRHybridRingdownWave(rdwave1l3,rdwave2l3,deltaT,mass1,mass2,inspWaveR,inspWaveI,modefreqs,matchrange);
2639  for (idx=intLen;idx<count;idx++) hLMtmp->data->data[idx]=hL3->data->data[7*idx+(l+m)];
2640  for (idx=count;idx<nRDWave;idx++) hLMtmp->data->data[idx]=rdwave1l3->data[idx-count]+I*rdwave2l3->data[idx-count];
2641  XLALSimAddMode(hPtmp,hCtmp,hLMtmp,iota,phiN,l,m,0);
2642  }
2643  XLALDestroyREAL8Vector(rdwave1l3);
2644  XLALDestroyREAL8Vector(rdwave2l3);
2645  }
2646  XLALDestroyCOMPLEX16Vector(hL3tmp);
2647  XLALDestroyCOMPLEX16TimeSeries(hL3);
2648 
2649  if ( ( modesChoice & LAL_SIM_INSPIRAL_MODES_CHOICE_4L) == LAL_SIM_INSPIRAL_MODES_CHOICE_4L ) {
2650  printf("Aggiunto modo l=4\n");
2651  REAL8Vector *rdwave1l4 = XLALCreateREAL8Vector( nRDWave );
2652  REAL8Vector *rdwave2l4 = XLALCreateREAL8Vector( nRDWave );
2653  l=4;
2654  for (m=-l;m<=l;m++) {
2655  for (idx=0;idx<nPtsComb+2;idx++) {
2656  kdx=upcnt+idx-nPtsComb-1;
2657  errcode =XLALSimSpinInspiralFillL4Modes(hL4tmp,velMat->data[kdx],eta,dm,PsiMat->data[kdx],alpMat->data[kdx],&trigAngle);
2658  waveR->data[idx]=creal(hL4tmp->data[m+l])*amp44ini*pow(velMat->data[idx],4);
2659  waveI->data[idx]=cimag(hL4tmp->data[m+l])*amp44ini*pow(velMat->data[idx],4);
2660  }
2661  errcode =XLALGenerateWaveDerivative(waveR,dwaveR,dtMat);
2662  errcode+=XLALGenerateWaveDerivative(waveI,dwaveI,dtMat);
2663  for (idx=0;idx<nPtsComb;idx++) {
2664  inspWaveR->data[idx] = waveR->data[idx+1];
2665  inspWaveR->data[idx+ nPtsComb] = dwaveR->data[idx+1];
2666  inspWaveI->data[idx] = waveI->data[idx+1];
2667  inspWaveI->data[idx+ nPtsComb] = dwaveI->data[idx+1];
2668  }
2669  errcode+= XLALSimIMRPhenSpinGenerateQNMFreq(modefreqs,l,abs(m), finalMass, dsign(m)*finalSpin, Mass);
2670  errcode+= XLALSimIMRHybridRingdownWave(rdwave1l4,rdwave2l4,deltaT,mass1,mass2,inspWaveR,
2671  inspWaveI,modefreqs,matchrange);
2672  for (idx=intLen;idx<count;idx++) hLMtmp->data->data[idx-intLen]=hL4->data->data[9*idx+(l+m)];
2673  for (idx=0;idx<nRDWave;idx++) hLMtmp->data->data[count-intLen+idx]=rdwave1l4->data[idx]+I*rdwave2l4->data[idx];
2674  XLALSimAddMode(hPtmp,hCtmp,hLMtmp,iota,phiN,l,m,0);
2675  }
2676  XLALDestroyREAL8Vector(rdwave1l4);
2677  XLALDestroyREAL8Vector(rdwave2l4);
2678  }
2679  XLALDestroyCOMPLEX16Vector(hL4tmp);
2680  XLALDestroyCOMPLEX16TimeSeries(hL4);
2681 
2682  XLALDestroyCOMPLEX16TimeSeries(hLMtmp);
2683  XLALDestroyREAL8Vector(matchrange);
2684  XLALDestroyREAL8Vector(PsiMat);
2685  XLALDestroyREAL8Vector(velMat);
2686  XLALDestroyREAL8Vector(alpMat);
2687 
2688  XLALDestroyCOMPLEX16Vector( modefreqs );
2689  XLALDestroyREAL8Vector( waveR );
2690  XLALDestroyREAL8Vector( dwaveR );
2691  XLALDestroyREAL8Vector( waveI );
2692  XLALDestroyREAL8Vector( dwaveI );
2693  XLALDestroyREAL8VectorSequence( inspWaveR );
2694  XLALDestroyREAL8VectorSequence( inspWaveI );
2695  if (errcode) XLAL_ERROR( XLAL_EFUNC );
2696  count+=nRDWave-1;
2697 
2698  } /* End of: if phen part not sane*/
2699 
2700  } /*End of if errcodeInt==LALSIMINSPIRAL_PST4_TEST_OMEGAMATCH*/
2701  else {
2702  if (omega) XLALDestroyREAL8TimeSeries(omega);
2703  if (Phi) XLALDestroyREAL8TimeSeries(Phi);
2704  if (LNhatx) XLALDestroyREAL8TimeSeries(LNhatx);
2705  if (LNhaty) XLALDestroyREAL8TimeSeries(LNhaty);
2706  if (LNhatz) XLALDestroyREAL8TimeSeries(LNhatz);
2707  if (S1x) XLALDestroyREAL8TimeSeries(S1x);
2708  if (S1y) XLALDestroyREAL8TimeSeries(S1y);
2709  if (S1z) XLALDestroyREAL8TimeSeries(S1z);
2710  if (S2x) XLALDestroyREAL8TimeSeries(S2x);
2711  if (S2y) XLALDestroyREAL8TimeSeries(S2y);
2712  if (S2z) XLALDestroyREAL8TimeSeries(S2z);
2713  if (Energy) XLALDestroyREAL8TimeSeries(Energy);
2714  }
2715 
2716  if ((*hPlus) && (*hCross)) {
2717  if ((*hPlus)->data->length!=(*hCross)->data->length) {
2718  XLALPrintError("*** LALSimIMRPSpinInspiralRD ERROR: h+ and hx differ in length: %d vs. %d\n",(*hPlus)->data->length,(*hCross)->data->length);
2719  XLAL_ERROR(XLAL_EFAILED);
2720  }
2721  else {
2722  if ((int)(*hPlus)->data->length<count) {
2723  XLALPrintError("*** LALSimIMRPSpinInspiralRD ERROR: h+ and hx too short: %d vs. %d\n",(*hPlus)->data->length,count);
2724  XLAL_ERROR(XLAL_EFAILED);
2725  }
2726  else {
2727  XLALGPSAdd(&((*hPlus)->epoch),-tPeak);
2728  XLALGPSAdd(&((*hCross)->epoch),-tPeak);
2729  }
2730  }
2731  }
2732  else {
2733  XLALGPSAdd(&tStart,-tPeak);
2734  INT4 wfLen=1;
2735  while (wfLen<count) wfLen*=2;
2736  *hPlus = XLALCreateREAL8TimeSeries("H+", &tStart, 0.0, deltaT, &lalDimensionlessUnit, wfLen);
2737  *hCross = XLALCreateREAL8TimeSeries("Hx", &tStart, 0.0, deltaT, &lalDimensionlessUnit, wfLen);
2738  if(*hPlus == NULL || *hCross == NULL)
2739  XLAL_ERROR(XLAL_ENOMEM);
2740  }
2741 
2742  INT4 minLen=hPtmp->data->length < (*hPlus)->data->length ? hPtmp->data->length : (*hPlus)->data->length;
2743  for (idx=0;idx<minLen;idx++) {
2744  (*hPlus)->data->data[idx] =hPtmp->data->data[idx];
2745  (*hCross)->data->data[idx]=hCtmp->data->data[idx];
2746  }
2747  for (idx=minLen;idx<(int)(*hPlus)->data->length;idx++) {
2748  (*hPlus)->data->data[idx] =0.;
2749  (*hCross)->data->data[idx]=0.;
2750  }
2751 
2752  XLALDestroyREAL8TimeSeries(hPtmp);
2753  XLALDestroyREAL8TimeSeries(hCtmp);
2754 
2755  return count;
2756 } /* End of XLALSimIMRPhenSpinInspiralRDGenerator*/
2757 
2758 /** @} */
LALMalloc
#define LALMalloc(n)
LALFree
#define LALFree(p)
c1
const double c1
Definition: LALSimIMRLackeyTidal2013.c:51
c2
const double c2
Definition: LALSimIMRLackeyTidal2013.c:52
XLALSpinInspiralDerivatives
static INT4 XLALSpinInspiralDerivatives(UNUSED double t, const double values[], double dvalues[], void *mparams)
Definition: LALSimIMRPSpinInspiralRD.c:379
LAL_PST4_RELATIVE_TOLERANCE
#define LAL_PST4_RELATIVE_TOLERANCE
Definition: LALSimIMRPSpinInspiralRD.c:55
XLALSimIMRPhenSpinParamsSetup
static INT4 XLALSimIMRPhenSpinParamsSetup(LALSimInspiralPhenSpinTaylorT4Coeffs *params, REAL8 dt, REAL8 fStart, REAL8 fEnd, REAL8 mass1, REAL8 mass2, REAL8 lambda1, REAL8 lambda2, REAL8 quadparam1, REAL8 quadparam2, int spinO, int tideO, LALDict *LALparams, UINT4 order)
Convenience function to set up XLALSimInspiralSpinTaylotT4Coeffs struct.
Definition: LALSimIMRPSpinInspiralRD.c:160
XLALSimIMRPhenSpinInitialize
static INT4 XLALSimIMRPhenSpinInitialize(REAL8 mass1, REAL8 mass2, REAL8 lambda1, REAL8 lambda2, REAL8 quadparam1, REAL8 quadparam2, REAL8 *yinit, REAL8 fStart, REAL8 fEnd, REAL8 deltaT, INT4 phaseO, LALSimInspiralPhenSpinTaylorT4Coeffs *params, LALDict *LALparams, Approximant approx)
PhenSpin Initialization.
Definition: LALSimIMRPSpinInspiralRD.c:1306
LALSIMINSPIRAL_PST4_TEST_OMEGAMATCH
#define LALSIMINSPIRAL_PST4_TEST_OMEGAMATCH
Definition: LALSimIMRPSpinInspiralRD.c:50
dsign
static REAL8 dsign(INT4 i)
Definition: LALSimIMRPSpinInspiralRD.c:65
XLALSimInspiralComputeAlpha
static INT4 XLALSimInspiralComputeAlpha(LALSimInspiralPhenSpinTaylorT4Coeffs params, REAL8 LNhx, REAL8 LNhy, REAL8 S1x, REAL8 S1y, REAL8 S2x, REAL8 S2y, REAL8 *alpha)
The following lines are necessary in the case L is initially parallel to N so that alpha is undefined...
Definition: LALSimIMRPSpinInspiralRD.c:1079
LAL_NUM_PST4_VARIABLES
#define LAL_NUM_PST4_VARIABLES
Definition: LALSimIMRPSpinInspiralRD.c:52
XLALSimInspiralComputeInclAngle
static INT4 XLALSimInspiralComputeInclAngle(REAL8 ciota, LALSimInspiralInclAngle *angle)
Definition: LALSimIMRPSpinInspiralRD.c:1043
LALSIMINSPIRAL_PST4_TEST_OMEGADOT
#define LALSIMINSPIRAL_PST4_TEST_OMEGADOT
Definition: LALSimIMRPSpinInspiralRD.c:45
LALSIMINSPIRAL_PST4_TEST_FREQBOUND
#define LALSIMINSPIRAL_PST4_TEST_FREQBOUND
Definition: LALSimIMRPSpinInspiralRD.c:48
appendTSandFree
static REAL8TimeSeries * appendTSandFree(REAL8TimeSeries *start, REAL8TimeSeries *end)
Definition: LALSimIMRPSpinInspiralRD.c:1369
fracRD
static REAL8 fracRD(REAL8 LNhS1, REAL8 LNhS2, REAL8 S1S1, REAL8 S1S2, REAL8 S2S2)
Definition: LALSimIMRPSpinInspiralRD.c:140
OmMatch
static REAL8 OmMatch(REAL8 LNhS1, REAL8 LNhS2, REAL8 S1S1, REAL8 S1S2, REAL8 S2S2)
Definition: LALSimIMRPSpinInspiralRD.c:115
nModes
#define nModes
Definition: LALSimIMRPSpinInspiralRD.c:62
XLALGenerateWaveDerivative
static INT4 XLALGenerateWaveDerivative(REAL8Vector *dwave, REAL8Vector *wave, REAL8 dt)
Definition: LALSimIMRPSpinInspiralRD.c:566
rotateY
static void rotateY(REAL8 phi, REAL8 *vx, REAL8 *vy, REAL8 *vz)
Here we use the following convention: the coordinates of the spin vectors spin1,2 and the inclination...
Definition: LALSimIMRPSpinInspiralRD.c:1142
XLALSimSpinInspiralFillL4Modes
static INT4 XLALSimSpinInspiralFillL4Modes(COMPLEX16Vector *hL4, UNUSED REAL8 v, REAL8 eta, UNUSED REAL8 dm, REAL8 Psi, REAL8 alpha, LALSimInspiralInclAngle *an)
Definition: LALSimIMRPSpinInspiralRD.c:870
XLALUpSampling
static INT4 XLALUpSampling(REAL8Vector *vHi, REAL8 dtHi, REAL8Vector *v, REAL8 dt)
Definition: LALSimIMRPSpinInspiralRD.c:1605
LAL_PST4_ABSOLUTE_TOLERANCE
#define LAL_PST4_ABSOLUTE_TOLERANCE
Definition: LALSimIMRPSpinInspiralRD.c:54
LALSIMINSPIRAL_PST4_DERIVATIVE_OMEGANONPOS
#define LALSIMINSPIRAL_PST4_DERIVATIVE_OMEGANONPOS
Definition: LALSimIMRPSpinInspiralRD.c:49
XLALSimIMRPhenSpinTest
static INT4 XLALSimIMRPhenSpinTest(UNUSED double t, const double values[], double dvalues[], void *mparams)
Definition: LALSimIMRPSpinInspiralRD.c:690
XLALSimSpinInspiralFillL2Modes
static INT4 XLALSimSpinInspiralFillL2Modes(COMPLEX16Vector *hL2, REAL8 v, REAL8 eta, REAL8 dm, REAL8 Psi, REAL8 alpha, LALSimInspiralInclAngle *an)
Definition: LALSimIMRPSpinInspiralRD.c:762
rotateZ
static void rotateZ(REAL8 phi, REAL8 *vx, REAL8 *vy, REAL8 *vz)
Definition: LALSimIMRPSpinInspiralRD.c:1153
XLALSimIMRHybridRingdownWave
static INT4 XLALSimIMRHybridRingdownWave(REAL8Vector *rdwave1, REAL8Vector *rdwave2, const REAL8 dt, const REAL8 mass1, const REAL8 mass2, REAL8VectorSequence *inspwave1, REAL8VectorSequence *inspwave2, COMPLEX16Vector *modefreqs, REAL8Vector *matchrange)
Definition: LALSimIMRPSpinInspiralRD.c:1385
LALSIMINSPIRAL_PST4_TEST_ENERGY
#define LALSIMINSPIRAL_PST4_TEST_ENERGY
Definition: LALSimIMRPSpinInspiralRD.c:44
minIntLen
#define minIntLen
Definition: LALSimIMRPSpinInspiralRD.c:58
XLALSimInspiralSpinTaylorT4Engine
static INT4 XLALSimInspiralSpinTaylorT4Engine(REAL8TimeSeries **omega, REAL8TimeSeries **Phi, REAL8TimeSeries **LNhatx, REAL8TimeSeries **LNhaty, REAL8TimeSeries **LNhatz, REAL8TimeSeries **S1x, REAL8TimeSeries **S1y, REAL8TimeSeries **S1z, REAL8TimeSeries **S2x, REAL8TimeSeries **S2y, REAL8TimeSeries **S2z, REAL8TimeSeries **Energy, const REAL8 yinit[], const INT4 lengthH, const Approximant approx, LALSimInspiralPhenSpinTaylorT4Coeffs *params)
Definition: LALSimIMRPSpinInspiralRD.c:922
XLALSimSpinInspiralTest
static INT4 XLALSimSpinInspiralTest(UNUSED double t, const double values[], double dvalues[], void *mparams)
Definition: LALSimIMRPSpinInspiralRD.c:651
RD_EFOLDS
#define RD_EFOLDS
Definition: LALSimIMRPSpinInspiralRD.c:63
XLALSimSpinInspiralFillL3Modes
static INT4 XLALSimSpinInspiralFillL3Modes(COMPLEX16Vector *hL3, REAL8 v, REAL8 eta, REAL8 dm, REAL8 Psi, REAL8 alpha, LALSimInspiralInclAngle *an)
Definition: LALSimIMRPSpinInspiralRD.c:813
XLALSimIMRPhenSpinInspiralSetAxis
static INT4 XLALSimIMRPhenSpinInspiralSetAxis(REAL8 **yinitOut, REAL8 *iota, REAL8 *phiN, REAL8 *yinitIn, const REAL8 inclination, const REAL8 mass1, const REAL8 mass2, LALSimInspiralFrameAxis axisChoice)
Definition: LALSimIMRPSpinInspiralRD.c:1165
LALSIMINSPIRAL_PST4_TEST_OMEGANAN
#define LALSIMINSPIRAL_PST4_TEST_OMEGANAN
Definition: LALSimIMRPSpinInspiralRD.c:47
XLALGetApproximantFromString
int XLALGetApproximantFromString(const char *waveform)
Definition: LALSimInspiral.c:4206
XLALSimInspiralTaylorLength
REAL8 XLALSimInspiralTaylorLength(REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 f_min, int O)
Definition: LALSimInspiralTaylorLength.c:131
XLALSimInspiralSpinDot_4PNS2OCoeffAvg
static const REAL8 UNUSED XLALSimInspiralSpinDot_4PNS2OCoeffAvg
Definition: LALSimInspiralPNCoefficients.c:2347
XLALSimInspiralTaylorT4wdot_12PNTidalCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_12PNTidalCoeff(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:2162
XLALSimInspiralTaylorT4wdot_3PNSOCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_3PNSOCoeff(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:1841
XLALSimInspiralPNEnergy_5PNSOCoeff
static REAL8 UNUSED XLALSimInspiralPNEnergy_5PNSOCoeff(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:174
XLALSimInspiralPNEnergy_4PNCoeff
static REAL8 UNUSED XLALSimInspiralPNEnergy_4PNCoeff(REAL8 eta)
Definition: LALSimInspiralPNCoefficients.c:63
XLALSimInspiralTaylorT4wdot_7PNCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_7PNCoeff(REAL8 eta)
Definition: LALSimInspiralPNCoefficients.c:2120
XLALSimInspiralTaylorT4wdot_4PNCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_4PNCoeff(REAL8 eta)
Definition: LALSimInspiralPNCoefficients.c:1848
XLALSimInspiralPNEnergy_10PNTidalCoeff
static REAL8 UNUSED XLALSimInspiralPNEnergy_10PNTidalCoeff(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:322
XLALSimInspiralSpinDot_4PNS2CoeffAvg
static const REAL8 UNUSED XLALSimInspiralSpinDot_4PNS2CoeffAvg
Definition: LALSimInspiralPNCoefficients.c:2344
XLALSimInspiralPNEnergy_2PNCoeff
static REAL8 UNUSED XLALSimInspiralPNEnergy_2PNCoeff(REAL8 eta)
Definition: LALSimInspiralPNCoefficients.c:56
XLALSimInspiralTaylorT4wdot_0PNCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_0PNCoeff(REAL8 eta)
Computes the PN Coefficients for using in the TaylorT4 frequency equation.
Definition: LALSimInspiralPNCoefficients.c:1820
XLALSimInspiralTaylorT4wdot_2PNCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_2PNCoeff(REAL8 eta)
Definition: LALSimInspiralPNCoefficients.c:1827
XLALSimInspiralSpinDot_4PNQMSOCoeffAvg
static REAL8 UNUSED XLALSimInspiralSpinDot_4PNQMSOCoeffAvg(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:2353
XLALSimInspiralPNEnergy_4PNS1S2CoeffAvg
static REAL8 UNUSED XLALSimInspiralPNEnergy_4PNS1S2CoeffAvg(REAL8 eta)
Definition: LALSimInspiralPNCoefficients.c:110
XLALSimInspiralTaylorT4wdot_6PNCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_6PNCoeff(REAL8 eta)
Definition: LALSimInspiralPNCoefficients.c:1976
XLALSimInspiralPNEnergy_0PNCoeff
static REAL8 UNUSED XLALSimInspiralPNEnergy_0PNCoeff(REAL8 eta)
Computes the PN Coefficients for using in the PN energy equation.
Definition: LALSimInspiralPNCoefficients.c:49
XLALSimInspiralTaylorT4wdot_4PNS1OS2OCoeffAvg
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_4PNS1OS2OCoeffAvg(REAL8 eta)
Definition: LALSimInspiralPNCoefficients.c:1883
XLALSimInspiralTaylorT4wdot_4PNQMS1OS1OCoeffAvg
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_4PNQMS1OS1OCoeffAvg(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:1946
XLALSimInspiralPNEnergy_3PNSOCoeff
static REAL8 UNUSED XLALSimInspiralPNEnergy_3PNSOCoeff(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:101
XLALSimInspiralTaylorT4wdot_10PNTidalCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_10PNTidalCoeff(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:2155
XLALSimInspiralTaylorT4wdot_6PNLogCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_6PNLogCoeff(REAL8 UNUSED eta)
Definition: LALSimInspiralPNCoefficients.c:1987
XLALSimInspiralTaylorT4wdot_4PNS1S2CoeffAvg
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_4PNS1S2CoeffAvg(REAL8 eta)
Definition: LALSimInspiralPNCoefficients.c:1876
XLALSimInspiralPNEnergy_4PNQMS1S1CoeffAvg
static REAL8 UNUSED XLALSimInspiralPNEnergy_4PNQMS1S1CoeffAvg(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:142
XLALSimInspiralPNEnergy_12PNTidalCoeff
static REAL8 UNUSED XLALSimInspiralPNEnergy_12PNTidalCoeff(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:329
XLALSimInspiralSpinDot_3PNCoeff
static REAL8 UNUSED XLALSimInspiralSpinDot_3PNCoeff(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:2328
XLALSimInspiralPNEnergy_4PNS1OS2OCoeffAvg
static REAL8 UNUSED XLALSimInspiralPNEnergy_4PNS1OS2OCoeffAvg(REAL8 eta)
Definition: LALSimInspiralPNCoefficients.c:126
XLALSimInspiralTaylorT4wdot_4PNS1OS1OCoeffAvg
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_4PNS1OS1OCoeffAvg(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:1911
XLALSimInspiralPNEnergy_4PNQMS1OS1OCoeffAvg
static REAL8 UNUSED XLALSimInspiralPNEnergy_4PNQMS1OS1OCoeffAvg(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:151
XLALSimInspiralTaylorT4wdot_4PNQMS1S1CoeffAvg
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_4PNQMS1S1CoeffAvg(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:1939
XLALSimInspiralSpinDot_5PNCoeff
static REAL8 UNUSED XLALSimInspiralSpinDot_5PNCoeff(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:2370
XLALSimInspiralTaylorT4wdot_5PNCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_5PNCoeff(REAL8 eta)
Definition: LALSimInspiralPNCoefficients.c:1953
XLALSimInspiralPNEnergy_6PNCoeff
static REAL8 UNUSED XLALSimInspiralPNEnergy_6PNCoeff(REAL8 eta)
Definition: LALSimInspiralPNCoefficients.c:70
XLALSimInspiralTaylorT4wdot_3PNCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_3PNCoeff(REAL8 UNUSED eta)
Definition: LALSimInspiralPNCoefficients.c:1834
XLALSimInspiralTaylorT4wdot_6PNSOCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_6PNSOCoeff(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:1994
XLALSimInspiralTaylorT4wdot_5PNSOCoeff
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_5PNSOCoeff(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:1969
XLALSimInspiralTaylorT4wdot_4PNS1S1CoeffAvg
static REAL8 UNUSED XLALSimInspiralTaylorT4wdot_4PNS1S1CoeffAvg(REAL8 mByM)
Definition: LALSimInspiralPNCoefficients.c:1904
XLALSimInspiralWaveformParamsLookupFrameAxis
INT4 XLALSimInspiralWaveformParamsLookupFrameAxis(LALDict *params)
XLALSimInspiralWaveformParamsLookupModesChoice
INT4 XLALSimInspiralWaveformParamsLookupModesChoice(LALDict *params)
XLALSimInspiralWaveformParamsLookupPNTidalOrder
INT4 XLALSimInspiralWaveformParamsLookupPNTidalOrder(LALDict *params)
XLALSimInspiralWaveformParamsLookupNonGRPhi1
REAL8 XLALSimInspiralWaveformParamsLookupNonGRPhi1(LALDict *params)
XLALSimInspiralWaveformParamsLookupPNSpinOrder
INT4 XLALSimInspiralWaveformParamsLookupPNSpinOrder(LALDict *params)
XLALSimIMRPhenSpinGenerateQNMFreq
static int UNUSED XLALSimIMRPhenSpinGenerateQNMFreq(COMPLEX16Vector *modefreqs, UINT4 l, INT4 m, REAL8 finalMass, REAL8 finalSpin, REAL8 totalMass)
Definition: LALSimPhenSpinRingDown.c:34
fprintf
#define fprintf
XLAL_CALLGSL
#define XLAL_CALLGSL(statement)
s
int s
Definition: bh_qnmode.c:137
l
int l
Definition: bh_qnmode.c:135
dt
double dt
Definition: bh_ringdown.c:113
i
double i
Definition: bh_ringdown.c:118
vy
const double vy
Definition: exact_derivatives-HrealHeader.c:53
qq
const double qq
Definition: exact_derivatives-HrealHeader.c:95
a2
const double a2
Definition: exact_derivatives-HrealHeader.c:33
vz
const double vz
Definition: exact_derivatives-HrealHeader.c:54
vx
const double vx
Definition: exact_derivatives-HrealHeader.c:52
__attribute__
#define __attribute__(x)
XLALDestroyREAL8Array
void XLALDestroyREAL8Array(REAL8Array *)
XLALCreateREAL8VectorSequence
REAL8VectorSequence * XLALCreateREAL8VectorSequence(UINT4 length, UINT4 veclen)
XLALDestroyREAL8VectorSequence
void XLALDestroyREAL8VectorSequence(REAL8VectorSequence *vecseq)
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_G_SI
#define LAL_G_SI
LIGOTIMEGPSZERO
#define LIGOTIMEGPSZERO
REAL8
double REAL8
UINT4
uint32_t UINT4
INT4
int32_t INT4
XLALSimIMRPhenSpinInspiralRDGenerator
INT4 XLALSimIMRPhenSpinInspiralRDGenerator(REAL8TimeSeries **hPlus, REAL8TimeSeries **hCross, REAL8 phi_start, REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 f_start, REAL8 f_ref, REAL8 r, REAL8 iota, REAL8 s1x, REAL8 s1y, REAL8 s1z, REAL8 s2x, REAL8 s2y, REAL8 s2z, INT4 phaseO, INT4 UNUSED ampO, REAL8 lambda1, REAL8 lambda2, REAL8 quadparam1, REAL8 quadparam2, LALDict *LALparams)
Driver routine for generating PhenSpinRD waveforms.
Definition: LALSimIMRPSpinInspiralRD.c:2036
XLALSimIMRPhenSpinFinalMassSpin
INT4 XLALSimIMRPhenSpinFinalMassSpin(REAL8 *finalMass, REAL8 *finalSpin, REAL8 mass1, REAL8 mass2, REAL8 s1s1, REAL8 s2s2, REAL8 s1L, REAL8 s2L, REAL8 s1s2, REAL8 energy)
Definition: LALSimIMRPSpinInspiralRD.c:1956
XLALSimSpinInspiralGenerator
INT4 XLALSimSpinInspiralGenerator(REAL8TimeSeries **hPlus, REAL8TimeSeries **hCross, REAL8 phi_start, REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 f_start, REAL8 f_ref, REAL8 r, REAL8 iota, REAL8 s1x, REAL8 s1y, REAL8 s1z, REAL8 s2x, REAL8 s2y, REAL8 s2z, INT4 phaseO, INT4 UNUSED ampO, REAL8 lambda1, REAL8 lambda2, REAL8 quadparam1, REAL8 quadparam2, LALDict *LALparams)
Driver routine to compute the PhenSpin Inspiral waveform without ring-down.
Definition: LALSimIMRPSpinInspiralRD.c:1673
LAL_MAX_PN_ORDER
#define LAL_MAX_PN_ORDER
Definition: LALSimInspiral.h:220
LALSimInspiralFrameAxis
LALSimInspiralFrameAxis
Enumerator for choosing the reference frame associated with PSpinInspiralRD waveforms.
Definition: LALSimInspiralWaveformFlags.h:110
Approximant
Approximant
Enum that specifies the PN approximant to be used in computing the waveform.
Definition: LALSimInspiral.h:229
LAL_SIM_INSPIRAL_MODES_CHOICE_3L
@ LAL_SIM_INSPIRAL_MODES_CHOICE_3L
Inlude only l=3 modes.
Definition: LALSimInspiralWaveformFlags.h:126
LAL_SIM_INSPIRAL_MODES_CHOICE_4L
@ LAL_SIM_INSPIRAL_MODES_CHOICE_4L
Inlude only l=4 modes.
Definition: LALSimInspiralWaveformFlags.h:128
LAL_SIM_INSPIRAL_SPIN_ORDER_0PN
@ LAL_SIM_INSPIRAL_SPIN_ORDER_0PN
Definition: LALSimInspiralWaveformFlags.h:59
LAL_SIM_INSPIRAL_SPIN_ORDER_25PN
@ LAL_SIM_INSPIRAL_SPIN_ORDER_25PN
Definition: LALSimInspiralWaveformFlags.h:64
LAL_SIM_INSPIRAL_SPIN_ORDER_ALL
@ LAL_SIM_INSPIRAL_SPIN_ORDER_ALL
Definition: LALSimInspiralWaveformFlags.h:67
LAL_SIM_INSPIRAL_SPIN_ORDER_2PN
@ LAL_SIM_INSPIRAL_SPIN_ORDER_2PN
Definition: LALSimInspiralWaveformFlags.h:63
LAL_SIM_INSPIRAL_SPIN_ORDER_15PN
@ LAL_SIM_INSPIRAL_SPIN_ORDER_15PN
Definition: LALSimInspiralWaveformFlags.h:62
LAL_SIM_INSPIRAL_SPIN_ORDER_1PN
@ LAL_SIM_INSPIRAL_SPIN_ORDER_1PN
Definition: LALSimInspiralWaveformFlags.h:61
LAL_SIM_INSPIRAL_SPIN_ORDER_05PN
@ LAL_SIM_INSPIRAL_SPIN_ORDER_05PN
Definition: LALSimInspiralWaveformFlags.h:60
LAL_SIM_INSPIRAL_SPIN_ORDER_3PN
@ LAL_SIM_INSPIRAL_SPIN_ORDER_3PN
Definition: LALSimInspiralWaveformFlags.h:65
LAL_SIM_INSPIRAL_FRAME_AXIS_VIEW
@ LAL_SIM_INSPIRAL_FRAME_AXIS_VIEW
Set z-axis along direction of GW propagation (line of sight)
Definition: LALSimInspiralWaveformFlags.h:111
LAL_SIM_INSPIRAL_FRAME_AXIS_TOTAL_J
@ LAL_SIM_INSPIRAL_FRAME_AXIS_TOTAL_J
Set z-axis along the initial total angular momentum.
Definition: LALSimInspiralWaveformFlags.h:112
LAL_SIM_INSPIRAL_FRAME_AXIS_ORBITAL_L
@ LAL_SIM_INSPIRAL_FRAME_AXIS_ORBITAL_L
Set z-axis along the initial orbital angular momentum.
Definition: LALSimInspiralWaveformFlags.h:113
LAL_SIM_INSPIRAL_TIDAL_ORDER_5PN
@ LAL_SIM_INSPIRAL_TIDAL_ORDER_5PN
Definition: LALSimInspiralWaveformFlags.h:82
LAL_SIM_INSPIRAL_TIDAL_ORDER_6PN
@ LAL_SIM_INSPIRAL_TIDAL_ORDER_6PN
Definition: LALSimInspiralWaveformFlags.h:83
LAL_SIM_INSPIRAL_TIDAL_ORDER_ALL
@ LAL_SIM_INSPIRAL_TIDAL_ORDER_ALL
Definition: LALSimInspiralWaveformFlags.h:87
LAL_SIM_INSPIRAL_TIDAL_ORDER_0PN
@ LAL_SIM_INSPIRAL_TIDAL_ORDER_0PN
Definition: LALSimInspiralWaveformFlags.h:81
PhenSpinTaylor
@ PhenSpinTaylor
Inspiral part of the PhenSpinTaylorRD.
Definition: LALSimInspiral.h:288
PhenSpinTaylorRD
@ PhenSpinTaylorRD
Phenomenological waveforms, interpolating between a T4 spin-inspiral and the ringdown.
Definition: LALSimInspiral.h:290
XLALSimAddMode
int XLALSimAddMode(REAL8TimeSeries *hplus, REAL8TimeSeries *hcross, COMPLEX16TimeSeries *hmode, REAL8 theta, REAL8 phi, int l, int m, int sym)
Multiplies a mode h(l,m) by a spin-2 weighted spherical harmonic to obtain hplus - i hcross,...
Definition: LALSimSphHarmMode.c:44
r
static const INT4 r
m
static const INT4 m
XLALResizeREAL8TimeSeries
REAL8TimeSeries * XLALResizeREAL8TimeSeries(REAL8TimeSeries *series, int first, size_t length)
XLALCreateREAL8TimeSeries
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
XLALDestroyCOMPLEX16TimeSeries
void XLALDestroyCOMPLEX16TimeSeries(COMPLEX16TimeSeries *series)
XLALDestroyREAL8TimeSeries
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
XLALCreateCOMPLEX16TimeSeries
COMPLEX16TimeSeries * XLALCreateCOMPLEX16TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
lalDimensionlessUnit
const LALUnit lalDimensionlessUnit
XLALCreateCOMPLEX16Vector
COMPLEX16Vector * XLALCreateCOMPLEX16Vector(UINT4 length)
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLALDestroyCOMPLEX16Vector
void XLALDestroyCOMPLEX16Vector(COMPLEX16Vector *vector)
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
XLALClearErrno
int XLALClearErrno(void)
XLAL_EBADLEN
XLAL_EBADLEN
XLAL_ENOMEM
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_ERANGE
XLAL_ERANGE
XLAL_EFUNC
XLAL_EFUNC
XLAL_EDOM
XLAL_EDOM
XLAL_EINVAL
XLAL_EINVAL
XLAL_EFAILED
XLAL_EFAILED
XLAL_FAILURE
XLAL_FAILURE
XLALGPSAdd
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
x
list x
p
list p
y
list y
end
end
N
int N
alpha
double alpha
Definition: sgwb.c:106
COMPLEX16TimeSeries::data
COMPLEX16Sequence * data
COMPLEX16Vector::data
COMPLEX16 * data
LALSimInspiralInclAngle
Utility type and function to compute trigonometric functions from the cosine of an angle.
Definition: LALSimIMRPSpinInspiralRD.c:739
LALSimInspiralPhenSpinTaylorT4Coeffs::E4QMS2
REAL8 E4QMS2
non-dynamical S2^2 2PN quadrupole-monopole correction
Definition: LALSimIMRPSpinInspiralRD.c:96
LALSimInspiralPhenSpinTaylorT4Coeffs::fStart
REAL8 fStart
starting GW frequency of integration
Definition: LALSimIMRPSpinInspiralRD.c:111
LALSimInspiralPhenSpinTaylorT4Coeffs::wdotlogcoeff
REAL8 wdotlogcoeff
coefficient of log term in wdot
Definition: LALSimIMRPSpinInspiralRD.c:78
LALSimInspiralPhenSpinTaylorT4Coeffs::LNhat4SS
REAL8 LNhat4SS
non-dynamical 2PN SS correction
Definition: LALSimIMRPSpinInspiralRD.c:101
LALSimInspiralPhenSpinTaylorT4Coeffs::wdot4QMS1
REAL8 wdot4QMS1
non-dynamical S1^2 2PN quadrupole-monopole correction
Definition: LALSimIMRPSpinInspiralRD.c:86
LALSimInspiralPhenSpinTaylorT4Coeffs::Etidal5pn
REAL8 Etidal5pn
leading order tidal correction to energy
Definition: LALSimIMRPSpinInspiralRD.c:104
LALSimInspiralPhenSpinTaylorT4Coeffs::wdottidal5pn
REAL8 wdottidal5pn
leading order tidal correction
Definition: LALSimIMRPSpinInspiralRD.c:102
LALSimInspiralPhenSpinTaylorT4Coeffs::E4QMS1
REAL8 E4QMS1
non-dynamical S1^2 2PN quadrupole-monopole correction
Definition: LALSimIMRPSpinInspiralRD.c:94
LALSimInspiralPhenSpinTaylorT4Coeffs::Etidal6pn
REAL8 Etidal6pn
next to leading order tidal correction to energy
Definition: LALSimIMRPSpinInspiralRD.c:105
LALSimInspiralPhenSpinTaylorT4Coeffs::E4S1OS2O
REAL8 E4S1OS2O
non-dynamical 2PN SS correction
Definition: LALSimIMRPSpinInspiralRD.c:93
LALSimInspiralPhenSpinTaylorT4Coeffs::Enewt
REAL8 Enewt
coeffs. of PN corrections to energy
Definition: LALSimIMRPSpinInspiralRD.c:79
LALSimInspiralPhenSpinTaylorT4Coeffs::E4QMS1O
REAL8 E4QMS1O
non-dynamical (S1.L)^2 2PN quadrupole-monopole correction
Definition: LALSimIMRPSpinInspiralRD.c:95
LALSimInspiralPhenSpinTaylorT4Coeffs::wdot4QMS1O
REAL8 wdot4QMS1O
non-dynamical (S1.L)^2 2PN quadrupole-monopole correction
Definition: LALSimIMRPSpinInspiralRD.c:87
LALSimInspiralPhenSpinTaylorT4Coeffs::wdot4S1OS2O
REAL8 wdot4S1OS2O
non-dynamical 2PN SS correction
Definition: LALSimIMRPSpinInspiralRD.c:83
LALSimInspiralPhenSpinTaylorT4Coeffs::fEnd
REAL8 fEnd
ending GW frequency of integration
Definition: LALSimIMRPSpinInspiralRD.c:112
LALSimInspiralPhenSpinTaylorT4Coeffs::wdottidal6pn
REAL8 wdottidal6pn
next to leading order tidal correction
Definition: LALSimIMRPSpinInspiralRD.c:103
LALSimInspiralPhenSpinTaylorT4Coeffs::wdotnewt
REAL8 wdotnewt
leading order coefficient of wdot =
Definition: LALSimIMRPSpinInspiralRD.c:76
LALSimInspiralPhenSpinTaylorT4Coeffs::E4QMS2O
REAL8 E4QMS2O
non-dynamical (S2.L)^2 2PN quadrupole-monopole correction
Definition: LALSimIMRPSpinInspiralRD.c:97
LALSimInspiralPhenSpinTaylorT4Coeffs::wdot4S1S2
REAL8 wdot4S1S2
non-dynamical 2PN SS correction
Definition: LALSimIMRPSpinInspiralRD.c:82
LALSimInspiralPhenSpinTaylorT4Coeffs::wdot4QMS2
REAL8 wdot4QMS2
non-dynamical S2^2 2PN quadrupole-monopole correction
Definition: LALSimIMRPSpinInspiralRD.c:88
LALSimInspiralPhenSpinTaylorT4Coeffs::wdot4QMS2O
REAL8 wdot4QMS2O
non-dynamical (S2.L)^2 2PN quadrupole-monopole correction
Definition: LALSimIMRPSpinInspiralRD.c:89
REAL8Array::data
REAL8 * data
REAL8TimeSeries::data
REAL8Sequence * data
REAL8TimeSeries::epoch
LIGOTimeGPS epoch
REAL8Vector::data
REAL8 * data
params
Definition: burst.c:245
deltaT
double deltaT
Definition: unicorn.c:24