Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-b246709
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
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
65static REAL8 dsign(INT4 i){
66 if (i>=0) return 1.;
67 else return -1.;
68}
69
70typedef 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
115static 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
140static 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 */
160static 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
379static 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
566static 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
651static 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
690static 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
739typedef 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
762static 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
813static 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
870static 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
922static 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
1043static 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
1079static 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 }*/
1142static 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}
1153static 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
1165static 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
1306static 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. */
1369static 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
1385static 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
1605static 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 */
1673INT4 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
1956INT4 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 */
2036INT4 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:4232
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)
XLALAdaptiveRungeKutta4Init
LALAdaptiveRungeKuttaIntegrator * XLALAdaptiveRungeKutta4Init(int dim, int(*dydt)(double t, const double y[], double dydt[], void *params), int(*stop)(double t, const double y[], double dydt[], void *params), double eps_abs, double eps_rel)
XLALAdaptiveRungeKutta4Hermite
int XLALAdaptiveRungeKutta4Hermite(LALAdaptiveRungeKuttaIntegrator *integrator, void *params, REAL8 *yinit, REAL8 tinit, REAL8 tend_in, REAL8 deltat, REAL8Array **yout)
LAL_C_SI
#define LAL_C_SI
LAL_MSUN_SI
#define LAL_MSUN_SI
LAL_PI
#define LAL_PI
LAL_MTSUN_SI
#define LAL_MTSUN_SI
LAL_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
XLALDestroyCOMPLEX16TimeSeries
void XLALDestroyCOMPLEX16TimeSeries(COMPLEX16TimeSeries *series)
XLALCreateCOMPLEX16TimeSeries
COMPLEX16TimeSeries * XLALCreateCOMPLEX16TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
XLALDestroyREAL8TimeSeries
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
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)
lalDimensionlessUnit
const LALUnit lalDimensionlessUnit
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALCreateCOMPLEX16Vector
COMPLEX16Vector * XLALCreateCOMPLEX16Vector(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)
y
list y
end
end
N
int N
p
p
x
x
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