Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInspiral 5.0.3.1-b246709
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALEOBWaveform.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2007 Stas Babak, David Churches, Duncan Brown, David Chin, Jolien Creighton,
3* B.S. Sathyaprakash, Craig Robinson , Thomas Cokelaer, Evan Ochsner
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 * \author Sathyaprakash, B. S., Cokelaer T.
23 * \file
24 * \ingroup LALInspiral_h
25 *
26 * \brief Module to generate effective-one-body waveforms.
27 *
28 * ### Prototypes ###
29 *
30 * <tt>LALEOBWaveform()</tt>
31 * <ul>
32 * <li> \c signalvec: Output containing the inspiral waveform.</li>
33 * <li> \c params: Input containing binary chirp parameters.</li>
34 * </ul>
35 *
36 * <tt>LALEOBWaveformTemplates()</tt>
37 * <ul>
38 * <li> \c signalvec1: Output containing the 0-phase inspiral waveform.</li>
39 * <li> \c signalvec2: Output containing the \f$\pi/2\f$-phase inspiral waveform.</li>
40 * <li> \c params: Input containing binary chirp parameters.</li>
41 * </ul>
42 *
43 * <tt>LALEOBWaveformForInjection()</tt>
44 * <ul>
45 * <li> \c inject_hc: Output containing the 0-phase inspiral waveform.</li>
46 * <li> \c inject_hp: Output containing the \f$\pi/2\f$-phase inspiral waveform.</li>
47 * <li> \c inject_phase: Output containing the phase of inspiral waveform.</li>
48 * <li> \c inject_freq: Output containing the frequency of inspiral waveform.</li>
49 * <li> \c params: Input containing binary chirp parameters.</li>
50 * </ul>
51 *
52 * ### Description ###
53 *
54 * By solving four coupled ordinary differential equations in
55 * \eqref{eq_3_28}--\eqref{eq_3_31} this module computes the
56 * waveform in \eqref{eq_4_1} (see discussion in sec_EOB
57 * for details on how the initial conditions are chosen, when the
58 * waveform is terminated and so on).
59 * No quasi-normal mode oscillations are added to the plunge signal
60 * so the waveform is terminated around \f$2.8\,M\f$.
61 *
62 * ### 3PN vs 2PN ###
63 *
64 * At 3PN, two additional parameters exist namely OmegaS and Zeta2.
65 * The first parameters should be set to zero. If the second parameter
66 * is also set to zero then the waveform correponds to the standard
67 * waveforms.
68 *
69 * ### Algorithm ###
70 *
71 * A fourth order Runge-Kutta is used to solve the differential equations.
72 *
73 * ### Uses ###
74 *
75 * \code
76 * LALInspiralSetup
77 * LALInspiralChooseModel
78 * LALInspiralVelocity
79 * LALInspiralPhasing1
80 * LALDBisectionFindRoot
81 * LALRungeKutta4
82 * LALHCapDerivatives
83 * LALHCapDerivatives3PN
84 * LALHCapDerivativesP4PN
85 * LALlightRingRadius
86 * LALlightRingRadius3PN
87 * LALlightRingRadiusP4PN
88 * LALpphiInit
89 * LALpphiInit3PN
90 * LALpphiInitP4PN
91 * LALprInit
92 * LALprInit3PN
93 * LALprInitP4PN
94 * LALrOfOmega
95 * LALrOfOmega3PN
96 * LALrOfOmegaP4PN
97 * \endcode
98 *
99 * ### Notes ###
100 *
101 * The length of the waveform returned by \c LALInspiralWaveLength is
102 * occassionally smaller than what is required to hold an EOB waveform.
103 * This is because EOB goes beyond the last stable orbit up to the light
104 * ring while \c LALInspiralWaveLength assumes that the waveform terminates
105 * at the last stable orbit. It is recommended that a rather generous
106 * <tt>params->nEndPad</tt> be used to prevent the code from crashing.
107 *
108 */
109#include <lal/Units.h>
110#include <lal/LALInspiral.h>
111#include <lal/FindRoot.h>
112#include <lal/SeqFactories.h>
113#include <lal/NRWaveInject.h>
114#include <lal/TimeSeries.h>
115
116#ifdef __GNUC__
117#define UNUSED __attribute__ ((unused))
118#else
119#define UNUSED
120#endif
121
122#define ninty4by3etc 18.687902694437592603 /* (94/3 -41/31*pi*pi) */
123
124typedef struct tagrOfOmegaIn {
125 REAL8 eta, omega;
126} rOfOmegaIn;
127
128typedef struct tagPr3In {
129 REAL8 eta, zeta2, omegaS, omega, vr,r,q;
130 InspiralDerivativesIn in3copy;
131} pr3In;
132
133static REAL8
134XLALOmegaOfR2PN (
135 REAL8 r,
136 REAL8 eta) ;
137
138static REAL8
139XLALOmegaOfR3PN (
140 REAL8 r,
141 void *params) ;
142
143static REAL8
144XLALOmegaOfRP4PN (
145 REAL8 r,
146 void *params) ;
147
148static
149void LALHCapDerivatives( REAL8Vector *values,
150 REAL8Vector *dvalues,
151 void *funcParams);
152
153static
154REAL8 XLALprInit( REAL8 r,
155 InspiralDerivativesIn *ak);
156
157static
158REAL8 XLALpphiInit( REAL8 r,
159 REAL8 eta);
160
161static
162REAL8 XLALlightRingRadius( REAL8 r,
163 void *params);
164
165static REAL8 XLALrOfOmega( REAL8 r,
166 void *params);
167
168static
169void LALHCapDerivatives3PN( REAL8Vector *values,
170 REAL8Vector *dvalues,
171 void *funcParams);
172
173static
174REAL8 XLALprInit3PN( REAL8 , void *params);
175
176static
177REAL8 XLALpphiInit3PN( REAL8 r, REAL8 eta, REAL8 omegaS);
178
179static
180REAL8 XLALlightRingRadius3PN( REAL8 r, void *params);
181
182static
183REAL8 XLALrOfOmega3PN ( REAL8 r, void *params);
184
185static
186REAL8 XLALvr3PN(void *params);
187
188static
189void LALHCapDerivativesP4PN( REAL8Vector *values,
190 REAL8Vector *dvalues,
191 void *funcParams);
192
193static
194REAL8 XLALprInitP4PN( REAL8 , void *params);
195
196static
197REAL8 XLALpphiInitP4PN( REAL8 r, REAL8 eta, REAL8 omegaS);
198
199static
200REAL8 XLALlightRingRadiusP4PN( REAL8 r, void *params);
201
202static
203REAL8 XLALrOfOmegaP4PN ( REAL8 r, void *params);
204
205static
206REAL8 XLALvrP4PN(void *params);
207
208static int
209XLALEOBWaveformEngine (
210 REAL4Vector *signalvec1,
211 REAL4Vector *signalvec2,
212 REAL4Vector *h,
213 REAL4Vector *a,
214 REAL4Vector *ff,
215 REAL8Vector *phi,
216 UINT4 *countback,
217 InspiralTemplate *params,
218 InspiralInit *paramsInit
219 );
220
221/*--------------------------------------------------------------------*/
222
223static REAL8
224XLALprInit(
225 REAL8 r,
226 InspiralDerivativesIn *ak
227 )
228{
229
230/*
231 This combines Eq. (4.13), (4.14) (4.16) of BD2 to get
232 the initial value of pr
233*/
234
235 REAL8 eta,
236 z,
237 omega,
238 jadiab,
239 djadiabdr,
240 H0cap,
241 v,
242 FDIS,
243 A,
244 r2, /* temp variable */
245 r3, /* temp variable */
246 cr,
247 pr;
248
249 /* Temp variable to avoid the use of pow */
250 REAL8 tmpVar;
251
252 eta = ak->coeffs->eta;
253 r2 = r*r;
254 r3 = r2*r;
255
256 A = 1. - 2./r + 2.*eta/r3;
257 z = r3 * A * A/(r3 - 3.*r2 + 5.*eta);
258 omega = sqrt((1.-3*eta/r2)/(r3 * (1. + 2*eta*(sqrt(z)-1.))));
259 jadiab = sqrt (r2 * (r2 - 3.*eta)/(r3 - 3.*r2 + 5.*eta));
260
261 tmpVar = r3 - 3.*r2 + 5.*eta;
262 djadiabdr = (r3*r3 - 6.*r3*r2 + 3.*eta*r2*r2 +20.*eta*r3-30.*eta*eta*r)/
263 (tmpVar * tmpVar * 2. * jadiab);
264 H0cap = sqrt(1. + 2.*eta*(-1. + sqrt(z)))/eta;
265 cr = A*A/((1.-6.*eta/r2) * eta * H0cap * sqrt(z));
266 v = cbrt(omega);
267 FDIS = -ak->flux(v, ak->coeffs)/(eta * omega);
268 pr = FDIS/(djadiabdr*cr);
269 return pr;
270}
271
272/*--------------------------------------------------------------------*/
273static REAL8
274XLALpphiInit(
275 REAL8 r,
276 REAL8 eta
277 )
278{
279 REAL8 phase, r2, r3;
280 r2 = r*r;
281 r3 = r2*r;
282 phase = sqrt(r2 * (r2 - 3.*eta) / (r3 - 3.* r2 + 5.*eta));
283 return phase;
284}
285
286/*--------------------------------------------------------------------*/
287static REAL8
288XLALOmegaOfR2PN (
289 REAL8 r,
290 REAL8 eta
291 )
292{
293
294 REAL8 x, r2, r3, A, z;
295
296 r2 = r*r;
297 r3 = r2*r;
298
299 A = 1. - 2./r + 2.*eta/r3;
300 z = r3 * A * A/(r3 - 3.*r2 + 5.*eta);
301 x = sqrt((1.-3*eta/r2)/(r3 * (1. + 2*eta*(sqrt(z)-1.))));
302 return x;
303}
304
305/*--------------------------------------------------------------------*/
306
307
308
309static REAL8
310XLALrOfOmega (
311 REAL8 r,
312 void *params
313 )
314{
315
316 REAL8 x, r2, r3, A, z, eta, omega;
317 rOfOmegaIn *rofomegain;
318
319 rofomegain = (rOfOmegaIn *) params;
320 eta = rofomegain->eta;
321 omega = rofomegain->omega;
322 r2 = r*r;
323 r3 = r2*r;
324
325 A = 1. - 2./r + 2.*eta/r3;
326 z = r3 * A * A/(r3 - 3.*r2 + 5.*eta);
327 x = omega - sqrt((1.-3*eta/r2)/(r3 * (1. + 2*eta*(sqrt(z)-1.))));
328
329 return x;
330}
331
332/*--------------------------------------------------------------------*/
333
334static REAL8
335XLALlightRingRadius(
336 REAL8 r,
337 void *params
338 )
339{
340 REAL8 x, eta, r2;
341 rOfOmegaIn *rofomegain;
342
343 rofomegain = (rOfOmegaIn *) params;
344 eta = rofomegain->eta;
345 r2 = r*r;
346 x = r2*r - 3.*r2 + 5.* eta;
347 return x;
348}
349
350
351static void
352LALHCapDerivatives(
353 REAL8Vector *values,
354 REAL8Vector *dvalues,
355 void *funcParams
356 )
357{
358 REAL8 r, p, q, r2, r3, p2, q2, A, B, dA, dB, hcap, Hcap, etahH;
359 REAL8 omega, v, eta;
360 InspiralDerivativesIn *ak;
361
362 ak = (InspiralDerivativesIn *) funcParams;
363
364 eta = ak->coeffs->eta;
365
366 r = values->data[0];
367 p = values->data[2];
368 q = values->data[3];
369
370 r2 = r*r;
371 r3 = r2*r;
372 p2 = p*p;
373 q2 = q*q;
374 A = 1. - 2./r + 2.*eta/r3;
375 B = (1. - 6.*eta/r2)/A;
376 dA = 2./r2 - 6.*eta/(r3*r);
377 dB = (-dA * B + 12.*eta/r3)/A;
378 hcap = sqrt(A*(1. + p2/B + q2/r2));
379 Hcap = sqrt(1. + 2.*eta*(hcap - 1.)) / eta;
380 etahH = eta*hcap*Hcap;
381
382 dvalues->data[0] = A*A*p/((1. - 6.*eta/r2) * etahH);
383 dvalues->data[1] = omega = A * q / (r2 * etahH);
384
385 v = cbrt(omega);
386
387 dvalues->data[2] = -0.5 * (dA * hcap * hcap/A - p2 * A * dB/(B*B)
388 - 2. * A * q2/r3) / etahH;
389 dvalues->data[3] = -ak->flux(v, ak->coeffs)/(eta * omega);
390 /*
391 printf("%e %e %e %e %e %e %e %e\n", r, s, p, q, A, B, hcap, Hcap);
392 */
393}
394
395
396
397/*--------------------------------------------------------------------*/
398
399static REAL8
400XLALpphiInit3PN(
401 REAL8 r,
402 REAL8 eta,
403 REAL8 omegaS
404 )
405{
406 REAL8 phase, u, u2, u3, a4, a4p4eta, a4peta2, NA, DA, A, dA;
407
408
409 u = 1./r;
410 u2 = u*u;
411 u3 = u2*u;
412 a4 = (ninty4by3etc - 2. * omegaS) * eta;
413 a4p4eta = a4 + 4. * eta;
414 a4peta2 = a4 + eta * eta;
415 NA = 2.*(4.-eta) + (a4 - 16. + 8. * eta) * u;
416 DA = 2.*(4.-eta) + a4p4eta * u + 2. * a4p4eta * u2 + 4.*a4peta2 * u3;
417 A = NA/DA;
418 dA = ( (a4 - 16. + 8. * eta) * DA - NA * (a4p4eta + 4. * a4p4eta * u
419 + 12. * a4peta2 * u2))/(DA*DA);
420
421 phase = sqrt(-dA/(2.*u*A + u2 * dA));
422
423 return phase;
424}
425/*---------------------------------------------------------------*/
426
427static REAL8
428XLALprInit3PN(
429 REAL8 p,
430 void *params
431 )
432{
433 REAL8 pr, u, u2, u3, p2, p3, p4, A, DA, NA;
434 REAL8 onebyD, AbyD, Heff, HReal, etahH;
435 REAL8 omegaS, eta, a4, a4p4eta, a4peta2, z3, r, vr, q;
436 pr3In *ak;
437
438 ak = (pr3In *) params;
439
440 eta = ak->eta;
441 vr = ak->vr;
442 r = ak->r;
443 q = ak->q;
444 omegaS = ak->omegaS;
445
446
447 p2 = p*p;
448 p3 = p2*p;
449 p4 = p2*p2;
450 u = 1./ r;
451 u2 = u*u;
452 u3 = u2 * u;
453 z3 = 2. * (4. - 3. * eta) * eta;
454 a4 = (ninty4by3etc - 2. * omegaS) * eta;
455 a4p4eta = a4 + 4. * eta;
456 a4peta2 = a4 + eta * eta;
457
458/* From DJS 14, 15 */
459 NA = 2.*(4.-eta) + (a4 - 16. + 8. * eta) * u;
460 DA = 2.*(4.-eta) + a4p4eta * u + 2. * a4p4eta * u2 + 4.*a4peta2 * u3;
461 A = NA/DA;
462 onebyD = 1. + 6.*eta*u2 + 2. * ( 26. - 3. * eta) * eta * u3;
463 AbyD = A * onebyD;
464
465 Heff = sqrt(A*(1. + AbyD * p2 + q*q * u2 + z3 * p4 * u2));
466 HReal = sqrt(1. + 2.*eta*(Heff - 1.)) / eta;
467 etahH = eta*Heff*HReal;
468
469 pr = -vr + A*(AbyD*p + 2. * z3 * u2 * p3)/etahH;
470
471 return pr;
472}
473
474static REAL8
475XLALOmegaOfR3PN (
476 REAL8 r,
477 void *params)
478{
479 REAL8 x, u, u2, u3, a4, a4p4eta, a4peta2, eta, NA, DA, A, dA;
480 REAL8 omegaS;
481
482 /*include a status here ?*/
483 pr3In *ak;
484 ak = (pr3In *) params;
485 omegaS = ak->omegaS;
486 eta = ak->eta;
487
488 u = 1./r;
489 u2 = u*u;
490 u3 = u2*u;
491 a4 = (ninty4by3etc - 2. * omegaS) * eta;
492
493 a4p4eta = a4 + 4. * eta;
494 a4peta2 = a4 + eta * eta;
495 NA = 2.*(4.-eta) + (a4 - 16. + 8. * eta) * u;
496 DA = 2.*(4.-eta) + a4p4eta * u + 2. * a4p4eta * u2 + 4.*a4peta2 * u3;
497 A = NA/DA;
498 dA = ( (a4 - 16. + 8. * eta) * DA - NA
499 * (a4p4eta + 4. * a4p4eta * u + 12. * a4peta2 * u2))/(DA*DA);
500 x = sqrt ( -u3 * 0.5 * dA /(1. + 2.*eta * (A/sqrt(A+0.5 * u*dA)-1.)));
501 return x;
502
503}
504
505static REAL8
506XLALrOfOmega3PN(
507 REAL8 r,
508 void *params)
509{
510 REAL8 x, omega1,omega2;
511 pr3In *pr3in;
512
513 pr3in = (pr3In *) params;
514
515 omega1 = pr3in->omega;
516 omega2 = XLALOmegaOfR3PN(r, params);
517 x = -omega1 + omega2;
518
519 return x;
520}
521/*--------------------------------------------------------------------*/
522static REAL8
523XLALlightRingRadius3PN(
524 REAL8 r,
525 void *params
526 )
527{
528 REAL8 x, eta, u, u2, u3, a4, a4p4eta,a4peta2, NA, DA, A, dA;
529 rOfOmegaIn *rofomegain;
530 rofomegain = (rOfOmegaIn *) params;
531 eta = rofomegain->eta;
532
533
534 u = 1./r;
535 u2 = u*u;
536 u3 = u2*u;
537 a4 = ninty4by3etc * eta;
538 a4p4eta = a4 + 4. * eta;
539 a4peta2 = a4 + eta * eta;
540 NA = 2.*(4.-eta) + (a4 - 16. + 8. * eta) * u;
541 DA = 2.*(4.-eta) + a4p4eta * u + 2. * a4p4eta * u2 + 4.*a4peta2 * u3;
542 A = NA/DA;
543 dA = ( (a4 - 16. + 8. * eta) * DA - NA * (a4p4eta + 4.
544 * a4p4eta * u + 12. * a4peta2 * u2))/(DA*DA);
545 x = 2 * A + dA * u;
546 return x;
547}
548/*--------------------------------------------------------------------*/
549
550static void
551LALHCapDerivatives3PN(
552 REAL8Vector *values,
553 REAL8Vector *dvalues,
554 void *funcParams
555 )
556{
557 REAL8 r, p, q, u, u2, u3, u4, p2, p3, p4, q2, Apot, DA, NA;
558 REAL8 dA, onebyD, DonebyD, AbyD, Heff, HReal, etahH;
559 REAL8 omega, v, eta, a4, a4p4eta, a4peta2, z2, z30, zeta2;
560 REAL8 n1, c1, d1, d2, d3, oneby4meta;
561 REAL8 flexNonAdiab = 0;
562 REAL8 flexNonCirc = 0;
563
564 InspiralDerivativesIn *ak;
565
566 ak = (InspiralDerivativesIn *) funcParams;
567 eta = ak->coeffs->eta;
568 zeta2 = ak->coeffs->zeta2;
569
570 r = values->data[0];
571 p = values->data[2];
572 q = values->data[3];
573
574 p2 = p*p;
575 p3 = p2*p;
576 p4 = p2*p2;
577 q2 = q*q;
578 u = 1./r;
579 u2 = u*u;
580 u3 = u2 * u;
581 u4 = u2 * u2;
582 z30 = 2.L * (4.L - 3.L * eta) * eta;
583 z2 = 0.75L * z30 * zeta2,
584
585 a4 = ninty4by3etc * eta;
586 a4p4eta = a4 + 4. * eta;
587 a4peta2 = a4 + eta * eta;
588
589 /* From DJS 3PN Hamiltonian */
590 oneby4meta = 1./(4.-eta);
591 n1 = 0.5 * (a4 - 16. + 8.*eta) * oneby4meta;
592 d1 = 0.5 * a4p4eta * oneby4meta;
593 d2 = a4p4eta * oneby4meta;
594 d3 = 2. * a4peta2 * oneby4meta;
595 NA = 1. + n1 * u;
596 DA = 1 + d1*u + d2*u2 + d3*u3;
597 Apot = NA/DA;
598
599 onebyD = 1. + 6.*eta*u2 + (2. * ( 26. - 3. * eta) * eta - z2)* u3;
600 AbyD = Apot * onebyD;
601 Heff = sqrt(Apot*(1. + AbyD * p2 + q*q * u2 + z30 * (p4 + zeta2*(-0.25*p4
602 + 0.75 * p2 * q2 * u2 )) * u2));
603 HReal = sqrt(1. + 2.*eta*(Heff - 1.)) / eta;
604
605 dA = -u2/(DA*DA) * (n1*DA - NA * (d1 + 2.*d2*u + 3.*d3*u2));
606
607 DonebyD = -12.*eta*u3 - (6.*(26. - 3.*eta)*eta - z2)*u4;
608 etahH = eta*Heff*HReal;
609
610 dvalues->data[0] = Apot*(AbyD*p + z30 * u2 *(2.* p3
611 + zeta2*(-0.5*p3 + 0.75*p*q2*u2)))/etahH;
612 dvalues->data[1] = omega = Apot * q * u2 * (1. + 0.75*z30*zeta2*p2*u2)/ etahH;
613 v = cbrt(omega);
614
615 dvalues->data[2] = -0.5 * Apot * (dA*Heff*Heff/(Apot*Apot) - 2.*q2*u3
616 + (dA * onebyD + Apot * DonebyD) * p2
617 + z30 * u3 *(-2.* p4+zeta2*(0.5*p4 - 3.0*p2*q2*u2))) / etahH;
618
619 c1 = 1.+(u2 - 2.*u3*Apot/dA) * q2;
620 dvalues->data[3] = -(1. - flexNonAdiab*c1) * (1. + flexNonCirc*p2/(q2*u2))
621 * ak->flux(v,ak->coeffs)/(eta * omega);
622}
623
624
625
626/*----------------------------------------------------------------------*/
627static REAL8 XLALvr3PN( void *params )
628{
629 REAL8 vr, A, dA, d2A, NA, DA, dDA, dNA, d2DA, DA2;
630 REAL8 u, u2, u3, v, x1;
631 REAL8 eta,a4, a4p4eta, a4peta2, FDIS;
632
633 /* Temp variable to avoid use of pow */
634 REAL8 x1Bit;
635
636 pr3In *pr3in;
637 pr3in = (pr3In *)params;
638
639 eta = pr3in->eta;
640 u = 1./ pr3in->r;
641
642 u2 = u*u;
643 u3 = u2*u;
644
645
646 a4 = (ninty4by3etc - 2. * pr3in->omegaS) * eta;
647 a4p4eta = a4 + 4. * eta;
648 a4peta2 = a4 + eta * eta;
649 NA = 2.*(4.-eta) + (a4 - 16. + 8. * eta) * u;
650 DA = 2.*(4.-eta) + a4p4eta * u + 2. * a4p4eta * u2 + 4.*a4peta2 * u3;
651 DA2 = DA * DA;
652 A = NA/DA;
653 dNA = (a4 - 16. + 8. * eta);
654 dDA = (a4p4eta + 4. * a4p4eta * u + 12. * a4peta2 * u2);
655 d2DA = 4. * a4p4eta + 24. * a4peta2 * u;
656
657 dA = (dNA * DA - NA * dDA)/ DA2;
658 d2A = (-NA * DA * d2DA - 2. * dNA * DA * dDA + 2. * NA * dDA * dDA)/(DA2*DA);
659 v = cbrt(pr3in->omega);
660 FDIS = -pr3in->in3copy.flux(v, pr3in->in3copy.coeffs)/(eta* pr3in->omega);
661
662 x1Bit = 2.* u * A + u2 * dA;
663 x1 = -1./u2 * sqrt (-dA * x1Bit*x1Bit*x1Bit )
664 / (2.* u * dA * dA + A*dA - u * A * d2A);
665 vr = FDIS * x1;
666 return vr;
667}
668
669/*-------------------------------------------------------------------*/
670/* pseudo-4PN functions */
671/*-------------------------------------------------------------------*/
672
673static REAL8
674XLALpphiInitP4PN(
675 REAL8 r,
676 REAL8 eta,
677 REAL8 omegaS
678 )
679{
680 REAL8 phase, u, u2, u3, u4, a4, a5, eta2, NA, DA, A, dA;
681
682 eta2 = eta*eta;
683 u = 1./r;
684 u2 = u*u;
685 u3 = u2*u;
686 u4 = u2*u2;
687 a4 = (ninty4by3etc - 2. * omegaS) * eta;
688 a5 = 60.;
689 NA = (32. - 24.*eta - 4.*a4 - a5*eta)*u + (a4 - 16. + 8.*eta);
690 DA = a4 - 16. + 8.*eta - (2.*a4 + a5*eta + 8.*eta)*u - (4.*a4 + 2.*a5*eta + 16.*eta)*u2
691 - (8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u3
692 + (-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u4;
693 A = NA/DA;
694 dA = ( (32. - 24.*eta - 4.*a4 - a5*eta) * DA - NA *
695 ( -(2.*a4 + a5*eta + 8.*eta) - 2.*(4.*a4 + 2.*a5*eta + 16.*eta)*u
696 - 3.*(8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u2
697 + 4.*(-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u3))/(DA*DA);
698 phase = sqrt(-dA/(2.*u*A + u2 * dA));
699/* why is it called phase? This is initial j!? */
700 return phase;
701}
702
703/*-------------------------------------------------------------------*/
704static REAL8
705XLALprInitP4PN(
706 REAL8 p,
707 void *params
708 )
709{
710 REAL8 pr, u, u2, u3, u4, p2, p3, p4, A, DA, NA;
711 REAL8 onebyD, AbyD, Heff, HReal, etahH;
712 REAL8 eta, eta2, a4, a5, z3, r, vr, q;
713 pr3In *ak;
714
715 ak = (pr3In *) params;
716
717 eta = ak->eta;
718 vr = ak->vr;
719 r = ak->r;
720 q = ak->q;
721 eta2 = eta*eta;
722
723
724 p2 = p*p;
725 p3 = p2*p;
726 p4 = p2*p2;
727 u = 1./ r;
728 u2 = u*u;
729 u3 = u2 * u;
730 u4 = u2 * u2;
731 z3 = 2. * (4. - 3. * eta) * eta;
732 a4 = ninty4by3etc * eta;
733 a5 = 60.;
734
735 NA = (32. - 24.*eta - 4.*a4 - a5*eta)*u + (a4 - 16. + 8.*eta);
736 DA = a4 - 16. + 8.*eta - (2.*a4 + a5*eta + 8.*eta)*u - (4.*a4 + 2.*a5*eta + 16.*eta)*u2
737 - (8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u3
738 + (-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u4;
739 A = NA/DA;
740 onebyD = 1. + 6.*eta*u2 + 2. * ( 26. - 3. * eta) * eta * u3 + 36.*eta2*u4;
741 AbyD = A * onebyD;
742
743 Heff = sqrt(A*(1. + AbyD * p2 + q*q * u2 + z3 * p4 * u2));
744 HReal = sqrt(1. + 2.*eta*(Heff - 1.)) / eta;
745 etahH = eta*Heff*HReal;
746
747 pr = -vr + A*(AbyD*p + 2. * z3 * u2 * p3)/etahH;
748/* This sets pr = dH/dpr - vr, calls rootfinder,
749 gets value of pr s.t. dH/pr = vr */
750 return pr;
751}
752
753
754/*-------------------------------------------------------------------*/
755static REAL8
756XLALOmegaOfRP4PN (
757 REAL8 r,
758 void *params)
759{
760 REAL8 x, u, u2, u3, u4, a4, a5, eta, eta2, NA, DA, A, dA;
761
762 /*include a status here ?*/
763 pr3In *ak;
764 ak = (pr3In *) params;
765 eta = ak->eta;
766 eta2 = eta*eta;
767
768 u = 1./r;
769 u2 = u*u;
770 u3 = u2*u;
771 u4 = u2*u2;
772 a4 = ninty4by3etc * eta;
773 a5 = 60.;
774 NA = (32. - 24.*eta - 4.*a4 - a5*eta)*u + (a4 - 16. + 8.*eta);
775 DA = a4 - 16. + 8.*eta - (2.*a4 + a5*eta + 8.*eta)*u - (4.*a4 + 2.*a5*eta + 16.*eta)*u2
776 - (8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u3
777 + (-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u4;
778 A = NA/DA;
779 dA = ( (32. - 24.*eta - 4.*a4 - a5*eta) * DA - NA *
780 ( -(2.*a4 + a5*eta + 8.*eta) - 2.*(4.*a4 + 2.*a5*eta + 16.*eta)*u
781 - 3.*(8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u2
782 + 4.*(-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u3))/(DA*DA);
783
784 x = sqrt ( - u3 * 0.5 * dA /(1. + 2.*eta * (A/sqrt(A+0.5 * u*dA)-1.)));
785 return x;
786}
787
788
789/*-------------------------------------------------------------------*/
790static REAL8
791XLALrOfOmegaP4PN(
792 REAL8 r,
793 void *params)
794{
795 REAL8 x, omega1,omega2;
796 pr3In *pr3in;
797
798 pr3in = (pr3In *) params;
799
800 omega1 = pr3in->omega;
801 omega2 = XLALOmegaOfRP4PN(r, params);
802 x = -omega1 + omega2;
803
804 return x;
805}
806
807
808/*-------------------------------------------------------------------*/
809static REAL8
810XLALlightRingRadiusP4PN(
811 REAL8 r,
812 void *params
813 )
814{
815 REAL8 x, eta, eta2, u, u2, u3, u4, a4, a5, NA, DA, A, dA;
816 rOfOmegaIn *rofomegain;
817 rofomegain = (rOfOmegaIn *) params;
818 eta = rofomegain->eta;
819 eta2 = eta*eta;
820
821
822 u = 1./r;
823 u2 = u*u;
824 u3 = u2*u;
825 u4 = u2*u2;
826 a4 = ninty4by3etc * eta;
827 a5 = 60.;
828 NA = (32. - 24.*eta - 4.*a4 - a5*eta)*u + (a4 - 16. + 8.*eta);
829 DA = a4 - 16. + 8.*eta - (2.*a4 + a5*eta + 8.*eta)*u - (4.*a4 + 2.*a5*eta + 16.*eta)*u2
830 - (8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u3
831 + (-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u4;
832 A = NA/DA;
833 dA = ( (32. - 24.*eta - 4.*a4 - a5*eta) * DA - NA *
834 ( -(2.*a4 + a5*eta + 8.*eta) - 2.*(4.*a4 + 2.*a5*eta + 16.*eta)*u
835 - 3.*(8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u2
836 + 4.*(-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u3))/(DA*DA);
837 x = 2 * A + dA * u;
838
839 return x;
840}
841
842/*-------------------------------------------------------------------*/
843static void
844LALHCapDerivativesP4PN(
845 REAL8Vector *values,
846 REAL8Vector *dvalues,
847 void *funcParams
848 )
849{
850 REAL8 r, p, q, u, u2, u3, u4, u5, p2, p3, p4, q2, Apot, DA, NA;
851 REAL8 dA, onebyD, DonebyD, AbyD, Heff, HReal, etahH;
852 REAL8 omega, v, eta, eta2, a4, z2, z30, z3, zeta2;
853 REAL8 a5;
854 double UNUSED dr, UNUSED ds, UNUSED dp, UNUSED dq;
855
856 InspiralDerivativesIn *ak;
857
858 ak = (InspiralDerivativesIn *) funcParams;
859 eta = ak->coeffs->eta;
860 zeta2 = ak->coeffs->zeta2;
861
862 r = values->data[0];
863 p = values->data[2];
864 q = values->data[3];
865
866 p2 = p*p;
867 p3 = p2*p;
868 p4 = p2*p2;
869 q2 = q*q;
870 u = 1./r;
871 u2 = u*u;
872 u3 = u2 * u;
873 u4 = u2 * u2;
874 u5 = u*u4;
875 z30 = 2.L * (4.L - 3.L * eta) * eta;
876 z2 = 0.75L * z30 * zeta2,
877 z3 = z30 * (1.L - zeta2);
878 eta2 = eta*eta;
879
880 a4 = ninty4by3etc * eta;
881 a5 = 60.;
882
883 NA = (32. - 24.*eta - 4.*a4 - a5*eta)*u + (a4 - 16. + 8.*eta);
884 DA = a4 - 16. + 8.*eta - (2.*a4 + a5*eta + 8.*eta)*u - (4.*a4 + 2.*a5*eta + 16.*eta)*u2
885 - (8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u3
886 + (-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u4;
887 Apot = NA/DA; /* This A(u) assume zeta2=0 (the default value) */
888
889 onebyD = 1. + 6.*eta*u2 + (2.*eta * ( 26. - 3.*eta ) - z2)* u3 + 36.*eta2*u4;
890 AbyD = Apot * onebyD;
891 Heff = sqrt(Apot*(1. + AbyD * p2 + q*q * u2 + z3 * p4 * u2));
892 HReal = sqrt(1. + 2.*eta*(Heff - 1.)) / eta;
893 dA = -u2 * ( (32. - 24.*eta - 4.*a4 - a5*eta) * DA - NA *
894 ( -(2.*a4 + a5*eta + 8.*eta) - 2.*(4.*a4 + 2.*a5*eta + 16.*eta)*u
895 - 3.*(8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u2
896 + 4.*(-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u3))/(DA*DA);
897
898 DonebyD = -12.*eta*u3 - (6.*eta*(26. - 3.*eta) - z2)*u4 - 144.*eta2*u5;
899 etahH = eta*Heff*HReal;
900
901 dr = dvalues->data[0] = Apot*(AbyD*p + z30 * u2 *(2.* p3
902 + zeta2*(-0.5*p3 + 0.75*p*q2*u2)))/etahH;
903 ds = dvalues->data[1] = omega = Apot * q * u2 * (1. + 0.75*z30*zeta2*p2*u2)/ etahH;
904 v = cbrt(omega);
905
906 dp = dvalues->data[2] = -0.5 * Apot * (dA*Heff*Heff/(Apot*Apot) - 2.*q2*u3
907 + (dA * onebyD + Apot * DonebyD) * p2
908 + z30 * u3 *(-2.* p4+zeta2*(0.5*p4 - 3.0*p2*q2*u2))) / etahH;
909 dq = dvalues->data[3] = - ak->flux(v,ak->coeffs)/(eta * omega);
910 /*
911 fprintf(stdout, "%e %e %e %e %e %e %e %e %e %e %e\n", r, s, p, q, Heff, v, Apot, dr, ds, dp, dq);
912 */
913}
914
915
916/*-------------------------------------------------------------------*/
917static REAL8 XLALvrP4PN( void *params )
918{
919 REAL8 vr, A, dA, d2A, NA, DA, DA2, dDA, dNA, d2DA;
920 REAL8 u, u2, u3, u4, v, x1;
921 REAL8 eta, eta2, a4, a5, FDIS;
922
923 /* Temporary variable to avoid the use of pow */
924 REAL8 x1Bit;
925
926 pr3In *pr3in;
927 pr3in = (pr3In *)params;
928
929 eta = pr3in->eta;
930 u = 1./ pr3in->r;
931
932 u2 = u*u;
933 u3 = u2*u;
934 u4 = u2*u2;
935 eta2 = eta*eta;
936
937
938 a4 = (ninty4by3etc - 2. * pr3in->omegaS) * eta;
939 a5 = 60.;
940 NA = (32. - 24.*eta - 4.*a4 - a5*eta)*u + (a4 - 16. + 8.*eta);
941 DA = a4 - 16. + 8.*eta - (2.*a4 + a5*eta + 8.*eta)*u - (4.*a4 + 2.*a5*eta + 16.*eta)*u2
942 - (8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u3
943 + (-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u4;
944 DA2 = DA * DA;
945 A = NA/DA;
946 dNA = (32. - 24.*eta - 4.*a4 - a5*eta);
947 dDA = - (2.*a4 + a5*eta + 8.*eta) - 2.*(4.*a4 + 2.*a5*eta + 16.*eta)*u
948 - 3.*(8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u2
949 + 4.*(-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u3;
950 d2DA = - 2.*(4.*a4 + 2.*a5*eta + 16.*eta)
951 - 6.*(8.*a4 + 4.*a5*eta + 2.*a4*eta + 16.*eta2)*u
952 + 12.*(-a4*a4 - 8.*a5*eta - 8.*a4*eta + 2.*a5*eta2 - 16.*eta2)*u2;
953
954 dA = (dNA * DA - NA * dDA)/ DA2;
955 d2A = (-NA * DA * d2DA - 2. * dNA * DA * dDA + 2. * NA * dDA * dDA)/(DA2 * DA);
956
957 v = cbrt(pr3in->omega);
958 FDIS = -pr3in->in3copy.flux(v, pr3in->in3copy.coeffs)/(eta* pr3in->omega);
959
960 x1Bit = 2.* u * A + u2 * dA;
961 x1 = -1./u2 * sqrt (-dA * x1Bit * x1Bit * x1Bit )
962 / (2.* u * dA * dA + A*dA - u * A * d2A);
963 vr = FDIS * x1;
964 return vr;
965}
966
967
968/*-------------------------------------------------------------------*/
969
970
971int
972XLALEOBWaveform(
973 REAL4Vector *signalvec,
974 InspiralTemplate *params
975 )
976{
977 UINT4 count;
978 InspiralInit paramsInit;
979
980 if ( !signalvec )
981 {
982 XLAL_ERROR( XLAL_EFAULT );
983 }
984 if ( !signalvec->data )
985 {
986 XLAL_ERROR( XLAL_EFAULT );
987 }
988 if ( !params )
989 {
990 XLAL_ERROR( XLAL_EFAULT );
991 }
992 if ( params->nStartPad < 0 || params->nEndPad < 0 )
993 {
994 XLALPrintError( "nStartPad and nEndPad must be >= 0.\n");
995 XLAL_ERROR( XLAL_EDOM );
996 }
997 if ( params->fLower <= 0. )
998 {
999 XLALPrintError( "fLower must be > 0.\n");
1000 XLAL_ERROR( XLAL_EDOM );
1001 }
1002 if ( params->tSampling <= 0. )
1003 {
1004 XLALPrintError( "tSampling must be > 0.\n");
1005 XLAL_ERROR( XLAL_EDOM );
1006 }
1007 if ( params->totalMass <= 0. )
1008 {
1009 XLALPrintError( "totalMass must be > 0.\n");
1010 XLAL_ERROR( XLAL_EDOM );
1011 }
1012
1013 if ( XLALInspiralSetup (&(paramsInit.ak), params) == XLAL_FAILURE )
1014 {
1015 XLAL_ERROR( XLAL_EFUNC );
1016 }
1017
1018 if ( XLALInspiralChooseModel( &(paramsInit.func),
1019 &(paramsInit.ak), params) == XLAL_FAILURE )
1020 {
1021 XLAL_ERROR( XLAL_EFUNC );
1022 }
1023
1024 memset(signalvec->data, 0, signalvec->length * sizeof( REAL4 ));
1025
1026 /* Call the engine function */
1027 if ( XLALEOBWaveformEngine(signalvec, NULL, NULL, NULL,
1028 NULL, NULL, &count, params, &paramsInit) == XLAL_FAILURE )
1029 {
1030 XLAL_ERROR( XLAL_EFUNC );
1031 }
1032
1033 return XLAL_SUCCESS;
1034}
1035
1036int
1037XLALEOBWaveformTemplates(
1038 REAL4Vector *signalvec1,
1039 REAL4Vector *signalvec2,
1040 InspiralTemplate *params
1041 )
1042{
1043 UINT4 count;
1044
1045 InspiralInit paramsInit;
1046
1047 if ( !signalvec1 )
1048 {
1049 XLAL_ERROR( XLAL_EFAULT );
1050 }
1051 if ( !signalvec2 )
1052 {
1053 XLAL_ERROR( XLAL_EFAULT );
1054 }
1055 if ( !signalvec1->data )
1056 {
1057 XLAL_ERROR( XLAL_EFAULT );
1058 }
1059 if ( !signalvec2->data )
1060 {
1061 XLAL_ERROR( XLAL_EFAULT );
1062 }
1063 if ( !params )
1064 {
1065 XLAL_ERROR( XLAL_EFAULT );
1066 }
1067 if ( params->nStartPad < 0 || params->nEndPad < 0 )
1068 {
1069 XLALPrintError( "nStartPad and nEndPad must be >= 0.\n");
1070 XLAL_ERROR( XLAL_EDOM );
1071 }
1072 if ( params->fLower <= 0. )
1073 {
1074 XLALPrintError( "fLower must be > 0.\n");
1075 XLAL_ERROR( XLAL_EDOM );
1076 }
1077 if ( params->tSampling <= 0. )
1078 {
1079 XLALPrintError( "tSampling must be > 0.\n");
1080 XLAL_ERROR( XLAL_EDOM );
1081 }
1082 if ( params->totalMass <= 0. )
1083 {
1084 XLALPrintError( "totalMass must be > 0.\n");
1085 XLAL_ERROR( XLAL_EDOM );
1086 }
1087
1088 if ( XLALInspiralSetup (&(paramsInit.ak), params) == XLAL_FAILURE )
1089 {
1090 XLAL_ERROR( XLAL_EFUNC );
1091 }
1092 if ( XLALInspiralChooseModel(&(paramsInit.func),
1093 &(paramsInit.ak), params) == XLAL_FAILURE )
1094 {
1095 XLAL_ERROR( XLAL_EFUNC );
1096 }
1097
1098 memset(signalvec1->data, 0, signalvec1->length * sizeof( REAL4 ));
1099 memset(signalvec2->data, 0, signalvec2->length * sizeof( REAL4 ));
1100
1101 /* Call the engine function */
1102 if ( XLALEOBWaveformEngine(signalvec1, signalvec2, NULL, NULL,
1103 NULL, NULL, &count, params, &paramsInit) == XLAL_FAILURE )
1104 {
1105 XLAL_ERROR( XLAL_EFUNC );
1106 }
1107
1108 return XLAL_SUCCESS;
1109}
1110
1111
1112/*=========================================================*/
1113/*======INJECTION =========================================*/
1114/*=========================================================*/
1115
1116
1117int XLALEOBWaveformForInjection(
1118 CoherentGW *waveform,
1119 InspiralTemplate *params,
1120 PPNParamStruc *ppnParams
1121 )
1122{
1123 UINT4 count, i;
1124
1125 REAL4Vector *a=NULL;/* pointers to generated amplitude data */
1126 REAL4Vector *h=NULL;/* pointers to generated polarization data */
1127 REAL4Vector *ff=NULL ;/* pointers to generated frequency data */
1128 REAL8Vector *phi=NULL;/* pointer to generated phase data */
1129
1130 REAL8 s;
1131
1132 REAL8 phiC;/* phase at coalescence */
1133 InspiralInit paramsInit;
1134
1135 REAL8 deltaT;
1136
1137 /* We need a blank LIGOTimeGPS for creating time series */
1138 LIGOTimeGPS epoch = {0, 0};
1139
1140 /* Make sure parameter and waveform structures exist. */
1141 if ( !params )
1142 XLAL_ERROR( XLAL_EFAULT );
1143
1144 if ( !waveform )
1145 XLAL_ERROR( XLAL_EFAULT );
1146
1147 /* Make sure waveform fields don't exist. */
1148 if ( waveform->a )
1149 {
1150 XLALPrintError( "Pointer for waveform->a exists. Was expecting NULL.\n" );
1151 XLAL_ERROR( XLAL_EFAULT );
1152 }
1153 if ( waveform->h )
1154 {
1155 XLALPrintError( "Pointer for waveform->h exists. Was expecting NULL.\n" );
1156 XLAL_ERROR( XLAL_EFAULT );
1157 }
1158 if ( waveform->f )
1159 {
1160 XLALPrintError( "Pointer for waveform->f exists. Was expecting NULL.\n" );
1161 XLAL_ERROR( XLAL_EFAULT );
1162 }
1163 if ( waveform->phi )
1164 {
1165 XLALPrintError( "Pointer for waveform->phi exists. Was expecting NULL.\n" );
1166 XLAL_ERROR( XLAL_EFAULT );
1167 }
1168 if ( waveform->shift )
1169 {
1170 XLALPrintError( "Pointer for waveform->shift exists. Was expecting NULL.\n" );
1171 XLAL_ERROR( XLAL_EFAULT );
1172 }
1173
1174 params->ampOrder = (LALPNOrder) 0;
1175 XLALPrintInfo( "WARNING: Amp Order has been reset to %d\n", params->ampOrder);
1176
1177 /* Compute some parameters*/
1178 if ( XLALInspiralInit( params, &paramsInit ) == XLAL_FAILURE )
1179 {
1180 XLAL_ERROR( XLAL_EFUNC );
1181 }
1182
1183 if (paramsInit.nbins==0)
1184 {
1185 XLALPrintWarning( "Waveform is of zero length!\n" );
1186 return XLAL_SUCCESS;
1187 }
1188
1189 deltaT = 1. / params->tSampling;
1190
1191 /* Now we can allocate memory and vector for coherentGW structure*/
1192 ff = XLALCreateREAL4Vector(paramsInit.nbins);
1193 a = XLALCreateREAL4Vector(2*paramsInit.nbins);
1194 phi = XLALCreateREAL8Vector(paramsInit.nbins);
1195 if ( !ff || !a || !phi )
1196 {
1197 if (ff) XLALDestroyREAL4Vector( ff );
1198 if (a) XLALDestroyREAL4Vector( a );
1199 if (phi) XLALDestroyREAL8Vector( phi );
1200 XLAL_ERROR( XLAL_EFUNC );
1201 }
1202
1203 /* By default the waveform is empty */
1204 memset(ff->data, 0, paramsInit.nbins * sizeof(REAL4));
1205 memset(a->data, 0, 2 * paramsInit.nbins * sizeof(REAL4));
1206 memset(phi->data, 0, paramsInit.nbins * sizeof(REAL8));
1207
1208 if( params->approximant == EOBNR )
1209 {
1210 h = XLALCreateREAL4Vector(2*paramsInit.nbins);
1211 if ( !h )
1212 {
1213 XLALDestroyREAL4Vector( ff );
1214 XLALDestroyREAL4Vector( a );
1215 XLALDestroyREAL8Vector( phi );
1216 XLAL_ERROR( XLAL_EFUNC );
1217 }
1218 memset(h->data, 0, 2 * paramsInit.nbins * sizeof(REAL4));
1219 }
1220
1221 /* Call the engine function */
1222 params->startPhase = ppnParams->phi;
1223 if ( XLALEOBWaveformEngine(NULL, NULL, h, a, ff,
1224 phi, &count, params, &paramsInit) == XLAL_FAILURE )
1225 {
1226 XLALDestroyREAL4Vector(ff);
1227 XLALDestroyREAL4Vector(a);
1228 XLALDestroyREAL8Vector(phi);
1229 if( params->approximant == EOBNR )
1230 {
1231 XLALDestroyREAL4Vector(h);
1232 }
1233 XLAL_ERROR( XLAL_EFUNC );
1234 }
1235
1236 /* Check an empty waveform hasn't been returned */
1237 for (i = 0; i < phi->length; i++)
1238 {
1239 if (phi->data[i] != 0.0) break;
1240
1241 if (i == phi->length - 1)
1242 {
1243 XLALDestroyREAL4Vector(ff);
1244 XLALDestroyREAL4Vector(a);
1245 XLALDestroyREAL8Vector(phi);
1246 if( params->approximant == EOBNR )
1247 {
1248 XLALDestroyREAL4Vector(h);
1249 }
1250 XLALPrintWarning( "An empty waveform has been returned!\n" );
1251 return XLAL_SUCCESS;
1252 }
1253 }
1254
1255 s = 0.5 * phi->data[count - 1];
1256
1257 XLALPrintInfo("fFinal = %f\n", params->fFinal);
1258
1259 XLALPrintInfo("cycles = %f\n", s/3.14159);
1260
1261 XLALPrintInfo("final coalescence phase with respect to actual data =%f\n",
1262 (ff->data[count]-ff->data[count-1])/2/3.14159);
1263
1264
1265
1266 if ( (s/LAL_PI) < 2 )
1267 {
1268 XLALPrintWarning("The waveform has only %f cycles; we don't keep waveform with less than 2 cycles.\n",
1269 (double) s/ (double)LAL_PI );
1270 }
1271 else
1272 {
1273 phiC = phi->data[count-1] ;
1274
1275 for (i=0; i<count;i++)
1276 {
1277 phi->data[i] = -phiC + phi->data[i] + ppnParams->phi;
1278 }
1279
1280 /* Allocate the waveform structures. */
1281 if ( ( waveform->a = (REAL4TimeVectorSeries *)
1282 LALMalloc( sizeof(REAL4TimeVectorSeries) ) ) == NULL )
1283 {
1284 XLALDestroyREAL4Vector(ff);
1285 XLALDestroyREAL4Vector(a);
1286 XLALDestroyREAL8Vector(phi);
1287 if( params->approximant == EOBNR )
1288 {
1289 XLALDestroyREAL4Vector(h);
1290 }
1291 XLAL_ERROR( XLAL_ENOMEM );
1292 }
1293
1294 memset( waveform->a, 0, sizeof(REAL4TimeVectorSeries) );
1295
1296 waveform->a->data = XLALCreateREAL4VectorSequence( count, 2 );
1297 if ( !waveform->a->data )
1298 {
1299 XLALDestroyREAL4Vector(ff);
1300 XLALDestroyREAL4Vector(a);
1301 XLALDestroyREAL8Vector(phi);
1302 if( params->approximant == EOBNR )
1303 {
1304 XLALDestroyREAL4Vector(h);
1305 }
1306 XLAL_ERROR( XLAL_ENOMEM );
1307 }
1308
1309 waveform->a->deltaT = deltaT;
1310 waveform->a->sampleUnits = lalStrainUnit;
1311
1312 snprintf( waveform->a->name,
1313 LALNameLength, "EOB inspiral amplitudes");
1314
1315 waveform->f = XLALCreateREAL4TimeSeries( "EOB inspiral frequency",
1316 &epoch, 0, deltaT, &lalHertzUnit, count );
1317 if ( !(waveform->f) )
1318 {
1319 XLALDestroyREAL4VectorSequence( waveform->a->data );
1320 LALFree( waveform->a );
1321 waveform->a = NULL;
1322 XLALDestroyREAL4Vector(ff);
1323 XLALDestroyREAL4Vector(a);
1324 XLALDestroyREAL8Vector(phi);
1325 if( params->approximant == EOBNR )
1326 {
1327 XLALDestroyREAL4Vector(h);
1328 }
1329 XLAL_ERROR( XLAL_EFUNC );
1330 }
1331
1332 waveform->phi = XLALCreateREAL8TimeSeries("EOB inspiral phase",
1333 &epoch, 0, deltaT, &lalDimensionlessUnit, count );
1334 if ( !(waveform->phi) )
1335 {
1336 XLALDestroyREAL4TimeSeries( waveform->f );
1337 waveform->f = NULL;
1338 XLALDestroyREAL4VectorSequence( waveform->a->data );
1339 LALFree( waveform->a );
1340 waveform->a = NULL;
1341 XLALDestroyREAL4Vector(ff);
1342 XLALDestroyREAL4Vector(a);
1343 XLALDestroyREAL8Vector(phi);
1344 if( params->approximant == EOBNR )
1345 {
1346 XLALDestroyREAL4Vector(h);
1347 }
1348 XLAL_ERROR( XLAL_EFUNC );
1349 }
1350
1351 memcpy(waveform->a->data->data , a->data, 2*count*(sizeof(REAL4)));
1352 memcpy(waveform->f->data->data , ff->data, count*(sizeof(REAL4)));
1353 memcpy(waveform->phi->data->data ,phi->data, count*(sizeof(REAL8)));
1354
1355 waveform->position = ppnParams->position;
1356 waveform->psi = ppnParams->psi;
1357
1358 /* --- fill some output ---*/
1359 ppnParams->tc = (double)(count-1) / params->tSampling ;
1360 ppnParams->length = count;
1361 ppnParams->dfdt = ((REAL4)(waveform->f->data->data[count-1]
1362 - waveform->f->data->data[count-2])) * ppnParams->deltaT;
1363 ppnParams->fStop = params->fFinal;
1364 ppnParams->termCode = GENERATEPPNINSPIRALH_EFSTOP;
1365 ppnParams->termDescription = GENERATEPPNINSPIRALH_MSGEFSTOP;
1366
1367 ppnParams->fStart = ppnParams->fStartIn;
1368
1369 if( params->approximant == EOBNR )
1370 {
1371 if ( ( waveform->h = (REAL4TimeVectorSeries *)
1372 LALMalloc( sizeof(REAL4TimeVectorSeries) ) ) == NULL )
1373 {
1374 XLALDestroyREAL8TimeSeries( waveform->phi );
1375 waveform->phi = NULL;
1376 XLALDestroyREAL4TimeSeries( waveform->f );
1377 waveform->f = NULL;
1378 XLALDestroyREAL4VectorSequence( waveform->a->data );
1379 LALFree( waveform->a );
1380 waveform->a = NULL;
1381 XLALDestroyREAL4Vector(ff);
1382 XLALDestroyREAL4Vector(a);
1383 XLALDestroyREAL8Vector(phi);
1384 XLALDestroyREAL4Vector(h);
1385 XLAL_ERROR( XLAL_EFUNC );
1386 }
1387 memset( waveform->h, 0, sizeof(REAL4TimeVectorSeries) );
1388
1389 waveform->h->data = XLALCreateREAL4VectorSequence( count, 2 );
1390 if ( !(waveform->h->data) )
1391 {
1392 LALFree( waveform->h );
1393 waveform->h = NULL;
1394 XLALDestroyREAL8TimeSeries( waveform->phi );
1395 waveform->phi = NULL;
1396 XLALDestroyREAL4TimeSeries( waveform->f );
1397 waveform->f = NULL;
1398 XLALDestroyREAL4VectorSequence( waveform->a->data );
1399 LALFree( waveform->a );
1400 waveform->a = NULL;
1401 XLALDestroyREAL4Vector(ff);
1402 XLALDestroyREAL4Vector(a);
1403 XLALDestroyREAL8Vector(phi);
1404 XLALDestroyREAL4Vector(h);
1405 XLAL_ERROR( XLAL_EFUNC );
1406 }
1407
1408 memcpy(waveform->h->data->data , h->data, 2*count*(sizeof(REAL4)));
1409 waveform->h->deltaT = 1./params->tSampling;
1410 waveform->h->sampleUnits = lalStrainUnit;
1411 snprintf( waveform->h->name,
1412 LALNameLength, "EOB inspiral polarizations");
1413 XLALDestroyREAL4Vector(h);
1414 }
1415 } /* end phase condition*/
1416
1417 /* --- free memory --- */
1418
1419 XLALDestroyREAL4Vector( a );
1420 XLALDestroyREAL4Vector( ff );
1421 XLALDestroyREAL8Vector( phi );
1422
1423 /*on peut utiliser tSampling pour dfdt*/
1424
1425 return XLAL_SUCCESS;
1426
1427}
1428
1429/* Engine function for generating waveform
1430 Craig Robinson 15/07/05 */
1431static int
1432XLALEOBWaveformEngine (
1433 REAL4Vector *signalvec1,
1434 REAL4Vector *signalvec2,
1435 REAL4Vector *h,
1436 REAL4Vector *a,
1437 REAL4Vector *ff,
1438 REAL8Vector *phi,
1439 UINT4 *countback,
1440 InspiralTemplate *params,
1441 InspiralInit *paramsInit
1442 )
1443{
1444
1445
1446 UINT4 count, nn=4, length = 0, hiSRndx=0, ndx=0, higherSR=0;
1447 REAL4Vector *sig1, *sig2, *ampl, *freq;
1448 REAL8Vector *phse;
1449
1450
1451 REAL8 v2, eta, m, rn, r, rOld, s, p, q, dt, t, v, omega, f, ampl0;
1452 REAL8 omegamatch;
1453 /* Track change to help step back two points in integraion */
1454 REAL8 rpr1=0, rpr2=0, spr1=0, spr2=0, ppr1=0, ppr2=0, qpr1=0, qpr2=0;
1455
1456 void *funcParams1, *funcParams2, *funcParams3;
1457
1458 REAL8Vector *values, *dvalues, *newvalues, *yt, *dym, *dyt;
1459 TofVIn in1;
1460 InspiralPhaseIn in2;
1461 InspiralDerivativesIn in3;
1462 rk4In in4;
1463 rk4GSLIntegrator *integrator = NULL;
1464 pr3In pr3in;
1465 expnCoeffs ak;
1466 expnFunc func;
1467
1468 rOfOmegaIn rofomegain;
1469
1470 /* Functions, tolerances and initial ranges of various bisection root finding calls */
1471 static const REAL8 xacc = 1.0e-12;
1472
1473 REAL8 (*lightRingRadiusFunc)(REAL8, void *) = NULL;
1474 static const REAL8 lightRingMin = 4.;
1475 static const REAL8 lightRingMax = 1.;
1476
1477 REAL8 (*rOfOmegaFunc)(REAL8, void *) = NULL;
1478 static const REAL8 rInitMin = 3.;
1479 static const REAL8 rInitMax = 1000.;
1480
1481 REAL8 (*prInitFunc)(REAL8, void *) = NULL;
1482 static const REAL8 prInitMin = -10.;
1483 static const REAL8 prInitMax = 5.;
1484
1485 /* Variables to allow the waveform to be generated */
1486 /* from a specific fLower */
1487 REAL8 fCurrent; /* The current frequency of the waveform */
1488 BOOLEAN writeToWaveform = 0; /* Set to true when the current frequency
1489 * crosses fLower */
1490 REAL8 sInit, s0 = 0.0; /* Initial phase, and phase to subtract */
1491 REAL8 rmin = 20; /* Smallest value of r at which to generate the waveform */
1492 COMPLEX16 MultSphHarmP; /* Spin-weighted spherical harmonics */
1493 COMPLEX16 MultSphHarmM; /* Spin-weighted spherical harmonics */
1494 REAL4 x1, x2;
1495 UINT4 i, j, k, modeL;
1496 INT4 modeM; /* number of modes required */
1497 REAL4 inclination; /* binary inclination */
1498 REAL4 coa_phase; /* binary coalescence phase */
1499 REAL8 y_1, y_2, z1, z2; /* (2,2) and (2,-2) spherical harmonics needed in (h+,hx) */
1500
1501 /* Used for EOBNR */
1502 COMPLEX8Vector *modefreqs;
1503 UINT4 resampFac = 8;
1504
1505 /* Variables used in injection */
1506 REAL8 unitHz;
1507 REAL8 cosI;/* cosine of system inclination */
1508 REAL8 apFac, acFac;/* extra factor in plus and cross amplitudes */
1509
1510 ak = paramsInit->ak;
1511 func = paramsInit->func;
1512
1513 /* Check order is consistent if using EOBNR and EOB */
1514 if ( params->approximant == EOBNR && params->order != LAL_PNORDER_PSEUDO_FOUR )
1515 {
1516 XLALPrintError( "Order must be LAL_PNORDER_PSEUDO_FOUR for approximant EOBNR.\n" );
1517 XLAL_ERROR( XLAL_EINVAL );
1518 }
1519 else if ( params->approximant == EOB && params->order < LAL_PNORDER_TWO )
1520 {
1521 XLALPrintError( "Order must be LAL_PNORDER_TWO or greater for approximant EOB.\n" );
1522 XLAL_ERROR( XLAL_EINVAL );
1523 }
1524
1525
1526 if (signalvec1) length = signalvec1->length; else if (ff) length = ff->length;
1527
1528 values = XLALCreateREAL8Vector( nn );
1529 dvalues = XLALCreateREAL8Vector( nn );
1530 newvalues = XLALCreateREAL8Vector( nn );
1531 yt = XLALCreateREAL8Vector( nn );
1532 dym = XLALCreateREAL8Vector( nn );
1533 dyt = XLALCreateREAL8Vector( nn );
1534
1535 if ( !values || !dvalues || !newvalues || !yt || !dym || !dyt )
1536 {
1537 XLALDestroyREAL8Vector( values );
1538 XLALDestroyREAL8Vector( dvalues );
1539 XLALDestroyREAL8Vector( newvalues );
1540 XLALDestroyREAL8Vector( yt );
1541 XLALDestroyREAL8Vector( dym );
1542 XLALDestroyREAL8Vector( dyt );
1543 XLAL_ERROR( XLAL_EFUNC );
1544 }
1545
1546 /* Set dt to sampling interval specified by user */
1547 /* But this is changed to 1/16 kHz if the approximant is EOBNR */
1548 dt = 1./params->tSampling;
1549 eta = ak.eta;
1550 m = ak.totalmass;
1551
1552 /* only used in injection case */
1553 unitHz = m*(REAL8)LAL_PI;
1554 cosI = cos( params->inclination );
1555 apFac = -2.0 * (1.0 + cosI*cosI) * eta * params->totalMass * LAL_MRSUN_SI/params->distance;
1556 acFac = -4.0 * cosI * eta * params->totalMass * LAL_MRSUN_SI/params->distance;
1557
1558 /* Set the amplitude depending on whether the distance is given */
1559 if ( params->distance > 0.0 )
1560 ampl0 = -sqrt(64.0*LAL_PI/5.) * eta * params->totalMass * LAL_MRSUN_SI/params->distance;
1561 else
1562 ampl0 = 2.0 * sqrt( LAL_PI / 5.0) * params->signalAmplitude;
1563
1564 /* Check we get a sensible answer */
1565 if ( ampl0 == 0.0 )
1566 {
1567 XLALPrintWarning( "Generating waveform of zero amplitude!!\n" );
1568 }
1569
1570 /* For EOBNR, Check that the 220 QNM freq. is less than the Nyquist freq. */
1571 if ( params->approximant == EOBNR )
1572 {
1573 /* Get QNM frequencies */
1574 modefreqs = XLALCreateCOMPLEX8Vector( 3 );
1575 if ( XLALGenerateQNMFreq( modefreqs, params, 2, 2, 3 ) == XLAL_FAILURE )
1576 {
1577 XLALDestroyCOMPLEX8Vector( modefreqs );
1578 XLALDestroyREAL8Vector( values );
1579 XLALDestroyREAL8Vector( dvalues );
1580 XLALDestroyREAL8Vector( newvalues );
1581 XLALDestroyREAL8Vector( yt );
1582 XLALDestroyREAL8Vector( dym );
1583 XLALDestroyREAL8Vector( dyt );
1584 XLAL_ERROR( XLAL_EFUNC );
1585 }
1586
1587 /* If 220 QNM freq. > Nyquist freq., print warning but continue */
1588 /* Note that we cancelled a factor of 2 occuring on both sides */
1589 if ( params->tSampling < crealf(modefreqs->data[0]) / LAL_PI )
1590 {
1591 XLALPrintWarning( "Ringdown freq. greater than Nyquist freq. "
1592 "Beware of aliasing! Consider increasing the sample rate.\n" );
1593 }
1594 XLALDestroyCOMPLEX8Vector( modefreqs );
1595 }
1596
1597 /* Find the initial velocity given the lower frequency */
1598 t = 0.0;
1599 in1.t = t;
1600 in1.t0 = ak.t0;
1601 in1.v0 = ak.v0;
1602 in1.vlso = ak.vlso;
1603 in1.totalmass = ak.totalmass;
1604 in1.dEnergy = func.dEnergy;
1605 in1.flux = func.flux;
1606 in1.coeffs = &ak;
1607
1608 v = XLALInspiralVelocity( &in1 );
1609 if ( XLAL_IS_REAL8_FAIL_NAN( v ) )
1610 {
1611 XLALDestroyREAL8Vector( values );
1612 XLALDestroyREAL8Vector( dvalues );
1613 XLALDestroyREAL8Vector( newvalues );
1614 XLALDestroyREAL8Vector( yt );
1615 XLALDestroyREAL8Vector( dym );
1616 XLALDestroyREAL8Vector( dyt );
1617 XLAL_ERROR( XLAL_EFUNC );
1618 }
1619
1620 omega = v*v*v;
1621 f = omega/(LAL_PI*m);
1622
1623 /* Then the initial phase */
1624 in2.v0 = ak.v0;
1625 in2.phi0 = params->startPhase;
1626 in2.dEnergy = func.dEnergy;
1627 in2.flux = func.flux;
1628 in2.coeffs = &ak;
1629 s = XLALInspiralPhasing1( v, &in2 );
1630 if ( XLAL_IS_REAL8_FAIL_NAN( s ) )
1631 {
1632 XLALDestroyREAL8Vector( values );
1633 XLALDestroyREAL8Vector( dvalues );
1634 XLALDestroyREAL8Vector( newvalues );
1635 XLALDestroyREAL8Vector( yt );
1636 XLALDestroyREAL8Vector( dym );
1637 XLALDestroyREAL8Vector( dyt );
1638 XLAL_ERROR( XLAL_EFUNC );
1639 }
1640/*
1641 LALInspiralPhasing1(v) gives the GW phase (= twice the orbital phase).
1642 The ODEs we solve give the orbital phase. Therefore, set the
1643 initial phase to be half the GW pahse.
1644*/
1645 s = s/2.;
1646 sInit = s;
1647
1648 /* light ring value - where to stop evolution */
1649 rofomegain.eta = eta;
1650 rofomegain.omega = omega;
1651 pr3in.eta = eta;
1652 pr3in.omegaS = params->OmegaS;
1653 pr3in.zeta2 = params->Zeta2;
1654
1655 /* We will be changing the starting r if it is less than rmin */
1656 /* Therefore, we should reset pr3in.omega later if necessary. */
1657 /* For now we need it so that we can see what the initial r is. */
1658
1659 pr3in.omega = omega;
1660 in3.totalmass = ak.totalmass;
1661 in3.dEnergy = func.dEnergy;
1662 in3.flux = func.flux;
1663 in3.coeffs = &ak;
1664 funcParams3 = (void *) &in3;
1665
1666 funcParams1 = (void *) &rofomegain;
1667
1668 switch (params->order)
1669 {
1670 case LAL_PNORDER_TWO:
1671 case LAL_PNORDER_TWO_POINT_FIVE:
1672 lightRingRadiusFunc = XLALlightRingRadius;
1673 rOfOmegaFunc = XLALrOfOmega;
1674 funcParams2 = (void *) &rofomegain;
1675 break;
1676 case LAL_PNORDER_THREE:
1677 case LAL_PNORDER_THREE_POINT_FIVE:
1678 lightRingRadiusFunc = XLALlightRingRadius3PN;
1679 rOfOmegaFunc = XLALrOfOmega3PN;
1680 funcParams2 = (void *) &pr3in;
1681 break;
1682 case LAL_PNORDER_PSEUDO_FOUR:
1683 lightRingRadiusFunc = XLALlightRingRadiusP4PN;
1684 rOfOmegaFunc = XLALrOfOmegaP4PN;
1685 funcParams2 = (void *) &pr3in;
1686 break;
1687 default:
1688 XLALPrintError( "There are no EOB/EOBNR waveforms implemented at order %d\n", params->order );
1689 XLALDestroyREAL8Vector( values );
1690 XLALDestroyREAL8Vector( dvalues );
1691 XLALDestroyREAL8Vector( newvalues );
1692 XLALDestroyREAL8Vector( yt );
1693 XLALDestroyREAL8Vector( dym );
1694 XLALDestroyREAL8Vector( dyt );
1695 XLAL_ERROR( XLAL_EINVAL );
1696 }
1697 rn = XLALDBisectionFindRoot( lightRingRadiusFunc, lightRingMin,
1698 lightRingMax, xacc, funcParams1);
1699 if ( XLAL_IS_REAL8_FAIL_NAN( rn ) )
1700 {
1701 XLALDestroyREAL8Vector( values );
1702 XLALDestroyREAL8Vector( dvalues );
1703 XLALDestroyREAL8Vector( newvalues );
1704 XLALDestroyREAL8Vector( yt );
1705 XLALDestroyREAL8Vector( dym );
1706 XLALDestroyREAL8Vector( dyt );
1707 XLAL_ERROR( XLAL_EFUNC );
1708 }
1709
1710 r = XLALDBisectionFindRoot( rOfOmegaFunc, rInitMin, rInitMax, xacc, funcParams2);
1711 if ( XLAL_IS_REAL8_FAIL_NAN( rn ) )
1712 {
1713 XLALDestroyREAL8Vector( values );
1714 XLALDestroyREAL8Vector( dvalues );
1715 XLALDestroyREAL8Vector( newvalues );
1716 XLALDestroyREAL8Vector( yt );
1717 XLALDestroyREAL8Vector( dym );
1718 XLALDestroyREAL8Vector( dyt );
1719 XLAL_ERROR( XLAL_EFUNC );
1720 }
1721
1722 /* Is the initial condition sensible? */
1723 /* For EOB this check should be done to prevent templates/injections being too short. */
1724 /* No need to do this for EOBNR, as this gets a ringdown attached. */
1725 if ( params->approximant == EOB && r < 6 )
1726 {
1727 XLALPrintError( "EOB:initialCondition:Initial r found = %f "
1728 "too small (below 6 no waveform is generated)\n", r );
1729 XLALDestroyREAL8Vector( values );
1730 XLALDestroyREAL8Vector( dvalues );
1731 XLALDestroyREAL8Vector( newvalues );
1732 XLALDestroyREAL8Vector( yt );
1733 XLALDestroyREAL8Vector( dym );
1734 XLALDestroyREAL8Vector( dyt );
1735 XLAL_ERROR( XLAL_ERANGE );
1736 }
1737
1738 /* We want the waveform to generate from a point which won't cause */
1739 /* problems with the initial conditions. Therefore we force the code */
1740 /* to start at least at r = rmin (in units of M). */
1741
1742 r = (r<rmin) ? rmin : r;
1743
1744 pr3in.in3copy = in3;
1745 pr3in.r = r;
1746
1747 /* Now that r is changed recompute omega corresponding */
1748 /* to that r and only then compute initial pr and pphi */
1749
1750 switch (params->order)
1751 {
1752 case LAL_PNORDER_TWO:
1753 case LAL_PNORDER_TWO_POINT_FIVE:
1754 omega = XLALOmegaOfR2PN (r, eta);
1755 pr3in.omega = omega;
1756 q = XLALpphiInit( r, eta);
1757 p = XLALprInit(r, &in3);
1758 in4.function = LALHCapDerivatives;
1759 break;
1760 case LAL_PNORDER_THREE:
1761 case LAL_PNORDER_THREE_POINT_FIVE:
1762 omega = XLALOmegaOfR3PN(r, &pr3in);
1763 pr3in.omega = omega;
1764 q = XLALpphiInit3PN( r, eta, params->OmegaS);
1765 prInitFunc = XLALprInit3PN;
1766 /* first we compute vr (we need coeef->Fp6) */
1767 pr3in.q = q;
1768 funcParams2 = (void *) &pr3in;
1769 pr3in.vr = XLALvr3PN(funcParams2);
1770 /* then we compute the initial value of p */
1771 p = XLALDBisectionFindRoot( prInitFunc, prInitMin, prInitMax, xacc, funcParams2);
1772 if ( XLAL_IS_REAL8_FAIL_NAN( p ) )
1773 {
1774 XLALDestroyREAL8Vector( values );
1775 XLALDestroyREAL8Vector( dvalues );
1776 XLALDestroyREAL8Vector( newvalues );
1777 XLALDestroyREAL8Vector( yt );
1778 XLALDestroyREAL8Vector( dym );
1779 XLALDestroyREAL8Vector( dyt );
1780 XLAL_ERROR( XLAL_EFUNC );
1781 }
1782 in4.function = LALHCapDerivatives3PN;
1783 break;
1784 case LAL_PNORDER_PSEUDO_FOUR:
1785 omega = XLALOmegaOfRP4PN (r, &pr3in);
1786 pr3in.omega = omega;
1787 q = XLALpphiInitP4PN( r, eta, params->OmegaS);
1788 prInitFunc = XLALprInitP4PN;
1789 /* first we compute vr (we need coeef->Fp6) */
1790 pr3in.q = q;
1791 funcParams2 = (void *) &pr3in;
1792 pr3in.vr = XLALvrP4PN(funcParams2);
1793 /* then we compute the initial value of p */
1794 p = XLALDBisectionFindRoot( prInitFunc, prInitMin, prInitMax, xacc, funcParams2);
1795 if ( XLAL_IS_REAL8_FAIL_NAN( p ) )
1796 {
1797 XLALDestroyREAL8Vector( values );
1798 XLALDestroyREAL8Vector( dvalues );
1799 XLALDestroyREAL8Vector( newvalues );
1800 XLALDestroyREAL8Vector( yt );
1801 XLALDestroyREAL8Vector( dym );
1802 XLALDestroyREAL8Vector( dyt );
1803 XLAL_ERROR( XLAL_EFUNC );
1804 }
1805 in4.function = LALHCapDerivativesP4PN;
1806 break;
1807 default:
1808 XLALPrintError( "There are no EOB/EOBNR waveforms implemented at order %d\n", params->order);
1809 XLALDestroyREAL8Vector( values );
1810 XLALDestroyREAL8Vector( dvalues );
1811 XLALDestroyREAL8Vector( newvalues );
1812 XLALDestroyREAL8Vector( yt );
1813 XLALDestroyREAL8Vector( dym );
1814 XLALDestroyREAL8Vector( dyt );
1815 XLAL_ERROR( XLAL_EFUNC );
1816 }
1817
1818 values->data[0] = r;
1819 values->data[1] = s;
1820 values->data[2] = p;
1821 values->data[3] = q;
1822#if 0
1823 sprintf(message, "In EOB Initial values of r=%10.5e p=%10.5e q=%10.5e\n", r, p, q);
1824 LALInfo(status, message);
1825#endif
1826
1827 in4.y = values;
1828 in4.h = dt/m;
1829 in4.n = nn;
1830 in4.yt = yt;
1831 in4.dym = dym;
1832 in4.dyt = dyt;
1833
1834 /* Allocate memory for temporary arrays */
1835 sig1 = XLALCreateREAL4Vector ( length );
1836 sig2 = XLALCreateREAL4Vector ( length );
1837 ampl = XLALCreateREAL4Vector ( length*2 );
1838 freq = XLALCreateREAL4Vector ( length );
1839 phse = XLALCreateREAL8Vector ( length );
1840
1841 if ( !sig1 || !sig2 || !ampl || !freq || !phse )
1842 {
1843 XLALDestroyREAL4Vector( sig1 );
1844 XLALDestroyREAL4Vector( sig2 );
1845 XLALDestroyREAL4Vector( ampl );
1846 XLALDestroyREAL4Vector( freq );
1847 XLALDestroyREAL8Vector( phse );
1848 XLALDestroyREAL8Vector( values );
1849 XLALDestroyREAL8Vector( dvalues );
1850 XLALDestroyREAL8Vector( newvalues );
1851 XLALDestroyREAL8Vector( yt );
1852 XLALDestroyREAL8Vector( dym );
1853 XLALDestroyREAL8Vector( dyt );
1854 XLAL_ERROR( XLAL_ENOMEM );
1855 }
1856
1857 memset(sig1->data, 0, sig1->length * sizeof( REAL4 ));
1858 memset(sig2->data, 0, sig2->length * sizeof( REAL4 ));
1859 memset(ampl->data, 0, ampl->length * sizeof( REAL4 ));
1860 memset(freq->data, 0, freq->length * sizeof( REAL4 ));
1861 memset(phse->data, 0, phse->length * sizeof( REAL8 ));
1862
1863 /* Initialize the GSL integrator */
1864 if (!(integrator = XLALRungeKutta4Init(nn, &in4)))
1865 {
1866 XLALDestroyREAL4Vector( sig1 );
1867 XLALDestroyREAL4Vector( sig2 );
1868 XLALDestroyREAL4Vector( ampl );
1869 XLALDestroyREAL4Vector( freq );
1870 XLALDestroyREAL8Vector( phse );
1871 XLALDestroyREAL8Vector( values );
1872 XLALDestroyREAL8Vector( dvalues );
1873 XLALDestroyREAL8Vector( newvalues );
1874 XLALDestroyREAL8Vector( yt );
1875 XLALDestroyREAL8Vector( dym );
1876 XLALDestroyREAL8Vector( dyt );
1877 XLAL_ERROR( XLAL_EFUNC );
1878 }
1879
1880 count = 0;
1881 if (a || signalvec2)
1882 params->nStartPad = 0; /* must be zero for templates and injection */
1883
1884 count = params->nStartPad;
1885
1886 /* Calculate the initial value of omega */
1887 in4.function(values, dvalues, funcParams3);
1888 omega = dvalues->data[1];
1889
1890 /* Begin integration loop here */
1891 t = 0.0;
1892 rOld = r+0.1;
1893
1894 omegamatch = -0.01 + 0.133 + 0.183 * params->eta + 0.161 * params->eta * params->eta;
1895
1896 while (r > rn && r < rOld)
1897 {
1898 if (count > length)
1899 {
1900 XLALRungeKutta4Free( integrator );
1901 XLALDestroyREAL4Vector( sig1 );
1902 XLALDestroyREAL4Vector( sig2 );
1903 XLALDestroyREAL4Vector( ampl );
1904 XLALDestroyREAL4Vector( freq );
1905 XLALDestroyREAL8Vector( phse );
1906 XLALDestroyREAL8Vector( values );
1907 XLALDestroyREAL8Vector( dvalues );
1908 XLALDestroyREAL8Vector( newvalues );
1909 XLALDestroyREAL8Vector( yt );
1910 XLALDestroyREAL8Vector( dym );
1911 XLALDestroyREAL8Vector( dyt );
1912 XLALPrintError( "Waveform evolution has run off the end of the vector" );
1913 XLAL_ERROR( XLAL_EBADLEN );
1914 }
1915
1916 rOld = r;
1917
1918 fCurrent = omega / (LAL_PI*m);
1919 if (!writeToWaveform)
1920 {
1921 s0 = s - sInit;
1922 if (r > rmin || fCurrent > f || fabs(fCurrent - f) < 1.0e-5)
1923 {
1924 writeToWaveform = 1;
1925 }
1926 }
1927
1928 v = cbrt(omega);
1929 v2 = v*v;
1930
1931 if (writeToWaveform)
1932 {
1933 double st, amp;
1934 i = count;
1935 j = 2*count;
1936 k = j+1;
1937 st = 2.*(s - s0);
1938 amp = ampl0 * v2;
1939 /*--------------------------------------------------------
1940 First we generate the real and imagninary parts of h22
1941 --------------------------------------------------------*/
1942 sig1->data[i] = (REAL4)( amp * cos(st) );
1943 sig2->data[i] = -(REAL4)( amp * sin(st) );
1944 /*----------------------------------------------------------
1945 ... then the frequency, amplitude of h+ and hx and phase
1946 ----------------------------------------------------------*/
1947 freq->data[i] = (REAL4)( omega );
1948 ampl->data[j] = (REAL4)( apFac * v2 );
1949 ampl->data[k] = (REAL4)( acFac * v2 );
1950 phse->data[i] = (REAL8)( st );
1951 }
1952
1953 /* Integrate one step forward */
1954 in4.dydx = dvalues;
1955 in4.x = t/m;
1956 if ( XLALRungeKutta4( newvalues, integrator, funcParams3) == XLAL_FAILURE )
1957 {
1958 XLALRungeKutta4Free( integrator );
1959 XLALDestroyREAL4Vector( sig1 );
1960 XLALDestroyREAL4Vector( sig2 );
1961 XLALDestroyREAL4Vector( ampl );
1962 XLALDestroyREAL4Vector( freq );
1963 XLALDestroyREAL8Vector( phse );
1964 XLALDestroyREAL8Vector( values );
1965 XLALDestroyREAL8Vector( dvalues );
1966 XLALDestroyREAL8Vector( newvalues );
1967 XLALDestroyREAL8Vector( yt );
1968 XLALDestroyREAL8Vector( dym );
1969 XLALDestroyREAL8Vector( dyt );
1970 XLAL_ERROR( XLAL_EFUNC );
1971 }
1972
1973 /* We need to track the dynamical variables prior to the current step */
1974 if(ndx>1)
1975 {
1976 rpr2 = rpr1;
1977 spr2=spr1;
1978 ppr2=ppr1;
1979 qpr2=qpr1;
1980 }
1981
1982 /* These are the current values of the dynamical variables */
1983 rpr1=r;
1984 spr1=s;
1985 ppr1=p;
1986 qpr1=q;
1987
1988 /* Update the values of the dynamical variables */
1989 r = values->data[0] = newvalues->data[0];
1990 s = values->data[1] = newvalues->data[1];
1991 p = values->data[2] = newvalues->data[2];
1992 q = values->data[3] = newvalues->data[3];
1993
1994 /* Compute the derivaties at the new location */
1995 in4.function(values, dvalues, funcParams3);
1996 omega = dvalues->data[1];
1997
1998 /*----------------------------------------------------------------------*/
1999 /* We are going to terminate waveform generation if omega is greater */
2000 /* than omegamatch - the frequency at which the ringdown is matched to */
2001 /* merger waveform */
2002 /*----------------------------------------------------------------------*/
2003 if ( params->approximant == EOBNR && (r < rn || omega > omegamatch ) && !higherSR)
2004 {
2005 /* We are now going to work with a higher sampling rate */
2006 /* Sometime in the future we might change code so that */
2007 /* a higher sampling rate is used only if required */
2008 higherSR = 1;
2009 /*-------------------------------------------------------------*/
2010 /* We are going to decrease the number of points by 2 */
2011 /* In reality, note that we are really using the previous */
2012 /* point from the current step; 2 is needed below instead of 1 */
2013 /* only because count is incremented before returning to the */
2014 /* continuing the integration; the same is true with dt */
2015 /*-------------------------------------------------------------*/
2016 count -= 2;
2017 hiSRndx = count+1;
2018 t -= dt;
2019 dt /= (double) resampFac;
2020 t -= dt;
2021 in4.h = dt/m;
2022
2023 r = values->data[0] = rpr2;
2024 s = values->data[1] = spr2;
2025 p = values->data[2] = ppr2;
2026 q = values->data[3] = qpr2;
2027
2028 /*----------------------------------------------------------------------*/
2029 /* Integration will stop if rOld is not reset to a value greater than r */
2030 /*----------------------------------------------------------------------*/
2031
2032 rOld = r+0.1;
2033
2034 in4.function(values, dvalues, funcParams3);
2035 omega = dvalues->data[1];
2036 fCurrent = omega/(LAL_PI*m);
2037 }
2038
2039 if (writeToWaveform)
2040 {
2041 if (!higherSR)
2042 t = (++count-params->nStartPad) * dt;
2043 else
2044 {
2045 t += dt;
2046 count++;
2047 }
2048
2049 }
2050 ndx++;
2051 }
2052
2053 /*----------------------------------------------------------------------*/
2054 /* Record the final cutoff frequency of BD Waveforms for record keeping */
2055 /* ---------------------------------------------------------------------*/
2056 params->vFinal = v;
2057 if (signalvec1 && !signalvec2) params->tC = t;
2058 if (params->approximant == EOB)
2059 {
2060 params->fFinal = v*v*v/(LAL_PI*m);
2061 }
2062 else
2063 {
2064 params->fFinal = params->tSampling/2.;
2065 }
2066
2067 /* ------------------------------------------------------------------*/
2068 /* This is the count for the inspiral part only. It is changed below */
2069 /* when the merger part is added; the phase is changed artificially */
2070 /* by a small increment for each data point added but be warned that */
2071 /* it is not the total accumuated phase. */
2072 /*-------------------------------------------------------------------*/
2073 *countback = count;
2074
2075 XLALRungeKutta4Free( integrator );
2076 XLALDestroyREAL8Vector( values );
2077 XLALDestroyREAL8Vector( dvalues );
2078 XLALDestroyREAL8Vector( newvalues );
2079 XLALDestroyREAL8Vector( yt );
2080 XLALDestroyREAL8Vector( dym );
2081 XLALDestroyREAL8Vector( dyt );
2082 /*--------------------------------------------------------------
2083 * Attach the ringdown waveform to the end of inspiral if
2084 * the approximant is EOBNR
2085 -------------------------------------------------------------*/
2086 if (params->approximant == EOBNR)
2087 {
2088 REAL8 tmpSamplingRate = params->tSampling;
2089 params->tSampling *= resampFac;
2090
2091 /* Check that we have enough points to perform the attachment */
2092 /* at the higher sample rate. If not, we bail out here */
2093 /* Hard coded for now... */
2094 if ( count - hiSRndx < 7 )
2095 {
2096 XLALPrintError( " We do not have enough points at a higher SR.\n" );
2097 XLALDestroyREAL4Vector( sig1 );
2098 XLALDestroyREAL4Vector( sig2 );
2099 XLALDestroyREAL4Vector( ampl );
2100 XLALDestroyREAL4Vector( freq );
2101 XLALDestroyREAL8Vector( phse );
2102 XLAL_ERROR( XLAL_EFAILED );
2103 }
2104
2105 if ( XLALInspiralAttachRingdownWave( freq, sig1, sig2, params ) == XLAL_FAILURE )
2106 {
2107 XLALDestroyREAL4Vector( sig1 );
2108 XLALDestroyREAL4Vector( sig2 );
2109 XLALDestroyREAL4Vector( ampl );
2110 XLALDestroyREAL4Vector( freq );
2111 XLALDestroyREAL8Vector( phse );
2112 XLAL_ERROR( XLAL_EFUNC );
2113 }
2114 params->tSampling = tmpSamplingRate;
2115 count = hiSRndx;
2116 for(j=hiSRndx; j<length; j+=resampFac)
2117 {
2118 sig1->data[count] = sig1->data[j];
2119 sig2->data[count] = sig2->data[j];
2120 freq->data[count] = freq->data[j];
2121 if (sig1->data[count] == 0)
2122 {
2123 k = count;
2124 while (++k<length)
2125 {
2126 sig1->data[k] = 0.;
2127 sig2->data[k] = 0.;
2128 freq->data[k] = 0.;
2129 }
2130 break;
2131 }
2132 count++;
2133 }
2134 }
2135 *countback = count;
2136
2137 /*-------------------------------------------------------------------
2138 * Compute the spherical harmonics required for constructing (h+,hx).
2139 * We are going to choose coa_phase to be zero. This perhaps should be
2140 * made compatible with the wave CoherentGW handles the phase at
2141 * coalecence. I have no idea how I (i.e., Sathya) might be able to
2142 * do this for EOBNR as there is no such thing as "phase at merger".
2143 -------------------------------------------------------------------*/
2144 inclination = (REAL4)params->inclination;
2145 coa_phase = 0.;
2146 /* -----------------------------------------------------------------
2147 * Attaching the (2,2) Spherical Harmonic
2148 * need some error checking
2149 *----------------------------------------*/
2150 modeL = 2;
2151 modeM = 2;
2152 if ( XLALSphHarm( &MultSphHarmP, modeL, modeM, inclination, coa_phase ) == XLAL_FAILURE )
2153 {
2154 XLALDestroyREAL4Vector( sig1 );
2155 XLALDestroyREAL4Vector( sig2 );
2156 XLALDestroyREAL4Vector( ampl );
2157 XLALDestroyREAL4Vector( freq );
2158 XLALDestroyREAL8Vector( phse );
2159 XLAL_ERROR( XLAL_EFUNC );
2160 }
2161
2162 modeM = -2;
2163 if ( XLALSphHarm( &MultSphHarmM, modeL, modeM, inclination, coa_phase ) == XLAL_FAILURE )
2164 {
2165 XLALDestroyREAL4Vector( sig1 );
2166 XLALDestroyREAL4Vector( sig2 );
2167 XLALDestroyREAL4Vector( ampl );
2168 XLALDestroyREAL4Vector( freq );
2169 XLALDestroyREAL8Vector( phse );
2170 XLAL_ERROR( XLAL_EFUNC );
2171 }
2172
2173 y_1 = creal(MultSphHarmP) + creal(MultSphHarmM);
2174 y_2 = cimag(MultSphHarmM) - cimag(MultSphHarmP);
2175 z1 = - cimag(MultSphHarmM) - cimag(MultSphHarmP);
2176 z2 = creal(MultSphHarmM) - creal(MultSphHarmP);
2177
2178#if 0
2179 sprintf(message, "MultSphHarm2,+2 re=%10.5e im=%10.5e\n", MultSphHarmP.re, MultSphHarmP.im);
2180 LALInfo(status, message);
2181 sprintf(message, "MultSphHarm2,-2 re=%10.5e im=%10.5e\n", MultSphHarmP.re, MultSphHarmM.im);
2182 LALInfo(status, message);
2183#endif
2184
2185 /* Next, compute h+ and hx from h22, h22*, Y22, Y2-2 */
2186 for ( i = 0; i < sig1->length; i++)
2187 {
2188 freq->data[i] /= unitHz;
2189 x1 = sig1->data[i];
2190 x2 = sig2->data[i];
2191 sig1->data[i] = (x1 * y_1) + (x2 * y_2);
2192 sig2->data[i] = (x1 * z1) + (x2 * z2);
2193 if (x1 || x2)
2194 {
2195 /*
2196 * If the ringdown modes were added then artificially increase the phasing so
2197 * that it is nonzero until the end of the ringdown. When ringdown modes are
2198 * added the phase information is not used in injetions, only hplus and hcross
2199 * and therefore it shouldn't matter what phasing is as long as it is nonzero.
2200 */
2201 if ( params->approximant == EOBNR && i >= hiSRndx )
2202 {
2203 phse->data[i] = phse->data[i-1]+LAL_PI/20.;
2204 }
2205 }
2206 }
2207
2208 /*------------------------------------------------------
2209 * If required by the user copy other data sets to the
2210 * relevant arrays
2211 ------------------------------------------------------*/
2212 if (h)
2213 {
2214 for(i = 0; i < length; i++)
2215 {
2216 j = 2*i;
2217 k = j+1;
2218 h->data[j] = sig1->data[i];
2219 h->data[k] = sig2->data[i];
2220 }
2221 }
2222 if (signalvec1) memcpy(signalvec1->data , sig1->data, length * (sizeof(REAL4)));
2223 if (signalvec2) memcpy(signalvec2->data , sig2->data, length * (sizeof(REAL4)));
2224 if (ff) memcpy(ff->data , freq->data, length * (sizeof(REAL4)));
2225 if (a) memcpy(a->data , ampl->data, 2*length*(sizeof(REAL4)));
2226 if (phi) memcpy(phi->data , phse->data, length * (sizeof(REAL8)));
2227
2228#if 0
2229 sprintf(message, "fFinal=%10.5e count=%d\n", params->fFinal, *countback);
2230 LALInfo(status, message);
2231#endif
2232
2233 /* Clean up */
2234 XLALDestroyREAL4Vector ( sig1 );
2235 XLALDestroyREAL4Vector ( sig2 );
2236 XLALDestroyREAL4Vector ( ampl );
2237 XLALDestroyREAL4Vector ( freq );
2238 XLALDestroyREAL8Vector ( phse );
2239
2240
2241 return XLAL_SUCCESS;
2242}
XLALprInitP4PN
static REAL8 XLALprInitP4PN(REAL8, void *params)
Definition: LALEOBWaveform.c:705
XLALpphiInit3PN
static REAL8 XLALpphiInit3PN(REAL8 r, REAL8 eta, REAL8 omegaS)
Definition: LALEOBWaveform.c:400
XLALpphiInitP4PN
static REAL8 XLALpphiInitP4PN(REAL8 r, REAL8 eta, REAL8 omegaS)
Definition: LALEOBWaveform.c:674
XLALvrP4PN
static REAL8 XLALvrP4PN(void *params)
Definition: LALEOBWaveform.c:917
XLALlightRingRadius
static REAL8 XLALlightRingRadius(REAL8 r, void *params)
Definition: LALEOBWaveform.c:335
XLALOmegaOfRP4PN
static REAL8 XLALOmegaOfRP4PN(REAL8 r, void *params)
Definition: LALEOBWaveform.c:756
XLALrOfOmega
static REAL8 XLALrOfOmega(REAL8 r, void *params)
Definition: LALEOBWaveform.c:310
LALHCapDerivatives
static void LALHCapDerivatives(REAL8Vector *values, REAL8Vector *dvalues, void *funcParams)
Definition: LALEOBWaveform.c:352
XLALprInit3PN
static REAL8 XLALprInit3PN(REAL8, void *params)
Definition: LALEOBWaveform.c:428
XLALEOBWaveformEngine
static int XLALEOBWaveformEngine(REAL4Vector *signalvec1, REAL4Vector *signalvec2, REAL4Vector *h, REAL4Vector *a, REAL4Vector *ff, REAL8Vector *phi, UINT4 *countback, InspiralTemplate *params, InspiralInit *paramsInit)
Definition: LALEOBWaveform.c:1432
XLALlightRingRadius3PN
static REAL8 XLALlightRingRadius3PN(REAL8 r, void *params)
Definition: LALEOBWaveform.c:523
XLALOmegaOfR3PN
static REAL8 XLALOmegaOfR3PN(REAL8 r, void *params)
Definition: LALEOBWaveform.c:475
XLALlightRingRadiusP4PN
static REAL8 XLALlightRingRadiusP4PN(REAL8 r, void *params)
Definition: LALEOBWaveform.c:810
XLALpphiInit
static REAL8 XLALpphiInit(REAL8 r, REAL8 eta)
Definition: LALEOBWaveform.c:274
LALHCapDerivativesP4PN
static void LALHCapDerivativesP4PN(REAL8Vector *values, REAL8Vector *dvalues, void *funcParams)
Definition: LALEOBWaveform.c:844
XLALEOBWaveform
int XLALEOBWaveform(REAL4Vector *signalvec, InspiralTemplate *params)
Definition: LALEOBWaveform.c:972
XLALEOBWaveformForInjection
int XLALEOBWaveformForInjection(CoherentGW *waveform, InspiralTemplate *params, PPNParamStruc *ppnParams)
Definition: LALEOBWaveform.c:1117
XLALEOBWaveformTemplates
int XLALEOBWaveformTemplates(REAL4Vector *signalvec1, REAL4Vector *signalvec2, InspiralTemplate *params)
Definition: LALEOBWaveform.c:1037
XLALrOfOmegaP4PN
static REAL8 XLALrOfOmegaP4PN(REAL8 r, void *params)
Definition: LALEOBWaveform.c:791
XLALOmegaOfR2PN
static REAL8 XLALOmegaOfR2PN(REAL8 r, REAL8 eta)
Definition: LALEOBWaveform.c:288
XLALprInit
static REAL8 XLALprInit(REAL8 r, InspiralDerivativesIn *ak)
Definition: LALEOBWaveform.c:224
LALHCapDerivatives3PN
static void LALHCapDerivatives3PN(REAL8Vector *values, REAL8Vector *dvalues, void *funcParams)
Definition: LALEOBWaveform.c:551
XLALvr3PN
static REAL8 XLALvr3PN(void *params)
Definition: LALEOBWaveform.c:627
ninty4by3etc
#define ninty4by3etc
Definition: LALEOBWaveform.c:122
XLALrOfOmega3PN
static REAL8 XLALrOfOmega3PN(REAL8 r, void *params)
Definition: LALEOBWaveform.c:506
XLALInspiralPhasing1
REAL8 XLALInspiralPhasing1(REAL8 v, void *params)
Definition: LALInspiralPhasing1.c:68
XLALRungeKutta4Free
void XLALRungeKutta4Free(rk4GSLIntegrator *integrator)
Definition: LALRungeKutta4.c:228
XLALInspiralAttachRingdownWave
INT4 XLALInspiralAttachRingdownWave(REAL4Vector *Omega, REAL4Vector *signalvec1, REAL4Vector *signalvec2, InspiralTemplate *params)
Definition: LALInspiralRingdownWave.c:1136
XLALInspiralSetup
int XLALInspiralSetup(expnCoeffs *ak, InspiralTemplate *params)
Definition: LALInspiralSetup.c:88
XLALInspiralInit
int XLALInspiralInit(InspiralTemplate *params, InspiralInit *paramsInit)
Definition: LALInspiralInit.c:87
XLALInspiralVelocity
REAL8 XLALInspiralVelocity(TofVIn *params)
Definition: LALInspiralVelocity.c:59
XLALGenerateQNMFreq
INT4 XLALGenerateQNMFreq(COMPLEX8Vector *modefreqs, InspiralTemplate *params, UINT4 l, UINT4 m, UINT4 nmodes)
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
s
int s
i
double i
pr
const double pr
u
const double u
a4
const double a4
u3
const double u3
r2
const double r2
u5
const double u5
u2
const double u2
u4
const double u4
B
const double B
XLALDestroyREAL4VectorSequence
void XLALDestroyREAL4VectorSequence(REAL4VectorSequence *vecseq)
XLALCreateREAL4VectorSequence
REAL4VectorSequence * XLALCreateREAL4VectorSequence(UINT4 length, UINT4 veclen)
XLALDBisectionFindRoot
REAL8 XLALDBisectionFindRoot(REAL8(*y)(REAL8, void *), REAL8 xmin, REAL8 xmax, REAL8 xacc, void *params)
GENERATEPPNINSPIRALH_EFSTOP
#define GENERATEPPNINSPIRALH_EFSTOP
Reached requested termination frequency.
Definition: GeneratePPNInspiral.h:202
LAL_PI
#define LAL_PI
LAL_MRSUN_SI
#define LAL_MRSUN_SI
BOOLEAN
unsigned char BOOLEAN
COMPLEX16
double complex COMPLEX16
REAL8
double REAL8
UINT4
uint32_t UINT4
INT4
int32_t INT4
REAL4
float REAL4
LALNameLength
LALNameLength
LALInfo
int LALInfo(LALStatus *status, const char *info)
LALPNOrder
LALPNOrder
EOB
EOB
EOBNR
EOBNR
LAL_PNORDER_TWO_POINT_FIVE
LAL_PNORDER_TWO_POINT_FIVE
LAL_PNORDER_THREE
LAL_PNORDER_THREE
LAL_PNORDER_TWO
LAL_PNORDER_TWO
LAL_PNORDER_PSEUDO_FOUR
LAL_PNORDER_PSEUDO_FOUR
LAL_PNORDER_THREE_POINT_FIVE
LAL_PNORDER_THREE_POINT_FIVE
r
static const INT4 r
m
static const INT4 m
q
static const INT4 q
a
static const INT4 a
XLALSphHarm
INT4 XLALSphHarm(COMPLEX16 *out, UINT4 L, INT4 M, REAL4 theta, REAL4 phi)
XLALCreateREAL4TimeSeries
REAL4TimeSeries * XLALCreateREAL4TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
XLALDestroyREAL8TimeSeries
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
XLALDestroyREAL4TimeSeries
void XLALDestroyREAL4TimeSeries(REAL4TimeSeries *series)
XLALCreateREAL8TimeSeries
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
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)
XLALPrintInfo
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_ERROR
#define XLAL_ERROR(...)
XLALPrintError
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLALPrintWarning
int int XLALPrintWarning(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_IS_REAL8_FAIL_NAN
#define XLAL_IS_REAL8_FAIL_NAN(val)
XLAL_EBADLEN
XLAL_EBADLEN
XLAL_ENOMEM
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFAULT
XLAL_ERANGE
XLAL_ERANGE
XLAL_EFUNC
XLAL_EFUNC
XLAL_EDOM
XLAL_EDOM
XLAL_EINVAL
XLAL_EINVAL
XLAL_EFAILED
XLAL_EFAILED
XLAL_FAILURE
XLAL_FAILURE
lalinspiral.git_version.status
string status
Definition: git_version.py:32
p
p
x
x
thinca_brute_force_coinc.dt
dt
Definition: thinca_brute_force_coinc.py:52
COMPLEX8Vector::data
COMPLEX8 * data
CoherentGW::shift
REAL4TimeSeries * shift
CoherentGW::position
SkyPosition position
CoherentGW::a
REAL4TimeVectorSeries * a
CoherentGW::phi
REAL8TimeSeries * phi
CoherentGW::h
REAL4TimeVectorSeries * h
CoherentGW::f
REAL4TimeSeries * f
InspiralDerivativesIn
Structure used as an input to compute the derivatives needed in solving the phasing formula when the ...
Definition: LALInspiral.h:607
InspiralDerivativesIn::dEnergy
EnergyFunction * dEnergy
Definition: LALInspiral.h:609
InspiralDerivativesIn::flux
FluxFunction * flux
Definition: LALInspiral.h:610
InspiralDerivativesIn::coeffs
expnCoeffs * coeffs
Definition: LALInspiral.h:611
InspiralInit::func
expnFunc func
Definition: LALInspiral.h:680
InspiralInit::ak
expnCoeffs ak
Definition: LALInspiral.h:679
InspiralPhaseIn
Structures used to compute the phase of the signal from the ‘beginning’, when the veolcity parameter ...
Definition: LALInspiral.h:654
InspiralPhaseIn::dEnergy
EnergyFunction * dEnergy
Definition: LALInspiral.h:657
InspiralPhaseIn::flux
FluxFunction * flux
Definition: LALInspiral.h:658
InspiralPhaseIn::coeffs
expnCoeffs * coeffs
Definition: LALInspiral.h:659
InspiralTemplate
The inspiral waveform parameter structure containing information about the waveform to be generated.
Definition: LALInspiral.h:205
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::phi
REAL4 phi
The phase at coalescence (or arbitrary reference phase for a post -Newtonian approximation),...
Definition: GeneratePPNInspiral.h:248
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::data
REAL4Sequence * data
REAL4TimeVectorSeries::name
CHAR name[LALNameLength]
REAL4TimeVectorSeries::data
REAL4VectorSequence * data
REAL4Vector::data
REAL4 * data
REAL8TimeSeries::data
REAL8Sequence * data
REAL8Vector::data
REAL8 * data
TofVIn::t0
REAL8 t0
Definition: LALInspiral.h:567
TofVIn::v0
REAL8 v0
Definition: LALInspiral.h:566
TofVIn::totalmass
REAL8 totalmass
Definition: LALInspiral.h:569
TofVIn::dEnergy
EnergyFunction dEnergy
TofVIn::coeffs
expnCoeffsdEnergyFlux * coeffs
TofVIn::t
REAL8 t
Definition: LALInspiral.h:565
TofVIn::flux
FluxFunction flux
TofVIn::vlso
REAL8 vlso
Definition: LALInspiral.h:568
expnCoeffs
This structure contains various post-Newtonian and P-approximant expansion coefficients; the meanings...
Definition: LALInspiral.h:399
expnCoeffs::totalmass
REAL8 totalmass
Definition: LALInspiral.h:474
expnCoeffs::vlso
REAL8 vlso
Definition: LALInspiral.h:487
expnCoeffs::zeta2
REAL8 zeta2
Definition: LALInspiral.h:479
expnFunc
Structure to hold the pointers to the generic functions defined above.
Definition: LALInspiral.h:546
expnFunc::dEnergy
EnergyFunction * dEnergy
Definition: LALInspiral.h:547
expnFunc::flux
FluxFunction * flux
Definition: LALInspiral.h:548
params::inclination
double inclination
params::distance
double distance
pr3In::vr
REAL8 vr
Definition: LALEOBPPWaveform.c:71
pr3In::eta
REAL8 eta
Definition: LALEOBPPWaveform.c:71
pr3In::omega
REAL8 omega
Definition: LALEOBPPWaveform.c:71
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
epoch
enum @1 epoch
deltaT
double deltaT