Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInspiral 5.0.3.1-5e288d3
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALPSpinInspiralRD.c
Go to the documentation of this file.
1
2/*
3* Copyright (C) 2011 Riccardo Sturani
4*
5* This program is free software; you can redistribute it and/or modify
6* it under the terms of the GNU General Public License as published by
7* the Free Software Foundation; either version 2 of the License, or
8* (at your option) any later version.
9*
10* This program is distributed in the hope that it will be useful,
11* but WITHOUT ANY WARRANTY; without even the implied warranty of
12* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13* GNU General Public License for more details.
14*
15* You should have received a copy of the GNU General Public License
16* along with with program; see the file COPYING. If not, write to the
17* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
18* MA 02110-1301 USA
19*/
20
21/**
22 * \defgroup LALPSpinInspiralRD_c Module LALPSpinInspiralRD.c
23 * \ingroup LALInspiral_h
24 *
25 * \brief Module to generate generic spinning binaries waveforms complete with ring-down
26 *
27 * ### Prototypes ###
28 *
29 * <tt>LALPSpinInspiralRD()</tt>
30 * <dl>
31 * <dt>status:</dt><dd> Input/Output
32 * <dt>signalvec:</dt><dd> Output containing the inspiral waveform.</dd>
33 * <dt>params:</dt><dd> Input containing binary chirp parameters.</dd>
34 * </dl>
35 *
36 * <tt>LALPSpinInspiralRDTemplates()</tt>
37 * <dl>
38 * <dt>status:</dt><dd>Input/Output
39 * <dt>signalvec1:</dt><dd>Output containing the \f$+\f$ inspiral waveform.</dd>
40 * <dt>signalvec2:</dt><dd>Output containing the \f$\times\f$ inspiral waveform.</dd>
41 * <dt>params:</dt><dd>Input containing binary chirp parameters.</dd>
42 * </dl>
43 *
44 * <tt>LALPSpinInspiralRDInjection()</tt>
45 * <dl>
46 * <dt>status:</dt><dd> Input/Output
47 * <dt>signalvec:</dt><dd>Output containing the inspiral waveform.</dd>
48 * <dt>params:</dt><dd>Input containing binary chirp parameters.</dd>
49 * </dl>
50 *
51 * <tt>LALPSpinInspiralRDFreqDom()</tt>
52 * <dl>
53 * <dt>status:</dt><dd> Input/Output
54 * <dt>signalvec:</dt><dd>Output containing the inspiral waveform in frequency domain.</dd>
55 * <dt>params:</dt><dd>Input containing binary chirp parameters.</dd>
56 * </dl>
57 *
58 * ### Description ###
59 *
60 * This codes provide complete waveforms for generically spinning binary systems.
61 * In order to construct the waveforms three phases are joined together:
62 * an initial inspiral phase, a phenomenological phase encompassing the description
63 * of the merger and the ring-down of the final black hole.
64 * During the inspiral phase the system is evolved according to the standard
65 * PN formulas, valid up to 3.5PN for the orbital motion,
66 * to 3PN level for spin-orbital momentum and to 2PN for spin-spin contributions
67 * to the orbital phase.
68 * Then a phenomenological phase is added during which the frequency of the
69 * waveform has a pole-like behaviour. The stitching is performed in order to
70 * ensure continuity of the phase, the frequency and its first and second
71 * derivatives. Finally a ring-down phase is attached.
72 *
73 * ### Algorithm ###
74 *
75 *
76 * ### Uses ###
77 *
78 * \code
79 * LALPSpinInspiralRDderivatives()
80 * LALInspiralSetup()
81 * LALInspiralChooseModel()
82 * LALRungeKutta4()
83 * OmMatch()
84 * fracRD()
85 * XLALPSpinInspiralRDSetParams()
86 * XLALSpinInspiralTest()
87 * XLALSpinInspiralDerivatives()
88 * LALSpinInspiralDerivatives()
89 * XLALSpinInspiralFillH2Modes()
90 * XLALSpinInspiralFillH3Modes()
91 * XLALSpinInspiralFillH4Modes()
92 * \endcode
93 *
94 * ### Notes ###
95 *
96 *//** @{ */
97
98#include <lal/LALAdaptiveRungeKuttaIntegrator.h>
99
100#include <lal/Units.h>
101#include <lal/LALInspiral.h>
102#include <lal/SeqFactories.h>
103#include <lal/NRWaveInject.h>
104#include <lal/RealFFT.h>
105
106/* use error codes above 1024 to avoid conflicts with GSL */
107#define LALPSIRDPN_TEST_ENERGY 1025
108#define LALPSIRDPN_TEST_OMEGADOT 1026
109#define LALPSIRDPN_TEST_OMEGANAN 1028
110#define LALPSIRDPN_TEST_OMEGAMATCH 1029
111#define LALPSIRDPN_TEST_OMEGANONPOS 1031
112#define LALPSIRDPN_TEST_OMEGACUT 1032
113
114#define sqrtOnePointFive 1.22474
115#define sqrtPoint15 0.387298
116#define sqrtFiveOver2 1.1183
117#define minIntLen 8
118
119#define UNUSED(expr) do { (void)(expr); } while (0)
120
121static REAL8 OmMatch(REAL8 LNhS1, REAL8 LNhS2, REAL8 S1S1, REAL8 S1S2, REAL8 S2S2) {
122
123 const REAL8 omM = 0.0555;
124 const REAL8 omMsz12 = 9.97e-4;
125 const REAL8 omMs1d2 = -2.032e-3;
126 const REAL8 omMssq = 5.629e-3;
127 const REAL8 omMsz1d2 = 8.646e-3;
128 const REAL8 omMszsq = -5.909e-3;
129 const REAL8 omM3s1d2 = 1.801e-3;
130 const REAL8 omM3ssq = -1.4059e-2;
131 const REAL8 omM3sz1d2 = 1.5483e-2;
132 const REAL8 omM3szsq = 8.922e-3;
133
134 return omM + /*6.05e-3 * sqrtOneMinus4Eta +*/
135 omMsz12 * (LNhS1 + LNhS2) +
136 omMs1d2 * (S1S2) +
137 omMssq * (S1S1 + S2S2) +
138 omMsz1d2 * (LNhS1 * LNhS2) +
139 omMszsq * (LNhS1 * LNhS1 + LNhS2 * LNhS2) +
140 omM3s1d2 * (LNhS1 + LNhS2) * (S1S2) +
141 omM3ssq * (LNhS1 + LNhS2) * (S1S1+S2S2) +
142 omM3sz1d2 * (LNhS1 + LNhS2) * (LNhS1*LNhS2) +
143 omM3szsq * (LNhS1 + LNhS2) * (LNhS1*LNhS1+LNhS2*LNhS2);
144}
145
146static REAL8 fracRD(REAL8 LNhS1, REAL8 LNhS2, REAL8 S1S1, REAL8 S1S2, REAL8 S2S2) {
147
148 const double frac0 = 0.840;
149 const double fracsz12 = -2.145e-2;
150 const double fracs1d2 = -4.421e-2;
151 const double fracssq = -2.643e-2;
152 const double fracsz1d2 = -5.876e-2;
153 const double fracszsq = -2.215e-2;
154
155 return frac0 +
156 fracsz12 * (LNhS1 + LNhS2) +
157 fracs1d2 * (S1S2) +
158 fracssq * (S1S1 + S2S2) +
159 fracsz1d2 * (LNhS1 * LNhS2) +
160 fracszsq * (LNhS1 * LNhS1 + LNhS2 * LNhS2);
161}
162
163typedef struct LALPSpinInspiralRDstructparams {
164 REAL8 dt;
165 REAL8 eta; ///< symmetric mass ratio
166 REAL8 dm; ///< \f$m_1-m_2\f$
167 REAL8 m1m2; ///< \f$m_1/m_2\f$
168 REAL8 m2m1; ///< \f$m_2/m_1\f$
169 REAL8 m1m;
170 REAL8 m2m;
171 REAL8 m1msq;
172 REAL8 m2msq;
173 REAL8 m;
174 REAL8 wdotorb[8]; ///< Coefficients of the analytic PN expansion of \f$\dot\omega_{orb}\f$
175 REAL8 wdotorblog; ///< Log coefficient of the PN expansion of of \f$\dot\omega_{orb}\f$
176 REAL8 wdotspin15S1LNh;
177 REAL8 wdotspin15S2LNh;
178 REAL8 wdotspin20S1S2;
179 REAL8 wdotspin20S1S1; ///< Coeff. of the \f$s_1s_1\f$ cntrb. to \f$\dot\omega\f$
180 REAL8 wdotspin20S1S2LNh;
181 REAL8 wdotspin25S1LNh;
182 REAL8 wdotspin25S2LNh; ///< Coeff. of the \f$s_2\cdot \hat L_N\f$ cntrb. to \f$\dot\omega\f$
183 REAL8 wdotspin30S1LNh;
184 REAL8 wdotspin30S2LNh;
185 REAL8 S1dot15;
186 REAL8 S2dot15;
187 REAL8 Sdot20;
188 REAL8 Sdot20S;
189 REAL8 S1dot25;
190 REAL8 S2dot25;
191 REAL8 LNhdot15;
192 REAL8 epnorb[4]; ///< Coefficients of the PN expansion of the energy
193 REAL8 epnspin15S1dotLNh; ///< Coeff. of the \f$S_1\cdot L\f$ term in energy
194 REAL8 epnspin15S2dotLNh; ///< Coeff. of the \f$S_2\cdot L\f$ term in energy
195 REAL8 epnspin20S1S2; ///< Coeff. of the \f$S_1\cdot S_2\f$ term in energy
196 REAL8 epnspin20S1S2dotLNh; ///< Coeff. of the \f$S_{1,2}\cdot L\f$ term in energy
197 REAL8 epnspin20S1S1; ///< Coeff. of the \f$S_1\cdot S_1\f$ term in energy
198 REAL8 epnspin20S1S1dotLNh;
199 REAL8 epnspin20S2S2; ///< Coeff. of the \f$S_2\cdot S_2\f$ term in energy
200 REAL8 epnspin20S2S2dotLNh;
201 REAL8 epnspin25S1dotLNh;
202 REAL8 epnspin25S2dotLNh;
203 REAL8 OmCutoff;
204 REAL8 lengths;
205 REAL8 omOffset;
206 UINT4 length;
207 UINT4 inspiralOnly;
208} LALPSpinInspiralRDparams;
209
210typedef struct LALPSpinInspiralPhenParsStruct {
211 REAL8 endtime;
212 REAL8 Psi;
213 REAL8 alpha;
214 REAL8 ci;
215 REAL8 omega;
216 REAL8 domega;
217 REAL8 ddomega;
218 REAL8 diota;
219 REAL8 ddiota;
220 REAL8 dalpha;
221 REAL8 ddalpha;
222 REAL8 energy;
223 REAL8 LNhS1;
224 REAL8 LNhS2;
225 REAL8 S1S1;
226 REAL8 S1S2;
227 REAL8 S2S2;
228 INT4 countback;
229 INT4 intreturn;
230} LALPSpinInspiralPhenPars;
231
232typedef struct LALSpinInspiralAngleStruct {
233 REAL8 ci2;
234 REAL8 si2;
235 REAL8 ci;
236 REAL8 si;
237 REAL8 c2i;
238 REAL8 s2i;
239 REAL8 c2i2;
240 REAL8 s2i2;
241 REAL8 c3i2;
242 REAL8 s3i2;
243 REAL8 c4i2;
244 REAL8 s4i2;
245 REAL8 c5i2;
246 REAL8 s5i2;
247 REAL8 c6i2;
248 REAL8 s6i2;
249 REAL8 c8i2;
250 REAL8 s8i2;
251 REAL8 cdi;
252 REAL8 sdi;
253} LALSpinInspiralAngle;
254
255static int XLALPSpinInspiralRDSetParams(LALPSpinInspiralRDparams *mparams,InspiralTemplate *params, InspiralInit *paramsInit) {
256
257 mparams->inspiralOnly = params->inspiralOnly;
258 mparams->dt = 1./params->tSampling;
259 mparams->OmCutoff = params->fCutoff*params->totalMass * LAL_MTSUN_SI * (REAL8) LAL_PI;
260 mparams->lengths = (5.0 / 256.0) / LAL_PI * pow(LAL_PI * params->chirpMass * LAL_MTSUN_SI * params->fLower,-5.0 / 3.0) / params->fLower;
261 mparams->omOffset = 0.006;
262
263 /* setup coefficients for PN equations */
264 mparams->m = params->totalMass;
265 mparams->m2m1 = params->mass2 / params->mass1;
266 mparams->m1m2 = params->mass1 / params->mass2;
267 mparams->m1m = params->mass1 / params->totalMass;
268 mparams->m2m = params->mass2 / params->totalMass;
269 mparams->m1msq = mparams->m1m * mparams->m1m;
270 mparams->m2msq = mparams->m2m * mparams->m2m;
271 mparams->dm = (params->mass1 - params->mass2) / params->totalMass;
272
273 /* params->eta might have been set up before but just for safety, we
274 * recompute it here below.*/
275 params->eta = (params->mass1 * params->mass2) / (params->mass1 + params->mass2) / (params->mass1 + params->mass2);
276 mparams->eta = params->eta;
277
278 switch (params->order) {
279
280 case LAL_PNORDER_THREE_POINT_FIVE:
281 mparams->wdotorb[7] = paramsInit->ak.ST[8];
282#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
283 __attribute__ ((fallthrough));
284#endif
285
286 case LAL_PNORDER_THREE:
287 mparams->epnorb[3] = paramsInit->ak.ETa3;
288 mparams->wdotorb[6] = paramsInit->ak.ST[6];
289 mparams->wdotorblog = paramsInit->ak.ST[7];
290 mparams->wdotspin30S1LNh = -LAL_PI/3. * ( 188. - 151./2./mparams->m1m);
291 mparams->wdotspin30S2LNh = -LAL_PI/3. * ( 188. + 151./2./mparams->m2m);
292#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
293 __attribute__ ((fallthrough));
294#endif
295
296 case LAL_PNORDER_TWO_POINT_FIVE:
297 mparams->wdotorb[5] = paramsInit->ak.ST[5];
298 mparams->epnspin25S1dotLNh = 8. - 31. / 9. * mparams->eta + (3. - 10. / 3. * mparams->eta) * mparams->m2m1;
299 mparams->epnspin25S2dotLNh = 8. - 31. / 9. * mparams->eta + (3. - 10. / 3. * mparams->eta) * mparams->m1m2;
300 mparams->wdotspin25S1LNh = -31319. / 1008. + 1159. / 24. * mparams->eta + (-809. / 84. + 281. / 8. * mparams->eta) * mparams->m2m1;
301 mparams->wdotspin25S2LNh = -31319. / 1008. + 1159. / 24. * mparams->eta + (-809. / 84. + 281. / 8. * mparams->eta) * mparams->m1m2;
302 mparams->S1dot25 = 0.5625 + 1.25 * mparams->eta - mparams->eta * mparams->eta / 24. + mparams->dm * (-0.5625 + 0.625 * mparams->eta);
303 mparams->S2dot25 = 0.5625 + 1.25 * mparams->eta - mparams->eta * mparams->eta / 24. - mparams->dm * (-0.5625 + 0.625 * mparams->eta);
304#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
305 __attribute__ ((fallthrough));
306#endif
307
308 case LAL_PNORDER_TWO:
309 mparams->epnorb[2] = paramsInit->ak.ETa2;
310 mparams->wdotorb[4] = paramsInit->ak.ST[4];
311 mparams->wdotspin20S1S2 = -(1.0 / 48.0) / mparams->eta;
312 mparams->wdotspin20S1S1 = 1. / 96.;
313 mparams->Sdot20 = 0.5;
314 mparams->Sdot20S = 0.5;
315 mparams->epnspin20S1S2 = 1. / mparams->eta;
316 mparams->epnspin20S1S2dotLNh = -3. / mparams->eta;
317 mparams->epnspin20S1S1 = (1. + mparams->m2m1) * (1. + mparams->m2m1) / 2.;
318 mparams->epnspin20S2S2 = (1. + mparams->m1m2) * (1. + mparams->m1m2) / 2.;
319 mparams->epnspin20S1S1dotLNh = -3. * (1. + mparams->m2m1) * (1. + mparams->m2m1) / 2.;
320 mparams->epnspin20S2S2dotLNh = -3. * (1. + mparams->m1m2) * (1. + mparams->m1m2) / 2.;
321#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
322 __attribute__ ((fallthrough));
323#endif
324
325 case LAL_PNORDER_ONE_POINT_FIVE:
326 mparams->wdotorb[3] = paramsInit->ak.ST[3];
327 mparams->epnspin15S1dotLNh = 8. / 3. + 2. * mparams->m2m1;
328 mparams->epnspin15S2dotLNh = 8. / 3. + 2. * mparams->m1m2;
329 mparams->wdotspin15S1LNh = -(113.0 + 75.0 * mparams->m2m1) / 12.0;
330 mparams->wdotspin15S2LNh = -(113.0 + 75.0 * mparams->m1m2) / 12.0;
331 mparams->LNhdot15 = 0.5;
332 mparams->S1dot15 = (4.0 + 3.0 * mparams->m2m1) / 2.0 * mparams->eta;
333 mparams->S2dot15 = (4.0 + 3.0 * mparams->m1m2) / 2.0 * mparams->eta;
334#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
335 __attribute__ ((fallthrough));
336#endif
337
338 case LAL_PNORDER_ONE:
339 mparams->epnorb[1] = paramsInit->ak.ETa1;
340 mparams->wdotorb[2] = paramsInit->ak.ST[2];
341#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
342 __attribute__ ((fallthrough));
343#endif
344
345 case LAL_PNORDER_HALF:
346 mparams->wdotorb[1] = paramsInit->ak.ST[1];
347#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
348 __attribute__ ((fallthrough));
349#endif
350
351 case LAL_PNORDER_NEWTONIAN:
352 mparams->epnorb[0] = paramsInit->ak.ETaN;
353 mparams->wdotorb[0] = paramsInit->ak.ST[0];
354 break;
355
356 case LAL_PNORDER_PSEUDO_FOUR:
357 XLALPrintError("*** LALPhenSpinInspiralRD ERROR: PhenSpin approximant not available at pseudo4PN order\n");
358 break;
359
360 case LAL_PNORDER_NUM_ORDER:
361 XLALPrintError("*** LALPhenSpinInspiralRD ERROR: NUM_ORDER not a valid PN order\n");
362 break;
363
364 default:
365 XLALPrintError("*** LALPhenSpinInspiralRD ERROR: Impossible to create waveform with %d order\n",params->order);
366 break;
367 }
368
369 switch (params->interaction) {
370
371 case LAL_INSPIRAL_INTERACTION_NONE:
372 /*This kills all spin effects in the phase. Still there are spin effects
373 in the waveform due to orbital plane precession*/
374 mparams->epnspin15S1dotLNh = 0.;
375 mparams->epnspin15S2dotLNh = 0.;
376 mparams->wdotspin15S1LNh = 0.;
377 mparams->wdotspin15S2LNh = 0.;
378 mparams->S1dot15 = 0.;
379 mparams->S2dot15 = 0.;
380#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
381 __attribute__ ((fallthrough));
382#endif
383
384 case LAL_INSPIRAL_INTERACTION_SPIN_ORBIT_15PN:
385 /* This keeps only the leading spin-orbit interactions*/
386 mparams->wdotspin20S1S2 = 0.;
387 mparams->epnspin20S1S2 = 0.;
388 mparams->epnspin20S1S2dotLNh = 0.;
389#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
390 __attribute__ ((fallthrough));
391#endif
392
393 case LAL_INSPIRAL_INTERACTION_SPIN_SPIN_2PN:
394 /* This keeps S1-S2 interactions but kill spin self-interactions*/
395 mparams->wdotspin20S1S1 = 0.;
396 mparams->epnspin20S1S1 = 0.;
397 mparams->epnspin20S2S2 = 0.;
398 mparams->Sdot20S = 0.;
399 mparams->epnspin20S1S1 = 0.;
400 mparams->epnspin20S2S2 = 0.;
401 mparams->epnspin20S1S1dotLNh = 0.;
402 mparams->epnspin20S2S2dotLNh = 0.;
403#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
404 __attribute__ ((fallthrough));
405#endif
406
407 case LAL_INSPIRAL_INTERACTION_SPIN_SPIN_SELF_2PN:
408 /* This kills all spin interaction intervening at 2.5PN order or higher*/
409 mparams->epnspin25S1dotLNh = 0.;
410 mparams->epnspin25S2dotLNh = 0.;
411 mparams->wdotspin25S1LNh = 0.;
412 mparams->wdotspin25S2LNh = 0.;
413 mparams->S1dot25 = 0.;
414 mparams->S2dot25 = 0.;
415#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
416 __attribute__ ((fallthrough));
417#endif
418
419 case LAL_INSPIRAL_INTERACTION_QUAD_MONO_2PN:
420
421 case LAL_INSPIRAL_INTERACTION_SPIN_ORBIT_25PN:
422 mparams->wdotspin30S1LNh = 0.;
423 mparams->wdotspin30S2LNh = 0.;
424#if __GNUC__ >= 7 && !defined __INTEL_COMPILER
425 __attribute__ ((fallthrough));
426#endif
427
428 case LAL_INSPIRAL_INTERACTION_SPIN_ORBIT_3PN:
429
430 case LAL_INSPIRAL_INTERACTION_ALL_SPIN:
431
432 case LAL_INSPIRAL_INTERACTION_TIDAL_5PN:
433
434 case LAL_INSPIRAL_INTERACTION_TIDAL_6PN:
435
436 case LAL_INSPIRAL_INTERACTION_ALL:
437
438 default:
439 break;
440 }
441
442 return XLAL_SUCCESS;
443}
444
445/*
446 *
447 * function to set derivatives: values and mparams input, dvalues output
448 *
449 */
450
451static int XLALSpinInspiralTest(double t, const double values[], double dvalues[], void *mparams) {
452
453 LALPSpinInspiralRDparams *params = (LALPSpinInspiralRDparams *) mparams;
454 REAL8 omega;
455 REAL8 energy;
456 REAL8 denergy;
457 INT4 returnint;
458
459 UNUSED(t);
460 omega = values[1];
461 energy = values[11];
462 denergy = dvalues[11];
463
464 REAL8 LNhS1=(values[2]*values[5]+values[3]*values[6]+values[4]*values[7])/params->m1msq;
465 REAL8 LNhS2=(values[2]*values[8]+values[3]*values[9]+values[4]*values[10])/params->m2msq;
466 REAL8 S1sq=(values[5]*values[5]+values[6]*values[6]+values[7]*values[7])/params->m1msq/params->m1msq;
467 REAL8 S2sq=(values[8]*values[8]+values[9]*values[9]+values[10]*values[10])/params->m2msq/params->m2msq;
468 REAL8 S1S2=(values[5]*values[8]+values[6]*values[9]+values[7]*values[10])/params->m1msq/params->m2msq;
469
470 REAL8 omegaMatch=OmMatch(LNhS1,LNhS2,S1sq,S1S2,S2sq)+params->omOffset;
471
472 if ( (energy > 0.0) || (( denergy > - 0.01*energy/params->dt*params->m*LAL_MTSUN_SI )&&(energy<0.) ) ) {
473 /*energy increase*/
474 /*XLALPrintWarning("** LALPSpinInspiralRD WARNING **: Energy increases dE %12.6e E %12.6e m/M:(%12.4e, %12.4e) om %12.6e vs. omM %12.6e\n",denergy, 0.01*energy, params->m1m, params->m2m, omega, omegaMatch);*/
475 returnint= LALPSIRDPN_TEST_ENERGY;
476 }
477 else if (omega < 0.0) {
478 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: Omega has become -ve, this should lead to nan's \n");
479 returnint= LALPSIRDPN_TEST_OMEGANONPOS;
480 }
481 else if (dvalues[1] < 0.0) {
482 /* omegadot < 0 */
483 returnint= LALPSIRDPN_TEST_OMEGADOT;
484 }
485 else if (isnan(omega)) {
486 /* omega is nan */
487 returnint= LALPSIRDPN_TEST_OMEGANAN;
488 }
489 else if ((params->inspiralOnly==1)&&(omega>params->OmCutoff)) {
490 returnint= LALPSIRDPN_TEST_OMEGACUT;
491 }
492 else if ((params->inspiralOnly!=1)&&(omega>omegaMatch)) {
493 returnint= LALPSIRDPN_TEST_OMEGAMATCH;
494 }
495 else
496 returnint= GSL_SUCCESS;
497
498 return returnint;
499
500}
501
502/**
503 * Function to compute detivative of dynamical variables
504 */
505static int XLALSpinInspiralDerivatives(double t, const double values[], double dvalues[], void *mparams) {
506
507 REAL8 omega; // time-derivative of the orbital phase
508 REAL8 LNhx, LNhy, LNhz; // orbital angolar momentum unit vector
509 REAL8 S1x, S1y, S1z; // dimension-less spin variable S/M^2
510 REAL8 S2x, S2y, S2z;
511 REAL8 LNhS1, LNhS2; // scalar products
512 REAL8 domega; // derivative of omega
513 REAL8 dLNhx, dLNhy, dLNhz; // derivatives of \f$\hat L_N\f$ components
514 REAL8 dS1x, dS1y, dS1z; // derivative of \f$S_i\f$
515 REAL8 dS2x, dS2y, dS2z;
516 REAL8 energy,energyold;
517
518 /* auxiliary variables*/
519 REAL8 S1S2, S1S1, S2S2; // Scalar products
520 REAL8 alphadotcosi; // alpha is the right ascension of L, i(iota) the angle between L and J
521 REAL8 v, v2, v3, v4, v5, v6, v7;
522 REAL8 tmpx, tmpy, tmpz, cross1x, cross1y, cross1z, cross2x, cross2y, cross2z, LNhxy;
523
524 LALPSpinInspiralRDparams *params= (LALPSpinInspiralRDparams *) mparams;
525
526 UNUSED(t);
527
528 /* --- computation start here --- */
529 omega = values[1];
530
531 LNhx = values[2];
532 LNhy = values[3];
533 LNhz = values[4];
534
535 S1x = values[5];
536 S1y = values[6];
537 S1z = values[7];
538
539 S2x = values[8];
540 S2y = values[9];
541 S2z = values[10];
542
543 energyold = values[11];
544
545 v = cbrt(omega);
546 v2 = v * v;
547 v3 = omega;
548 v4 = omega * v;
549 v5 = omega * v2;
550 v6 = omega * omega;
551 v7 = v6 * v;
552
553 // Omega derivative without spin effects up to 3.5 PN
554 // this formula does not include the 1.5PN shift mentioned in arXiv:0810.5336, five lines below (3.11)
555 domega = params->wdotorb[0]
556 + v * (params->wdotorb[1]
557 + v * (params->wdotorb[2]
558 + v * (params->wdotorb[3]
559 + v * (params->wdotorb[4]
560 + v * (params->wdotorb[5]
561 + v * (params->wdotorb[6] + params->wdotorblog * log(omega)
562 + v * params->wdotorb[7]))))));
563
564 energy = (1. + v2 * (params->epnorb[1] +
565 v2 * (params->epnorb[2] +
566 v2 * params->epnorb[3])));
567
568 // Adding spin effects
569 // L dot S1,2
570 LNhS1 = (LNhx * S1x + LNhy * S1y + LNhz * S1z);
571 LNhS2 = (LNhx * S2x + LNhy * S2y + LNhz * S2z);
572
573 // wdotspin15SiLNh = -1/12 (...)
574 domega += v3 * (params->wdotspin15S1LNh * LNhS1 + params->wdotspin15S2LNh * LNhS2); // see e.g. Buonanno et al. gr-qc/0211087
575
576 energy += v3 * (params->epnspin15S1dotLNh * LNhS1 + params->epnspin15S2dotLNh * LNhS2); // see e.g. Blanchet et al. gr-qc/0605140
577
578 // wdotspin20S1S1 = -1/48 eta
579 S1S1 = (S1x * S1x + S1y * S1y + S1z * S1z);
580 S2S2 = (S2x * S2x + S2y * S2y + S2z * S2z);
581 S1S2 = (S1x * S2x + S1y * S2y + S1z * S2z);
582 domega += params->wdotspin20S1S2 * v4 * (247.0 * S1S2 - 721.0 * LNhS1 * LNhS2); // see e.g. Buonanno et al. arXiv:0810.5336
583 domega += params->wdotspin20S1S1 * v4 * (719. * (LNhS1 * LNhS1 + LNhS2 * LNhS2) - 233. * (S1S1 + S2S2)); // see Racine et al. arXiv:0812.4413
584
585 energy += v4 * (params->epnspin20S1S2 * S1S2 + params->epnspin20S1S2dotLNh * LNhS1 * LNhS2); // see e.g. Buonanno et al. as above
586 energy += v4 * (params->epnspin20S1S1 * S1S1 + params->epnspin20S2S2 * S2S2 + params->epnspin20S1S1dotLNh * LNhS1 * LNhS1 + params->epnspin20S2S2 * LNhS2 * LNhS2); // see Racine et al. as above
587
588 // wdotspin25SiLNh = see below
589 domega += v5 * (params->wdotspin25S1LNh * LNhS1 + params->wdotspin25S2LNh * LNhS2); //see (8.3) of Blanchet et al.
590 energy += v5 * (params->epnspin25S1dotLNh * LNhS1 + params->epnspin25S2dotLNh * LNhS2); //see (7.9) of Blanchet et al.
591
592 domega += v6 * (params->wdotspin30S1LNh * LNhS1 + params->wdotspin30S2LNh * LNhS2); // see (6.5) of arXiv:1104.5659
593
594 // Setting the right pre-factor
595 domega *= 96. / 5. * params->eta * v5 * omega* omega;
596
597 energy *= params->epnorb[0] * v2;
598
599 /*Derivative of the angular momentum and spins */
600
601 /* dS1, 1.5PN */
602 /* S1dot15= (4+3m2/m1)/2 * eta */
603 cross1x = (LNhy * S1z - LNhz * S1y);
604 cross1y = (LNhz * S1x - LNhx * S1z);
605 cross1z = (LNhx * S1y - LNhy * S1x);
606
607 dS1x = params->S1dot15 * v5 * cross1x;
608 dS1y = params->S1dot15 * v5 * cross1y;
609 dS1z = params->S1dot15 * v5 * cross1z;
610
611 /* dS1, 2PN */
612 /* Sdot20= 0.5 */
613 tmpx = S1z * S2y - S1y * S2z;
614 tmpy = S1x * S2z - S1z * S2x;
615 tmpz = S1y * S2x - S1x * S2y;
616
617 // S1S2 contribution see. eq. 2.23 of arXiv:0812.4413
618 dS1x += params->Sdot20 * v6 * (tmpx - 3. * LNhS2 * cross1x);
619 dS1y += params->Sdot20 * v6 * (tmpy - 3. * LNhS2 * cross1y);
620 dS1z += params->Sdot20 * v6 * (tmpz - 3. * LNhS2 * cross1z);
621 // S1S1 contribution
622 dS1x -= 3. * params->Sdot20S * v6 * LNhS1 * cross1x * (1. + params->m2m1) * params->m2m;
623 dS1y -= 3. * params->Sdot20S * v6 * LNhS1 * cross1y * (1. + params->m2m1) * params->m2m;
624 dS1z -= 3. * params->Sdot20S * v6 * LNhS1 * cross1z * (1. + params->m2m1) * params->m2m;
625
626 // dS1, 2.5PN, eq. 7.8 of Blanchet et al. gr-qc/0605140
627 // S1dot25= 9/8-eta/2.+eta+mparams->eta*29./24.+mparams->m1m2*(-9./8.+5./4.*mparams->eta)
628 dS1x += params->S1dot25 * v7 * cross1x;
629 dS1y += params->S1dot25 * v7 * cross1y;
630 dS1z += params->S1dot25 * v7 * cross1z;
631
632 /* dS2, 1.5PN */
633 cross2x = (LNhy * S2z - LNhz * S2y);
634 cross2y = (LNhz * S2x - LNhx * S2z);
635 cross2z = (LNhx * S2y - LNhy * S2x);
636
637 dS2x = params->S2dot15 * v5 * cross2x;
638 dS2y = params->S2dot15 * v5 * cross2y;
639 dS2z = params->S2dot15 * v5 * cross2z;
640
641 /* dS2, 2PN */
642 dS2x += params->Sdot20 * v6 * (-tmpx - 3.0 * LNhS1 * cross2x);
643 dS2y += params->Sdot20 * v6 * (-tmpy - 3.0 * LNhS1 * cross2y);
644 dS2z += params->Sdot20 * v6 * (-tmpz - 3.0 * LNhS1 * cross2z);
645 // S2S2 contribution
646 dS2x -= 3. * params->Sdot20S * v6 * LNhS2 * cross2x * (1. + params->m1m2) * params->m1m;
647 dS2y -= 3. * params->Sdot20S * v6 * LNhS2 * cross2y * (1. + params->m1m2) * params->m1m;
648 dS2z -= 3. * params->Sdot20S * v6 * LNhS2 * cross2z * (1. + params->m1m2) * params->m1m;
649
650 // dS2, 2.5PN, eq. 7.8 of Blanchet et al. gr-qc/0605140
651 dS2x += params->S2dot25 * v7 * cross2x;
652 dS2y += params->S2dot25 * v7 * cross2y;
653 dS2z += params->S2dot25 * v7 * cross2z;
654
655 dLNhx = -(dS1x + dS2x) * v / params->eta;
656 dLNhy = -(dS1y + dS2y) * v / params->eta;
657 dLNhz = -(dS1z + dS2z) * v / params->eta;
658
659 /* dphi */
660 LNhxy = LNhx * LNhx + LNhy * LNhy;
661
662 if (LNhxy > 0.0)
663 alphadotcosi = LNhz * (LNhx * dLNhy - LNhy * dLNhx) / LNhxy;
664 else
665 {
666 //XLALPrintWarning("*** LALPSpinInspiralRD WARNING ***: alphadot set to 0 by hand LNh:(%12.4e %12.4e %12.4e)\n",LNhx,LNhy,LNhz);
667 alphadotcosi = 0.;
668 }
669
670 /* dvalues->data[0] is the phase derivative */
671 /* omega is the derivative of the orbital phase omega \neq dvalues->data[0] */
672 dvalues[0] = omega - alphadotcosi;
673 dvalues[1] = domega;
674
675 dvalues[2] = dLNhx;
676 dvalues[3] = dLNhy;
677 dvalues[4] = dLNhz;
678
679 dvalues[5] = dS1x;
680 dvalues[6] = dS1y;
681 dvalues[7] = dS1z;
682
683 dvalues[8] = dS2x;
684 dvalues[9] = dS2y;
685 dvalues[10] = dS2z;
686
687 dvalues[11] = (energy-energyold)/params->dt*params->m*LAL_MTSUN_SI;
688
689 return GSL_SUCCESS;
690
691} /* end of XLALSpinInspiralDerivatives */
692
693
694static void LALSpinInspiralDerivatives(REAL8Vector * values, REAL8Vector * dvalues, void *mparams)
695{
696 XLALSpinInspiralDerivatives(0.,values->data,dvalues->data,mparams);
697} /* end of LALSpinInspiralDerivatives */
698
699/**
700 * Main function to produce waveforms
701 */
702static int XLALPSpinInspiralRDEngine(
703 REAL8Vector * signalvec1,
704 REAL8Vector * signalvec2,
705 REAL8Vector * hh,
706 REAL8Vector * ff,
707 REAL8Vector * phi,
708 InspiralTemplate *params,
709 InspiralInit *paramsInit);
710
711
712int XLALPSpinInspiralRD(REAL4Vector * signalvec, InspiralTemplate * params)
713{
714 INT4 count=0;
715 InspiralInit paramsInit;
716
717 if(!signalvec) XLAL_ERROR(XLAL_EFAULT);
718 if(!signalvec->data) XLAL_ERROR(XLAL_EFAULT);
719 if(!params) XLAL_ERROR(XLAL_EFAULT);
720 if(params->nStartPad<0) XLAL_ERROR(XLAL_EBADLEN);
721 if(params->nEndPad<0) XLAL_ERROR(XLAL_EBADLEN);
722 if(params->fLower <=0) XLAL_ERROR(XLAL_EFREQ);
723 if(params->tSampling<=0) XLAL_ERROR(XLAL_ETIME);
724 if(params->totalMass<=0) XLAL_ERROR(XLAL_EDOM);
725
726 if (XLALInspiralSetup(&(paramsInit.ak), params))
727 XLAL_ERROR(XLAL_EFUNC);
728 if (XLALInspiralChooseModel(&(paramsInit.func),&(paramsInit.ak), params))
729 XLAL_ERROR(XLAL_EFUNC);
730
731 REAL8Vector *s=XLALCreateREAL8Vector(signalvec->length);
732 memset(s->data, 0, s->length * sizeof(REAL8));
733
734 /* Call the engine function */
735 count = XLALPSpinInspiralRDEngine(s, NULL, NULL, NULL, NULL, params, &paramsInit);
736 if (count == XLAL_FAILURE)
737 XLAL_ERROR(XLAL_EFUNC);
738
739 UINT4 i;
740 for (i=0;i<s->length;i++) {
741 signalvec->data[i]=(REAL4) s->data[i];
742 }
743
744 return XLAL_SUCCESS;
745}
746
747/**
748 * Function to produce waveform templates
749 */
750
751int XLALPSpinInspiralRDTemplates(
752 REAL4Vector * signalvec1,
753 REAL4Vector * signalvec2,
754 InspiralTemplate * params)
755{
756 INT4 count=0;
757 InspiralInit paramsInit;
758
759 if(!signalvec1) XLAL_ERROR(XLAL_EFAULT);
760 if(!signalvec1->data) XLAL_ERROR(XLAL_EFAULT);
761 if(!signalvec2) XLAL_ERROR(XLAL_EFAULT);
762 if(!signalvec2->data) XLAL_ERROR(XLAL_EFAULT);
763 if(!params) XLAL_ERROR(XLAL_EFAULT);
764 if(params->nStartPad<0) XLAL_ERROR(XLAL_EBADLEN);
765 if(params->nEndPad<0) XLAL_ERROR(XLAL_EBADLEN);
766 if(params->fLower <=0) XLAL_ERROR(XLAL_EFREQ);
767 if(params->tSampling<=0) XLAL_ERROR(XLAL_ETIME);
768 if(params->totalMass<=0) XLAL_ERROR(XLAL_EDOM);
769
770 if (XLALInspiralSetup(&(paramsInit.ak), params))
771 XLAL_ERROR(XLAL_EFUNC);
772 if (XLALInspiralChooseModel(&(paramsInit.func), &(paramsInit.ak), params))
773 XLAL_ERROR(XLAL_EFUNC);
774
775 REAL8Vector* s1=XLALCreateREAL8Vector(signalvec1->length);
776 REAL8Vector* s2=XLALCreateREAL8Vector(signalvec2->length);
777
778 memset(s1->data, 0, signalvec1->length * sizeof(REAL8));
779 memset(s2->data, 0, signalvec2->length * sizeof(REAL8));
780
781 /* Call the engine function*/
782 count = XLALPSpinInspiralRDEngine(s1, s2, NULL, NULL, NULL, params, &paramsInit);
783 if (count == XLAL_FAILURE)
784 XLAL_ERROR(XLAL_EFUNC);
785
786 UINT4 i;
787 for (i=0;i<s1->length;i++) {
788 signalvec1->data[i]=(REAL4) s1->data[i];
789 }
790 for (i=0;i<s2->length;i++) {
791 signalvec2->data[i]=(REAL4) s2->data[i];
792 }
793
794 XLALDestroyREAL8Vector(s1);
795 XLALDestroyREAL8Vector(s2);
796
797 return XLAL_SUCCESS;
798}
799
800/**
801 * Function Module to produce injection waveforms
802 */
803
804int XLALPSpinInspiralRDForInjection(
805 CoherentGW * waveform,
806 InspiralTemplate * params,
807 PPNParamStruc * ppnParams)
808{
809 REAL8Vector *h = NULL; /* pointers to generated amplitude data */
810 REAL8Vector *f = NULL; /* pointers to generated frequency data */
811 REAL8Vector *phi = NULL; /* pointer to generated phase data */
812
813 InspiralInit paramsInit;
814 INT4 nbins,count;
815 UINT4 i;
816
817 /* Make sure parameter and waveform structures exist. */
818 if(!params) XLAL_ERROR(XLAL_EFAULT);
819 if(!waveform) XLAL_ERROR(XLAL_EFAULT);
820 /* Make sure waveform fields don't exist. */
821 if(waveform->a) XLAL_ERROR(XLAL_EFAULT);
822 if(waveform->f) XLAL_ERROR(XLAL_EFAULT);
823 if(waveform->phi) XLAL_ERROR(XLAL_EFAULT);
824 if(waveform->h) XLAL_ERROR(XLAL_EFAULT);
825
826 /* Compute some parameters */
827 XLALInspiralInit(params, &paramsInit);
828 if (paramsInit.nbins == 0) {
829 return XLAL_SUCCESS;
830 }
831
832 nbins = paramsInit.nbins;
833
834 /* Now we can allocate memory and vector for coherentGW structure */
835 f = XLALCreateREAL8Vector(nbins);
836 h = XLALCreateREAL8Vector(2 * nbins);
837 phi = XLALCreateREAL8Vector(nbins);
838
839 /* By default the waveform is empty */
840 memset(f->data, 0, nbins * sizeof(REAL8));
841 memset(h->data, 0, 2 * nbins * sizeof(REAL8));
842 memset(phi->data, 0, nbins * sizeof(REAL8));
843
844 /* Call the engine function */
845 count = XLALPSpinInspiralRDEngine(NULL, NULL, h, f, phi, params, &paramsInit);
846 if (count == XLAL_FAILURE)
847 XLAL_ERROR(XLAL_EFUNC);
848
849 /* Check an empty waveform hasn't been returned */
850 for (i = 0; i < phi->length; i++) {
851 if (phi->data[i] != 0.0)
852 break;
853 if (i == phi->length - 1) {
854 XLALDestroyREAL8Vector(f);
855 XLALDestroyREAL8Vector(h);
856 XLALDestroyREAL8Vector(phi);
857 }
858 }
859
860 /* Allocate the waveform structures. */
861 if ((waveform->h = (REAL4TimeVectorSeries *)
862 LALMalloc(sizeof(REAL4TimeVectorSeries))) == NULL) {
863 XLAL_ERROR(XLAL_ENOMEM);
864 }
865 memset(waveform->h, 0, sizeof(REAL4TimeVectorSeries));
866
867 if ((waveform->f = (REAL4TimeSeries *)
868 LALMalloc(sizeof(REAL4TimeSeries))) == NULL) {
869 LALFree(waveform->h);
870 waveform->a = NULL;
871 XLAL_ERROR(XLAL_ENOMEM);
872 }
873 memset(waveform->f, 0, sizeof(REAL4TimeSeries));
874
875 if ((waveform->phi = (REAL8TimeSeries *)
876 LALMalloc(sizeof(REAL8TimeSeries))) == NULL) {
877 LALFree(waveform->h);
878 waveform->h = NULL;
879 LALFree(waveform->f);
880 waveform->f = NULL;
881 XLAL_ERROR(XLAL_ENOMEM);
882 }
883 memset(waveform->phi, 0, sizeof(REAL8TimeSeries));
884
885 waveform->h->data = XLALCreateREAL4VectorSequence(count, 2);
886 waveform->f->data = XLALCreateREAL4Vector(count);
887 waveform->phi->data = XLALCreateREAL8Vector(count);
888
889 for (i=0; i<(UINT4)count; i++) {
890 waveform->f->data->data[i] =(REAL4) f->data[i];
891 waveform->h->data->data[2*i] =(REAL4) h->data[2*i];
892 waveform->h->data->data[2*i+1] =(REAL4) h->data[2*i+1];
893 waveform->phi->data->data[i] = phi->data[i];
894 }
895
896 waveform->h->deltaT = waveform->f->deltaT = waveform->phi->deltaT = 1. / params->tSampling;
897
898 waveform->h->sampleUnits = lalStrainUnit;
899 waveform->f->sampleUnits = lalHertzUnit;
900 waveform->phi->sampleUnits = lalDimensionlessUnit;
901
902 waveform->position = ppnParams->position;
903 waveform->psi = ppnParams->psi;
904
905 snprintf(waveform->h->name, LALNameLength, "PSpinInspiralRD amplitudes");
906 snprintf(waveform->f->name, LALNameLength, "PSpinInspiralRD frequency");
907 snprintf(waveform->phi->name, LALNameLength,"PSpinInspiralRD main phase");
908
909 /* --- fill some output --- */
910 ppnParams->tc = params->tC;
911 ppnParams->length = count;
912 ppnParams->dfdt = ((REAL4) (waveform->f->data->data[count - 1]- waveform->f->data->data[count - 2])) * ppnParams->deltaT;
913 ppnParams->fStop = params->fFinal;
914 ppnParams->termCode = GENERATEPPNINSPIRALH_EFSTOP;
915 ppnParams->termDescription = GENERATEPPNINSPIRALH_MSGEFSTOP;
916
917 ppnParams->fStart = ppnParams->fStartIn;
918
919 /* --- free memory --- */
920
921 XLALDestroyREAL8Vector(f);
922 XLALDestroyREAL8Vector(h);
923 XLALDestroyREAL8Vector(phi);
924
925 return XLAL_SUCCESS;
926} /* End LALPSpinInspiralRDForInjection */
927
928int XLALPSpinInspiralRDFreqDom(
929 REAL4Vector * signalvec,
930 InspiralTemplate * params)
931{
932 REAL8Vector *tsignalvec = NULL;
933 REAL4Vector *fsignalvec = NULL;
934 REAL4FFTPlan *forwPlan = NULL;
935
936 InspiralInit paramsInit;
937
938 INT4 count;
939 UINT4 nbins,idx;
940 UINT4 length = signalvec->length;
941
942 if(!signalvec) XLAL_ERROR(XLAL_EFAULT);
943 if(!signalvec->data) XLAL_ERROR(XLAL_EFAULT);
944 if(!params) XLAL_ERROR(XLAL_EFAULT);
945 if(params->nStartPad<0) XLAL_ERROR(XLAL_EBADLEN);
946 if(params->nEndPad<0) XLAL_ERROR(XLAL_EBADLEN);
947 if(params->fLower <=0) XLAL_ERROR(XLAL_EFREQ);
948 if(params->tSampling<=0) XLAL_ERROR(XLAL_ETIME);
949 if(params->totalMass<=0) XLAL_ERROR(XLAL_EDOM);
950
951 XLALInspiralInit(params, &paramsInit);
952
953 if (paramsInit.nbins == 0) {
954 return XLAL_SUCCESS;
955 }
956 nbins = paramsInit.nbins;
957 if (nbins < length) nbins = length;
958 tsignalvec = XLALCreateREAL8Vector(nbins);
959 fsignalvec = XLALCreateREAL4Vector(nbins);
960
961 memset(signalvec->data, 0, length * sizeof(REAL4));
962 memset(tsignalvec->data, 0, nbins * sizeof(REAL8));
963 memset(fsignalvec->data, 0, nbins * sizeof(REAL4));
964
965 /* Call the engine function */
966 count = XLALPSpinInspiralRDEngine(tsignalvec, NULL, NULL, NULL, NULL, params, &paramsInit);
967 if (count == XLAL_FAILURE)
968 XLAL_ERROR(XLAL_EFUNC);
969
970 REAL4Vector *tsigR4=XLALCreateREAL4Vector(nbins);
971 for (idx=0;idx<nbins;idx++) tsigR4->data[idx]=(REAL4) tsignalvec->data[idx];
972 XLALDestroyREAL8Vector(tsignalvec);
973 XLALSimInspiralREAL4WaveTaper(tsigR4, (LALSimInspiralApplyTaper) 3);
974
975 forwPlan = XLALCreateForwardREAL4FFTPlan(nbins, 0);
976 if (forwPlan == NULL) {
977 XLALDestroyREAL4Vector(fsignalvec);
978 XLALDestroyREAL8Vector(tsignalvec);
979 XLALDestroyREAL4Vector(tsigR4);
980 XLALDestroyREAL4FFTPlan(forwPlan);
981 XLAL_ERROR(XLAL_ENOMEM);
982 }
983 XLALREAL4VectorFFT(fsignalvec, tsigR4, forwPlan);
984 XLALDestroyREAL4Vector(tsigR4);
985
986 for (idx = 0; idx < nbins; idx++) fsignalvec->data[idx] /= params->tSampling;
987
988 if (nbins > length) {
989 /*do interpolation*/
990 REAL8Vector *fsigRe = XLALCreateREAL8Vector(nbins/2);
991 REAL8Vector *fsigIm = XLALCreateREAL8Vector(nbins/2);
992 REAL8Vector *freq = XLALCreateREAL8Vector(length/2);
993 REAL8Vector *freqSup = XLALCreateREAL8Vector(nbins/2);
994
995 REAL8 dF = 1./params->tSampling/((REAL8) length);
996 REAL8 dFsup = 1./params->tSampling/((REAL8) nbins);
997
998 for (idx=0; idx< length/2; idx++) freq->data[idx]=((REAL4) idx) *dF;
999
1000 for (idx = 0; idx < nbins/2; idx++) {
1001 fsigRe->data[idx]=(REAL8)fsignalvec->data[idx];
1002 fsigIm->data[idx]=(REAL8)fsignalvec->data[nbins-idx];
1003 freqSup->data[idx]=((REAL8) idx) *dFsup;
1004 }
1005
1006 gsl_interp_accel* acc;
1007 gsl_spline* spline_real;
1008 gsl_spline* spline_imag;
1009 XLAL_CALLGSL( spline_imag = (gsl_spline*) gsl_spline_alloc(gsl_interp_cspline, nbins/2));
1010 XLAL_CALLGSL( spline_real = (gsl_spline*) gsl_spline_alloc(gsl_interp_cspline, nbins/2));
1011 XLAL_CALLGSL( acc = (gsl_interp_accel*) gsl_interp_accel_alloc());
1012 XLAL_CALLGSL( gsl_spline_init(spline_real, freqSup->data, fsigRe->data, nbins/2));
1013 XLAL_CALLGSL( gsl_spline_init(spline_imag, freqSup->data, fsigIm->data, nbins/2));
1014
1015 for (idx=0;idx<length/2;idx++){
1016 signalvec->data[idx] = gsl_spline_eval ( spline_real , freq->data[idx] , acc);
1017 signalvec->data[length-idx] = gsl_spline_eval ( spline_imag , freq->data[idx] , acc);
1018 }
1019 signalvec->data[0] = 0.;
1020 signalvec->data[length / 2] = 0.;
1021
1022 XLALDestroyREAL8Vector(fsigRe);
1023 XLALDestroyREAL8Vector(fsigIm);
1024 XLALDestroyREAL8Vector(freq);
1025 XLALDestroyREAL8Vector(freqSup);
1026 gsl_spline_free(spline_real);
1027 gsl_spline_free(spline_imag);
1028 gsl_interp_accel_free(acc);
1029
1030 }
1031 else {
1032 for (idx=0; idx<length; idx++)
1033 signalvec->data[idx]=fsignalvec->data[idx];
1034 }
1035
1036 XLALDestroyREAL4Vector(fsignalvec);
1037 XLALDestroyREAL4FFTPlan(forwPlan);
1038
1039 return XLAL_SUCCESS;
1040}
1041
1042/**
1043 * Function actually computing PSIRD waveforms
1044 */
1045
1046static int XLALSpinInspiralFillH2Modes(
1047 REAL8Vector* h2P2,
1048 REAL8Vector* h2M2,
1049 REAL8Vector* h2P1,
1050 REAL8Vector* h2M1,
1051 REAL8Vector* h20,
1052 UINT4 j,
1053 REAL4 amp,
1054 REAL4 v,
1055 REAL4 eta,
1056 REAL4 dm,
1057 REAL8 Psi,
1058 REAL8 alpha,
1059 LALSpinInspiralAngle *an
1060 ){
1061
1062 REAL8 amp20 = amp * sqrtOnePointFive;
1063 REAL8 v2 = v*v;
1064 REAL8 omega = v*v*v;
1065 const REAL8 omegaC = 0.05;
1066 const REAL8 Afac = .5;
1067 REAL8 damp = omega > omegaC ? 1. : exp(-Afac*(1.-omegaC/omega)*(1.-omegaC/omega));
1068
1069 h2P2->data[2 * j] = amp * ( ( 1. - damp * v2 / 42. * (107. - 55. * eta) )*
1070 ( cos(2. * (Psi + alpha)) * an->c4i2 + cos(2. * (Psi - alpha)) * an->s4i2) +
1071 v * dm / 3. * an->si *
1072 ( cos(Psi - 2. * alpha) * an->s2i2 + cos(Psi + 2. * alpha) * an->c2i2 ) );
1073
1074 h2M2->data[2 * j] = amp * ( ( 1. - damp * v2 / 42. * (107. - 55. * eta) )*
1075 ( cos(2. * (Psi + alpha)) * an->c4i2 + cos(2. * (Psi - alpha)) * an->s4i2) -
1076 v * dm / 3. * an->si *
1077 ( cos(Psi - 2. * alpha) * an->s2i2 + cos(Psi + 2. * alpha) * an->c2i2) );
1078
1079 h2P2->data[2 * j + 1] = amp * ( (1. - damp * v2 / 42. * (107. - 55. * eta) )*
1080 ( -sin(2. * (Psi + alpha)) * an->c4i2 + sin(2. * (Psi - alpha)) * an->s4i2) +
1081 v * dm / 3. * an->si *
1082 ( sin(Psi - 2. * alpha) * an->s2i2 - sin(Psi + 2. * alpha) * an->c2i2) );
1083
1084 h2M2->data[2 * j + 1] = amp * ( (1. - damp * v2 / 42. * (107. - 55. * eta) )*
1085 ( sin(2. * (Psi + alpha)) * an->c4i2 - sin(2. * (Psi - alpha)) * an->s4i2) +
1086 v * dm / 3. * an->si *
1087 ( sin(Psi - 2. * alpha) * an->s2i2 - sin(Psi + 2. * alpha) * an->c2i2) );
1088
1089 h2P1->data[2 * j] = amp * (an->si * (1. - damp * v2 / 42. * (107. - 55. * eta) ) *
1090 ( -cos(2. * Psi - alpha) * an->s2i2 + cos(2. * Psi + alpha) * an->c2i2) +
1091 v * dm / 3. *
1092 (-cos(Psi + alpha) * (an->ci + an->cdi)/2. -
1093 cos(Psi - alpha) * an->s2i2 * (1. + 2. * an->ci) ) );
1094
1095 h2M1->data[2 * j] = amp * (an->si * (1. - damp * v2 / 42. * (107. - 55. * eta)) *
1096 ( cos(2. * Psi - alpha) * an->s2i2 - cos(2. * Psi + alpha) * an->c2i2) +
1097 v * dm / 3. *
1098 (-cos(Psi + alpha) * (an->ci + an->cdi)/2. -
1099 cos(Psi - alpha) * an->s2i2 * (1. + 2. * an->ci) ) );
1100
1101 h2P1->data[2 * j + 1] = amp * (an->si * (1. - damp * v2 / 42. * (107. - 55. * eta) ) *
1102 ( -sin(2. * Psi - alpha ) * an->s2i2 - sin(2. * Psi + alpha) * an->c2i2) +
1103 v * dm / 3. *
1104 (sin(Psi + alpha) * (an->ci + an->cdi)/2. -
1105 sin(Psi - alpha) * an->s2i2 * (1. + 2. * an->ci) ) );
1106
1107 h2M1->data[2 * j + 1] = amp * (an->si * (1. - damp * v2 / 42. * (107. - 55. * eta)) *
1108 ( -sin(2. * Psi - alpha) * an->s2i2 - sin(2. * Psi + alpha) * an->c2i2) -
1109 v * dm / 3. *
1110 (sin(Psi + alpha) * (an->ci + an->cdi) / 2. -
1111 sin(Psi - alpha) * an->s2i2 * (1. + 2. * an->ci) ) );
1112
1113 h20->data[2 * j] = amp20 * ( an->s2i * (1.- damp * v2/42. * (107. - 55.*eta) ) * cos(2. * Psi) );
1114
1115 h20->data[2 * j + 1] = amp20 * ( v * dm / 3. * an->sdi * sin(Psi) );
1116
1117 return 0;
1118
1119}
1120
1121static int XLALSpinInspiralFillH3Modes(
1122 REAL8Vector* h3P3,
1123 REAL8Vector* h3M3,
1124 REAL8Vector* h3P2,
1125 REAL8Vector* h3M2,
1126 REAL8Vector* h3P1,
1127 REAL8Vector* h3M1,
1128 REAL8Vector* h30,
1129 UINT4 j,
1130 REAL4 amp,
1131 REAL4 v,
1132 REAL4 eta,
1133 REAL4 dm,
1134 REAL8 Psi,
1135 REAL8 alpha,
1136 LALSpinInspiralAngle *an
1137 ){
1138
1139 REAL8 amp32 = amp * sqrtOnePointFive;
1140 REAL8 amp31 = amp * sqrtPoint15;
1141 REAL8 amp30 = amp / sqrtFiveOver2;
1142 REAL8 v2 = v*v;
1143
1144 h3P3->data[2 * j] = amp * ( v * dm *
1145 (-9. * cos(3. * (Psi - alpha)) * an->s6i2 -
1146 cos( Psi - 3. * alpha) * an->s4i2 * an->c2i2 +
1147 cos( Psi + 3. * alpha) * an->s2i2 * an->c4i2 +
1148 9. * cos(3. * (Psi + alpha)) * an->c6i2) +
1149 v2 * 4. * an->si * (1. - 3. * eta) *
1150 (-cos(2. * Psi - 3. * alpha) * an->s4i2 +
1151 cos(2. * Psi + 3. * alpha) * an->c4i2 ) );
1152
1153 h3M3->data[2 * j] = amp * (-v * dm *
1154 (-9. * cos(3. * (Psi - alpha)) * an->s6i2 -
1155 cos( Psi - 3. * alpha) * an->s4i2 * an->c2i2 +
1156 cos( Psi + 3. * alpha) * an->s2i2 * an->c4i2 +
1157 9. * cos(3. * (Psi + alpha)) * an->c6i2) +
1158 v2 * 4. * an->si * (1. - 3. * eta) *
1159 (-cos(2. * Psi - 3. * alpha) * an->s4i2 +
1160 cos(2. * Psi + 3. * alpha) * an->c4i2 ) );
1161
1162 h3P3->data[2 * j + 1] = amp * ( v * dm *
1163 (-9. * sin(3. * (Psi - alpha)) * an->s6i2 -
1164 sin( Psi - 3. * alpha) * an->s4i2 * an->c2i2 -
1165 sin( Psi + 3. * alpha) * an->s2i2 * an->c4i2 -
1166 9. * sin(3. * (Psi + alpha)) * an->c6i2) +
1167 v2 * 4. * an->si * (1. - 3. * eta) *
1168 (-sin(2. * Psi - 3. * alpha) * an->s4i2
1169 -sin(2. * Psi + 3. * alpha) * an->c4i2 ) );
1170
1171 h3M3->data[2 * j + 1] = amp * ( v * dm *
1172 (-9. * sin(3. * (Psi - alpha)) * an->s6i2 -
1173 sin( Psi - 3. * alpha) * an->s4i2 * an->c2i2 -
1174 sin( Psi + 3. * alpha) * an->s2i2 * an->c4i2 -
1175 9. * sin(3. * (Psi + alpha)) * an->c6i2) -
1176 v2 * 4. * an->si * (1. - 3. * eta) *
1177 ( - sin(2. * Psi - 3. * alpha) * an->s4i2
1178 - sin(2. * Psi + 3. * alpha) * an->c4i2 ) );
1179
1180 h3P2->data[2 * j] = amp32 * ( v * dm / 3. *
1181 ( 27. * cos(3. * Psi - 2. * alpha) * an->si*an->s4i2 +
1182 27. * cos(3. * Psi + 2. * alpha) * an->si*an->c4i2 +
1183 cos( Psi + 2. * alpha) * an->c3i2 * (5.*an->si2-3.*an->si*an->ci2-3.*an->ci*an->si2) /2. +
1184 cos( Psi - 2. * alpha) * an->s3i2 * (5.*an->ci2+3.*an->ci*an->ci2-3.*an->si*an->si2) /2.) +
1185 v2*(1./3.-eta) *
1186 ( - 8.*an->c4i2*(3.*an->ci-2.)*cos(2.*(Psi+alpha)) +
1187 8.*an->s4i2*(3.*an->ci+2.)*cos(2.*(Psi-alpha)) ) );
1188
1189 h3M2->data[2 * j] = amp32 * ( v * dm / 3. *
1190 ( 27. * cos(3. * Psi - 2. * alpha) * an->si*an->s4i2 +
1191 27. * cos(3. * Psi + 2. * alpha) * an->si*an->c4i2 +
1192 cos( Psi + 2. * alpha) * an->c3i2 * (5.*an->si2-3.*an->si*an->ci2-3.*an->ci*an->si2) /2. +
1193 cos( Psi - 2. * alpha) * an->s3i2 * (5.*an->ci2+3.*an->ci*an->ci2-3.*an->si*an->si2) /2.) -
1194 v2*(1./3.-eta) *
1195 ( 8.*an->c4i2*(3.*an->ci-2.)*cos(2.*(Psi+alpha)) -
1196 8.*an->s4i2*(3.*an->ci+2.)*cos(2.*(Psi-alpha)) ) );
1197
1198 h3P2->data[2 * j + 1 ] = amp32 * ( v * dm / 3. *
1199 ( 27. * sin(3. * Psi - 2. * alpha) * an->si*an->s4i2 -
1200 27. * cos(3. * Psi + 2. * alpha) * an->si*an->c4i2 -
1201 sin( Psi + 2. * alpha) * an->c3i2 * (5.*an->si2-3.*an->si*an->ci2-3.*an->ci*an->si2) /2. +
1202 sin( Psi - 2. * alpha) * an->s3i2 * (5.*an->ci2+3.*an->ci*an->ci2-3.*an->si*an->si2) /2.) +
1203 v2*(1./3.-eta) *
1204 ( 8.*an->c4i2*(3.*an->ci-2.)*sin(2.*(Psi+alpha)) +
1205 8.*an->s4i2*(3.*an->ci+2.)*sin(2.*(Psi-alpha)) ) );
1206
1207 h3M2->data[2 * j + 1 ] = amp32 * ( -v * dm / 3. *
1208 ( 27. * sin(3. * Psi - 2. * alpha) * an->si*an->s4i2 -
1209 27. * cos(3. * Psi + 2. * alpha) * an->si*an->c4i2 -
1210 sin( Psi + 2. * alpha) * an->c3i2 * (5.*an->si2-3.*an->si*an->ci2-3.*an->ci*an->si2) /2.+
1211 sin( Psi - 2. * alpha) * an->s3i2 * (5.*an->ci2+3.*an->ci*an->ci2-3.*an->si*an->si2) /2.)+
1212 v2*(1./3.-eta) *
1213 ( 8.*an->c4i2*(3.*an->ci-2.)*sin(2.*(Psi+alpha)) +
1214 8.*an->s4i2*(3.*an->ci+2.)*sin(2.*(Psi-alpha)) ) );
1215
1216 h3P1->data[2 * j] = amp31 * ( v * dm / 6. *
1217 ( -135. * cos(3.*Psi - alpha) * an->s2i*an->s2i2 +
1218 135. * cos(3.*Psi + alpha) * an->s2i*an->c2i2 +
1219 cos(Psi+alpha) * an->c2i2*(15.*an->cdi-20.*an->ci+13.)/2.-
1220 cos(Psi-alpha) * an->s2i2*(15.*an->cdi+20.*an->ci+13.)/2. )
1221 -v2 * (1./3.-eta)*
1222 ( 20.*an->c3i2*cos(2.*Psi+alpha)*(3.*(an->si2*an->ci+an->ci2*an->si)-5.*an->si2) +
1223 20.*an->s3i2*cos(2.*Psi-alpha)*(3.*(an->ci2*an->ci-an->si2*an->si)+5.*an->ci2) ) );
1224
1225 h3M1->data[2 * j] = amp31 * (-v * dm / 6. *
1226 ( -135. * cos(3.*Psi - alpha) * an->s2i*an->s2i2 +
1227 135. * cos(3.*Psi + alpha) * an->s2i*an->c2i2 +
1228 cos(Psi+alpha) * an->c2i2*(15.*an->cdi-20.*an->ci+13.)/2.-
1229 cos(Psi-alpha) * an->s2i2*(15.*an->cdi+20.*an->ci+13.)/2. )
1230 -v2 * (1./3.-eta)*
1231 ( 20.*an->c3i2*cos(2.*Psi+alpha)*(3.*(an->si2*an->ci+an->ci2*an->si)-5.*an->si2) +
1232 20.*an->s3i2*cos(2.*Psi-alpha)*(3.*(an->ci2*an->ci-an->si2*an->si)+5.*an->ci2) ) );
1233
1234 h3P1->data[2 * j + 1] = amp31 * ( v * dm / 6. *
1235 ( -135. * sin(3.*Psi - alpha) * an->s2i*an->s2i2 -
1236 135.* sin(3.*Psi + alpha) * an->s2i*an->c2i2 -
1237 sin(Psi+alpha) * an->c2i2*(15.*an->cdi-20.*an->ci+13.)/2.-
1238 sin(Psi-alpha) * an->s2i2*(15.*an->cdi+20.*an->ci+13.)/2. )
1239 +v2 * (1./3.-eta)*
1240 ( 20.*an->c3i2*sin(2.*Psi+alpha)*(3.*(an->si2*an->ci+an->ci2*an->si)-5.*an->si2)
1241 -20.*an->s3i2*sin(2.*Psi-alpha)*(3.*(an->ci2*an->ci-an->si2*an->si)+5.*an->ci2) ) );
1242
1243 h3M1->data[2 * j + 1] = amp31 * ( v * dm / 6. *
1244 ( -135. * sin(3.*Psi - alpha) *an->s2i*an->s2i2 -
1245 135. * sin(3.*Psi + alpha) *an->s2i*an->c2i2 -
1246 sin(Psi+alpha) * an->c2i2*(15.*an->cdi-20.*an->ci+13.)/2.-
1247 sin(Psi-alpha) * an->s2i2*(15.*an->cdi+20.*an->ci+13.)/2. )
1248 -v2 * (1./3.-eta)*
1249 ( 20.*an->c3i2*sin(2.*Psi+alpha)*(3.*(an->si2*an->ci+an->ci2*an->si)-5.*an->si2)
1250 -20.*an->s3i2*sin(2.*Psi-alpha)*(3.*(an->ci2*an->ci-an->si2*an->si)+5.*an->ci2) ) );
1251
1252 h30->data[2 * j] = 0.;
1253
1254 h30->data[2 * j + 1] = amp30 * ( v * dm *
1255 ( cos(Psi) * an->si*(cos(2.*Psi)*(45.*an->s2i)-(25.*an->cdi-21.) ) )
1256 +v2 * (1.-3.*eta) *
1257 (80. * an->s2i*an->c2i*sin(2.*Psi) ) );
1258
1259 return 0;
1260
1261}
1262
1263static int XLALSpinInspiralFillH4Modes(
1264 REAL8Vector* h4P4,
1265 REAL8Vector* h4M4,
1266 REAL8Vector* h4P3,
1267 REAL8Vector* h4M3,
1268 REAL8Vector* h4P2,
1269 REAL8Vector* h4M2,
1270 REAL8Vector* h4P1,
1271 REAL8Vector* h4M1,
1272 REAL8Vector* h40,
1273 INT4 j,
1274 REAL8 amp44,
1275 REAL8 v,
1276 REAL8 eta,
1277 REAL8 dm,
1278 REAL8 Psi,
1279 REAL8 alpha,
1280 LALSpinInspiralAngle *an
1281 ){
1282
1283 UNUSED(v);
1284 UNUSED(dm);
1285
1286 REAL8 amp43 = - amp44 * sqrt(2.);
1287 REAL8 amp42 = amp44 * sqrt(7.)/2.;
1288 REAL8 amp41 = amp44 * sqrt(3.5)/4.;
1289 REAL8 amp40 = amp44 * sqrt(17.5)/16.;
1290
1291 h4P4->data[2 * j] = amp44 * (1. - 3.*eta) *
1292 ( 4.* an->s8i2 * cos(4.*(Psi-alpha)) + cos(2.*Psi-4.*alpha) *an->s6i2*an->c2i2
1293 + an->s2i2*an->c6i2* cos(2.*Psi+4.*alpha) + 4.*an->c8i2* cos(4.*(Psi+alpha)) );
1294
1295 h4M4->data[2 * j] = amp44 * (1. - 3.*eta) *
1296 ( 4.* an->s8i2 * cos(4.*(Psi-alpha)) + cos(2.*Psi-4.*alpha) *an->s6i2*an->c2i2
1297 + an->s2i2*an->c6i2* cos(2.*Psi+4.*alpha) + 4.*an->c8i2* cos(4.*(Psi+alpha)) );
1298
1299 h4P4->data[2 * j + 1] = amp44 * (1. - 3.*eta) *
1300 ( 4.* an->s8i2 * sin(4.*(Psi-alpha)) + sin(2.*Psi-4.*alpha) *an->s6i2*an->c2i2
1301 - an->s2i2*an->c6i2* sin(2.*Psi+4.*alpha) - 4.*an->c8i2* sin(4.*(Psi+alpha)) );
1302
1303 h4M4->data[2 * j + 1] = - amp44 * (1. - 3.*eta) *
1304 ( 4.* an->s8i2 * sin(4.*(Psi-alpha)) + sin(2*Psi-4.*alpha) *an->s6i2*an->c2i2
1305 - an->s2i2*an->c6i2* sin(2.*Psi+4.*alpha) - 4.*an->c8i2* sin(4.*(Psi+alpha)) );
1306
1307 h4P3->data[2 * j] = amp43 * (1. - 3.*eta) * an->si *
1308 ( 4.*an->s6i2* cos(4.*Psi-3.*alpha) - 4.*an->c6i2* cos(4.*Psi+3.*alpha) -
1309 an->s4i2*(an->ci+0.5)/2. * cos(2.*Psi-3.*alpha) - an->c4i2*(an->ci-0.5) * cos(2.*Psi+3.*alpha) );
1310
1311 h4M3->data[2 * j] = - amp43 * (1. - 3.*eta) * an->si *
1312 ( 4.*an->s6i2* cos(4.*Psi-3.*alpha) - 4.*an->c6i2* cos(4.*Psi+3.*alpha) -
1313 an->s4i2*(an->ci+0.5)/2. * cos(2.*Psi-3.*alpha) - an->c4i2*(an->ci-0.5) * cos(2.*Psi+3.*alpha) );
1314
1315 h4P3->data[2 * j + 1] = amp43 * (1. - 3.*eta) * an->si *
1316 ( 4.*an->s6i2* sin(4.*Psi-3.*alpha) + 4.*an->c6i2* sin(4.*Psi+3.*alpha) -
1317 an->s4i2*(an->ci+0.5)/2. * sin(2.*Psi-3.*alpha) + an->c4i2*(an->ci-0.5) * sin(2.*Psi+3.*alpha) );
1318
1319 h4M3->data[2 * j + 1] = amp43 * (1. - 3.*eta) * an->si *
1320 ( 4.*an->s6i2* sin(4.*Psi-3.*alpha) + 4.*an->c6i2* sin(4.*Psi+3.*alpha) -
1321 an->s4i2*(an->ci+0.5)/2. * sin(2.*Psi-3.*alpha) + an->c4i2*(an->ci-0.5) * sin(2.*Psi+3.*alpha) );
1322
1323 h4P2->data[2 * j] = amp42 * (1. - 3.*eta) *
1324 ( 16.*an->s6i2*an->c2i2 * cos(4.*Psi-2.*alpha) + 16.*an->c6i2*an->s2i2 * cos(4.*Psi+2.*alpha)
1325 - an->c4i2 * cos(2.*(Psi+alpha))*(an->cdi-2.*an->ci+9./7.)/2. - an->s4i2 * cos(2.*(Psi-alpha))*(an->cdi+2.*an->ci+9./7.)/2. );
1326
1327 h4M2->data[2 * j] = amp42 * (1. - 3.*eta) *
1328 ( 16.*an->s6i2*an->c2i2 * cos(4.*Psi-2.*alpha) + 16.*an->c6i2*an->s2i2 * cos(4.*Psi+2.*alpha)
1329 - an->c4i2 * cos(2.*(Psi+alpha))*(an->cdi-2.*an->ci+9./7.)/2. - an->s4i2 * cos(2.*(Psi-alpha))*(an->cdi+2.*an->ci+9./7.)/2. );
1330
1331 h4P2->data[2 * j + 1] = amp42 * (1. - 3.*eta) *
1332 ( 16.*an->s6i2*an->c2i2 * sin(4.*Psi-2.*alpha) - 16.*an->c6i2*an->s2i2 * sin(4.*Psi+2.*alpha)
1333 + an->c4i2 * sin(2.*(Psi+alpha))*(an->cdi-2.*an->ci+9./7.)/2. - an->s4i2 * sin(2.*(Psi-alpha))*(an->cdi+2.*an->ci+9./7.)/2. );
1334
1335 h4M2->data[2 * j + 1] = -amp42 * (1. - 3.*eta) *
1336 ( 16.*an->s6i2*an->c2i2 * sin(4.*Psi-2.*alpha) - 16.*an->c6i2*an->s2i2 * sin(4.*Psi+2.*alpha)
1337 + an->c4i2 * sin(2.*(Psi+alpha))*(an->cdi-2.*an->ci+9./7.)/2. - an->s4i2 * sin(2.*(Psi-alpha))*(an->cdi+2.*an->ci+9./7.)/2. );
1338
1339 h4P1->data[2 * j] = amp41 * (1. - 3.*eta) *
1340 ( -64.*an->s5i2*an->c3i2 * cos(4.*Psi-alpha) +64.*an->s3i2*an->c5i2 * cos(4.*Psi+alpha) -
1341 an->s3i2*cos(2.*Psi-alpha)*((an->cdi*an->ci2-an->sdi*an->si2)+2.*(an->ci2*an->ci-an->si2*an->si)+19./7.*an->ci2) +
1342 an->c3i2*cos(2.*Psi+alpha)*((an->cdi*an->si2+an->sdi*an->ci2)-2.*(an->si*an->ci2+an->ci*an->si2)+19./7.*an->ci2) );
1343
1344 h4M1->data[2 * j] = -amp41 * (1. - 3.*eta) *
1345 ( -64*an->s5i2*an->c3i2 * cos(4.*Psi-alpha) +64.*an->s3i2*an->c5i2 * cos(4.*Psi+alpha) -
1346 an->s3i2*cos(2.*Psi-alpha)*((an->cdi*an->ci2-an->sdi*an->si2)+2.*(an->ci2*an->ci-an->si2*an->si)+19./7.*an->ci2) +
1347 an->c3i2*cos(2.*Psi+alpha)*((an->cdi*an->si2+an->sdi*an->ci2)-2.*(an->si*an->ci2+an->ci*an->si2)+19./7.*an->ci2) );
1348
1349 h4P1->data[2 * j + 1] = amp41 * (1. - 3.*eta) *
1350 ( -64.*an->s5i2*an->c3i2 * sin(4.*Psi-alpha) - 64.*an->s3i2*an->c5i2 * sin(4.*Psi+alpha) -
1351 an->s3i2*sin(2.*Psi-alpha)*((an->cdi*an->ci2-an->sdi*an->si2)+2.*(an->ci2*an->ci-an->si2*an->si)+19./7.*an->ci2) -
1352 an->c3i2*sin(2.*Psi+alpha)*((an->cdi*an->si2+an->sdi*an->ci2)-2.*(an->si*an->ci2+an->ci*an->si2)+19./7.*an->ci2) );
1353
1354 h4M1->data[2 * j + 1] = amp41 * (1. - 3.*eta) *
1355 ( -64.*an->s5i2*an->c3i2 * sin(4.*Psi-alpha) - 64.*an->s3i2*an->c5i2 * sin(4.*Psi+alpha) -
1356 an->s3i2*sin(2.*Psi-alpha)*((an->cdi*an->ci2-an->sdi*an->si2)+2.*(an->ci2*an->ci-an->si2*an->si)+19./7.*an->ci2) -
1357 an->c3i2*sin(2.*Psi+alpha)*((an->cdi*an->si2+an->sdi*an->ci2)-2.*(an->si*an->ci2+an->ci*an->si2)+19./7.*an->ci2) );
1358
1359 h40->data[2 * j] = amp40 * (1.-3.*eta) * an->s2i * (8.*an->s2i*cos(4.*Psi) +
1360 cos(2.*Psi)*(an->cdi+5./7.) );
1361 h40->data[2 * j +1] = 0.;
1362
1363 return 0;
1364}
1365
1366static int XLALSpinInspiralEngine(
1367 UINT4 neqs,
1368 const REAL8 yinit[],
1369 REAL8 amp22ini,
1370 LALPSpinInspiralRDparams *mparams,
1371 REAL8Vector* h2P2,
1372 REAL8Vector* h2M2,
1373 REAL8Vector* h2P1,
1374 REAL8Vector* h2M1,
1375 REAL8Vector* h20,
1376 REAL8Vector* h3P3,
1377 REAL8Vector* h3M3,
1378 REAL8Vector* h3P2,
1379 REAL8Vector* h3M2,
1380 REAL8Vector* h3P1,
1381 REAL8Vector* h3M1,
1382 REAL8Vector* h30,
1383 REAL8Vector* h4P4,
1384 REAL8Vector* h4M4,
1385 REAL8Vector* h4P3,
1386 REAL8Vector* h4M3,
1387 REAL8Vector* h4P2,
1388 REAL8Vector* h4M2,
1389 REAL8Vector* h4P1,
1390 REAL8Vector* h4M1,
1391 REAL8Vector* h40,
1392 REAL8Vector* freq,
1393 REAL8Vector* phase,
1394 LALPSpinInspiralPhenPars *phenPars
1395 )
1396{
1397 INT4 intreturn;
1398 UINT4 count=0;
1399 UINT4 write=0;
1400 UINT4 modcount=0;
1401 UINT4 j;
1402
1403 REAL8 dt,tm,timewrite;
1404 REAL8 Mass,dm;
1405 REAL8 v,v2;
1406 REAL8 Phi=0.;
1407 REAL8 omega;
1408
1409 REAL8 LNhx,LNhy,LNhz;
1410 REAL8 S1x,S1y,S1z;
1411 REAL8 S2x,S2y,S2z;
1412 REAL8 LNhxy;
1413 REAL8 energy = 0.;
1414 REAL8 energywrite = 0.;
1415
1416 REAL8 LNhS1=0.;
1417 REAL8 LNhS2=0.;
1418 REAL8 S1S1=0.;
1419 REAL8 S1S2=0.;
1420 REAL8 S2S2=0.;
1421 REAL8 dLNhx,dLNhy,dLNhz;
1422 REAL8 LNhS1w = 0.;
1423 REAL8 LNhS2w = 0.;
1424 REAL8 S1S1w = 0.;
1425 REAL8 S2S2w = 0.;
1426 REAL8 S1S2w = 0.;
1427
1428 //REAL8 Phiold;
1429 REAL8 alpha,alphaold;
1430
1431 REAL8 Phiwrite = 0.;
1432 REAL8 alphawrite = 0.;
1433 REAL8 omegawrite = 0.;
1434
1435 REAL8 amp22,amp33,amp44;
1436 REAL8 unitHz;
1437 LALSpinInspiralAngle trigAngle;
1438 UINT4 subsampling=1;
1439
1440 UINT4 Npoints = 20;
1441 INT4 errcode;
1442
1443 rk4In in4; // used to setup the Runge-Kutta integration
1444 rk4GSLIntegrator *integrator;
1445
1446 REAL8Vector dummy, values, dvalues, newvalues, yt, dym, dyt;
1447 // support variables
1448
1449 dummy.length = neqs * 6;
1450
1451 values.length = dvalues.length = newvalues.length = yt.length = dym.length = dyt.length = neqs;
1452
1453 if (!(dummy.data = (REAL8 *) LALMalloc(sizeof(REAL8) * neqs * 6))) {
1454 XLAL_ERROR(XLAL_ENOMEM);
1455 }
1456
1457 dt=mparams->dt;
1458 Mass=mparams->m * LAL_MTSUN_SI;
1459 dm=mparams->dm;
1460 unitHz= Mass * (REAL8) LAL_PI;
1461
1462 values.data = &dummy.data[0];
1463 dvalues.data = &dummy.data[neqs];
1464 newvalues.data = &dummy.data[2 * neqs];
1465 yt.data = &dummy.data[3 * neqs];
1466 dym.data = &dummy.data[4 * neqs];
1467 dyt.data = &dummy.data[5 * neqs];
1468
1469 REAL8 S1x0=values.data[5];
1470 REAL8 S1y0=values.data[6];
1471 REAL8 S1z0=values.data[7];
1472 REAL8 S2x0=values.data[8];
1473 REAL8 S2y0=values.data[9];
1474 REAL8 S2z0=values.data[10];
1475
1476 /* Variables initializations */
1477 for (j=0;j<neqs;j++) values.data[j]=yinit[j];
1478
1479 omega = values.data[1];
1480 //Phiold = Phi;
1481 Phi = values.data[0];
1482 v = cbrt(omega);
1483 v2 = v*v;
1484 alpha = atan2(values.data[3],values.data[2]);
1485 trigAngle.ci = LNhz = values.data[4];
1486
1487 LNhS1 = (values.data[2]*values.data[5]+values.data[3]*values.data[6]+values.data[4]*values.data[7])/mparams->m1msq;
1488 LNhS2 = (values.data[2]*values.data[8]+values.data[3]*values.data[9]+values.data[4]*values.data[10])/mparams->m2msq;
1489 S1S1 = (values.data[5]*values.data[5]+values.data[6]*values.data[6]+values.data[7]*values.data[7])/mparams->m1msq/mparams->m1msq;
1490 S2S2 = (values.data[8]*values.data[8]+values.data[9]*values.data[9]+values.data[10]*values.data[10])/mparams->m2msq/mparams->m2msq;
1491 S1S2= (values.data[5]*values.data[8]+values.data[6]*values.data[9]+values.data[7]*values.data[10])/mparams->m1msq/mparams->m2msq;;
1492
1493 while ( (OmMatch(LNhS1,LNhS1,S1S1,S1S2,S2S2) * 16. / unitHz) > (REAL4) (subsampling) / dt ) {
1494 subsampling *= 2;
1495 dt /= (REAL8) (subsampling);
1496 }
1497
1498 in4.function = LALSpinInspiralDerivatives;
1499 in4.y = &values;
1500 in4.dydx = &dvalues;
1501 in4.h = dt / Mass;
1502 in4.n = neqs;
1503 in4.yt = &yt;
1504 in4.dym = &dym;
1505 in4.dyt = &dyt;
1506
1507 /* Start of the calculation of the inspiral part via the fixed step integration method */
1508
1509 /* Initialize GSL integrator */
1510 if (!(integrator = XLALRungeKutta4Init(neqs, &in4))) {
1511 XLALFree(dummy.data);
1512 INT4 errNum = XLALClearErrno();
1513 if (errNum == XLAL_ENOMEM)
1514 XLAL_ERROR(XLAL_ENOMEM);
1515 else
1516 XLAL_ERROR(XLAL_EDOM);
1517 }
1518
1519 count = write= 0;
1520 tm = timewrite = 0.;
1521
1522 LALSpinInspiralDerivatives(&values, &dvalues, (void *) &mparams);
1523
1524 omega = values.data[1];
1525 //Phiold = Phi;
1526 Phi = values.data[0];
1527 v = cbrt(omega);
1528 v2 = v*v;
1529 alpha = atan2(values.data[3],values.data[2]);
1530
1531 REAL8Vector *domega = XLALCreateREAL8Vector(Npoints);
1532 REAL8Vector *diota = XLALCreateREAL8Vector(Npoints);
1533 REAL8Vector *dalpha = XLALCreateREAL8Vector(Npoints);
1534 REAL8Vector *ddomega = XLALCreateREAL8Vector(Npoints);
1535 REAL8Vector *ddiota = XLALCreateREAL8Vector(Npoints);
1536 REAL8Vector *ddalpha = XLALCreateREAL8Vector(Npoints);
1537
1538 do {
1539
1540 if (count%subsampling==0) {
1541
1542 modcount=0;
1543
1544 if (write >= mparams->length) {
1545 XLALRungeKutta4Free(integrator);
1546 XLALFree(dummy.data);
1547 XLAL_ERROR(XLAL_EDOM);
1548 }
1549
1550 amp22 = amp22ini * v2;
1551 amp33 = -amp22 / 4. * sqrt(5. / 42.);
1552 amp44 = amp22 * v2 * 2.*sqrt(5./7.)/9.;
1553
1554 Phiwrite = Phi;
1555 omegawrite = omega;
1556 alphawrite = alpha;
1557 energywrite = energy;
1558 LNhS1w = LNhS1;
1559 LNhS2w = LNhS2;
1560 S1S1w = S1S1;
1561 S1S2w = S1S2;
1562 S2S2w = S2S2;
1563
1564 trigAngle.ci = LNhz;
1565 trigAngle.cdi = 2. * trigAngle.ci * trigAngle.ci - 1.;
1566 trigAngle.c2i = trigAngle.ci * trigAngle.ci;
1567 trigAngle.s2i = 1. - trigAngle.ci * trigAngle.ci;
1568 trigAngle.si = sqrt(trigAngle.s2i);
1569 trigAngle.sdi = 2. * trigAngle.ci * trigAngle.si;
1570 trigAngle.c2i2 = (1. + trigAngle.ci) / 2.;
1571 trigAngle.s2i2 = (1. - trigAngle.ci) / 2.;
1572 trigAngle.ci2 = sqrt(trigAngle.c2i2);
1573 trigAngle.si2 = sqrt(trigAngle.s2i2);
1574 trigAngle.c3i2 = trigAngle.c2i2 * trigAngle.ci2;
1575 trigAngle.s3i2 = trigAngle.s2i2 * trigAngle.si2;
1576 trigAngle.c4i2 = trigAngle.c2i2 * trigAngle.c2i2;
1577 trigAngle.s4i2 = trigAngle.s2i2 * trigAngle.s2i2;
1578 trigAngle.c5i2 = trigAngle.c4i2 * trigAngle.ci2;
1579 trigAngle.s5i2 = trigAngle.s4i2 * trigAngle.si2;
1580 trigAngle.c6i2 = trigAngle.c4i2 * trigAngle.c2i2;
1581 trigAngle.s6i2 = trigAngle.s4i2 * trigAngle.s2i2;
1582 trigAngle.c8i2 = trigAngle.c4i2 * trigAngle.c4i2;
1583 trigAngle.s8i2 = trigAngle.s4i2 * trigAngle.s4i2;
1584
1585 XLALSpinInspiralFillH2Modes(h2P2,h2M2,h2P1,h2M1,h20,write,amp22,v,mparams->eta,dm,Phiwrite,alphawrite,&trigAngle);
1586
1587 XLALSpinInspiralFillH3Modes(h3P3,h3M3,h3P2,h3M2,h3P1,h3M1,h30,write,amp33,v,mparams->eta,dm,Phiwrite,alphawrite,&trigAngle);
1588
1589 XLALSpinInspiralFillH4Modes(h4P4,h4M4,h4P3,h4M3,h4P2,h4M2,h4P1,h4M1,h40,write,amp44,v,mparams->eta,dm,Phiwrite,alphawrite,&trigAngle);
1590
1591 freq->data[write]=omega;
1592 phase->data[write]=Phi;
1593
1594 write++;
1595
1596 timewrite+=mparams->dt;
1597
1598 }
1599
1600 in4.x = tm / mparams->m;
1601
1602 if (XLALRungeKutta4(&newvalues, integrator,(void *) mparams) == XLAL_FAILURE)
1603 XLAL_ERROR(XLAL_EFUNC);
1604 /* updating values of dynamical variables */
1605
1606 Phi = values.data[0] = newvalues.data[0];
1607 omega = values.data[1] = newvalues.data[1];
1608
1609 LNhx = values.data[2] = newvalues.data[2];
1610 LNhy = values.data[3] = newvalues.data[3];
1611 LNhz = values.data[4] = newvalues.data[4];
1612
1613 S1x = values.data[5] = newvalues.data[5];
1614 S1y = values.data[6] = newvalues.data[6];
1615 S1z = values.data[7] = newvalues.data[7];
1616
1617 S2x = values.data[8] = newvalues.data[8];
1618 S2y = values.data[9] = newvalues.data[9];
1619 S2z = values.data[10] = newvalues.data[10];
1620
1621 energy = values.data[11] = newvalues.data[11];
1622
1623 alphaold = alpha;
1624 LNhxy = sqrt(LNhx * LNhx + LNhy * LNhy);
1625 if (LNhxy>0.)
1626 alpha = atan2(LNhy, LNhx);
1627 else
1628 alpha = alphaold;
1629
1630 /*if (count>1) {
1631 if ((alpha*alphaold)<0.) {
1632 if (fabs(cos(2.*(Phi+alpha))-cos(2.*(Phiold+alphaold)))>0.2) {
1633 fprintf(stdout,"*** LALPSpinInspiralRD WARNING ***: Possible problem with coordinate singularity:\n Step %d LNhy: %12.6e LNhx: %12.6e Psi+alpha: %12.6e\n Step %d Psiold+alphaold %12.6e\n",write,LNhy,LNhx,(Phi+alpha)/LAL_PI,write-1,(Phiold+alphaold)/LAL_PI);
1634 fprintf(stdout," m: (%12.6e,%12.6e)\n", mparams->m1m*mparams->m, mparams->m2m*mparams->m);
1635 fprintf(stdout," S1: (%9.6f,%9.6f,%9.6f)\n",yinit[5]/mparams->m1msq,yinit[6]/mparams->m1msq,yinit[7]/mparams->m1msq);
1636 fprintf(stdout," S2: (%9.6f,%9.6f,%9.6f)\n",yinit[8]/mparams->m2msq,yinit[9]/mparams->m2msq,yinit[10]/mparams->m2msq);
1637 }
1638 }
1639 }*/
1640
1641 LNhS1 = (S1x * LNhx + S1y * LNhy + S1z * LNhz) / mparams->m1msq;
1642 LNhS2 = (S2x * LNhx + S2y * LNhy + S2z * LNhz) / mparams->m2msq;
1643 S1S1 = (S1x * S1x + S1y * S1y + S1z * S1z)/mparams->m1msq/mparams->m1msq;
1644 S2S2 = (S2x * S2x + S2y * S2y + S2z * S2z)/mparams->m2msq/mparams->m2msq;
1645 S1S2 = (S1x * S2x + S1y * S2y + S1z * S2z)/mparams->m1msq/mparams->m2msq;
1646
1647 LALSpinInspiralDerivatives(&values, &dvalues, (void *) mparams);
1648
1649 dLNhx = dvalues.data[2];
1650 dLNhy = dvalues.data[3];
1651 dLNhz = dvalues.data[4];
1652 if (LNhxy > 0.) {
1653 diota->data[Npoints-1] = -dLNhz / LNhxy;
1654 dalpha->data[Npoints-1] = (LNhx * dLNhy - LNhy * dLNhx) / LNhxy;
1655 } else {
1656 diota->data[Npoints-1] = 0.;
1657 dalpha->data[Npoints-1] = 0.;
1658 }
1659
1660 v = cbrt(omega);
1661 v2 = v*v;
1662
1663 for (j=0;j<Npoints-1;j++) {
1664 domega->data[j] = domega->data[j+1];
1665 diota->data[j] = diota->data[j+1];
1666 dalpha->data[j] = dalpha->data[j+1];
1667 }
1668 domega->data[Npoints-1] = dvalues.data[1];
1669
1670 tm += dt;
1671
1672 count++;
1673 modcount++;
1674
1675 intreturn=XLALSpinInspiralTest(0.,values.data,dvalues.data,mparams);
1676
1677 } while (intreturn==GSL_SUCCESS);
1678
1679 XLALRungeKutta4Free(integrator);
1680 XLALFree(dummy.data);
1681
1682 if (count<Npoints) {
1683 XLALPrintError("*** LALPSpinInspiralRD WARNING: inspiral integration vey short: %12.f sec\n",tm);
1684 XLAL_ERROR(XLAL_EFAILED);
1685 }
1686
1687 errcode = XLALGenerateWaveDerivative(ddomega,domega,dt);
1688 errcode += XLALGenerateWaveDerivative(ddalpha,dalpha,dt);
1689 errcode += XLALGenerateWaveDerivative(ddiota,diota,dt);
1690
1691 if (errcode != 0) {
1692 XLALPrintError("**** LALPSpinInspiralRD ERROR ****: error generating derivatives\n");
1693 XLALPrintError(" m: : %12.5f %12.5f\n",mparams->m1m*mparams->m,mparams->m2m*mparams->m);
1694 XLALPrintError(" S1: : %12.5f %12.5f %12.5f\n",S1x0,S1y0,S1z0);
1695 XLALPrintError(" S2: : %12.5f %12.5f %12.5f\n",S2x0,S2y0,S2z0);
1696 XLAL_ERROR(XLAL_EFUNC);
1697 }
1698
1699 phenPars->endtime = timewrite-mparams->dt;
1700 phenPars->intreturn = intreturn;
1701 phenPars->Psi = Phiwrite;
1702 phenPars->alpha = alphawrite;
1703 phenPars->energy = energywrite;
1704 phenPars->omega = omegawrite;
1705 phenPars->domega = domega->data[Npoints-1-modcount]/Mass;
1706 phenPars->ddomega = ddomega->data[Npoints-1-modcount]/Mass;
1707 phenPars->diota = diota->data[Npoints-1-modcount];
1708 phenPars->ddiota = ddiota->data[Npoints-1-modcount];
1709 phenPars->dalpha = dalpha->data[Npoints-1-modcount] ;
1710 phenPars->ddalpha = ddalpha->data[Npoints-1-modcount];
1711 phenPars->ci = trigAngle.ci;
1712 phenPars->countback = write-1;
1713 phenPars->LNhS1 = LNhS1w;
1714 phenPars->LNhS2 = LNhS2w;
1715 phenPars->S1S2 = S1S2w;
1716 phenPars->S1S1 = S1S1w;
1717 phenPars->S2S2 = S2S2w;
1718
1719 XLALDestroyREAL8Vector(domega);
1720 XLALDestroyREAL8Vector(diota);
1721 XLALDestroyREAL8Vector(dalpha);
1722 XLALDestroyREAL8Vector(ddomega);
1723 XLALDestroyREAL8Vector(ddiota);
1724 XLALDestroyREAL8Vector(ddalpha);
1725
1726 return XLAL_SUCCESS;
1727} /* End of XLALSpinInspiralEngine*/
1728
1729static int XLALSpinInspiralAdaptiveEngine(
1730 const UINT4 neqs,
1731 const REAL8 yinit[],
1732 REAL8 amp22ini,
1733 LALPSpinInspiralRDparams *mparams,
1734 REAL8Vector* h2P2,
1735 REAL8Vector* h2M2,
1736 REAL8Vector* h2P1,
1737 REAL8Vector* h2M1,
1738 REAL8Vector* h20,
1739 REAL8Vector* h3P3,
1740 REAL8Vector* h3M3,
1741 REAL8Vector* h3P2,
1742 REAL8Vector* h3M2,
1743 REAL8Vector* h3P1,
1744 REAL8Vector* h3M1,
1745 REAL8Vector* h30,
1746 REAL8Vector* h4P4,
1747 REAL8Vector* h4M4,
1748 REAL8Vector* h4P3,
1749 REAL8Vector* h4M3,
1750 REAL8Vector* h4P2,
1751 REAL8Vector* h4M2,
1752 REAL8Vector* h4P1,
1753 REAL8Vector* h4M1,
1754 REAL8Vector* h40,
1755 REAL8Vector* freq,
1756 REAL8Vector* phase,
1757 LALPSpinInspiralPhenPars *phenPars
1758 )
1759{
1760
1761 UINT4 j;
1762 UINT4 k;
1763 UINT4 kMatch=0;
1764 UINT4 jMatch=0;
1765 UINT4 Npoints=10;
1766 UINT4 intlen;
1767 UINT4 intreturn;
1768
1769 LALSpinInspiralAngle trigAngle;
1770
1771 REAL8Array *yout=NULL;
1772 //yout=malloc(sizeof(REAL8Array));
1773 LALAdaptiveRungeKuttaIntegrator *integrator=NULL;
1774 //integrator=malloc(sizeof(LALAdaptiveRungeKutta4Integrator));
1775 //memset(&integrator,0,sizeof(LALAdaptiveRungeKutta4Integrator)+1);
1776
1777 REAL8 Psi;
1778 REAL8 alpha=0.;
1779 REAL8 alphaold;
1780 REAL8 v,v2;
1781 REAL8 dt;
1782 REAL8 Mass;
1783 REAL8 amp22;
1784 REAL8 amp33;
1785 REAL8 amp44;
1786
1787 REAL8 LNhxy;
1788 REAL8 LNhS1;
1789 REAL8 LNhS2;
1790 REAL8 S1S1;
1791 REAL8 S1S2;
1792 REAL8 S2S2;
1793 REAL8 omegaMatch;
1794 REAL8 c1,c2;
1795
1796 INT4 errcode;
1797
1798 REAL8 *yin = XLALMalloc(sizeof(REAL8) * neqs);
1799 if (!yin) XLAL_ERROR(XLAL_ENOMEM);
1800
1801 /* allocate the integrator */
1802 integrator = XLALAdaptiveRungeKutta4Init(neqs,XLALSpinInspiralDerivatives,XLALSpinInspiralTest,1.0e-6,1.0e-6);
1803 if (!integrator) {
1804 XLALPrintError("**** LALPSpinInspiralRD ERROR ****: Cannot allocate adaptive integrator.\n");
1805 if (XLALClearErrno() == XLAL_ENOMEM)
1806 XLAL_ERROR( XLAL_ENOMEM );
1807 else
1808 XLAL_ERROR( XLAL_EDOM );
1809 }
1810
1811 /* stop the integration only when the test is true */
1812 integrator->stopontestonly = 1;
1813
1814 /* run the integration; note: time is measured in units of total mass */
1815
1816 Mass = mparams->m * LAL_MTSUN_SI;
1817 dt = mparams->dt;
1818
1819 for (j=0; j<neqs; j++) yin[j]=yinit[j];
1820
1821 REAL8 S1x0=yinit[5];
1822 REAL8 S1y0=yinit[6];
1823 REAL8 S1z0=yinit[7];
1824 REAL8 S2x0=yinit[8];
1825 REAL8 S2y0=yinit[9];
1826 REAL8 S2z0=yinit[10];
1827
1828 intlen = XLALAdaptiveRungeKutta4Hermite(integrator,(void *)mparams,yin,0.0,mparams->lengths/Mass,dt/Mass,&yout);
1829
1830 intreturn = integrator->returncode;
1831 XLALAdaptiveRungeKuttaFree(integrator);
1832
1833 /* End integration*/
1834
1835 /* Start of the integration checks*/
1836 if (!intlen) {
1837 if (XLALClearErrno() == XLAL_ENOMEM) {
1838 XLAL_ERROR( XLAL_ENOMEM);
1839 } else {
1840 XLALPrintError("**** LALPSpinInspiralRD ERROR ****: integration failed with errorcode %d, integration length %d\n",intreturn,intlen);
1841 XLAL_ERROR( XLAL_EFAILED);
1842 }
1843 }
1844
1845 /* if we have enough space, compute the waveform components; otherwise abort */
1846 if ( intlen >= mparams->length ) {
1847 XLALPrintError("**** LALPSpinInspiralRD ERROR ****: no space to write in waveforms: %d vs. %d\n",intlen,mparams->length);
1848 XLAL_ERROR(XLAL_ESIZE);
1849 }
1850
1851 if ( intlen < minIntLen ) {
1852 XLALPrintError("**** LALPSpinInspiralRD ERROR ****: incorrect integration with length %d\n",intlen);
1853 XLAL_ERROR(XLAL_ESIZE);
1854 }
1855 /* End of integration checks*/
1856
1857 REAL8 *Phi = &yout->data[1*intlen];
1858 REAL8 *omega = &yout->data[2*intlen];
1859 REAL8 *LNhx = &yout->data[3*intlen];
1860 REAL8 *LNhy = &yout->data[4*intlen];
1861 REAL8 *LNhz = &yout->data[5*intlen];
1862 REAL8 *S1x = &yout->data[6*intlen];
1863 REAL8 *S1y = &yout->data[7*intlen];
1864 REAL8 *S1z = &yout->data[8*intlen];
1865 REAL8 *S2x = &yout->data[9*intlen];
1866 REAL8 *S2y = &yout->data[10*intlen];
1867 REAL8 *S2z = &yout->data[11*intlen];
1868 REAL8 *energy = &yout->data[12*intlen];
1869
1870 if (mparams->inspiralOnly!=1) {
1871
1872 j=intlen;
1873
1874 do {
1875 j--;
1876 LNhS1=(LNhx[j]*S1x[j]+LNhy[j]*S1y[j]+LNhz[j]*S1z[j])/mparams->m1msq;
1877 LNhS2=(LNhx[j]*S2x[j]+LNhy[j]*S2y[j]+LNhz[j]*S2z[j])/mparams->m2msq;
1878 S1S1=(S1x[j]*S1x[j]+S1y[j]*S1y[j]+S1z[j]*S1z[j])/mparams->m1msq/mparams->m1msq;
1879 S1S2=(S1x[j]*S2x[j]+S1y[j]*S2y[j]+S1z[j]*S2z[j])/mparams->m1msq/mparams->m2msq;
1880 S2S2=(S2x[j]*S2x[j]+S2y[j]*S2y[j]+S2z[j]*S2z[j])/mparams->m2msq/mparams->m2msq;
1881 omegaMatch=OmMatch(LNhS1,LNhS2,S1S1,S1S2,S2S2);
1882 if (omegaMatch>omega[j]) {
1883 if (omega[j-1]<omega[j]) jMatch=j;
1884 // The numerical integrator sometimes stops and stores twice the last
1885 // omega value, this 'if' instruction avoids keeping two identical
1886 // values of omega at the end of the integration.
1887 }
1888 } while ((j>0)&&(jMatch==0));
1889
1890 if (omegaMatch<omega[jMatch]) {
1891 XLALPrintError("*** LALPSpinInspiralRD ERROR ***: Impossible to attach phenom. part\n");
1892 XLAL_ERROR(XLAL_EFAILED);
1893 }
1894
1895 // Data structure are copied into Npoints-long
1896 // REAL8Array for interpolation and derivative computation
1897 if (Npoints > intlen) Npoints = intlen;
1898
1899 if ( (omega[jMatch+1]>omega[jMatch]) && ((jMatch+1)<intlen) )
1900 kMatch=Npoints-2;
1901 else
1902 kMatch=Npoints-1;
1903 //We keep until the point where omega > omegaMatch for better derivative
1904 //computation, but do the matching at the last point at which
1905 // omega < omegaMatch
1906
1907 REAL8Vector *omega_s = XLALCreateREAL8Vector(Npoints);
1908 REAL8Vector *LNhx_s = XLALCreateREAL8Vector(Npoints);
1909 REAL8Vector *LNhy_s = XLALCreateREAL8Vector(Npoints);
1910 REAL8Vector *LNhz_s = XLALCreateREAL8Vector(Npoints);
1911 REAL8Vector *alpha_s = XLALCreateREAL8Vector(Npoints);
1912
1913 REAL8Vector *domega = XLALCreateREAL8Vector(Npoints);
1914 REAL8Vector *dLNhx = XLALCreateREAL8Vector(Npoints);
1915 REAL8Vector *dLNhy = XLALCreateREAL8Vector(Npoints);
1916 REAL8Vector *dLNhz = XLALCreateREAL8Vector(Npoints);
1917 REAL8Vector *diota = XLALCreateREAL8Vector(Npoints);
1918 REAL8Vector *dalpha = XLALCreateREAL8Vector(Npoints);
1919
1920 REAL8Vector *ddomega = XLALCreateREAL8Vector(Npoints);
1921 REAL8Vector *ddiota = XLALCreateREAL8Vector(Npoints);
1922 REAL8Vector *ddalpha = XLALCreateREAL8Vector(Npoints);
1923
1924 for (k=0;k<Npoints;k++) {
1925 j=k+jMatch-kMatch;
1926 omega_s->data[k] = omega[j];
1927 LNhx_s->data[k] = LNhx[j];
1928 LNhy_s->data[k] = LNhy[j];
1929 LNhz_s->data[k] = LNhz[j];
1930 }
1931
1932 errcode = XLALGenerateWaveDerivative(domega,omega_s,dt);
1933 errcode += XLALGenerateWaveDerivative(dLNhx,LNhx_s,dt);
1934 errcode += XLALGenerateWaveDerivative(dLNhy,LNhy_s,dt);
1935 errcode += XLALGenerateWaveDerivative(dLNhz,LNhz_s,dt);
1936 if (errcode != XLAL_SUCCESS) {
1937 XLALPrintError("**** LALPSpinInspiralRD ERROR ****: error generating first derivatives: #points %d\n",Npoints);
1938 XLALPrintError(" m: : %12.5f %12.5f\n",mparams->m1m*mparams->m,mparams->m2m*mparams->m);
1939 XLALPrintError(" S1: : %12.5f %12.5f %12.5f\n",S1x0,S1y0,S1z0);
1940 XLALPrintError(" S2: : %12.5f %12.5f %12.5f\n",S2x0,S2y0,S2z0);
1941 XLALPrintError(" omM %12.5f om[%d] %12.5f\n",omegaMatch,jMatch,*omega);
1942 XLAL_ERROR(XLAL_EFAILED);
1943 }
1944
1945 for (k=0;k<Npoints;k++) {
1946 LNhxy = sqrt(LNhx_s->data[k] * LNhx_s->data[k] + LNhy_s->data[k] * LNhy_s->data[k]);
1947 if (LNhxy > 0.) {
1948 diota->data[k] = -dLNhz->data[k] / LNhxy;
1949 dalpha->data[k] = (LNhx_s->data[k] * dLNhy->data[k] - LNhy_s->data[k] * dLNhx->data[k]) / LNhxy;
1950 } else {
1951 diota->data[k] = 0.;
1952 dalpha->data[k] = 0.;
1953 }
1954 }
1955
1956 errcode = XLALGenerateWaveDerivative(ddiota,diota,dt);
1957 errcode += XLALGenerateWaveDerivative(ddalpha,dalpha,dt);
1958 errcode += XLALGenerateWaveDerivative(ddomega,domega,dt);
1959 if (errcode != XLAL_SUCCESS) {
1960 XLALPrintError("**** LALPSpinInspiralRD ERROR ****: error generating second derivatives\n");
1961 XLALPrintError(" m: : %12.5f %12.5f\n",mparams->m1m*mparams->m,mparams->m2m*mparams->m);
1962 XLALPrintError(" S1: : %12.5f %12.5f %12.5f\n",S1x0,S1y0,S1z0);
1963 XLALPrintError(" S2: : %12.5f %12.5f %12.5f\n",S2x0,S2y0,S2z0);
1964 XLALPrintError(" omM %12.5f om[%d] %12.5f\n",omegaMatch,jMatch,*omega);
1965 XLAL_ERROR(XLAL_EFAILED);
1966 }
1967
1968 if (ddomega->data[kMatch]<0.) {
1969 XLALPrintWarning("*** LALPSpinInspiralRD WARNING: the attach of the phenom. phase has been shifted back: m1 %12.6f m2 %12.6f\n",mparams->m1m*mparams->m,mparams->m2m*mparams->m);
1970 XLALPrintWarning(" Integration returned %d\n 1025: Energy increases\n 1026: Omegadot -ve\n 1028: Omega NAN\n 1029: Omega > Omegamatch\n 1031: Omega -ve\n 1032: Omega > OmegaCut %12.6e\n",intreturn,mparams->OmCutoff);
1971 while ((kMatch>0)&&(ddomega->data[kMatch]<0.)) {
1972 kMatch--;
1973 jMatch--;
1974 }
1975 }
1976
1977 phenPars->intreturn = intreturn;
1978 phenPars->energy = energy[jMatch];
1979 phenPars->omega = omega_s->data[kMatch];
1980 phenPars->domega = domega->data[kMatch];
1981 phenPars->ddomega = ddomega->data[kMatch];
1982 phenPars->diota = diota->data[kMatch];
1983 phenPars->ddiota = ddiota->data[kMatch];
1984 phenPars->dalpha = dalpha->data[kMatch];
1985 phenPars->ddalpha = ddalpha->data[kMatch];
1986 phenPars->countback = jMatch;
1987 phenPars->Psi = Phi[jMatch];
1988 phenPars->endtime = ((REAL8) jMatch)*dt;
1989 phenPars->ci = LNhz[jMatch];
1990 phenPars->LNhS1 = LNhS1;
1991 phenPars->LNhS2 = LNhS2;
1992 phenPars->S1S2 = S1S2;
1993 phenPars->S1S1 = S1S1;
1994 phenPars->S2S2 = S2S2;
1995
1996 XLALDestroyREAL8Vector(omega_s);
1997 XLALDestroyREAL8Vector(LNhx_s);
1998 XLALDestroyREAL8Vector(LNhy_s);
1999 XLALDestroyREAL8Vector(LNhz_s);
2000 XLALDestroyREAL8Vector(alpha_s);
2001 XLALDestroyREAL8Vector(dLNhx);
2002 XLALDestroyREAL8Vector(dLNhy);
2003 XLALDestroyREAL8Vector(dLNhz);
2004 XLALDestroyREAL8Vector(diota);
2005 XLALDestroyREAL8Vector(dalpha);
2006 XLALDestroyREAL8Vector(domega);
2007 XLALDestroyREAL8Vector(ddomega);
2008 XLALDestroyREAL8Vector(ddiota);
2009 XLALDestroyREAL8Vector(ddalpha);
2010 }
2011 else {
2012 jMatch=intlen-1;
2013 phenPars->intreturn = intreturn;
2014 phenPars->energy = 0.;
2015 phenPars->omega = 0.;
2016 phenPars->domega = 0.;
2017 phenPars->ddomega = 0.;
2018 phenPars->diota = 0.;
2019 phenPars->ddiota = 0.;
2020 phenPars->dalpha = 0.;
2021 phenPars->ddalpha = 0.;
2022 phenPars->countback = intlen-1;
2023 phenPars->Psi = 0.;
2024 phenPars->endtime = 0.;
2025 phenPars->ci = 0.;
2026 phenPars->LNhS1 = 0.;
2027 phenPars->LNhS2 = 0.;
2028 phenPars->S1S2 = 0.;
2029 phenPars->S1S1 = 0.;
2030 phenPars->S2S2 = 0.;
2031 }
2032
2033 /* Now fill the Hlm waveform structures*/
2034
2035 //REAL8 alphaoold = 0.;
2036 alphaold=alpha;
2037 if ((LNhy[0]*LNhy[0]+LNhx[0]*LNhx[0])>0.)
2038 alpha=atan2(LNhy[0],LNhx[0]);
2039 else {
2040 if ((S1x[0]*S1x[0]+S1y[0]*S1y[0]+S2x[0]*S2x[0]+S2y[0]*S2y[0])>0.) {
2041 c1=0.75+mparams->eta/2-0.75*mparams->dm;
2042 c2=0.75+mparams->eta/2+0.75*mparams->dm;
2043 alpha=atan2(-c1*S1x[0]-c2*S2x[0],c1*S1y[0]+c2*S2y[0]);
2044 }
2045 else
2046 alpha=0.;
2047 }
2048
2049 for (j=0;j<=jMatch;j++) {
2050
2051 freq->data[j]=omega[j];
2052 v=cbrt(omega[j]);
2053 v2=v*v;
2054
2055 // amp22= -2.0 * params->mu * LAL_MRSUN_SI/(params->distance) * sqrt( 16.*LAL_PI/5.)*v2;
2056 // amp20 = amp22*sqrt(3/2)
2057 // Y22 \pm Y2-2= sqrt(5/PI) ((1+cos^2 t)/4, (cos t)/2)
2058 // Y21 \pm Y2-1= sqrt(5/PI) ((sin t)/2, (sin 2t)/4)
2059 // Y20 = sqrt(15/2 PI) (sin^2 t)/4
2060
2061 amp22 = amp22ini * v2;
2062 amp33 = -amp22 / 4. * sqrt(5./42.);
2063 amp44 = amp22 * sqrt(5./7.) * 2./9.* v2;
2064
2065 Psi = phase->data[j] = Phi[j];// - 2. * omega[j] * log(omega[j]);
2066
2067 trigAngle.ci = (LNhz[j]);
2068 trigAngle.cdi = 2. * trigAngle.ci * trigAngle.ci - 1.;
2069 trigAngle.c2i = trigAngle.ci * trigAngle.ci;
2070 trigAngle.s2i = 1. - trigAngle.ci * trigAngle.ci;
2071 trigAngle.si = sqrt(trigAngle.s2i);
2072 trigAngle.sdi = 2. * trigAngle.ci * trigAngle.si;
2073 trigAngle.c2i2 = (1. + trigAngle.ci) / 2.;
2074 trigAngle.s2i2 = (1. - trigAngle.ci) / 2.;
2075 trigAngle.ci2 = sqrt(trigAngle.c2i2);
2076 trigAngle.si2 = sqrt(trigAngle.s2i2);
2077 trigAngle.c3i2 = trigAngle.c2i2 * trigAngle.ci2;
2078 trigAngle.s3i2 = trigAngle.s2i2 * trigAngle.si2;
2079 trigAngle.c4i2 = trigAngle.c2i2 * trigAngle.c2i2;
2080 trigAngle.s4i2 = trigAngle.s2i2 * trigAngle.s2i2;
2081 trigAngle.c5i2 = trigAngle.c4i2 * trigAngle.ci2;
2082 trigAngle.s5i2 = trigAngle.s4i2 * trigAngle.si2;
2083 trigAngle.c6i2 = trigAngle.c4i2 * trigAngle.c2i2;
2084 trigAngle.s6i2 = trigAngle.s4i2 * trigAngle.s2i2;
2085 trigAngle.c8i2 = trigAngle.c4i2 * trigAngle.c4i2;
2086 trigAngle.s8i2 = trigAngle.s4i2 * trigAngle.s4i2;
2087
2088 //alphaoold = alphaold;
2089 alphaold = alpha;
2090 if ((LNhy[j]*LNhy[j]+LNhx[j]*LNhx[j])>0.) {
2091 alpha = atan2(LNhy[j], LNhx[j]);
2092 }
2093 else alpha = alphaold;
2094
2095 errcode = XLALSpinInspiralFillH2Modes(h2P2,h2M2,h2P1,h2M1,h20,j,amp22,v,mparams->eta,mparams->dm,Psi,alpha,&trigAngle);
2096
2097 /*if (j>2) {
2098 if ((alphaold*alphaoold)<0.) {
2099 if ( fabs(cos(2.*(Phi[j-1]+alphaold))-cos(2.*(Phi[j-2]+alphaoold)))>0.2) {
2100 fprintf(stdout,"*** LALPSpinInspiralRD WARNING ***: Possible problem with coordinate singularity:\n Step %d LNhy: %12.6e LNhx: %12.6e Psi+alpha: %12.6e alpha %12.6e\n Step %d LNhy: %12.6e LNhx: %12.6e Psi+alpha: %12.6e alpha %12.6e\n Step %d LNhy: %12.6e LNhx: %12.6e Psi+alpha: %12.6e alpha %12.6e\n",j,LNhy[j],LNhx[j],(Phi[j]+alpha)/LAL_PI,alpha/LAL_PI,j-1,LNhy[j-1],LNhx[j-1],(Phi[j-1]+alphaold)/LAL_PI,alphaold/LAL_PI,j-2,LNhy[j-2],LNhx[j-2],(Phi[j-2]+alphaoold)/LAL_PI,alphaoold/LAL_PI);
2101 fprintf(stdout," m: (%12.6e,%12.6e)\n", mparams->m1m*mparams->m, mparams->m2m*mparams->m);
2102 fprintf(stdout," S1: (%9.6f,%9.6f,%9.6f)\n",yinit[5]/mparams->m1msq,yinit[6]/mparams->m1msq,yinit[7]/mparams->m1msq);
2103 fprintf(stdout," S2: (%9.6f,%9.6f,%9.6f)\n",yinit[8]/mparams->m2msq,yinit[9]/mparams->m2msq,yinit[10]/mparams->m2msq);
2104 }
2105 }
2106 }*/
2107
2108 errcode += XLALSpinInspiralFillH3Modes(h3P3,h3M3,h3P2,h3M2,h3P1,h3M1,h30,j,amp33,v,mparams->eta,mparams->dm,Psi,alpha,&trigAngle);
2109
2110 errcode += XLALSpinInspiralFillH4Modes(h4P4,h4M4,h4P3,h4M3,h4P2,h4M2,h4P1,h4M1,h40,j,amp44,v,mparams->eta,mparams->dm,Psi,alpha,&trigAngle);
2111
2112 if (errcode != XLAL_SUCCESS)
2113 XLAL_ERROR(XLAL_EFUNC);
2114 }
2115
2116 phenPars->alpha=alpha;
2117
2118 if (yin) XLALFree(yin);
2119 if (yout) XLALDestroyREAL8Array(yout);
2120
2121 return XLAL_SUCCESS;
2122
2123} /* End of the inspiral part created via the adaptive integration method */
2124
2125
2126static int XLALPSpinInspiralRDEngine(
2127 REAL8Vector * signalvec1,
2128 REAL8Vector * signalvec2,
2129 REAL8Vector * hh,
2130 REAL8Vector * ff,
2131 REAL8Vector * phi,
2132 InspiralTemplate *params,
2133 InspiralInit *paramsInit)
2134{
2135
2136 /* declare code parameters and variables */
2137 const INT4 neqs = 11+1; // number of dynamical variables plus the energy function
2138 UINT4 origcount,apcount; // integration steps performed
2139 UINT4 count = 0; // integration steps performed
2140 UINT4 length; // signal vector length
2141 UINT4 i, j, k, l; // counters
2142
2143 REAL8 v = 0.;
2144 REAL8 v2 = 0.;
2145 REAL8 v2old;
2146 REAL8 mass; // Total mass in SI units
2147 REAL8 tim; // time (units of total mass)
2148 REAL8 unitHz;
2149 REAL8 initomega,initphi;
2150 REAL8 inc;
2151 REAL8 LNhmag,initJmag;
2152 REAL8 initS1[3],initS2[3],initLNh[3],initJ[3];
2153 REAL8 iS1[3],iS2[3];
2154 REAL8 phiJ,thetaJ;
2155 REAL8 ry[3][3],rz[3][3];
2156 REAL8 dt;
2157
2158 INT4 intreturn;
2159 REAL8 yinit[neqs];
2160
2161 REAL8Vector* h2P2;
2162 REAL8Vector* h2M2;
2163 REAL8Vector* h2P1;
2164 REAL8Vector* h2M1;
2165 REAL8Vector* h20;
2166 REAL8Vector* h3P3;
2167 REAL8Vector* h3M3;
2168 REAL8Vector* h3P2;
2169 REAL8Vector* h3M2;
2170 REAL8Vector* h3P1;
2171 REAL8Vector* h3M1;
2172 REAL8Vector* h30;
2173 REAL8Vector* h4P4;
2174 REAL8Vector* h4M4;
2175 REAL8Vector* h4P3;
2176 REAL8Vector* h4M3;
2177 REAL8Vector* h4P2;
2178 REAL8Vector* h4M2;
2179 REAL8Vector* h4P1;
2180 REAL8Vector* h4M1;
2181 REAL8Vector* h40;
2182
2183 REAL8Vector* sigp;
2184 REAL8Vector* sigc;
2185 REAL8Vector* fap;
2186 REAL8Vector* hap;
2187 REAL8Vector* phap;
2188
2189 LALPSpinInspiralRDparams mparams;
2190 LALPSpinInspiralPhenPars phenPars;
2191 LALSpinInspiralAngle trigAngle;
2192
2193 REAL8 Psi=0.;
2194 REAL8 amp22ini,amp22,amp33,amp44;
2195
2196 REAL8 alpha=0.;
2197
2198 REAL8 t0,tAs;
2199 REAL8 om0,om1,om;
2200 REAL8 Psi0,alpha0;
2201 REAL8 dalpha0,dalpha1;
2202 //REAL8 omold,iota0,diota0,diota1;
2203 REAL8 LNhS1,LNhS2;
2204 REAL8 S1S1,S1S2,S2S2;
2205
2206 COMPLEX8Vector *modefreqs;
2207 COMPLEX16 MultSphHarmP; // Generic spin-weighted spherical harmonics
2208 COMPLEX16 MultSphHarmM; // Generic spin-weighted spherical harmonics
2209 REAL8 x0, x1, x2, x3;
2210
2211 /* The number of Ring Down modes is hard-coded here */
2212 const UINT4 nmodes=2;
2213 /* Nmodes should be restricted to either 1 or 2*/
2214
2215 UINT4 errcode;
2216
2217 REAL8 finalMass,finalSpin;
2218 REAL8 energy=0.;
2219 REAL8 omegaMatch;
2220 REAL8 frOmRD,omegaRD;
2221
2222
2223 if(!params)
2224 XLAL_ERROR(XLAL_EFAULT);
2225
2226 if ((params->fCutoff<=0.)&&(params->inspiralOnly==1)) {
2227 XLALPrintError("*** LALPSIRD ERROR ***: fCutoff %12.6e, with inspiral flag on it is mandatory to specify a positive cutoff frequency\n",params->fCutoff);
2228 XLAL_ERROR(XLAL_EDOM);
2229 }
2230
2231 mass = params->totalMass * LAL_MTSUN_SI;
2232 unitHz = params->totalMass * LAL_MTSUN_SI * (REAL8) LAL_PI;
2233
2234 if ((signalvec2)||(hh))
2235 params->nStartPad = 0; /* must be zero for templates and injection */
2236 /* -- length in seconds from Newtonian formula; */
2237
2238 dt = 1. / params->tSampling;
2239
2240 /* setup coefficients for PN equations */
2241 XLALPSpinInspiralRDSetParams(&mparams,params,paramsInit);
2242
2243 /* Check that initial frequency is smaller than omegamatch ~ xxyy for m=100 Msun */
2244 initphi = params->startPhase/2.;
2245 initomega = params->fLower*unitHz;
2246
2247 /* Check that initial frequency is smaller than omegamatch ~ xxyy for m=100 Msun */
2248
2249 LNhS1=params->spin1[2];
2250 LNhS2=params->spin2[2];
2251 S1S1=params->spin1[0]*params->spin1[0]+params->spin1[1]*params->spin1[1]+params->spin1[2]*params->spin1[2];
2252 S1S2=params->spin1[0]*params->spin2[0]+params->spin1[1]*params->spin2[1]+params->spin1[2]*params->spin2[2];
2253 S2S2=params->spin2[0]*params->spin2[0]+params->spin2[1]*params->spin2[1]+params->spin2[2]*params->spin2[2];
2254
2255 omegaMatch = OmMatch(LNhS1,LNhS2,S1S1,S1S2,S2S2);
2256
2257 if ( initomega > omegaMatch ) {
2258 /*if ((params->spin1[0]==params->spin1[1])&&(params->spin1[1]==params->spin2[0])&&(params->spin2[0]==params->spin2[1])&&(params->spin2[1]==0.)) {
2259 //Beware, this correspond to a shift of the initial phase!
2260 initomega = 0.95*omegaMatch;
2261 fprintf(stdout,"*** LALPSpinInspiralRD WARNING ***: Initial frequency reset from %12.6e to %12.6e Hz, m:(%12.4e,%12.4e)\n",params->fLower,initomega/unitHz,params->mass1,params->mass2);
2262 }*/
2263 /*else {*/
2264 XLALPrintError("**** LALPSpinInspiralRD ERROR ****: Initial frequency too high: %11.5e for omM ~ %11.5e and m:(%8.3f, %8.3f)\n",params->fLower,omegaMatch/unitHz,params->mass1,params->mass2);
2265 XLAL_ERROR(XLAL_EFAILED);
2266 /*}*/
2267 }
2268
2269 /* Here we use the following convention:
2270 the coordinates of the spin vectors params->spin1,2 and the params->inclination
2271 variable refers to different physical parameters according to the value of
2272 params->axisChoice:
2273
2274 * OrbitalL: params->inclination denotes the angle between the view direction
2275 N and the initial L (initial L//z, N in the x-z plane) and the spin
2276 coordinates are given with respect to initial L.
2277 * TotalJ: params->inclination denotes the angle between the view directoin
2278 and J (J is constant during the evolution, J//z, both N and initial
2279 L are in the x-z plane) and the spin coordinates are given wrt
2280 initial ** L **.
2281
2282 * View: params->inclination denotes the angle between the initial L and N
2283 (N//z, initial L in the x-z plane) and the spin coordinates
2284 are given with respect to N.
2285
2286 In order to reproduce the results of the SpinTaylor code View must be chosen.
2287 The spin magnitude are normalized to the individual mass^2, i.e.
2288 they are dimension-less.
2289 The modulus of the initial angular momentum is fixed by m1,m2 and
2290 initial frequency.
2291 The polarization angle is not used here, it enters the pattern
2292 functions along with the angles marking the sky position of the
2293 source. */
2294
2295 // Physical magnitude of the orbital angular momentum
2296 LNhmag = params->eta * params->totalMass * params->totalMass / cbrt(initomega);
2297
2298 // Physical values of the spins
2299 for (i = 0; i < 3; i++) {
2300 initS1[i] = params->spin1[i] * params->mass1 * params->mass1;
2301 initS2[i] = params->spin2[i] * params->mass2 * params->mass2;
2302 }
2303
2304 switch (params->axisChoice) {
2305
2306 case LAL_SIM_INSPIRAL_FRAME_AXIS_ORBITAL_L:
2307 //printf("*** OrbitalL ***\n");
2308 initLNh[0] = 0.;
2309 initLNh[1] = 0.;
2310 initLNh[2] = 1.;
2311 inc = params->inclination;
2312 break;
2313
2314 case LAL_SIM_INSPIRAL_FRAME_AXIS_TOTAL_J:
2315 //printf("*** TotalJ ***\n");
2316 for (j=0;j<3;j++) {
2317 iS1[j] = initS1[j];
2318 iS2[j] = initS2[j];
2319 initJ[j] = iS1[j] + iS2[j];
2320 initS1[j] = initS2[j]=0.;
2321 initLNh[j] = 0.;
2322 }
2323 initJ[2] += LNhmag;
2324 initJmag = sqrt(initJ[0] * initJ[0] + initJ[1] * initJ[1] + initJ[2] * initJ[2]);
2325 if (initJ[1])
2326 phiJ = atan2(initJ[1], initJ[0]);
2327 else
2328 phiJ = 0.;
2329 thetaJ = acos(initJ[2]/initJmag);
2330 rz[0][0] = -cos(phiJ);
2331 rz[0][1] = -sin(phiJ);
2332 rz[0][2] = 0.;
2333 rz[1][0] = sin(phiJ);
2334 rz[1][1] = -cos(phiJ);
2335 rz[1][2] = 0.;
2336 rz[2][0] = 0.;
2337 rz[2][1] = 0.;
2338 rz[2][2] = 1.;
2339 ry[0][0] = cos(thetaJ);
2340 ry[0][1] = 0;
2341 ry[0][2] = sin(thetaJ);
2342 ry[1][0] = 0.;
2343 ry[1][1] = 1.;
2344 ry[1][2] = 0.;
2345 ry[2][0] = -sin(thetaJ);
2346 ry[2][1] = 0.;
2347 ry[2][2] = cos(thetaJ);
2348 for (j = 0; j < 3; j++) {
2349 for (k = 0; k < 3; k++) {
2350 initLNh[j] += ry[j][k] * rz[k][2];
2351 for (l = 0; l < 3; l++) {
2352 initS1[j] += ry[j][k] * rz[k][l] * iS1[l];
2353 initS2[j] += ry[j][k] * rz[k][l] * iS2[l];
2354 }
2355 }
2356 }
2357 inc = params->inclination;
2358 break;
2359
2360 case LAL_SIM_INSPIRAL_FRAME_AXIS_VIEW:
2361 default:
2362 //printf("*** View ***\n");
2363 initLNh[0] = sin(params->inclination);
2364 initLNh[1] = 0.;
2365 initLNh[2] = cos(params->inclination);
2366 inc = 0.;
2367 break;
2368 }
2369
2370 /*All the PN formulas used in the differential equation integration
2371 assume that the spin variables are the physical ones divided by
2372 totalmasss^2, here we introduce the correct normalization, changing the
2373 input one, where spin components were normalized on individual mass. */
2374 for (j = 0; j < 3; j++) {
2375 initS1[j] /= params->totalMass * params->totalMass;
2376 initS2[j] /= params->totalMass * params->totalMass;
2377 }
2378
2379 if (signalvec1) {
2380 length = signalvec1->length;
2381 } else {
2382 if (ff)
2383 length = ff->length;
2384 else
2385 length = 0;
2386 }
2387 mparams.length = length;
2388
2389 /* Allocate memory for temporary arrays */
2390
2391 h2P2 = XLALCreateREAL8Vector(length * 2);
2392 h2M2 = XLALCreateREAL8Vector(length * 2);
2393 h2P1 = XLALCreateREAL8Vector(length * 2);
2394 h2M1 = XLALCreateREAL8Vector(length * 2);
2395 h20 = XLALCreateREAL8Vector(length * 2);
2396 h3P3 = XLALCreateREAL8Vector(length * 2);
2397 h3M3 = XLALCreateREAL8Vector(length * 2);
2398 h3P2 = XLALCreateREAL8Vector(length * 2);
2399 h3M2 = XLALCreateREAL8Vector(length * 2);
2400 h3P1 = XLALCreateREAL8Vector(length * 2);
2401 h3M1 = XLALCreateREAL8Vector(length * 2);
2402 h30 = XLALCreateREAL8Vector(length * 2);
2403 h4P4 = XLALCreateREAL8Vector(length * 2);
2404 h4M4 = XLALCreateREAL8Vector(length * 2);
2405 h4P3 = XLALCreateREAL8Vector(length * 2);
2406 h4M3 = XLALCreateREAL8Vector(length * 2);
2407 h4P2 = XLALCreateREAL8Vector(length * 2);
2408 h4M2 = XLALCreateREAL8Vector(length * 2);
2409 h4P1 = XLALCreateREAL8Vector(length * 2);
2410 h4M1 = XLALCreateREAL8Vector(length * 2);
2411 h40 = XLALCreateREAL8Vector(length * 2);
2412 sigp = XLALCreateREAL8Vector(length);
2413 sigc = XLALCreateREAL8Vector(length);
2414 hap = XLALCreateREAL8Vector(length * 2);
2415 fap = XLALCreateREAL8Vector(length);
2416 phap = XLALCreateREAL8Vector(length);
2417
2418 if (!h2P2 || !h2M2 || !h2P1 || !h2M1 || !h20 || !sigp || !sigc || !fap || !phap || !hap || !h3P3 || !h3M3 || !h3P2 || !h3M2 || !h3P1 || !h3M1 || !h30 || !h4P4 || !h4M4 || !h4P3 || !h4M3 || !h4P2 || !h4M2 || !h4P1 || !h4M1 || !h40 ) {
2419 if (h2P2)
2420 XLALDestroyREAL8Vector(h2P2);
2421 if (h2M2)
2422 XLALDestroyREAL8Vector(h2M2);
2423 if (h2P1)
2424 XLALDestroyREAL8Vector(h2P1);
2425 if (h2M2)
2426 XLALDestroyREAL8Vector(h2M1);
2427 if (h20)
2428 XLALDestroyREAL8Vector(h20);
2429 if (h3P3)
2430 XLALDestroyREAL8Vector(h3P3);
2431 if (h3M3)
2432 XLALDestroyREAL8Vector(h3M3);
2433 if (h3P2)
2434 XLALDestroyREAL8Vector(h3P2);
2435 if (h3M2)
2436 XLALDestroyREAL8Vector(h3M2);
2437 if (h3P1)
2438 XLALDestroyREAL8Vector(h3P1);
2439 if (h3M1)
2440 XLALDestroyREAL8Vector(h3M1);
2441 if (h30)
2442 XLALDestroyREAL8Vector(h30);
2443 if (h4P4)
2444 XLALDestroyREAL8Vector(h4P4);
2445 if (h4M4)
2446 XLALDestroyREAL8Vector(h4M4);
2447 if (h4P3)
2448 XLALDestroyREAL8Vector(h4P3);
2449 if (h4M3)
2450 XLALDestroyREAL8Vector(h4M3);
2451 if (h4P2)
2452 XLALDestroyREAL8Vector(h4P2);
2453 if (h4M2)
2454 XLALDestroyREAL8Vector(h4M2);
2455 if (h4P1)
2456 XLALDestroyREAL8Vector(h4P1);
2457 if (h4M1)
2458 XLALDestroyREAL8Vector(h4M1);
2459 if (h40)
2460 XLALDestroyREAL8Vector(h40);
2461 if (sigp)
2462 XLALDestroyREAL8Vector(sigp);
2463 if (sigc)
2464 XLALDestroyREAL8Vector(sigc);
2465 if (fap)
2466 XLALDestroyREAL8Vector(fap);
2467 if (hap)
2468 XLALDestroyREAL8Vector(hap);
2469 if (phap)
2470 XLALDestroyREAL8Vector(phap);
2471 XLAL_ERROR(XLAL_ENOMEM);
2472 }
2473
2474 memset(h2P2->data, 0, h2P2->length * sizeof(REAL8));
2475 memset(h2M2->data, 0, h2M2->length * sizeof(REAL8));
2476 memset(h2P1->data, 0, h2P1->length * sizeof(REAL8));
2477 memset(h2M1->data, 0, h2P1->length * sizeof(REAL8));
2478 memset(h20->data, 0, h20->length * sizeof(REAL8));
2479 memset(h3P3->data, 0, h3P3->length * sizeof(REAL8));
2480 memset(h3M3->data, 0, h3M3->length * sizeof(REAL8));
2481 memset(h3P2->data, 0, h3P2->length * sizeof(REAL8));
2482 memset(h3M2->data, 0, h3M2->length * sizeof(REAL8));
2483 memset(h3P1->data, 0, h3P1->length * sizeof(REAL8));
2484 memset(h3M1->data, 0, h3M1->length * sizeof(REAL8));
2485 memset(h30->data, 0, h30->length * sizeof(REAL8));
2486 memset(h4P4->data, 0, h3P3->length * sizeof(REAL8));
2487 memset(h4M4->data, 0, h3M3->length * sizeof(REAL8));
2488 memset(h4P3->data, 0, h3P3->length * sizeof(REAL8));
2489 memset(h4M3->data, 0, h3M3->length * sizeof(REAL8));
2490 memset(h4P2->data, 0, h3P2->length * sizeof(REAL8));
2491 memset(h4M2->data, 0, h3M2->length * sizeof(REAL8));
2492 memset(h4P1->data, 0, h3P1->length * sizeof(REAL8));
2493 memset(h4M1->data, 0, h3M1->length * sizeof(REAL8));
2494 memset(h40->data, 0, h30->length * sizeof(REAL8));
2495 memset(sigp->data, 0, sigp->length * sizeof(REAL8));
2496 memset(sigc->data, 0, sigc->length * sizeof(REAL8));
2497 memset(hap->data, 0, hap->length * sizeof(REAL8));
2498 memset(fap->data, 0, fap->length * sizeof(REAL8));
2499 memset(phap->data, 0, phap->length * sizeof(REAL8));
2500
2501 /* Here there used to be a check that OmegaRD is smaller than Nyquist, it
2502 has been taken out */
2503
2504 params->ampOrder = (LALPNOrder) 1;
2505 if (params->distance > 0.)
2506 amp22ini = -2.0 * params->mu * LAL_MRSUN_SI / params->distance * sqrt(16. * LAL_PI / 5.);
2507 else
2508 amp22ini = 2. * sqrt(LAL_PI / 5.0) * params->signalAmplitude;
2509
2510 /* initialize the coordinates */
2511 yinit[0] = initphi; /* phi */
2512 yinit[1] = initomega; /* omega (really pi M f) */
2513 yinit[2] = initLNh[0]; /* LNh(x,y,z) */
2514 yinit[3] = initLNh[1];
2515 yinit[4] = initLNh[2];
2516
2517 yinit[5] = initS1[0]; /* Spin1(x,y,z) */
2518 yinit[6] = initS1[1];
2519 yinit[7] = initS1[2];
2520
2521 yinit[8] = initS2[0]; /* Spin2(x,y,z) */
2522 yinit[9] = initS2[1];
2523 yinit[10]= initS2[2];
2524
2525 yinit[11]= 0.;
2526
2527 phenPars.intreturn = 0;
2528 phenPars.energy = 0.;
2529 phenPars.omega = 0.;
2530 phenPars.domega = 0.;
2531 phenPars.ddomega = 0.;
2532 phenPars.diota = 0.;
2533 phenPars.ddiota = 0.;
2534 phenPars.dalpha = 0.;
2535 phenPars.ddalpha = 0.;
2536 phenPars.countback = 0;
2537 phenPars.endtime = 0.;
2538 phenPars.Psi = 0.;
2539 phenPars.alpha = 0.;
2540 phenPars.ci = 0.;
2541 phenPars.LNhS1 = 0.;
2542 phenPars.LNhS2 = 0.;
2543 phenPars.S1S1 = 0.;
2544 phenPars.S1S2 = 0.;
2545 phenPars.S2S2 = 0.;
2546
2547 if (params->fixedStep == 1) {
2548 if (XLALSpinInspiralEngine(neqs,yinit,amp22ini,&mparams,h2P2,h2M2,h2P1,h2M1,h20,h3P3,h3M3,h3P2,h3M2,h3P1,h3M1,h30,h4P4,h4M4,h4P3,h4M3,h4P2,h4M2,h4P1,h4M1,h40,fap,phap,&phenPars) == XLAL_FAILURE) {
2549 XLALDestroyREAL8Vector(h2P2);
2550 XLALDestroyREAL8Vector(h2M2);
2551 XLALDestroyREAL8Vector(h2P1);
2552 XLALDestroyREAL8Vector(h2M1);
2553 XLALDestroyREAL8Vector(h20);
2554 XLALDestroyREAL8Vector(h3P3);
2555 XLALDestroyREAL8Vector(h3M3);
2556 XLALDestroyREAL8Vector(h3P2);
2557 XLALDestroyREAL8Vector(h3M2);
2558 XLALDestroyREAL8Vector(h3P1);
2559 XLALDestroyREAL8Vector(h3M1);
2560 XLALDestroyREAL8Vector(h30);
2561 XLALDestroyREAL8Vector(h4P4);
2562 XLALDestroyREAL8Vector(h4M4);
2563 XLALDestroyREAL8Vector(h4P3);
2564 XLALDestroyREAL8Vector(h4M3);
2565 XLALDestroyREAL8Vector(h4P2);
2566 XLALDestroyREAL8Vector(h4M2);
2567 XLALDestroyREAL8Vector(h4P1);
2568 XLALDestroyREAL8Vector(h4M1);
2569 XLALDestroyREAL8Vector(h40);
2570 XLALDestroyREAL8Vector(sigp);
2571 XLALDestroyREAL8Vector(sigc);
2572 XLALDestroyREAL8Vector(fap);
2573 XLALDestroyREAL8Vector(hap);
2574 XLALDestroyREAL8Vector(phap);
2575 XLAL_ERROR(XLAL_EFUNC);
2576 }
2577 } else {
2578 if (XLALSpinInspiralAdaptiveEngine(neqs,yinit,amp22ini,&mparams,h2P2,h2M2,h2P1,h2M1,h20,h3P3,h3M3,h3P2,h3M2,h3P1,h3M1,h30,h4P4,h4M4,h4P3,h4M3,h4P2,h4M2,h4P1,h4M1,h40,fap,phap,&phenPars) == XLAL_FAILURE) {
2579 XLALDestroyREAL8Vector(h2P2);
2580 XLALDestroyREAL8Vector(h2M2);
2581 XLALDestroyREAL8Vector(h2P1);
2582 XLALDestroyREAL8Vector(h2M1);
2583 XLALDestroyREAL8Vector(h20);
2584 XLALDestroyREAL8Vector(h3P3);
2585 XLALDestroyREAL8Vector(h3M3);
2586 XLALDestroyREAL8Vector(h3P2);
2587 XLALDestroyREAL8Vector(h3M2);
2588 XLALDestroyREAL8Vector(h3P1);
2589 XLALDestroyREAL8Vector(h3M1);
2590 XLALDestroyREAL8Vector(h30);
2591 XLALDestroyREAL8Vector(h4P4);
2592 XLALDestroyREAL8Vector(h4M4);
2593 XLALDestroyREAL8Vector(h4P3);
2594 XLALDestroyREAL8Vector(h4M3);
2595 XLALDestroyREAL8Vector(h4P2);
2596 XLALDestroyREAL8Vector(h4M2);
2597 XLALDestroyREAL8Vector(h4P1);
2598 XLALDestroyREAL8Vector(h4M1);
2599 XLALDestroyREAL8Vector(h40);
2600 XLALDestroyREAL8Vector(sigp);
2601 XLALDestroyREAL8Vector(sigc);
2602 XLALDestroyREAL8Vector(fap);
2603 XLALDestroyREAL8Vector(hap);
2604 XLALDestroyREAL8Vector(phap);
2605 XLAL_ERROR(XLAL_EFUNC);
2606 }
2607 }
2608 intreturn=phenPars.intreturn;
2609 /* report on abnormal termination:
2610 Termination is fine if omegamatch is passed or if energy starts
2611 increasing */
2612
2613 if ( (intreturn!=LALPSIRDPN_TEST_OMEGACUT) && (intreturn != LALPSIRDPN_TEST_OMEGAMATCH) && (intreturn != LALPSIRDPN_TEST_ENERGY) )
2614 {
2615 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: integration terminated with code %d.\n",intreturn);
2616 XLALPrintWarning(" 1025: Energy increases\n 1026: Omegadot -ve\n 1028: Omega NAN\n 1029: Omega > Omegamatch\n 1031: Omega -ve\n 1032: Omega > OmegaCut %12.6e\n",mparams.OmCutoff);
2617 XLALPrintWarning(" Waveform parameters were m1 = %14.6e, m2 = %14.6e, inc = %10.6f,\n", params->mass1, params->mass2, params->inclination);
2618 XLALPrintWarning(" S1 = (%10.6f,%10.6f,%10.6f)\n", params->spin1[0], params->spin1[1], params->spin1[2]);
2619 XLALPrintWarning(" S2 = (%10.6f,%10.6f,%10.6f)\n", params->spin2[0], params->spin2[1], params->spin2[2]);
2620 }
2621
2622 count = phenPars.countback;
2623 params->tC = ((REAL8) count) * dt;
2624
2625 if ((params->inspiralOnly != 1)&&(intreturn==LALPSIRDPN_TEST_OMEGAMATCH)) {
2626
2627 tim = t0 = phenPars.endtime;
2628 tAs = t0 + 2. * phenPars.domega / phenPars.ddomega;
2629 om1 = phenPars.domega * tAs * (1. - t0 / tAs) * (1. - t0 / tAs);
2630 om0 = phenPars.omega - om1 / (1. - t0 / tAs);
2631 om = phenPars.omega;
2632
2633 //diota1 = phenPars.ddiota * tAs * (1. - t0 / tAs) * (1. - t0 / tAs);
2634 //diota0 = phenPars.diota - diota1 / (1. - t0 / tAs);
2635
2636 dalpha1 = phenPars.ddalpha * tAs * (1. - t0 / tAs) * (1. - t0 / tAs);
2637 dalpha0 = phenPars.dalpha - dalpha1 / (1. - t0 / tAs);
2638
2639 //printf("time %12.6e count %d\n",tim,phenPars.countback);
2640
2641 if ((tAs < t0) || (om1 < 0.)) {
2642 XLALPrintError("**** LALPSpinInspiralRD ERROR ****: Could not attach phen part for:\n");
2643 XLALPrintError(" tAs %12.6e t0 %12.6e om1 %12.6e\n",tAs,t0,om1);
2644 XLALPrintError(" m1 = %14.6e, m2 = %14.6e, inc = %10.6f,\n", params->mass1, params->mass2, params->inclination);
2645 XLALPrintError(" S1 = (%10.6f,%10.6f,%10.6f)\n", params->spin1[0], params->spin1[1], params->spin1[2]);
2646 XLALPrintError(" S2 = (%10.6f,%10.6f,%10.6f)\n", params->spin2[0], params->spin2[1], params->spin2[2]);
2647 XLALDestroyREAL8Vector(h2P2);
2648 XLALDestroyREAL8Vector(h2M2);
2649 XLALDestroyREAL8Vector(h2P1);
2650 XLALDestroyREAL8Vector(h2M1);
2651 XLALDestroyREAL8Vector(h20);
2652 XLALDestroyREAL8Vector(h3P3);
2653 XLALDestroyREAL8Vector(h3M3);
2654 XLALDestroyREAL8Vector(h3P2);
2655 XLALDestroyREAL8Vector(h3M2);
2656 XLALDestroyREAL8Vector(h3P1);
2657 XLALDestroyREAL8Vector(h3M1);
2658 XLALDestroyREAL8Vector(h30);
2659 XLALDestroyREAL8Vector(h4P4);
2660 XLALDestroyREAL8Vector(h4M4);
2661 XLALDestroyREAL8Vector(h4P3);
2662 XLALDestroyREAL8Vector(h4M3);
2663 XLALDestroyREAL8Vector(h4P2);
2664 XLALDestroyREAL8Vector(h4M2);
2665 XLALDestroyREAL8Vector(h4P1);
2666 XLALDestroyREAL8Vector(h4M1);
2667 XLALDestroyREAL8Vector(h40);
2668 XLALDestroyREAL8Vector(hap);
2669 XLALDestroyREAL8Vector(fap);
2670 XLALDestroyREAL8Vector(phap);
2671 XLALDestroyREAL8Vector(sigp);
2672 XLALDestroyREAL8Vector(sigc);
2673 XLAL_ERROR(XLAL_EFAILED);
2674 }
2675 else {
2676 trigAngle.ci = phenPars.ci;
2677 trigAngle.cdi = 2. * trigAngle.ci * trigAngle.ci - 1.;
2678 trigAngle.c2i = trigAngle.ci * trigAngle.ci;
2679 trigAngle.s2i = 1. - trigAngle.ci * trigAngle.ci;
2680 trigAngle.si = sqrt(trigAngle.s2i);
2681 trigAngle.sdi = 2. * trigAngle.ci * trigAngle.si;
2682 trigAngle.c2i2 = (1. + trigAngle.ci) / 2.;
2683 trigAngle.s2i2 = (1. - trigAngle.ci) / 2.;
2684 trigAngle.ci2 = sqrt(trigAngle.c2i2);
2685 trigAngle.si2 = sqrt(trigAngle.s2i2);
2686 trigAngle.c3i2 = trigAngle.c2i2 * trigAngle.ci2;
2687 trigAngle.s3i2 = trigAngle.s2i2 * trigAngle.si2;
2688 trigAngle.c4i2 = trigAngle.c2i2 * trigAngle.c2i2;
2689 trigAngle.s4i2 = trigAngle.s2i2 * trigAngle.s2i2;
2690 trigAngle.c5i2 = trigAngle.c4i2 * trigAngle.ci2;
2691 trigAngle.s5i2 = trigAngle.s4i2 * trigAngle.si2;
2692 trigAngle.c6i2 = trigAngle.c4i2 * trigAngle.c2i2;
2693 trigAngle.s6i2 = trigAngle.s4i2 * trigAngle.s2i2;
2694 trigAngle.c8i2 = trigAngle.c4i2 * trigAngle.c4i2;
2695 trigAngle.s8i2 = trigAngle.s4i2 * trigAngle.s4i2;
2696
2697 Psi = phenPars.Psi;// - 2. * om * log(om);
2698 Psi0 = Psi + tAs * (om1/mass -dalpha1*trigAngle.ci) * log(1. - t0 / tAs);
2699 alpha0 = phenPars.alpha + tAs * dalpha1 * log(1. - t0 / tAs);
2700 //iota0 = acos(phenPars.ci) + diota1 * tAs * log(1. - t0 / tAs);
2701 energy = phenPars.energy;
2702
2703 /* Get QNM frequencies */
2704 errcode = XLALPSpinFinalMassSpin(&finalMass, &finalSpin, params, energy, initLNh);
2705 modefreqs=XLALCreateCOMPLEX8Vector(nmodes);
2706 errcode+=XLALPSpinGenerateQNMFreq(modefreqs, params, 2, 2, nmodes, finalMass, finalSpin);
2707 if (errcode != XLAL_SUCCESS) {
2708 XLALPrintError("**** LALPhenSpinInspiralRD ERROR ****: impossible to generate RingDown frequency\n");
2709 XLALPrintError( " m (%11.4e %11.4e) f0 %11.4e\n",params->mass1, params->mass2, params->fLower);
2710 XLALPrintError( " S1 (%8.4f %8.4f %8.4f)\n", initS1[0],initS1[1], initS1[2]);
2711 XLALPrintError( " S2 (%8.4f %8.4f %8.4f)\n", initS2[0],initS2[1], initS2[2]);
2712 XLALDestroyREAL8Vector(h2P2);
2713 XLALDestroyREAL8Vector(h2M2);
2714 XLALDestroyREAL8Vector(h2P1);
2715 XLALDestroyREAL8Vector(h2M1);
2716 XLALDestroyREAL8Vector(h20);
2717 XLALDestroyREAL8Vector(h3P3);
2718 XLALDestroyREAL8Vector(h3M3);
2719 XLALDestroyREAL8Vector(h3P2);
2720 XLALDestroyREAL8Vector(h3M2);
2721 XLALDestroyREAL8Vector(h3P1);
2722 XLALDestroyREAL8Vector(h3M1);
2723 XLALDestroyREAL8Vector(h30);
2724 XLALDestroyREAL8Vector(h4P4);
2725 XLALDestroyREAL8Vector(h4M4);
2726 XLALDestroyREAL8Vector(h4P3);
2727 XLALDestroyREAL8Vector(h4M3);
2728 XLALDestroyREAL8Vector(h4P2);
2729 XLALDestroyREAL8Vector(h4M2);
2730 XLALDestroyREAL8Vector(h4P1);
2731 XLALDestroyREAL8Vector(h4M1);
2732 XLALDestroyREAL8Vector(h40);
2733 XLALDestroyREAL8Vector(fap);
2734 XLALDestroyREAL8Vector(phap);
2735 XLALDestroyREAL8Vector(hap);
2736 XLALDestroyREAL8Vector(sigp);
2737 XLALDestroyREAL8Vector(sigc);
2738 XLALDestroyCOMPLEX8Vector(modefreqs);
2739 XLAL_ERROR(XLAL_EFAILED);
2740 }
2741
2742 omegaRD = crealf(modefreqs->data[0]) * unitHz / LAL_PI / 2.;
2743 frOmRD = fracRD(phenPars.LNhS1,phenPars.LNhS2,phenPars.S1S1,phenPars.S1S2,phenPars.S2S2)*omegaRD;
2744
2745 v = cbrt(om);
2746 v2 = v*v;
2747 amp22 = amp22ini*v2;
2748
2749 do {
2750
2751 count++;
2752 if (count >= length) {
2753 XLALPrintError("**** LALPhenSpinInspiralRD ERROR ****: phen. part exceeds array length");
2754 XLALPrintError( " m (%11.4e %11.4e) f0 %11.4e\n",params->mass1, params->mass2, params->fLower);
2755 XLALPrintError( " S1 (%8.4f %8.4f %8.4f)\n", initS1[0],initS1[1], initS1[2]);
2756 XLALPrintError( " S2 (%8.4f %8.4f %8.4f)\n", initS2[0],initS2[1], initS2[2]);
2757 XLALDestroyREAL8Vector(h2P2);
2758 XLALDestroyREAL8Vector(h2M2);
2759 XLALDestroyREAL8Vector(h2P1);
2760 XLALDestroyREAL8Vector(h2M1);
2761 XLALDestroyREAL8Vector(h20);
2762 XLALDestroyREAL8Vector(h3P3);
2763 XLALDestroyREAL8Vector(h3M3);
2764 XLALDestroyREAL8Vector(h3P2);
2765 XLALDestroyREAL8Vector(h3M2);
2766 XLALDestroyREAL8Vector(h3P1);
2767 XLALDestroyREAL8Vector(h3M1);
2768 XLALDestroyREAL8Vector(h30);
2769 XLALDestroyREAL8Vector(h4P4);
2770 XLALDestroyREAL8Vector(h4M4);
2771 XLALDestroyREAL8Vector(h4P3);
2772 XLALDestroyREAL8Vector(h4M3);
2773 XLALDestroyREAL8Vector(h4P2);
2774 XLALDestroyREAL8Vector(h4M2);
2775 XLALDestroyREAL8Vector(h4P1);
2776 XLALDestroyREAL8Vector(h4M1);
2777 XLALDestroyREAL8Vector(h40);
2778 XLALDestroyREAL8Vector(hap);
2779 XLALDestroyREAL8Vector(fap);
2780 XLALDestroyREAL8Vector(phap);
2781 XLALDestroyREAL8Vector(sigp);
2782 XLALDestroyREAL8Vector(sigc);
2783 XLALDestroyCOMPLEX8Vector(modefreqs);
2784 XLAL_ERROR(XLAL_ENOMEM);
2785 }
2786
2787 tim += dt;
2788 v2old = v2;
2789 //omold = om;
2790 om = om1 / (1. - tim / tAs) + om0;
2791 fap->data[count] = om;
2792 Psi = Psi0 + (- tAs * (om1/mass-dalpha1*trigAngle.ci) * log(1. - tim / tAs) + (om0/mass-dalpha0*trigAngle.ci) * (tim - t0) );
2793 //trigAngle.ci = cos(diota0 * (tim - t0) - diota1 * tAs * log(1. - tim / tAs) + iota0);
2794 alpha = alpha0 + ( dalpha0 * (tim - t0) - dalpha1 * tAs * log(1. - tim / tAs) );
2795
2796 v = cbrt(om);
2797 v2 = v*v;
2798 amp22 *= v2 / v2old;
2799
2800 amp33 = -amp22 / 4. * sqrt(5. / 42.);
2801 amp44 = amp22 * sqrt(5./7.) * 2./9. * v2;
2802
2803 errcode=XLALSpinInspiralFillH2Modes(h2P2,h2M2,h2P1,h2M1,h20,count,amp22,v,mparams.eta,mparams.dm,Psi,alpha,&trigAngle);
2804
2805 errcode += XLALSpinInspiralFillH3Modes(h3P3,h3M3,h3P2,h3M2,h3P1,h3M1,h30,count,amp33,v,mparams.eta,mparams.dm,Psi,alpha,&trigAngle);
2806
2807 errcode += XLALSpinInspiralFillH4Modes(h4P4,h4M4,h4P3,h4M3,h4P2,h4M2,h4P1,h4M1,h40,count,amp44,v,mparams.eta,mparams.dm,Psi,alpha,&trigAngle);
2808
2809 if (errcode != XLAL_SUCCESS)
2810 XLAL_ERROR(XLAL_EFUNC);
2811
2812 fap->data[count] = om;
2813 phap->data[count] = Psi;
2814
2815 } while ( (om < frOmRD) && (tim < tAs) );
2816
2817 XLALDestroyCOMPLEX8Vector(modefreqs);
2818 origcount=count;
2819 params->tC = ((REAL8) count) * dt;
2820
2821 /*--------------------------------------------------------------
2822 * Attach the ringdown waveform to the end of inspiral
2823 -------------------------------------------------------------*/
2824
2825 //printf("time %12.6e count %d\n",tim,count);
2826
2827 apcount = origcount;
2828 errcode = XLALPSpinInspiralAttachRingdownWave(h2P2, params, &apcount, nmodes, 2, 2, finalMass, finalSpin);
2829 for (i = 2 * apcount; i < 2 * length; i++) h2P2->data[i] = 0.;
2830 if (apcount > count) count = apcount;
2831
2832 apcount = origcount;
2833 errcode += XLALPSpinInspiralAttachRingdownWave(h2M2, params, &apcount, nmodes, 2, -2, finalMass, finalSpin);
2834 for (i = 2 * apcount; i < 2 * length; i++) h2M2->data[i] = 0.;
2835 if (apcount > count) count = apcount;
2836
2837 apcount = origcount;
2838 errcode += XLALPSpinInspiralAttachRingdownWave(h2P1, params, &apcount, nmodes, 2, 1, finalMass, finalSpin);
2839 for (i = 2 * apcount; i < 2 * length; i++) h2P1->data[i] = 0.;
2840 if (apcount > count) count = apcount;
2841
2842 apcount = origcount;
2843 errcode += XLALPSpinInspiralAttachRingdownWave(h2M1, params, &apcount, nmodes, 2, -1, finalMass, finalSpin);
2844 for (i = 2 * apcount; i < 2 * length; i++) h2M1->data[i] = 0.;
2845 if (apcount > count) count = apcount;
2846
2847 apcount = origcount;
2848 errcode += XLALPSpinInspiralAttachRingdownWave(h20, params, &apcount, nmodes, 2, 0, finalMass, finalSpin);
2849 for (i = 2 * apcount; i < 2 * length; i++) h20->data[i] = 0.;
2850 if (apcount > count) count = apcount;
2851
2852 apcount = origcount;
2853 errcode += XLALPSpinInspiralAttachRingdownWave(h3P3, params, &apcount, nmodes, 3, 3, finalMass, finalSpin);
2854 for (i = 2 * apcount; i < 2 * length; i++) h3P3->data[i] = 0.;
2855 if (apcount > count) count = apcount;
2856
2857 apcount = origcount;
2858 errcode += XLALPSpinInspiralAttachRingdownWave(h3M3, params, &apcount, nmodes, 3, -3, finalMass, finalSpin);
2859 for (i = 2 * apcount; i < 2 * length; i++) h3M3->data[i] = 0.;
2860 if (apcount > count) count = apcount;
2861
2862 apcount = origcount;
2863 errcode += XLALPSpinInspiralAttachRingdownWave(h3P2, params, &apcount, nmodes, 3, 2, finalMass, finalSpin);
2864 for (i = 2 * apcount; i < 2 * length; i++) h3P2->data[i] = 0.;
2865 if (apcount > count) count = apcount;
2866
2867 apcount = origcount;
2868 errcode += XLALPSpinInspiralAttachRingdownWave(h3M2, params, &apcount, nmodes, 3, -2, finalMass, finalSpin);
2869 for (i = 2 * apcount; i < 2 * length; i++) h3P2->data[i] = 0.;
2870 if (apcount > count) count = apcount;
2871
2872 apcount = origcount;
2873 errcode += XLALPSpinInspiralAttachRingdownWave(h3P1, params, &apcount, nmodes, 3, 1, finalMass, finalSpin);
2874 for (i = 2 * apcount; i < 2 * length; i++) h3P1->data[i] = 0.;
2875 if (apcount > count) count = apcount;
2876
2877 apcount = origcount;
2878 errcode += XLALPSpinInspiralAttachRingdownWave(h3M1, params, &apcount, nmodes, 3, -1, finalMass, finalSpin);
2879 for (i = 2 * apcount; i < 2 * length; i++) h3M1->data[i] = 0.;
2880 if (apcount > count) count = apcount;
2881
2882 apcount = origcount;
2883 errcode += XLALPSpinInspiralAttachRingdownWave(h30, params, &apcount, nmodes, 3, 0, finalMass, finalSpin);
2884 for (i = 2 * apcount; i < 2 * length; i++) h30->data[i] = 0.;
2885 if (apcount > count) count = apcount;
2886
2887 apcount = origcount;
2888 errcode += XLALPSpinInspiralAttachRingdownWave(h4P4, params, &apcount, nmodes, 4, 4, finalMass, finalSpin);
2889 for (i = 2 * apcount; i < 2 * length; i++) h4P4->data[i] = 0.;
2890 if (apcount > count) count = apcount;
2891
2892 apcount = origcount;
2893 errcode += XLALPSpinInspiralAttachRingdownWave(h4M4, params, &apcount, nmodes, 4, -4, finalMass, finalSpin);
2894 for (i = 2 * apcount; i < 2 * length; i++) h4M4->data[i] = 0.;
2895 if (apcount > count) count = apcount;
2896
2897 apcount = origcount;
2898 errcode += XLALPSpinInspiralAttachRingdownWave(h4P3, params, &apcount, nmodes, 4, 3, finalMass, finalSpin);
2899 for (i = 2 * apcount; i < 2 * length; i++) h4P3->data[i] = 0.;
2900 if (apcount > count) count = apcount;
2901
2902 apcount = origcount;
2903 errcode += XLALPSpinInspiralAttachRingdownWave(h4M3, params, &apcount, nmodes, 4, -3, finalMass, finalSpin);
2904 for (i = 2 * apcount; i < 2 * length; i++) h4M3->data[i] = 0.;
2905 if (apcount > count) count = apcount;
2906
2907 apcount = origcount;
2908 errcode += XLALPSpinInspiralAttachRingdownWave(h4P2, params, &apcount, nmodes, 4, 2, finalMass, finalSpin);
2909 for (i = 2 * apcount; i < 2 * length; i++) h4P4->data[i] = 0.;
2910 if (apcount > count) count = apcount;
2911
2912 apcount = origcount;
2913 errcode += XLALPSpinInspiralAttachRingdownWave(h4M2, params, &apcount, nmodes, 4, -2, finalMass, finalSpin);
2914 for (i = 2 * apcount; i < 2 * length; i++) h4M4->data[i] = 0.;
2915 if (apcount > count) count = apcount;
2916
2917 apcount = origcount;
2918 errcode += XLALPSpinInspiralAttachRingdownWave(h4P1, params, &apcount, nmodes, 4, 1, finalMass, finalSpin);
2919 for (i = 2 * apcount; i < 2 * length; i++) h4P3->data[i] = 0.;
2920 if (apcount > count) count = apcount;
2921
2922 apcount = origcount;
2923 errcode += XLALPSpinInspiralAttachRingdownWave(h4M1, params, &apcount, nmodes, 4, -1, finalMass, finalSpin);
2924 for (i = 2 * apcount; i < 2 * length; i++) h4M3->data[i] = 0.;
2925 if (apcount > count) count = apcount;
2926
2927 apcount = origcount;
2928 errcode += XLALPSpinInspiralAttachRingdownWave(h40, params, &apcount, nmodes, 4, 0, finalMass, finalSpin);
2929 for (i = 2 * apcount; i < 2 * length; i++) h40->data[i] = 0.;
2930 if (apcount > count) count = apcount;
2931
2932 if (errcode != XLAL_SUCCESS) {
2933 XLALPrintError("**** LALPSpinInspiralRD ERROR ****: impossible to create RingDownWave\n");
2934 XLALDestroyREAL8Vector(h2P2);
2935 XLALDestroyREAL8Vector(h2M2);
2936 XLALDestroyREAL8Vector(h2P1);
2937 XLALDestroyREAL8Vector(h2M1);
2938 XLALDestroyREAL8Vector(h20);
2939 XLALDestroyREAL8Vector(h3P3);
2940 XLALDestroyREAL8Vector(h3M3);
2941 XLALDestroyREAL8Vector(h3P2);
2942 XLALDestroyREAL8Vector(h3M2);
2943 XLALDestroyREAL8Vector(h3P1);
2944 XLALDestroyREAL8Vector(h3M1);
2945 XLALDestroyREAL8Vector(h30);
2946 XLALDestroyREAL8Vector(h4P4);
2947 XLALDestroyREAL8Vector(h4M4);
2948 XLALDestroyREAL8Vector(h4P3);
2949 XLALDestroyREAL8Vector(h4M3);
2950 XLALDestroyREAL8Vector(h4P2);
2951 XLALDestroyREAL8Vector(h4M2);
2952 XLALDestroyREAL8Vector(h4P1);
2953 XLALDestroyREAL8Vector(h4M1);
2954 XLALDestroyREAL8Vector(h40);
2955 XLALDestroyREAL8Vector(hap);
2956 XLALDestroyREAL8Vector(fap);
2957 XLALDestroyREAL8Vector(phap);
2958 XLAL_ERROR(XLAL_EFAILED);
2959 }
2960 }
2961
2962 } /*End of if not inspiralonly and test_omegamatch*/
2963
2964 /*-------------------------------------------------------------------
2965 * Compute the spherical harmonics required for constructing (h+,hx).
2966 -------------------------------------------------------------------*/
2967
2968 /* The angles theta for the spherical harmonics has been set according to
2969 the input inclination parameter and the axisChoice */
2970
2971 for (i = 0; i < length; i++) {
2972 fap->data[i] /= unitHz;
2973 sigp->data[i] = 0.;
2974 sigc->data[i] = 0.;
2975 }
2976
2977 errcode = XLALSphHarm(&MultSphHarmP, 2, 2, inc, 0.);
2978 errcode += XLALSphHarm(&MultSphHarmM, 2, -2, inc, 0.);
2979 if (errcode != XLAL_SUCCESS) {
2980 XLALDestroyREAL8Vector(h2P2);
2981 XLALDestroyREAL8Vector(h2M2);
2982 XLALDestroyREAL8Vector(h2P1);
2983 XLALDestroyREAL8Vector(h2M1);
2984 XLALDestroyREAL8Vector(h20);
2985 XLALDestroyREAL8Vector(h3P3);
2986 XLALDestroyREAL8Vector(h3M3);
2987 XLALDestroyREAL8Vector(h3P2);
2988 XLALDestroyREAL8Vector(h3M2);
2989 XLALDestroyREAL8Vector(h3P1);
2990 XLALDestroyREAL8Vector(h3M1);
2991 XLALDestroyREAL8Vector(h30);
2992 XLALDestroyREAL8Vector(h4P4);
2993 XLALDestroyREAL8Vector(h4M4);
2994 XLALDestroyREAL8Vector(h4P3);
2995 XLALDestroyREAL8Vector(h4M3);
2996 XLALDestroyREAL8Vector(h4P2);
2997 XLALDestroyREAL8Vector(h4M2);
2998 XLALDestroyREAL8Vector(h4P1);
2999 XLALDestroyREAL8Vector(h4M1);
3000 XLALDestroyREAL8Vector(h40);
3001 XLALDestroyREAL8Vector(hap);
3002 XLALDestroyREAL8Vector(fap);
3003 XLALDestroyREAL8Vector(phap);
3004 XLALPrintError("**** LALPSpinInspiralRD ERROR ****: impossible to create Y22 or Y2-2\n");
3005 XLAL_ERROR(XLAL_EFAILED);
3006 }
3007 for (i = 0; i < length; i++) {
3008 x0 = h2P2->data[2 * i];
3009 x1 = h2P2->data[2 * i + 1];
3010 x2 = h2M2->data[2 * i];
3011 x3 = h2M2->data[2 * i + 1];
3012 sigp->data[i] += x0 * creal(MultSphHarmP) - x1 * cimag(MultSphHarmP) + x2 * creal(MultSphHarmM) - x3 * cimag(MultSphHarmM);
3013 sigc->data[i] += - x0 * cimag(MultSphHarmP) - x1 * creal(MultSphHarmP) - x2 * cimag(MultSphHarmM) - x3 * creal(MultSphHarmM);
3014 }
3015
3016 errcode = XLALSphHarm(&MultSphHarmP, 2, 1, inc, 0.);
3017 errcode += XLALSphHarm(&MultSphHarmM, 2, -1, inc, 0.);
3018 if (errcode != XLAL_SUCCESS){
3019 XLALDestroyREAL8Vector(h2P1);
3020 XLALDestroyREAL8Vector(h2M1);
3021 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: impossible to create Y21\n");
3022 } else {
3023 for (i = 0; i < length; i++) {
3024 x0 = h2P1->data[2 * i];
3025 x1 = h2P1->data[2 * i + 1];
3026 x2 = h2M1->data[2 * i];
3027 x3 = h2M1->data[2 * i + 1];
3028 sigp->data[i] += x0 * creal(MultSphHarmP) - x1 * cimag(MultSphHarmP) + x2 * creal(MultSphHarmM) - x3 * cimag(MultSphHarmM);
3029 sigc->data[i] += - x0 * cimag(MultSphHarmP) - x1 * creal(MultSphHarmP) - x2 * cimag(MultSphHarmM) - x3 * creal(MultSphHarmM);
3030 }
3031 }
3032
3033 errcode = XLALSphHarm(&MultSphHarmP, 2, 0, inc, 0.);
3034 if (errcode != XLAL_SUCCESS) {
3035 XLALDestroyREAL8Vector(h20);
3036 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: impossible to create Y20\n");
3037 } else {
3038 for (i = 0; i < length; i++) {
3039 x0 = h20->data[2 * i];
3040 x1 = h20->data[2 * i + 1];
3041 sigp->data[i] += x1 * creal(MultSphHarmP) - x1 * cimag(MultSphHarmP);
3042 sigc->data[i] -= x1 * cimag(MultSphHarmP) + x1 * creal(MultSphHarmP);
3043 }
3044 }
3045
3046 errcode = XLALSphHarm(&MultSphHarmP, 3, 3, inc, 0.);
3047 errcode += XLALSphHarm(&MultSphHarmM, 3, -3, inc, 0.);
3048
3049 if (errcode != XLAL_SUCCESS) {
3050 XLALDestroyREAL8Vector(h3P3);
3051 XLALDestroyREAL8Vector(h3M3);
3052 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: impossible to create Y33,Y3-3\n");
3053 } else {
3054 for (i = 0; i < length; i++) {
3055 x0 = h3P3->data[2 * i];
3056 x1 = h3P3->data[2 * i + 1];
3057 x2 = h3M3->data[2 * i];
3058 x3 = h3M3->data[2 * i + 1];
3059 sigp->data[i] += x0 * creal(MultSphHarmP) - x1 * cimag(MultSphHarmP) + x2 * creal(MultSphHarmM) - x3 * cimag(MultSphHarmM);
3060 sigc->data[i] -= x0 * cimag(MultSphHarmP) + x1 * creal(MultSphHarmP) + x2 * cimag(MultSphHarmM) + x3 * creal(MultSphHarmM);
3061 }
3062 }
3063
3064 errcode = XLALSphHarm(&MultSphHarmP, 3, 2, inc, 0.);
3065 errcode += XLALSphHarm(&MultSphHarmM, 3, -2, inc, 0.);
3066 if (errcode != XLAL_SUCCESS) {
3067 XLALDestroyREAL8Vector(h3P2);
3068 XLALDestroyREAL8Vector(h3M2);
3069 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: impossible to create Y32,Y3-2\n");
3070 } else {
3071 for (i = 0; i < length; i++) {
3072 x0 = h3P2->data[2 * i];
3073 x1 = h3P2->data[2 * i + 1];
3074 x2 = h3M2->data[2 * i];
3075 x3 = h3M2->data[2 * i + 1];
3076 sigp->data[i] += x0 * creal(MultSphHarmP) - x1 * cimag(MultSphHarmP) + x2 * creal(MultSphHarmM) - x3 * cimag(MultSphHarmM);
3077 sigc->data[i] -= x0 * cimag(MultSphHarmP) + x1 * creal(MultSphHarmP) + x2 * cimag(MultSphHarmM) + x3 * creal(MultSphHarmM);
3078 }
3079 }
3080
3081 errcode = XLALSphHarm(&MultSphHarmP, 3, 1, inc, 0.);
3082 errcode += XLALSphHarm(&MultSphHarmM, 3, -1, inc, 0.);
3083 if (errcode != XLAL_SUCCESS) {
3084 XLALDestroyREAL8Vector(h3P1);
3085 XLALDestroyREAL8Vector(h3M1);
3086 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: impossible to create Y31,Y3-1\n");
3087 } else {
3088 for (i = 0; i < length; i++) {
3089 x0 = h3P1->data[2 * i];
3090 x1 = h3P1->data[2 * i + 1];
3091 x2 = h3M1->data[2 * i];
3092 x3 = h3M1->data[2 * i + 1];
3093 sigp->data[i] += x0 * creal(MultSphHarmP) - x1 * cimag(MultSphHarmP) + x2 * creal(MultSphHarmM) - x3 * cimag(MultSphHarmM);
3094 sigc->data[i] -= x0 * cimag(MultSphHarmP) + x1 * creal(MultSphHarmP) + x2 * cimag(MultSphHarmM) + x3 * creal(MultSphHarmM);
3095 }
3096 }
3097
3098 errcode = XLALSphHarm(&MultSphHarmP, 3, 0, inc, 0.);
3099 if (errcode != XLAL_SUCCESS) {
3100 XLALDestroyREAL8Vector(h30);
3101 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: impossible to create Y30\n");
3102 } else {
3103 for (i = 0; i < length; i++) {
3104 x0 = h30->data[2 * i];
3105 x1 = h30->data[2 * i + 1];
3106 sigp->data[i] += x0 * creal(MultSphHarmP) - x1 * cimag(MultSphHarmM);
3107 sigc->data[i] -= x0 * cimag(MultSphHarmP) + x1 * creal(MultSphHarmP);
3108 }
3109 }
3110
3111 errcode = XLALSphHarm(&MultSphHarmP, 4, 4, inc, 0.);
3112 errcode += XLALSphHarm(&MultSphHarmM, 4, -4, inc, 0.);
3113 if (errcode != XLAL_SUCCESS) {
3114 XLALDestroyREAL8Vector(h4P4);
3115 XLALDestroyREAL8Vector(h4M4);
3116 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: impossible to create Y44,Y4-4\n");
3117 } else {
3118 for (i = 0; i < length; i++) {
3119 x0 = h4P4->data[2 * i];
3120 x1 = h4P4->data[2 * i + 1];
3121 x2 = h4P4->data[2 * i];
3122 x3 = h4M4->data[2 * i + 1];
3123 sigp->data[i] += x0 * creal(MultSphHarmP) - x1 * cimag(MultSphHarmP) + x2 * creal(MultSphHarmM) - x3 * cimag(MultSphHarmM);
3124 sigc->data[i] -= x0 * cimag(MultSphHarmP) + x1 * creal(MultSphHarmP) + x2 * cimag(MultSphHarmM) + x3 * creal(MultSphHarmM);
3125 }
3126 }
3127
3128 errcode = XLALSphHarm(&MultSphHarmP, 4, 3, inc, 0.);
3129 errcode += XLALSphHarm(&MultSphHarmM, 4, -3, inc, 0.);
3130 if (errcode != XLAL_SUCCESS) {
3131 XLALDestroyREAL8Vector(h4P3);
3132 XLALDestroyREAL8Vector(h4M3);
3133 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: impossible to create Y43,Y4-3\n");
3134 } else {
3135 for (i = 0; i < length; i++) {
3136 x0 = h4P3->data[2 * i];
3137 x1 = h4P3->data[2 * i + 1];
3138 x2 = h4M3->data[2 * i];
3139 x3 = h4M3->data[2 * i + 1];
3140 sigp->data[i] += x0 * creal(MultSphHarmP) - x1 * cimag(MultSphHarmP) + x2 * creal(MultSphHarmM) - x3 * cimag(MultSphHarmM);
3141 sigc->data[i] -= x0 * cimag(MultSphHarmP) + x1 * creal(MultSphHarmP) + x2 * cimag(MultSphHarmM) + x3 * creal(MultSphHarmM);
3142 }
3143 }
3144
3145 errcode = XLALSphHarm(&MultSphHarmP, 4, 2, inc, 0.);
3146 errcode += XLALSphHarm(&MultSphHarmM, 4, -2, inc, 0.);
3147 if (errcode != XLAL_SUCCESS) {
3148 XLALDestroyREAL8Vector(h4P2);
3149 XLALDestroyREAL8Vector(h4M2);
3150 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: impossible to create Y42,Y4-2\n");
3151 } else {
3152 for (i = 0; i < length; i++) {
3153 x0 = h4P2->data[2 * i];
3154 x1 = h4P2->data[2 * i + 1];
3155 x2 = h4M2->data[2 * i];
3156 x3 = h4M2->data[2 * i + 1];
3157 sigp->data[i] += x0 * creal(MultSphHarmP) - x1 * cimag(MultSphHarmP) + x2 * creal(MultSphHarmM) - x3 * cimag(MultSphHarmM);
3158 sigc->data[i] -= x0 * cimag(MultSphHarmP) + x1 * creal(MultSphHarmP) + x2 * cimag(MultSphHarmM) + x3 * creal(MultSphHarmM);
3159 }
3160 }
3161
3162 errcode = XLALSphHarm(&MultSphHarmP, 4, 1, inc, 0.);
3163 errcode += XLALSphHarm(&MultSphHarmM, 4, -1, inc, 0.);
3164 if (errcode != XLAL_SUCCESS) {
3165 XLALDestroyREAL8Vector(h4P1);
3166 XLALDestroyREAL8Vector(h4M1);
3167 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: impossible to create Y41,Y4-1\n");
3168 } else {
3169 for (i = 0; i < length; i++) {
3170 x0 = h4P1->data[2 * i];
3171 x1 = h4P1->data[2 * i + 1];
3172 x2 = h4M1->data[2 * i];
3173 x3 = h4M1->data[2 * i + 1];
3174 sigp->data[i] += x0 * creal(MultSphHarmP) - x1 * cimag(MultSphHarmP) + x2 * creal(MultSphHarmM) - x3 * cimag(MultSphHarmM);
3175 sigc->data[i] -= x0 * cimag(MultSphHarmP) + x1 * creal(MultSphHarmP) + x2 * cimag(MultSphHarmM) + x3 * creal(MultSphHarmM);
3176 }
3177 }
3178
3179 errcode = XLALSphHarm(&MultSphHarmP, 4, 0, inc, 0.);
3180 if (errcode != XLAL_SUCCESS) {
3181 XLALDestroyREAL8Vector(h40);
3182 XLALPrintWarning("** LALPSpinInspiralRD WARNING **: impossible to create Y40\n");
3183 } else {
3184 for (i = 0; i < length; i++) {
3185 x0 = h40->data[2 * i];
3186 x1 = h40->data[2 * i + 1];
3187 sigp->data[i] += x0 * creal(MultSphHarmP) - x1 * cimag(MultSphHarmP);
3188 sigc->data[i] -= x0 * cimag(MultSphHarmP) + x1 * creal(MultSphHarmP);
3189 }
3190 }
3191
3192 params->fFinal = params->tSampling / 2.;
3193
3194 /*------------------------------------------------------
3195 * If required by the user copy other data sets to the
3196 * relevant arrays
3197 ------------------------------------------------------*/
3198
3199 if (hh) {
3200 for (i = 0; i < length; i++) {
3201 j = 2 * i;
3202 k = 2 * i + 1;
3203 hap->data[j] = sigp->data[i];
3204 hap->data[k] = sigc->data[i];
3205 }
3206 }
3207
3208 if (signalvec1)
3209 memcpy(signalvec1->data, sigp->data, length * (sizeof(REAL8)));
3210 if (signalvec2)
3211 memcpy(signalvec2->data, sigc->data, length * (sizeof(REAL8)));
3212 if (hh)
3213 memcpy(hh->data, hap->data, 2 * length * (sizeof(REAL8)));
3214 if (ff)
3215 memcpy(ff->data, fap->data, length * (sizeof(REAL8)));
3216 if (phi)
3217 memcpy(phi->data, phap->data, length * (sizeof(REAL8)));
3218
3219 /* Clean up */
3220 XLALDestroyREAL8Vector(h2P2);
3221 XLALDestroyREAL8Vector(h2M2);
3222 XLALDestroyREAL8Vector(h2P1);
3223 XLALDestroyREAL8Vector(h2M1);
3224 XLALDestroyREAL8Vector(h20);
3225 XLALDestroyREAL8Vector(h3P3);
3226 XLALDestroyREAL8Vector(h3M3);
3227 XLALDestroyREAL8Vector(h3P2);
3228 XLALDestroyREAL8Vector(h3M2);
3229 XLALDestroyREAL8Vector(h3P1);
3230 XLALDestroyREAL8Vector(h3M1);
3231 XLALDestroyREAL8Vector(h30);
3232 XLALDestroyREAL8Vector(h4P4);
3233 XLALDestroyREAL8Vector(h4M4);
3234 XLALDestroyREAL8Vector(h4P3);
3235 XLALDestroyREAL8Vector(h4M3);
3236 XLALDestroyREAL8Vector(h4P2);
3237 XLALDestroyREAL8Vector(h4M2);
3238 XLALDestroyREAL8Vector(h4P1);
3239 XLALDestroyREAL8Vector(h4M1);
3240 XLALDestroyREAL8Vector(h40);
3241 XLALDestroyREAL8Vector(fap);
3242 XLALDestroyREAL8Vector(phap);
3243 XLALDestroyREAL8Vector(hap);
3244 XLALDestroyREAL8Vector(sigp);
3245 XLALDestroyREAL8Vector(sigc);
3246
3247 return count;
3248
3249 /*End */
3250}
3251
3252/** @} */
XLALRungeKutta4Free
void XLALRungeKutta4Free(rk4GSLIntegrator *integrator)
Definition: LALRungeKutta4.c:228
XLALPSpinGenerateQNMFreq
INT4 XLALPSpinGenerateQNMFreq(COMPLEX8Vector *modefreqs, InspiralTemplate *params, UINT4 l, INT4 m, UINT4 nmodes, REAL8 finalMass, REAL8 finalSpin)
Definition: LALPSpinInspiralRingdownWave.c:356
XLALInspiralSetup
int XLALInspiralSetup(expnCoeffs *ak, InspiralTemplate *params)
Definition: LALInspiralSetup.c:88
XLALInspiralInit
int XLALInspiralInit(InspiralTemplate *params, InspiralInit *paramsInit)
Definition: LALInspiralInit.c:87
XLALPSpinFinalMassSpin
INT4 XLALPSpinFinalMassSpin(REAL8 *finalMass, REAL8 *finalSpin, InspiralTemplate *params, REAL8 energy, REAL8 *LNhvec)
Definition: LALPSpinInspiralRingdownWave.c:568
XLALGenerateWaveDerivative
INT4 XLALGenerateWaveDerivative(REAL8Vector *dwave, REAL8Vector *wave, REAL8 dt)
Definition: LALPSpinInspiralRingdownWave.c:268
XLALPSpinInspiralAttachRingdownWave
INT4 XLALPSpinInspiralAttachRingdownWave(REAL8Vector *signalvec, InspiralTemplate *params, UINT4 *attpos, UINT4 nmodes, UINT4 l, INT4 m, REAL8 finalMass, REAL8 finalSpin)
Definition: LALPSpinInspiralRingdownWave.c:644
XLALInspiralChooseModel
int XLALInspiralChooseModel(expnFunc *func, expnCoeffs *ak, InspiralTemplate *params)
Definition: LALInspiralChooseModel.c:506
XLALRungeKutta4Init
rk4GSLIntegrator * XLALRungeKutta4Init(INT4 n, rk4In *input)
Definition: LALRungeKutta4.c:81
XLALRungeKutta4
int XLALRungeKutta4(REAL8Vector *yout, rk4GSLIntegrator *integrator, void *params)
Definition: LALRungeKutta4.c:147
LALMalloc
#define LALMalloc(n)
LALFree
#define LALFree(p)
c1
const double c1
c2
const double c2
XLAL_CALLGSL
#define XLAL_CALLGSL(statement)
s
int s
l
int l
i
double i
__attribute__
#define __attribute__(x)
XLALDestroyREAL8Array
void XLALDestroyREAL8Array(REAL8Array *)
XLALCreateREAL4VectorSequence
REAL4VectorSequence * XLALCreateREAL4VectorSequence(UINT4 length, UINT4 veclen)
GENERATEPPNINSPIRALH_EFSTOP
#define GENERATEPPNINSPIRALH_EFSTOP
Reached requested termination frequency.
Definition: GeneratePPNInspiral.h:202
XLALAdaptiveRungeKuttaFree
void XLALAdaptiveRungeKuttaFree(LALAdaptiveRungeKuttaIntegrator *integrator)
XLALAdaptiveRungeKutta4Init
LALAdaptiveRungeKuttaIntegrator * XLALAdaptiveRungeKutta4Init(int dim, int(*dydt)(double t, const double y[], double dydt[], void *params), int(*stop)(double t, const double y[], double dydt[], void *params), double eps_abs, double eps_rel)
XLALAdaptiveRungeKutta4Hermite
int XLALAdaptiveRungeKutta4Hermite(LALAdaptiveRungeKuttaIntegrator *integrator, void *params, REAL8 *yinit, REAL8 tinit, REAL8 tend_in, REAL8 deltat, REAL8Array **yout)
LAL_PI
#define LAL_PI
LAL_MTSUN_SI
#define LAL_MTSUN_SI
LAL_MRSUN_SI
#define LAL_MRSUN_SI
COMPLEX16
double complex COMPLEX16
REAL8
double REAL8
UINT4
uint32_t UINT4
INT4
int32_t INT4
REAL4
float REAL4
LALNameLength
LALNameLength
LAL_INSPIRAL_INTERACTION_TIDAL_5PN
@ LAL_INSPIRAL_INTERACTION_TIDAL_5PN
Leading-order tidal interaction.
Definition: LALInspiral.h:124
LAL_INSPIRAL_INTERACTION_ALL
@ LAL_INSPIRAL_INTERACTION_ALL
all spin and tidal interactions
Definition: LALInspiral.h:127
LAL_INSPIRAL_INTERACTION_TIDAL_6PN
@ LAL_INSPIRAL_INTERACTION_TIDAL_6PN
Next-to-leading-order tidal interaction.
Definition: LALInspiral.h:125
LAL_INSPIRAL_INTERACTION_SPIN_ORBIT_25PN
@ LAL_INSPIRAL_INTERACTION_SPIN_ORBIT_25PN
Next-to-leading-order spin-orbit interaction.
Definition: LALInspiral.h:122
LAL_INSPIRAL_INTERACTION_QUAD_MONO_2PN
@ LAL_INSPIRAL_INTERACTION_QUAD_MONO_2PN
Quadrupole-monopole interaction.
Definition: LALInspiral.h:121
LAL_INSPIRAL_INTERACTION_SPIN_ORBIT_3PN
@ LAL_INSPIRAL_INTERACTION_SPIN_ORBIT_3PN
Spin-spin interaction.
Definition: LALInspiral.h:123
LAL_INSPIRAL_INTERACTION_NONE
@ LAL_INSPIRAL_INTERACTION_NONE
No spin, tidal or other interactions.
Definition: LALInspiral.h:117
LAL_INSPIRAL_INTERACTION_SPIN_SPIN_2PN
@ LAL_INSPIRAL_INTERACTION_SPIN_SPIN_2PN
Spin-spin interaction.
Definition: LALInspiral.h:119
LAL_INSPIRAL_INTERACTION_SPIN_ORBIT_15PN
@ LAL_INSPIRAL_INTERACTION_SPIN_ORBIT_15PN
Leading order spin-orbit interaction.
Definition: LALInspiral.h:118
LAL_INSPIRAL_INTERACTION_ALL_SPIN
@ LAL_INSPIRAL_INTERACTION_ALL_SPIN
all spin interactions, no tidal interactions
Definition: LALInspiral.h:126
LAL_INSPIRAL_INTERACTION_SPIN_SPIN_SELF_2PN
@ LAL_INSPIRAL_INTERACTION_SPIN_SPIN_SELF_2PN
Spin-spin-self interaction.
Definition: LALInspiral.h:120
XLALMalloc
void * XLALMalloc(size_t n)
XLALFree
void XLALFree(void *p)
XLALSpinInspiralFillH4Modes
static int XLALSpinInspiralFillH4Modes(REAL8Vector *h4P4, REAL8Vector *h4M4, REAL8Vector *h4P3, REAL8Vector *h4M3, REAL8Vector *h4P2, REAL8Vector *h4M2, REAL8Vector *h4P1, REAL8Vector *h4M1, REAL8Vector *h40, INT4 j, REAL8 amp44, REAL8 v, REAL8 eta, REAL8 dm, REAL8 Psi, REAL8 alpha, LALSpinInspiralAngle *an)
Definition: LALPSpinInspiralRD.c:1263
LALPSIRDPN_TEST_OMEGANAN
#define LALPSIRDPN_TEST_OMEGANAN
Definition: LALPSpinInspiralRD.c:109
LALPSIRDPN_TEST_OMEGACUT
#define LALPSIRDPN_TEST_OMEGACUT
Definition: LALPSpinInspiralRD.c:112
LALPSIRDPN_TEST_OMEGANONPOS
#define LALPSIRDPN_TEST_OMEGANONPOS
Definition: LALPSpinInspiralRD.c:111
fracRD
static REAL8 fracRD(REAL8 LNhS1, REAL8 LNhS2, REAL8 S1S1, REAL8 S1S2, REAL8 S2S2)
Definition: LALPSpinInspiralRD.c:146
OmMatch
static REAL8 OmMatch(REAL8 LNhS1, REAL8 LNhS2, REAL8 S1S1, REAL8 S1S2, REAL8 S2S2)
Definition: LALPSpinInspiralRD.c:121
XLALSpinInspiralAdaptiveEngine
static int XLALSpinInspiralAdaptiveEngine(const UINT4 neqs, const REAL8 yinit[], REAL8 amp22ini, LALPSpinInspiralRDparams *mparams, REAL8Vector *h2P2, REAL8Vector *h2M2, REAL8Vector *h2P1, REAL8Vector *h2M1, REAL8Vector *h20, REAL8Vector *h3P3, REAL8Vector *h3M3, REAL8Vector *h3P2, REAL8Vector *h3M2, REAL8Vector *h3P1, REAL8Vector *h3M1, REAL8Vector *h30, REAL8Vector *h4P4, REAL8Vector *h4M4, REAL8Vector *h4P3, REAL8Vector *h4M3, REAL8Vector *h4P2, REAL8Vector *h4M2, REAL8Vector *h4P1, REAL8Vector *h4M1, REAL8Vector *h40, REAL8Vector *freq, REAL8Vector *phase, LALPSpinInspiralPhenPars *phenPars)
Definition: LALPSpinInspiralRD.c:1729
XLALPSpinInspiralRDEngine
static int XLALPSpinInspiralRDEngine(REAL8Vector *signalvec1, REAL8Vector *signalvec2, REAL8Vector *hh, REAL8Vector *ff, REAL8Vector *phi, InspiralTemplate *params, InspiralInit *paramsInit)
Main function to produce waveforms.
Definition: LALPSpinInspiralRD.c:2126
sqrtOnePointFive
#define sqrtOnePointFive
Definition: LALPSpinInspiralRD.c:114
XLALPSpinInspiralRDSetParams
static int XLALPSpinInspiralRDSetParams(LALPSpinInspiralRDparams *mparams, InspiralTemplate *params, InspiralInit *paramsInit)
Definition: LALPSpinInspiralRD.c:255
XLALSpinInspiralEngine
static int XLALSpinInspiralEngine(UINT4 neqs, const REAL8 yinit[], REAL8 amp22ini, LALPSpinInspiralRDparams *mparams, REAL8Vector *h2P2, REAL8Vector *h2M2, REAL8Vector *h2P1, REAL8Vector *h2M1, REAL8Vector *h20, REAL8Vector *h3P3, REAL8Vector *h3M3, REAL8Vector *h3P2, REAL8Vector *h3M2, REAL8Vector *h3P1, REAL8Vector *h3M1, REAL8Vector *h30, REAL8Vector *h4P4, REAL8Vector *h4M4, REAL8Vector *h4P3, REAL8Vector *h4M3, REAL8Vector *h4P2, REAL8Vector *h4M2, REAL8Vector *h4P1, REAL8Vector *h4M1, REAL8Vector *h40, REAL8Vector *freq, REAL8Vector *phase, LALPSpinInspiralPhenPars *phenPars)
Definition: LALPSpinInspiralRD.c:1366
sqrtPoint15
#define sqrtPoint15
Definition: LALPSpinInspiralRD.c:115
XLALSpinInspiralTest
static int XLALSpinInspiralTest(double t, const double values[], double dvalues[], void *mparams)
Definition: LALPSpinInspiralRD.c:451
XLALSpinInspiralFillH3Modes
static int XLALSpinInspiralFillH3Modes(REAL8Vector *h3P3, REAL8Vector *h3M3, REAL8Vector *h3P2, REAL8Vector *h3M2, REAL8Vector *h3P1, REAL8Vector *h3M1, REAL8Vector *h30, UINT4 j, REAL4 amp, REAL4 v, REAL4 eta, REAL4 dm, REAL8 Psi, REAL8 alpha, LALSpinInspiralAngle *an)
Definition: LALPSpinInspiralRD.c:1121
sqrtFiveOver2
#define sqrtFiveOver2
Definition: LALPSpinInspiralRD.c:116
LALSpinInspiralDerivatives
static void LALSpinInspiralDerivatives(REAL8Vector *values, REAL8Vector *dvalues, void *mparams)
Definition: LALPSpinInspiralRD.c:694
XLALSpinInspiralDerivatives
static int XLALSpinInspiralDerivatives(double t, const double values[], double dvalues[], void *mparams)
Function to compute detivative of dynamical variables.
Definition: LALPSpinInspiralRD.c:505
XLALPSpinInspiralRDTemplates
int XLALPSpinInspiralRDTemplates(REAL4Vector *signalvec1, REAL4Vector *signalvec2, InspiralTemplate *params)
Function to produce waveform templates.
Definition: LALPSpinInspiralRD.c:751
LALPSIRDPN_TEST_ENERGY
#define LALPSIRDPN_TEST_ENERGY
Definition: LALPSpinInspiralRD.c:107
LALPSIRDPN_TEST_OMEGAMATCH
#define LALPSIRDPN_TEST_OMEGAMATCH
Definition: LALPSpinInspiralRD.c:110
minIntLen
#define minIntLen
Definition: LALPSpinInspiralRD.c:117
LALPSIRDPN_TEST_OMEGADOT
#define LALPSIRDPN_TEST_OMEGADOT
Definition: LALPSpinInspiralRD.c:108
XLALSpinInspiralFillH2Modes
static int XLALSpinInspiralFillH2Modes(REAL8Vector *h2P2, REAL8Vector *h2M2, REAL8Vector *h2P1, REAL8Vector *h2M1, REAL8Vector *h20, UINT4 j, REAL4 amp, REAL4 v, REAL4 eta, REAL4 dm, REAL8 Psi, REAL8 alpha, LALSpinInspiralAngle *an)
Function actually computing PSIRD waveforms.
Definition: LALPSpinInspiralRD.c:1046
XLALPSpinInspiralRD
int XLALPSpinInspiralRD(REAL4Vector *signalvec, InspiralTemplate *params)
Definition: LALPSpinInspiralRD.c:712
XLALPSpinInspiralRDFreqDom
int XLALPSpinInspiralRDFreqDom(REAL4Vector *signalvec, InspiralTemplate *params)
Definition: LALPSpinInspiralRD.c:928
XLALPSpinInspiralRDForInjection
int XLALPSpinInspiralRDForInjection(CoherentGW *waveform, InspiralTemplate *params, PPNParamStruc *ppnParams)
Function Module to produce injection waveforms.
Definition: LALPSpinInspiralRD.c:804
LALSimInspiralApplyTaper
LALSimInspiralApplyTaper
LALPNOrder
LALPNOrder
LAL_SIM_INSPIRAL_FRAME_AXIS_VIEW
LAL_SIM_INSPIRAL_FRAME_AXIS_VIEW
LAL_SIM_INSPIRAL_FRAME_AXIS_TOTAL_J
LAL_SIM_INSPIRAL_FRAME_AXIS_TOTAL_J
LAL_SIM_INSPIRAL_FRAME_AXIS_ORBITAL_L
LAL_SIM_INSPIRAL_FRAME_AXIS_ORBITAL_L
LAL_PNORDER_TWO_POINT_FIVE
LAL_PNORDER_TWO_POINT_FIVE
LAL_PNORDER_NUM_ORDER
LAL_PNORDER_NUM_ORDER
LAL_PNORDER_THREE
LAL_PNORDER_THREE
LAL_PNORDER_TWO
LAL_PNORDER_TWO
LAL_PNORDER_ONE
LAL_PNORDER_ONE
LAL_PNORDER_PSEUDO_FOUR
LAL_PNORDER_PSEUDO_FOUR
LAL_PNORDER_THREE_POINT_FIVE
LAL_PNORDER_THREE_POINT_FIVE
LAL_PNORDER_HALF
LAL_PNORDER_HALF
LAL_PNORDER_ONE_POINT_FIVE
LAL_PNORDER_ONE_POINT_FIVE
LAL_PNORDER_NEWTONIAN
LAL_PNORDER_NEWTONIAN
XLALSimInspiralREAL4WaveTaper
int XLALSimInspiralREAL4WaveTaper(REAL4Vector *signalvec, LALSimInspiralApplyTaper bookends)
XLALCreateForwardREAL4FFTPlan
REAL4FFTPlan * XLALCreateForwardREAL4FFTPlan(UINT4 size, int measurelvl)
XLALDestroyREAL4FFTPlan
void XLALDestroyREAL4FFTPlan(REAL4FFTPlan *plan)
XLALREAL4VectorFFT
int XLALREAL4VectorFFT(REAL4Vector *_LAL_RESTRICT_ output, const REAL4Vector *_LAL_RESTRICT_ input, const REAL4FFTPlan *plan)
XLALSphHarm
INT4 XLALSphHarm(COMPLEX16 *out, UINT4 L, INT4 M, REAL4 theta, REAL4 phi)
lalStrainUnit
const LALUnit lalStrainUnit
lalHertzUnit
const LALUnit lalHertzUnit
lalDimensionlessUnit
const LALUnit lalDimensionlessUnit
XLALCreateREAL4Vector
REAL4Vector * XLALCreateREAL4Vector(UINT4 length)
XLALDestroyREAL4Vector
void XLALDestroyREAL4Vector(REAL4Vector *vector)
XLALCreateCOMPLEX8Vector
COMPLEX8Vector * XLALCreateCOMPLEX8Vector(UINT4 length)
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLALDestroyCOMPLEX8Vector
void XLALDestroyCOMPLEX8Vector(COMPLEX8Vector *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_EFREQ
XLAL_EFREQ
XLAL_ENOMEM
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFAULT
XLAL_EFUNC
XLAL_EFUNC
XLAL_EDOM
XLAL_EDOM
XLAL_ESIZE
XLAL_ESIZE
XLAL_EFAILED
XLAL_EFAILED
XLAL_FAILURE
XLAL_FAILURE
XLAL_ETIME
XLAL_ETIME
thinca_brute_force_coinc.dt
dt
Definition: thinca_brute_force_coinc.py:52
alpha
double alpha
t0
double t0
COMPLEX8Vector::data
COMPLEX8 * data
CoherentGW::position
SkyPosition position
CoherentGW::a
REAL4TimeVectorSeries * a
CoherentGW::phi
REAL8TimeSeries * phi
CoherentGW::h
REAL4TimeVectorSeries * h
CoherentGW::f
REAL4TimeSeries * f
InspiralInit::func
expnFunc func
Definition: LALInspiral.h:680
InspiralInit::ak
expnCoeffs ak
Definition: LALInspiral.h:679
InspiralTemplate
The inspiral waveform parameter structure containing information about the waveform to be generated.
Definition: LALInspiral.h:205
LALPSpinInspiralRDparams::epnspin20S1S2
REAL8 epnspin20S1S2
Coeff. of the term in energy.
Definition: LALPSpinInspiralRD.c:195
LALPSpinInspiralRDparams::epnspin20S2S2
REAL8 epnspin20S2S2
Coeff. of the term in energy.
Definition: LALPSpinInspiralRD.c:199
LALPSpinInspiralRDparams::wdotspin25S2LNh
REAL8 wdotspin25S2LNh
Coeff. of the cntrb. to .
Definition: LALPSpinInspiralRD.c:182
LALPSpinInspiralRDparams::wdotorblog
REAL8 wdotorblog
Log coefficient of the PN expansion of of .
Definition: LALPSpinInspiralRD.c:175
LALPSpinInspiralRDparams::epnspin20S1S1
REAL8 epnspin20S1S1
Coeff. of the term in energy.
Definition: LALPSpinInspiralRD.c:197
LALPSpinInspiralRDparams::epnspin15S2dotLNh
REAL8 epnspin15S2dotLNh
Coeff. of the term in energy.
Definition: LALPSpinInspiralRD.c:194
LALPSpinInspiralRDparams::epnspin15S1dotLNh
REAL8 epnspin15S1dotLNh
Coeff. of the term in energy.
Definition: LALPSpinInspiralRD.c:193
LALPSpinInspiralRDparams::epnspin20S1S2dotLNh
REAL8 epnspin20S1S2dotLNh
Coeff. of the term in energy.
Definition: LALPSpinInspiralRD.c:196
LALPSpinInspiralRDparams::epnorb
REAL8 epnorb[4]
Coefficients of the PN expansion of the energy.
Definition: LALPSpinInspiralRD.c:192
LALPSpinInspiralRDparams::eta
REAL8 eta
symmetric mass ratio
Definition: LALPSpinInspiralRD.c:165
LALPSpinInspiralRDparams::wdotorb
REAL8 wdotorb[8]
Coefficients of the analytic PN expansion of .
Definition: LALPSpinInspiralRD.c:174
LALPSpinInspiralRDparams::wdotspin20S1S1
REAL8 wdotspin20S1S1
Coeff. of the cntrb. to .
Definition: LALPSpinInspiralRD.c:179
PPNParamStruc
This structure stores the parameters for constructing a restricted post-Newtonian waveform.
Definition: GeneratePPNInspiral.h:230
PPNParamStruc::termDescription
const CHAR * termDescription
The termination code description (above)
Definition: GeneratePPNInspiral.h:274
PPNParamStruc::termCode
INT4 termCode
The termination condition (above) that stopped computation of the waveform.
Definition: GeneratePPNInspiral.h:273
PPNParamStruc::dfdt
REAL4 dfdt
The maximum value of encountered over any timestep used in generating the waveform.
Definition: GeneratePPNInspiral.h:269
PPNParamStruc::psi
REAL4 psi
polarization angle (radians)
Definition: GeneratePPNInspiral.h:234
PPNParamStruc::fStart
REAL4 fStart
The actual starting frequency of the waveform, in Hz (normally close but not identical to fStartIn)
Definition: GeneratePPNInspiral.h:270
PPNParamStruc::position
SkyPosition position
location of source on sky
Definition: GeneratePPNInspiral.h:233
PPNParamStruc::length
UINT4 length
The length of the generated waveform.
Definition: GeneratePPNInspiral.h:272
PPNParamStruc::fStop
REAL4 fStop
The frequency at the termination of the waveform, in Hz.
Definition: GeneratePPNInspiral.h:271
PPNParamStruc::tc
REAL8 tc
The time from the start of the waveform to coalescence (in the point-mass approximation),...
Definition: GeneratePPNInspiral.h:268
PPNParamStruc::fStartIn
REAL4 fStartIn
The requested starting frequency of the waveform, in Hz.
Definition: GeneratePPNInspiral.h:252
PPNParamStruc::deltaT
REAL8 deltaT
The requested sampling interval of the waveform, in s.
Definition: GeneratePPNInspiral.h:251
REAL4TimeSeries::name
CHAR name[LALNameLength]
REAL4TimeSeries::data
REAL4Sequence * data
REAL4TimeSeries::sampleUnits
LALUnit sampleUnits
REAL4TimeVectorSeries::name
CHAR name[LALNameLength]
REAL4TimeVectorSeries::data
REAL4VectorSequence * data
REAL4Vector::data
REAL4 * data
REAL8Array::data
REAL8 * data
REAL8TimeSeries::data
REAL8Sequence * data
REAL8TimeSeries::sampleUnits
LALUnit sampleUnits
REAL8TimeSeries::name
CHAR name[LALNameLength]
REAL8Vector::data
REAL8 * data
expnCoeffs::ETaN
REAL8 ETaN
Definition: LALInspiral.h:414
expnCoeffs::ETa1
REAL8 ETa1
Definition: LALInspiral.h:414
expnCoeffs::ST
REAL8 ST[9]
Definition: LALInspiral.h:464
expnCoeffs::ETa3
REAL8 ETa3
Definition: LALInspiral.h:414
expnCoeffs::ETa2
REAL8 ETa2
Definition: LALInspiral.h:414
params::inclination
double inclination
params::distance
double distance
rk4GSLIntegrator
Structure containing steps and controls for the GSL Runge-Kutta solver.
Definition: LALInspiral.h:637
rk4In
Structure used as an input to Runge-Kutta solver.
Definition: LALInspiral.h:620
rk4In::dym
REAL8Vector * dym
Definition: LALInspiral.h:626
rk4In::x
REAL8 x
Definition: LALInspiral.h:622
rk4In::yt
REAL8Vector * yt
Definition: LALInspiral.h:625
rk4In::dydx
REAL8Vector * dydx
Definition: LALInspiral.h:624
rk4In::h
REAL8 h
Definition: LALInspiral.h:628
rk4In::dyt
REAL8Vector * dyt
Definition: LALInspiral.h:627
rk4In::n
INT4 n
Definition: LALInspiral.h:629
rk4In::y
REAL8Vector * y
Definition: LALInspiral.h:623
rk4In::function
TestFunction * function
Definition: LALInspiral.h:621