Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInspiral 5.0.3.1-ea7c608
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALInspiralStationaryPhaseApprox2.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2007 Jolien Creighton, B.S. Sathyaprakash, Thomas Cokelaer
3*
4* This program is free software; you can redistribute it and/or modify
5* it under the terms of the GNU General Public License as published by
6* the Free Software Foundation; either version 2 of the License, or
7* (at your option) any later version.
8*
9* This program is distributed in the hope that it will be useful,
10* but WITHOUT ANY WARRANTY; without even the implied warranty of
11* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12* GNU General Public License for more details.
13*
14* You should have received a copy of the GNU General Public License
15* along with with program; see the file COPYING. If not, write to the
16* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17* MA 02110-1301 USA
18*/
19
20/**
21 * \author B.S. Sathyaprakash
22 * \file
23 * \ingroup LALInspiral_h
24 *
25 * \brief This module computes the usual stationary phase approximation to the
26 * Fourier transform of a chirp waveform given by \eqref{eq_InspiralFourierPhase_f2}.
27 *
28 * ### Prototypes ###
29 *
30 * <tt>XLALInspiralStationaryPhaseApprox2()</tt>
31 * <ul>
32 * <li> \c signalvec: Output containing the inspiral waveform.
33 * </li><li> \c params: Input containing binary chirp parameters.
34 * </li></ul>
35 *
36 * ### Description ###
37 *
38 * Computes the Fourier transform of the chirp signal in the stationary
39 * phase approximation and returns the real and imaginary parts of the
40 * Fourier domain signal in the convention of fftw. For a signal vector
41 * of length <tt>n=signalvec->length</tt> (\c n even):
42 * <ul>
43 * <li> <tt>signalvec->data[0]</tt> is the \e real 0th frequency component of the Fourier transform.
44 * </li><li> <tt>signalvec->data[n/2]</tt> is the \e real Nyquist frequency component of the Fourier transform.
45 * </li><li> <tt>signalvec->data[k]</tt> and <tt>signalvec->data[n-k],</tt> for <tt>k=1,..., n/2-1,</tt> are
46 * the real and imaginary parts of the Fourier transform at a frequency \f$k\Delta f=k/T,\f$ \f$T\f$ being
47 * the duration of the signal and \f$\Delta f=1/T\f$ is the frequency resolution.
48 * </li></ul>
49 *
50 * ### Algorithm ###
51 *
52 * The standard SPA is given by \eqref{eq_InspiralFourierPhase_f2}.
53 * We define a variable function pointer \c LALInspiralTaylorF2Phasing and point
54 * it to one of the \c static functions defined within this function
55 * that explicitly calculates the Fourier phase at the PN order chosen by the user.
56 * The reference points are chosen so that on inverse Fourier transforming
57 * the time-domain waveform will
58 * <ul>
59 * <li> be padded with zeroes in the first <tt>params->nStartPad</tt> bins,
60 * </li><li> begin with a phase shift of <tt>params->nStartPhase</tt> radians,
61 * </li><li> have an amplitude of \f$n v^2.\f$
62 * </li></ul>
63 *
64 * ### Uses ###
65 *
66 * \code
67 * XLALInspiralSetup()
68 * XLALInspiralChooseModel()
69 * XLALInspiralTaylorF2Phasing[0234567]PN()
70 * \endcode
71 *
72 * ### Notes ###
73 *
74 * If it is required to compare the output of this module with a time domain
75 * signal one should use an inverse Fourier transform routine that packs data
76 * in the same way as fftw. Moreover, one should divide the resulting inverse
77 * Fourier transform by a factor \f$n/2\f$ to be consistent with the
78 * amplitude used in time-domain signal models.
79 *
80 */
81
82#include <math.h>
83#include "LALInspiral.h"
84
85#ifdef __GNUC__
86#define UNUSED __attribute__ ((unused))
87#else
88#define UNUSED
89#endif
90
98
99int
101 REAL4Vector *signalvec,
103 )
104{
105 REAL8 UNUSED h1, UNUSED h2, pimmc, f, v, df, shft, phi, amp0, amp, psif, psi;
106 INT4 n, nby2, i, f0, fn;
107 expnCoeffs ak;
108 expnFunc func;
109 REAL8 (*XLALInspiralTaylorF2Phasing)(REAL8, expnCoeffs *) = NULL;
110
111 /* Perform some initial checks */
112 if (signalvec == NULL)
114 if (signalvec->data == NULL)
116 if (params == NULL)
118 if (signalvec->length<=2)
120
121 /* chose the required phasing function */
122 switch (params->order)
123 {
125 case LAL_PNORDER_HALF:
126 XLALInspiralTaylorF2Phasing = XLALInspiralTaylorF2Phasing0PN;
127 break;
128 case LAL_PNORDER_ONE:
129 XLALInspiralTaylorF2Phasing = XLALInspiralTaylorF2Phasing2PN;
130 break;
132 XLALInspiralTaylorF2Phasing = XLALInspiralTaylorF2Phasing3PN;
133 break;
134 case LAL_PNORDER_TWO:
135 XLALInspiralTaylorF2Phasing = XLALInspiralTaylorF2Phasing4PN;
136 break;
138 XLALInspiralTaylorF2Phasing = XLALInspiralTaylorF2Phasing5PN;
139 break;
141 XLALInspiralTaylorF2Phasing = XLALInspiralTaylorF2Phasing6PN;
142 break;
144 XLALInspiralTaylorF2Phasing = XLALInspiralTaylorF2Phasing7PN;
145 break;
146 default:
148
149 }
150
151 /* Set up the coefficients in post-Newtonian expansion, vlso, etc. */
152 n = signalvec->length;
153 nby2 = n/2;
154 memset( &ak, 0, sizeof( ak ) );
157
158
159 /* Set up the functions required for the chosen signal
160 approximation scheme */
161 if ( XLALInspiralChooseModel(&func, &ak, params) == XLAL_FAILURE)
163
164 /* compute some basic variables */
165 df = params->tSampling/signalvec->length;
166 pimmc = LAL_PI * params->totalMass * LAL_MTSUN_SI;
167 f0 = params->fLower;
168 fn = (params->fCutoff < ak.fn) ? params->fCutoff : ak.fn;
169 v = pow(pimmc*f0, (1./3.));
170
171 /* If we want to pad with zeroes in the beginning then the instant of
172 * coalescence will be the chirp time + the duration for which padding
173 * is needed. Thus, in the equation below nStartPad occurs with a +ve sign.
174 * This code doesn't support non-zero start-time. i.e. params->startTime
175 * should be necessarily zero.
176 */
177 shft = 2.L*LAL_PI * (ak.tn + params->nStartPad/params->tSampling +
178 params->startTime);
179 phi = params->startPhase + LAL_PI/4.L;
180 amp0 = params->signalAmplitude * ak.totalmass * pow(LAL_PI/12.L, 0.5L) * df;
181
182 /*
183 Compute the standard stationary phase approximation.
184 */
185 h1 = signalvec->data[0] = 0.L;
186 h2 = signalvec->data[nby2] = 0.L;
187 for (i=1; i<nby2; i++) {
188 f = i * df;
189 if (f < f0 || f > fn)
190 {
191 /*
192 * All frequency components below f0 and above fn are set to zero
193 */
194 signalvec->data[i] = 0.;
195 signalvec->data[n-i] = 0.;
196 }
197 else
198 {
199 v = pow(pimmc * f, (1./3.));
200 psif = XLALInspiralTaylorF2Phasing(v, &ak);
201 psi = shft * f + phi + psif;
202 /*
203 dEnergybyFlux computes 1/(dv/dt) while what we need is 1/(dF/dt):
204 dF/dt=(dF/dv)(dv/dt)=-dEnergybyFlux/(dF/dv)=-dEnergybyFlux 3v^2/(pi m^2)
205 Note that our energy is defined as E=-eta v^2/2 and NOT -eta m v^2/2.
206 This is the reason why there is an extra m in the last equation above
207 amp = amp0 * pow(-dEnergybyFlux(v)/v^2, 0.5) * v^2;
208 = amp0 * pow(-dEnergybyFlux(v), 0.5) * v;
209 */
210 amp = amp0 * pow(-func.dEnergy(v,&ak)/func.flux(v,&ak),0.5L) * v;
211 signalvec->data[i] = (REAL4) (amp * cos(psi));
212 signalvec->data[n-i] = (REAL4) (-amp * sin(psi));
213
214 }
215
216 }
217 params->fFinal = fn;
218
219 return XLAL_SUCCESS;
220}
221
222
224{
225 return ak->pfaN/pow(v,5.);
226}
227
228
230{
231 REAL8 x;
232 x = v*v;
233 return ak->pfaN * (1. + ak->pfa2 * x)/pow(v,5.);
234}
235
236
238{
239 REAL8 x;
240 x = v*v;
241 return ak->pfaN * (1. + ak->pfa2 * x + ak->pfa3 * v*x)/pow(v,5.);
242}
243
244
246{
247 REAL8 x;
248 x = v*v;
249 return ak->pfaN * (1. + ak->pfa2 * x + ak->pfa3 * v*x + ak->pfa4 * x*x)/pow(v,5.);
250}
251
252
254{
255 REAL8 x, y;
256 x = v*v;
257 y = x*x;
258 return ak->pfaN * (1. + ak->pfa2 * x + ak->pfa3 * v*x + ak->pfa4 * y
259 + (ak->pfa5 + ak->pfl5 * log(v/ak->v0)) * y*v)/pow(v,5.);
260}
261
262
264{
265 REAL8 x, y, z;
266 x = v*v;
267 y = x*x;
268 z = y*x;
269
270 return ak->pfaN * (1. + ak->pfa2 * x + ak->pfa3 * v*x + ak->pfa4 * y
271 + (ak->pfa5 + ak->pfl5 * log(v/ak->v0)) * y*v
272 + (ak->pfa6 + ak->pfl6 * log(4.*v) ) * z)
273 /pow(v,5.);
274}
275
276
278
279 REAL8 x, y, z;
280 x = v*v;
281 y = x*x;
282 z = y*x;
283
284 return ak->pfaN * (1. + ak->pfa2 * x + ak->pfa3 * v*x + ak->pfa4 * y
285 + (ak->pfa5 + ak->pfl5 * log(v/ak->v0)) * y*v
286 + (ak->pfa6 + ak->pfl6 * log(4.*v) ) * z
287 + ak->pfa7 * z*v)
288 /pow(v,5.);
289}
int XLALInspiralSetup(expnCoeffs *ak, InspiralTemplate *params)
int XLALInspiralChooseModel(expnFunc *func, expnCoeffs *ak, InspiralTemplate *params)
static REAL8 XLALInspiralTaylorF2Phasing3PN(REAL8 v, expnCoeffs *ak)
static REAL8 XLALInspiralTaylorF2Phasing6PN(REAL8 v, expnCoeffs *ak)
static REAL8 XLALInspiralTaylorF2Phasing0PN(REAL8 v, expnCoeffs *ak)
int XLALInspiralStationaryPhaseApprox2(REAL4Vector *signalvec, InspiralTemplate *params)
static REAL8 XLALInspiralTaylorF2Phasing2PN(REAL8 v, expnCoeffs *ak)
static REAL8 XLALInspiralTaylorF2Phasing5PN(REAL8 v, expnCoeffs *ak)
static REAL8 XLALInspiralTaylorF2Phasing4PN(REAL8 v, expnCoeffs *ak)
static REAL8 XLALInspiralTaylorF2Phasing7PN(REAL8 v, expnCoeffs *ak)
double i
#define LAL_PI
#define LAL_MTSUN_SI
double REAL8
int32_t INT4
float REAL4
LAL_PNORDER_TWO_POINT_FIVE
LAL_PNORDER_THREE
LAL_PNORDER_TWO
LAL_PNORDER_ONE
LAL_PNORDER_THREE_POINT_FIVE
LAL_PNORDER_HALF
LAL_PNORDER_ONE_POINT_FIVE
LAL_PNORDER_NEWTONIAN
#define XLAL_ERROR(...)
XLAL_EBADLEN
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFUNC
XLAL_ETYPE
XLAL_FAILURE
list y
int n
x
The inspiral waveform parameter structure containing information about the waveform to be generated.
Definition: LALInspiral.h:205
REAL4 * data
This structure contains various post-Newtonian and P-approximant expansion coefficients; the meanings...
Definition: LALInspiral.h:399
REAL8 pfa7
Definition: LALInspiral.h:459
REAL8 totalmass
Definition: LALInspiral.h:474
REAL8 pfa5
Definition: LALInspiral.h:459
REAL8 pfl5
Definition: LALInspiral.h:459
REAL8 pfa4
Definition: LALInspiral.h:459
REAL8 pfa2
Definition: LALInspiral.h:459
REAL8 pfa3
Definition: LALInspiral.h:459
REAL8 pfl6
Definition: LALInspiral.h:459
REAL8 pfa6
Definition: LALInspiral.h:459
REAL8 pfaN
Definition: LALInspiral.h:459
Structure to hold the pointers to the generic functions defined above.
Definition: LALInspiral.h:546
EnergyFunction * dEnergy
Definition: LALInspiral.h:547
FluxFunction * flux
Definition: LALInspiral.h:548
double df