Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInspiral 5.0.3.1-6c6b863
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALInspiralComputeMetric.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2007 David Churches, Duncan Brown, Jolien Creighton, Peter Shawhan, B.S. Sathyaprakash, Anand Sengupta, 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 Created: 7.9.96.
22 Author: B.S.Sathyaprakash, Caltech, Cardiff University.
23 Purpose: To compute the metric and the template bank
24 parameters, corresponding to 2PN chirps.
25 Revision History: Updates 15.7.97; 31.8.97.; 18.3.99; 25.05.02
26 First C version October 2000.
27 Major revision in May 2002: Direct computation in (tau0, tau3)
28 space avoid (m, eta) space. This means the code won't create
29 template parameters in (tau0, tau2) space. An interface to the
30 old code will have to be provided, if that is necessary.
31
32 Dependencies: moments.f, inverse.f transform.f (not needed in the version of 02/05)
33 Outputs:
34 g00:
35 g11:
36 theta: Angle which the t0-axis makes with semi-major (dx0) axis.
37 srate: The minimal sampling rate required (computed but not outputted.
38 Notes: Owen and Sathyaprakash (Caltech collaboration notes).
39 Also Sathyaprakash, note from October 2000, May 24/25, 2002.
40*/
41
42#include <stdlib.h>
43#include <lal/LALInspiralBank.h>
44#include <lal/LALNoiseModels.h>
45
46#define METRIC_DIMENSION 2
47#define METRIC_ORDER 5
48
49static void
52 InspiralMomentsEtc *moments,
53 REAL8 fLower,
54 REAL8 t0,
55 REAL8 t3
56 );
57
58/**
59 * \ingroup LALInspiralBank_h
60 * \deprecated Use XLALInspiralComputeMetric() instead.
61 */
62void
64 LALStatus *status, /**< LAL status pointer */
65 InspiralMetric *metric, /**< [out] the metric at the lattice point defined by <tt>params</tt> */
66 InspiralTemplate *params, /**< [in] the parameters where metric must be computed */
67 InspiralMomentsEtc *moments /**< [in] moments \f$J(1), \ldots, J(17),\f$ of the PSD and other constants needed in the computation of the metric */
68 )
69{
71
72 XLAL_PRINT_DEPRECATION_WARNING("XLALInspiralComputeMetric");
73
74 if (XLALInspiralComputeMetric(metric, moments, params->fLower, params->order, params->t0, params->t3) == XLAL_FAILURE){
76 };
77
78 RETURN( status );
79}
80
81
82/**
83 * \ingroup LALInspiralBank_h
84 * \brief Function to compute the components of the metric which is used to describe distances on the signal manifold.
85 * \author Churches, D. K., Cokelaer, T., Sathyaprakash, B. S.
86 *
87 * We calculate the components of the metric using the procedure outlined
88 * in Owen \cite Owen_96 .
89 * This uses the moments of the noise curve,
90 * \f{equation}{
91 * I(q) \equiv S_{h}(f_{0}) \int^{f_{c}/f_{0}}_{f_{s}/f_{0}} \frac{x^{-q/3}}{S_{h}(x)}
92 * \, dx
93 * \f}
94 * and
95 * \f{equation}{
96 * J(q) \equiv \frac{I(q)}{I(7)} \,.
97 * \f}
98 * (Please note that the function \c LALInspiralMoments doesn't compute \f$I(q)\f$ defined
99 * here; the index \f$q\f$ is definted differently there and the normalisation is supplied by the user.
100 * For ease of writing, here we shall follow the standard notation.)
101 * Then the moment functional \f$\mathcal{J}\f$ is defined such that, for a function \f$a\f$,
102 * \f{equation}{
103 * \mathcal{J} [a] \equiv \frac{1}{I(7)} \int^{f_{c}/f_{0}}_{f_{s}/f_{0}} \frac{x^{-7/3}}{S_{h}(x)}
104 * a(x) \, dx
105 * \f}
106 * which gives us
107 * \f{equation}{
108 * \mathcal{J} = \left[ \sum_{n} a_{n} x^{n} \right] = \sum_{n} a_{n} J(7-3n).
109 * \f}
110 * The above equation is used to calculate the components of the metric using the following formula:
111 * \f{equation}{
112 * \gamma_{\alpha \beta} = \frac{1}{2} \left( \mathcal{J} [ \psi_{\alpha} \psi_{\beta} ] -
113 * \mathcal{J} [ \psi_{\alpha} ] \mathcal{J} [ \psi_{\beta} ] \right)
114 * \f}
115 * where \f$\psi_\alpha\f$ is the derivative of the Fourier phase of the inspiral waveform
116 * with respect to the parameter \f$\lambda^\alpha,\f$ that is
117 * \f$\psi_\alpha \equiv \Psi_{,\alpha}.\f$ Writing the derivative index as
118 * \f$\alpha=0,j,\f$ with \f$j=1,\ldots,n\f$ we have
119 * \f{equation}{
120 * \psi_{0} \equiv 2 \pi f, \ \ \
121 * \psi_{j} \equiv \frac{\partial \Delta \Psi}{\partial \Delta \lambda^{j}}.
122 * \f}
123 * The phase \f$\Psi\f$ is that which appears in the usual stationary-phase formula for
124 * the Fourier transform:
125 * \f{equation}{
126 * \tilde{h}(f) \propto f^{-7/6} e^{i[-\pi/4 - \Phi_{0} + 2 \pi f t_{0} + \Psi(f;\vec{\lambda}]}.
127 * \f}
128 * If we take the usual chirp times and multiply each by \f$(\pi f_{0})\f$ then we get
129 * <em>dimensionless
130 * chirp times</em> \f$\tau_{k}\f$, and then the phase \f$\Psi\f$ may be written in the form
131 * \f{equation}{
132 * \Psi = 2 \pi f t_{c} + \sum_{k} \Psi_{k}(f) \tau_{k}
133 * \f}
134 * where, defining \f$v_0 = (\pi m f_0)^{1/3}\f$ (\f$m\f$ being total mass and \f$f_0\f$ a
135 * fiducial starting frequency), the chirptimes \f$\tau_{k},\f$ up to 2nd PN order,
136 * are given by
137 * \f[
138 * \tau_{0} = \frac{5}{256 \eta v_{0}^{5}},\ \
139 * \tau_{2} = \frac{5}{192 \eta v_{0}^{3}} \left( \frac{743}{336} + \frac{11}{4} \eta \right),
140 * \f]
141 * \f{equation}{
142 * \tau_{3} = \frac{\pi}{8 \eta v_{0}^{2}},\ \
143 * \tau_{4} = \frac{5}{128 \eta v_{0}} \left( \frac{3\,058\,673}{1\,016\,064} + \frac{5429}{1008}
144 * \eta +
145 * \frac{617}{144} \eta^{2} \right).
146 * \f}
147 * Up to second post-Newtonian approximation the \f$\Psi_{k}\f$ are given by
148 * \f{equation}{
149 * \Psi_{0} = \frac{6}{5 \nu^{5/3}},\ \
150 * \Psi_{2} = \frac{2}{\nu},\ \
151 * \Psi_{3} = - \frac{3}{\nu^{2/3}},\ \
152 * \Psi_{4} = \frac{6}{\nu^{1/3}}.
153 * \f}
154 * where \f$\nu = f/f_{0}\f$.
155 *
156 * If we now make the substitution \f$f = v^{3}/\pi m\f$ we then the find that the phase may be
157 * expressed in a simpler form
158 * \f{equation}{
159 * \Psi(v) = 2 \pi f t_{c} + \sum_{k} \theta_{k} v^{k-5}
160 * \f}
161 * where the <em>chirp parameters</em> \f$\theta_{k}\f$ are given by
162 * \f[
163 * \theta_{0} = \frac{3}{128 \eta}, \ \
164 * \theta_{2} = \frac{5}{96 \eta} \left( \frac{743}{336} + \frac{11}{4} \eta \right),
165 * \f]
166 * \f{equation}{
167 * \theta_{3} = - \frac{3 \pi}{8 \eta},\ \
168 * \theta_{4} = \frac{15}{64 \eta} \left( \frac{3\,058\,673}{1\,016\,064} + \frac{5429}{1008} \eta
169 * + \frac{617}{144}
170 * \eta^{2} \right).
171 * \f}
172 *
173 * If we want to express \f$\Psi\f$ in terms of \f$f\f$ rather than \f$v\f$ we simply substitute \f$v = (\pi m
174 * f)^{1/3}\f$
175 * to obtain
176 * \f{equation}{
177 * \label{phaselabel}
178 * \Psi(f) = 2 \pi f t_{c} + \sum_{k} \theta^{\prime}_{k} f^{(k-5)/3}
179 * \f}
180 * where
181 * \f{equation}{
182 * \theta^{\prime}_{k} = (\pi m)^{(k-5)/3} \theta_{k}.
183 * \f}
184 *
185 * We are now in a position to start calculating components of \f$\gamma_{\alpha \beta}\f$. We had
186 * \f{equation}{
187 * \psi_{j} \equiv \frac{\partial \Delta \Psi}{\partial \Delta \lambda^{j}}
188 * \f}
189 * where \f$\Psi\f$ is given by \eqref{phaselabel}. Therefore we may write
190 * \f{equation}{
191 * \Delta \Psi = \Delta \theta^{\prime}_{0} f^{-5/3} + \Delta \theta^{\prime}_{2} f^{-1} + \Delta
192 * \theta^{\prime}_{3} f^{-2/3} + \Delta \theta^{\prime}_{4} f^{-1/3}
193 * \f}
194 * All we need to do now is specify the coordinates \f$\lambda^{j}\f$ with respect to which the
195 * derivatives
196 * will be taken. In general, the template placement algorithm works in \f$(\tau_{0},\tau_{3})\f$
197 * coordinates. It is simplest for us to calculate the components of \f$\gamma_{\alpha \beta}\f$ in
198 * the \f$(m,\eta)\f$ coordinate system and then perform a coordinate transformation to get the
199 * components in
200 * the \f$(\tau_{0},\tau_{3})\f$ system.
201 * So, we first of all calculate the components of \f$\gamma_{\alpha \beta}\f$ in the \f$(m,\eta)\f$
202 * system.
203 *
204 * This involves calculating the following:
205 * \f{equation}{
206 * \frac{\partial \Delta \Psi}{\partial \Delta m} = \frac{\Delta \theta^{\prime}_{0}}{\Delta m}
207 * f^{-5/3} +
208 * \frac{\Delta \theta^{\prime}_{2}}{\Delta m} f^{-1} - \frac{\Delta \theta^{\prime}_{3}}{\Delta m}
209 * f^{-2/3} + \frac{\delta \theta^{\prime}_{4}}{\Delta m} f^{-1/3}
210 * \f}
211 * and
212 * \f{equation}{
213 * \frac{\partial \Delta \Psi}{\partial \Delta \eta} = \frac{\Delta \theta^{\prime}_{0}}{\Delta
214 * \eta} f^{-5/3} +
215 * \frac{\Delta \theta^{\prime}_{2}}{\Delta \eta} f^{-1} - \frac{\Delta \theta^{\prime}_{3}}{\Delta
216 * \eta}
217 * f^{-2/3} + \frac{\delta \theta^{\prime}_{4}}{\Delta \eta} f^{-1/3}
218 * \f}
219 * where all of the derivatives are easily calculable. This gives us the terms \f$\psi_{j}\f$ as a
220 * power
221 * series in \f$f\f$. These are then used in the formula
222 * \f{equation}{
223 * \gamma_{\alpha \beta} = \frac{1}{2} \left( \mathcal{J} [ \psi_{\alpha} \psi_{\beta} ] -
224 * \mathcal{J} [
225 * \psi_{\alpha}] \mathcal{J} [\psi_{\beta}] \right)
226 * \f}
227 * to calculate the components of \f$\gamma_{\alpha \beta}\f$. The fact that each of the \f$\psi_{j}\f$ is
228 * in the
229 * form of a power series in \f$f\f$ allows us to calculate \f$\gamma_{\alpha \beta}\f$ using
230 * \f{equation}{
231 * \mathcal{J} \left[ \sum_{n} a_{n} x^{n} \right] = \sum_{n} J(7-3n).
232 * \f}
233 * i.e.\ we can express all the \f$\mathcal{J}[]\f$ in terms of the integral \f$J(q)\f$ which we calculate
234 * numerically at the outset for the required values of \f$q\f$.
235 *
236 * Once we have obtained \f$\gamma_{\alpha \beta}\f$ in this way, we take the inverse of this matrix to
237 * give us \f$\gamma^{\alpha \beta}\f$ in the \f$(m,\eta)\f$ system. Then
238 * we perform the following coordinate transformation to give us the components of
239 * \f$\gamma^{\alpha^{\prime}
240 * \beta^{\prime}}\f$ in our chosen system,
241 * \f{equation}{
242 * \gamma^{\alpha^{\prime} \beta^{\prime}} = \Lambda^{\alpha^{\prime}}_{\,\,\sigma}
243 * \Lambda^{\beta^{\prime}}_{\,\,\delta} \gamma^{\sigma \delta}
244 * \f}
245 * where the transformation matrix \f$\Lambda^{\alpha^{\prime}}_{\,\,\beta}\f$ is defined by
246 * \f{equation}{
247 * \Lambda^{\alpha^{\prime}}_{\,\,\beta} = \frac{\partial x^{\alpha^{\prime}}}{\partial x^{\beta}}
248 * \f}
249 * Finally, we take the inverse of this matrix to obtain \f$\gamma_{\alpha^{\prime} \beta^{\prime}}\f$
250 * in the
251 * chosen system.
252 * Since the unprimed system corresponds to \f$(t_{c},m,\eta)\f$ coordinates and the primed system to
253 * \f$(t_{c},\tau_{0},\tau_{3})\f$ coordinates, the matrix
254 * \f$\Lambda^{\alpha^{\prime}}_{\,\,\beta^{\prime}}\f$ has
255 * element
256 * \f{equation}{
257 * \Lambda^{\alpha^{\prime}}_{\,\,\beta} = \left(
258 * \begin{array}{ccc}
259 * 1 & 0 & 0 \\
260 * 0 & \frac{\partial \tau_{0}}{\partial m} & \frac{\partial \tau_{0}}{\partial \eta} \\
261 * 0 & \frac{\partial \tau_{3}}{\partial m} & \frac{\partial \tau_{3}}{\partial \eta}
262 * \end{array}
263 * \right) = \left(
264 * \begin{array}{ccc}
265 * 1 & 0 & 0 \\
266 * 0 & -\frac{5 \tau_{0}}{3m} & -\frac{\tau_{0}}{\eta} \\
267 * 0 & -\frac{2 \tau_{3}}{3m} & -\frac{\tau_{3}}{\eta}
268 * \end{array}
269 * \right)
270 * \f}
271 *
272 * Finally, what is needed in laying a lattice in the space of dynamical parameters
273 * (also referred to as intrinsic parameters) is the metric with the kinematical
274 * parameter (also called extrinsic parameter) projected out: In other words one defines
275 * the 2-dimensional metric \f$g_{mn}\f$ by
276 * \f{equation}{
277 * g_{mn} = \gamma_{mn} - \frac{\gamma_{0m} \gamma_{0n}}{\gamma_{00}}.
278 * \f}
279 *
280 * ### Metric computation in the \f$\tau_0-\tau_3\f$ space ###
281 *
282 * The metric cannot be directly computed in the \f$(\tau_0,\tau_2)\f$ space.
283 * Therefore, in the previous Section we first computed the metric
284 * in the \f$(m,\eta)\f$ space and then transformed to \f$(\tau_0,\tau_2)\f$ space.
285 * The same method can also be used to find the metric in the \f$(\tau_0,\tau_3)\f$ space.
286 * However, in \f$(\tau_0,\tau_3)\f$ space one can directly compute the
287 * metric without recourse to \f$(m,\eta)\f$ coordinates. It is of interest to see
288 * whether this yields the same results as the previous method.
289 *
290 * The starting point of our derivation is Eq.\ (3.7) of Owen and Sathyaprakash
291 * (Phys. Rev. D 60, 022002, 1999) for the Fourier domain phase which we shall
292 * write as:
293 * \f{eqnarray}{
294 * \Psi(f; \theta_1, \theta_2) & = & a_{01}\theta_1 v^{-5}
295 * + \left [a_{21} \frac {\theta_1}{\theta_2} + a_{22} \left ( \theta_1 \theta_2^2 \right )^{1/3} \right ] v^{-3}
296 * + a_{31} \theta_2 v^{-2} \\
297 * & + & \left [a_{41} \frac {\theta_1}{\theta_2^2} + a_{42} \left ( \frac {\theta_1}{\theta_2} \right )^{1/3}
298 * + a_{43} \left ( \frac{\theta_2^4}{\theta_1} \right )^{1/3} \right ] v^{-1},
299 * \f}
300 * to 2nd post-Newtonain order. Here \f$v=(f/f_0)^{1/3},\f$ \f$\theta_1\f$ and \f$\theta_2\f$ are
301 * identical to the \f$\theta^1\f$ and \f$\theta^2\f$ parameters
302 * of Owen and Sathyaprakash defined in Eq.\ (3.3) there and the \f$a\f$ coefficients are given by:
303 * \f{eqnarray}{
304 * a_{01} = \frac{3}{5}, \ \ a_{21} = \frac{11\pi}{12}, \ \
305 * a_{22} = \frac{743}{2016} \left ( \frac {25}{2\pi^2} \right )^{1/3}, \ \ a_{31} = -\frac{3}{2}, \\
306 * a_{41} = \frac {617}{384} \pi^2, \ \ a_{42} = \frac{5429}{5376} \left ( \frac{25 \pi}{2} \right )^{1/3},\ \
307 * a_{43} = \frac {15293365}{10838016} \left ( \frac{5}{4\pi^4} \right )^{1/3}.
308 * \f}
309 * Differentials of the phase with respect to the coordinates \f$\theta_1\f$ and \f$\theta_2\f$ appear in the
310 * metric which we write as:
311 * \f{equation}{
312 * \psi_m \equiv \frac{\partial \Psi}{\partial \theta_m} = \sum_0^N \Psi_{mk} v^{k-5}.
313 * \f}
314 * where \f$N\f$ is the post-Newtonian order up to which the phase is known, or the post-Newtonian
315 * order at which the metric is desired.
316 * Expansion coefficients \f$\Psi_{mn}\f$ can be considered be \f$(2\times N)\f$ matrix which to
317 * second post-Newtonian order is given by:
318 * \f{equation}{
319 * \Psi =
320 * \left [ \begin{matrix}
321 * a_{01}
322 * & 0
323 * & {a_{21}}/{\theta_2} + ({a_{22}}/{3}) \left ( {\theta_2}/{\theta_1} \right )^{2/3}
324 * & 0
325 * & {a_{41}}/{\theta_2^2} + {a_{42}}/\left ({3 \left ( \theta_1^2\theta_2 \right )^{1/3} } \right )
326 * - ({a_{43}}/{3}) \left ( {\theta_2}/{\theta_1} \right )^{4/3} \cr
327 * 0
328 * & 0
329 * & - {a_{21}\theta_1}/{\theta_2^2} + (2 {a_{22}}/{3}) \left ( {\theta_1}/{\theta_2} \right )^{1/3}
330 * & a_{31}
331 * & - {2a_{41} \theta_1}/{\theta_2^3} - ({a_{42}}/{3}) \left ( {\theta_1}/{\theta_2^4} \right )^{1/3}
332 * + ({4a_{43}}/{3}) \left ( {\theta_2}/{\theta_1} \right )^{1/3}
333 * \end{matrix}
334 * \right ].
335 * \f}
336 *
337 * Using the definition of the
338 * metric introduced earlier and projecting out the \f$t_c\f$ coordinate, one finds that
339 * \f{eqnarray}{
340 * g_{mn} & = & \frac{1}{2}\sum_{k,l=0}^N \Psi_{mk} \Psi_{nl}
341 * \biggl [ J(17-k-l) - J(12-k) J(12-l) \biggr . \\
342 * & - & \biggl . \frac { \left ( J(9-k) - J(4)J(12-k) \right )
343 * \left ( J(9-l) - J(4)J(12-l) \right )} {\left (J(1) - J(4)^2 \right)}
344 * \biggr ]
345 * \f}
346 * where \f$J\f$'s are the moments introduced earlier.
347 */
348int
350 InspiralMetric *metric,
351 InspiralMomentsEtc *moments,
352 REAL8 fLower,
353 LALPNOrder order,
354 REAL8 t0,
355 REAL8 t3
356 )
357
358{
361
362 REAL8 a, b, c, q;
363 UINT4 PNorder, m, n;
364
365 const REAL8 two_pi_flower = LAL_TWOPI * fLower;
366 const REAL8 two_pi_flower_sq = two_pi_flower * two_pi_flower;
367
368 /* check inputs */
369 if (!metric || !moments){
370 XLALPrintError(LALINSPIRALH_MSGENULL);
372 };
373
374 if (t0 <= 0.L || t3 <= 0.L){
375 XLALPrintError("t0 and t3 must be positive");
377 };
378
379 if (fLower <= 0){
380 XLALPrintError("fLower must be positive");
382 }
383
384 /* use the order of the waveform to compute the metric */
385 /* summation below will be carried out up to PNorder */
386 if ( order != LAL_PNORDER_ONE &&
388 order != LAL_PNORDER_TWO )
389 {
390 /* Let us force the order to be twoPN because that is the only order
391 * available for the template bank anyway. */
392 XLALPrintWarning("forcing LAL_PNORDER_TWO");
393 PNorder = LAL_PNORDER_TWO;
394 }
395 else
396 {
397 PNorder = (UINT4) order;
398 }
399
400
401 /* Setting up \Psi_{mn} coefficients */
402 InspiralComputeMetricGetPsiCoefficients( Psi, moments, fLower, t0, t3 );
403
404 for ( m = 0; m < METRIC_DIMENSION; m++ )
405 {
406 for ( n = m; n < METRIC_DIMENSION; n++ )
407 {
408 UINT4 k, l;
409 g[m][n] = 0.L;
410
411 for ( k = 0 ; k <= PNorder; k++ )
412 {
413 for ( l = 0; l <= PNorder; l++ )
414 {
415 g[m][n] += Psi[m][k] * Psi[n][l] * (
416 moments->j[17-k-l] - moments->j[12-k] * moments->j[12-l]
417 - ( moments->j[9-k] - moments->j[4] * moments->j[12-k] )
418 * ( moments->j[9-l] - moments->j[4] * moments->j[12-l] )
419 / ( moments->j[1] - moments->j[4] * moments->j[4] )
420 );
421 }
422 }
423 g[m][n] /= 2.;
424 g[n][m] = g[m][n];
425 }
426 }
427
428#if 0
429 The minimum sampling rate for given MM is
430 srate =
431 2 * LAL_PI * f0 sqrt( (moments.j[1] - moments.j[4]*moments.j[4]) /
432 (2.0*(1.-MM)));
433#endif
434
435 /* The calculation above gives the metric in coordinates */
436 /* (t0=2\pi f_0 \tau0, t3=2\pi f_0 \tau3). Re-scale metric */
437 /* coefficients to get metric in (tau0, tau3) coordinates */
438 a = g[0][0] * two_pi_flower_sq;
439 b = g[0][1] * two_pi_flower_sq;
440 c = g[1][1] * two_pi_flower_sq;
441
442
443 /* The metric in tau0-tau2,3 space. */
444 metric->G00 = a;
445 metric->G01 = b;
446 metric->G11 = c;
447
448 /* Diagonalize the metric. */
449 q = sqrt( (a-c)*(a-c) + 4. * b*b );
450
451 metric->g00 = 0.5 * (a + c - q);
452 metric->g11 = 0.5 * (a + c + q);
453
454 if ( a == c )
455 {
456 metric->theta = LAL_PI/2.;
457 }
458 else
459 {
460 /* metric->theta = 0.5 * atan(2.*b/(a-c)); */
461 /* We want to always measure the angle from the */
462 /* semi-major axis to the tau0 axis which is given by */
463 /* the following line as opposed to the line above */
464 metric->theta = atan( b / (metric->g00 - c) );
465 }
466
467 /* memset the metric->Gamma array to zero before populating them with correct
468 * values. This prevents junk getting stored in unused fields */
469 memset (metric->Gamma, 0, 10*sizeof(REAL4));
470
471 /* Now we compute the 3d metric in tc,\tau_0,\tau_3 co-ordinates */
472 /* We only need metric->Gamma[0,...,5]. */
473 metric->Gamma[0] = 0.5*two_pi_flower_sq*
474 ( moments->j[1] - (moments->j[4]*moments->j[4]) );
475
476 metric->Gamma[1] = 0.5*two_pi_flower_sq*
477 ( Psi[0][0]*(moments->j[9] - (moments->j[4]*moments->j[12]) )
478 + Psi[0][2]*(moments->j[7] - (moments->j[4]*moments->j[10]) )
479 + Psi[0][4]*(moments->j[5] - (moments->j[4]*moments->j[8]) ));
480
481 metric->Gamma[2] = 0.5*two_pi_flower_sq*
482 ( Psi[1][2]*(moments->j[7] - (moments->j[4]*moments->j[10]) )
483 + Psi[1][3]*(moments->j[6] - (moments->j[4]*moments->j[9]) )
484 + Psi[1][4]*(moments->j[5] - (moments->j[4]*moments->j[8]) ));
485
486
487 metric->Gamma[3] = metric->G00 + metric->Gamma[1]*metric->Gamma[1]/metric->Gamma[0];
488 metric->Gamma[4] = metric->G01 + metric->Gamma[1]*metric->Gamma[2]/metric->Gamma[0];
489 metric->Gamma[5] = metric->G11 + metric->Gamma[2]*metric->Gamma[2]/metric->Gamma[0];
490
491
492 return XLAL_SUCCESS;
493}
494
495static void
498 InspiralMomentsEtc *moments,
499 REAL8 fLower,
500 REAL8 t0,
501 REAL8 t3
502 )
503{
504 REAL8 t1 = LAL_TWOPI * fLower * t0;
505 REAL8 t2 = LAL_TWOPI * fLower * t3;
506
507 Psi[0][0] = moments->a01;
508 Psi[0][1] = 0.L;
509 Psi[0][2] = moments->a21/t2 + moments->a22/3.L * cbrt(t2 * t2 / (t1 * t1));
510 Psi[0][3] = 0.L;
511 Psi[0][4] = moments->a41/(t2*t2) + moments->a42/(3.L* cbrt(t1*t1*t2))
512 - moments->a43/3.L * t2 / t1 * cbrt(t2 / t1);
513
514 Psi[1][0] = 0.L;
515 Psi[1][1] = 0.L;
516 Psi[1][2] = -moments->a21*t1/(t2*t2) + 2.L *
517 moments->a22/3.L * cbrt(t1/t2);
518 Psi[1][3] = moments->a31;
519 Psi[1][4] = - 2.L * moments->a41*t1 / (t2*t2*t2) -
520 moments->a42/3.L * cbrt(t1/(t2*t2*t2*t2)) +
521 4.L * moments->a43/3.L * cbrt(t2/t1);
522}
523
524/**
525 * UNDOCUMENTED
526 * \see See LALInspiralComputeMetric() for documentation
527 */
528void
531 InspiralMetric *metric,
534 )
535{
537 InspiralMomentsEtcBCV moments;
538 REAL8 num;
539 REAL8 a, b, c, q;
540
543
544 moments.alpha = params->alpha;
545 moments.n0 = 5.L/3.L;
546 moments.n15 = 2.L/3.L;
547
548 LALGetInspiralMomentsBCV( status->statusPtr, &moments, psd, params );
550
551 num = moments.M3[0][0] *moments.M3[1][1]
552 - moments.M3[0][1] * moments.M3[1][0];
553
554 g[0][0] =moments.M2[0][0]*(moments.M3[1][1]*moments.M2[0][0]
555 -moments.M3[0][1]*moments.M2[0][1])
556 +moments.M2[0][1]*(-moments.M3[0][1]*moments.M2[0][0]
557 +moments.M3[0][0]*moments.M2[0][1]);
558 g[0][0] /= num;
559
560
561 g[1][1] =moments.M2[0][1]*(moments.M3[1][1]*moments.M2[0][1]
562 -moments.M3[0][1]*moments.M2[1][1])
563 +moments.M2[1][1]*(-moments.M3[0][1]*moments.M2[0][1]
564 +moments.M3[0][0]*moments.M2[1][1]);
565 g[1][1] /= num;
566
567
568 g[0][1] = moments.M2[0][0]*(moments.M3[1][1]*moments.M2[0][1]
569 -moments.M3[0][1]*moments.M2[1][1])
570 +moments.M2[0][1]*(-moments.M3[0][1]*moments.M2[0][1]
571 +moments.M3[0][0]*moments.M2[1][1]);
572 g[0][1] /= num ;
573
574 metric->G00 = .5 *(moments.M1[0][0] - g[0][0] );
575 metric->G01 = .5 *(moments.M1[0][1] - g[0][1] );
576 metric->G11 = .5 *(moments.M1[1][1] - g[1][1] );
577
578 a = metric->G00;
579 b = metric->G01;
580 c = metric->G11;
581 q = sqrt( (a-c)*(a-c) + 4. * b*b );
582 metric->g00 = 0.5 * (a + c - q);
583 metric->g11 = 0.5 * (a + c + q);
584 if ( a == c )
585 {
586 metric->theta = LAL_PI/2.;
587 }
588 else
589 {
590 /* metric->theta = 0.5 * atan(2.*b/(a-c)); */
591 /* We want to always measure the angle from the */
592 /* semi-major axis to the tau0 axis which is given by */
593 /* the following line as opposed to the line above */
594 metric->theta = atan(b/(metric->g00 - c));
595 }
596
598 RETURN( status );
599}
600
601#undef METRIC_ORDER
602#undef METRIC_DIMENSION
static void InspiralComputeMetricGetPsiCoefficients(REAL8 Psi[METRIC_DIMENSION][METRIC_ORDER], InspiralMomentsEtc *moments, REAL8 fLower, REAL8 t0, REAL8 t3)
#define METRIC_ORDER
#define METRIC_DIMENSION
#define CHECKSTATUSPTR(statusptr)
#define ATTATCHSTATUSPTR(statusptr)
#define DETATCHSTATUSPTR(statusptr)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
#define ABORTXLAL(sp)
int l
double srate
#define LAL_PI
#define LAL_TWOPI
double REAL8
uint32_t UINT4
float REAL4
void LALInspiralComputeMetricBCV(LALStatus *status, InspiralMetric *metric, REAL8FrequencySeries *psd, InspiralTemplate *params)
UNDOCUMENTED.
int XLALInspiralComputeMetric(InspiralMetric *metric, InspiralMomentsEtc *moments, REAL8 fLower, LALPNOrder order, REAL8 t0, REAL8 t3)
Function to compute the components of the metric which is used to describe distances on the signal ma...
void LALGetInspiralMomentsBCV(LALStatus *status, InspiralMomentsEtcBCV *moments, REAL8FrequencySeries *psd, InspiralTemplate *params)
void LALInspiralComputeMetric(LALStatus *status, InspiralMetric *metric, InspiralTemplate *params, InspiralMomentsEtc *moments)
LALPNOrder
LAL_PNORDER_TWO
LAL_PNORDER_ONE
LAL_PNORDER_ONE_POINT_FIVE
static const INT4 m
static const INT4 q
static const INT4 a
#define XLAL_ERROR(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
int int XLALPrintWarning(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_PRINT_DEPRECATION_WARNING(replacement)
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EDOM
XLAL_FAILURE
c
int n
double t0
Structure to store metric at various points the signal manifold.
REAL4 Gamma[10]
3d metric co-efficients in coordinates; Gamma[6] is a vector that stores the upper triangular part o...
REAL8 g11
11-component of the diagonalised metric
REAL8 G01
01-component of the metric in coordinates
REAL8 g00
00-component of the diagonalised metric
REAL8 theta
Angle from tau0 to semi-major axis of the ellipse.
REAL8 G00
00-component of the metric in coordinates
REAL8 G11
11-component of the metric in coordinates
Parameter structure that holds the moments of the PSD and other useful constants required in the comp...
REAL8 j[18]
The required moments are all computed once and stored in this array; The required moments are from J(...
The inspiral waveform parameter structure containing information about the waveform to be generated.
Definition: LALInspiral.h:205