LALSimulation  5.4.0.1-fe68b98
LALSimInspiralTaylorF2ReducedSpinMetric.c
Go to the documentation of this file.
1 /*
2 * Copyright (C) 2011 P. Ajith
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with with program; see the file COPYING. If not, write to the
16 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17 * MA 02110-1301 USA
18 */
19 
20 #include <lal/AVFactories.h>
21 #include <stdlib.h>
22 #include <math.h>
23 #include <lal/LALConstants.h>
24 #include <lal/LALDatatypes.h>
25 #include <lal/LALSimInspiral.h>
26 #include <lal/Units.h>
27 #include <lal/XLALError.h>
28 #include <gsl/gsl_blas.h>
29 #include <gsl/gsl_linalg.h>
30 #include <gsl/gsl_permutation.h>
31 
32 /* Computed in Mathematica with N[\[Pi]^2, 35] */
33 #define Pi_p2 9.8696044010893586188344909998761511
34 #define Pi_p3 31.006276680299820175476315067101395
35 #define Pi_p4 97.409091034002437236440332688705111
36 #define Pi_p5 306.01968478528145326274131004343561
37 #define Pi_p4by3 4.6011511144704899609836928301589280
38 #define Pi_p5by3 6.7388085956981412852178979102191822
39 #define Pi_p10by3 45.411541289455155012188835024560824
40 #define Pi_p7by3 14.454942539276981063196177096150014
41 #define Pi_p1by3 1.4645918875615232630201425272637904
42 #define Pi_p2by3 2.1450293971110256000774441009412356
43 #define Pi_11by3 66.509374974201175524807943232081902
44 #define tenbyPi_p2by3 2.1638812222639821850478371484217674
45 #define twoPi_p2by3 3.4050219214767547368607032470366207
46 
47 /*
48  * Function to compute the metric elements using waveform derivatives
49  */
50 static REAL8 MetricCoeffs(REAL8Vector *A, REAL8Vector *dPsii, REAL8Vector *dPsij,
51  REAL8Vector *dAi, REAL8Vector*dAj, REAL8Vector *Sh, REAL8 hSqr, REAL8 df) {
52  size_t k = A->length;
53  REAL8 gij = 0.;
54  for (;k--;) {
55  gij += (A->data[k]*A->data[k]*dPsii->data[k]*dPsij->data[k]
56  + dAi->data[k]*dAj->data[k])/Sh->data[k];
57  }
58  return 4.*df*gij/(2.*hSqr);
59 }
60 
61 /*
62  * Frequency domain amplitude of the TaylorF2 Reduced Spin waveforms
63  */
64 static REAL8 XLALSimInspiralTaylorF2RedSpinAofF(
65  const REAL8 mc_msun, /**< chirp mass (solar masses) */
66  const REAL8 eta, /**< symmetric mass ratio */
67  const REAL8 chi, /**< reduced-spin parameter */
68  const REAL8 f /**< Fourier frequency (Hz) */
69 ) {
70  const REAL8 mc = mc_msun * LAL_MTSUN_SI;
71  const REAL8 etap2 = eta*eta;
72  const REAL8 etap3 = etap2*eta;
73  const REAL8 chip2 = chi*chi;
74  const REAL8 mcp2 = mc*mc;
75 
76  REAL8 eta_three_fifths = pow(eta, 0.6);
77  REAL8 f_mc_eta = (f*mc)/eta_three_fifths;
78  REAL8 f_mc_eta_one_third = cbrt(f_mc_eta);
79 
80  return (sqrt(0.8333333333333334)*pow(mc/eta_three_fifths,0.8333333333333334)*
81  sqrt(eta)*(1 + (Pi_p4by3*
82  f_mc_eta_one_third * f_mc_eta_one_third *
83  (743 + 924*eta))/672. -
84  (Pi_p10by3 *
85  f_mc_eta * f_mc_eta * f_mc_eta_one_third*
86  (5111593 + 8088752*eta + 151088*etap2))/2.709504e6 +
87  (f*mc*LAL_PI*(-2*LAL_PI + (113*chi)/24.))/eta_three_fifths +
88  (Pi_p5by3 *
89  f_mc_eta * f_mc_eta / f_mc_eta_one_third *
90  (12*LAL_PI*(-4757 + 4788*eta) -
91  (113*(502429 - 591368*eta + 1680*etap2)*chi)/
92  (-113 + 76*eta)))/16128. +
93  (Pi_p5 * Pi_p1by3 *
94  f_mc_eta * f_mc_eta_one_third *
95  (7266251 + 9532152*eta + 9730224*etap2 +
96  (3243530304*(-81 + 4*eta)*chip2)/
97  pow(113 - 76*eta,2)))/8.128512e6 +
98  (f*f*mcp2*Pi_p2*
99  (-58.601030974347324 + (10*Pi_p2)/3. +
100  (3526813753*eta)/2.7869184e7 -
101  (451*Pi_p2*eta)/96. -
102  (1041557*etap2)/258048. +
103  (67999*etap3)/82944. +
104  (856*(3*LAL_GAMMA + log((64*f*mc*LAL_PI)/eta_three_fifths)))/315.))/(eta_three_fifths * eta_three_fifths)))/
105  (2.*pow(f,1.1666666666666667)*Pi_p4by3);
106 }
107 
108 /*
109  * Derivative of the amplitude with respect to \f$\chi\f$ (reduced-spin parameter)
110  */
111 static REAL8 XLALSimInspiralTaylorF2RedSpinDerivAChi(
112  const REAL8 mc_msun, /**< chirp mass (solar masses) */
113  const REAL8 eta, /**< symmetric mass ratio */
114  const REAL8 chi, /**< reduced-spin parameter */
115  const REAL8 f /**< Fourier frequency (Hz) */
116 ) {
117  const REAL8 mc = mc_msun * LAL_MTSUN_SI;
118  const REAL8 etap2 = eta*eta;
119  REAL8 etap35 = pow(eta, 0.6);
120 
121  return (113*sqrt(0.8333333333333334)*Pi_p1by3*
122  pow(mc/etap35,0.8333333333333334)*sqrt(eta)*
123  ((672*f*mc)/etap35 -
124  (Pi_p1by3*Pi_p1by3*
125  pow((f*mc)/etap35,1.6666666666666667)*
126  (502429 - 591368*eta + 1680*etap2))/
127  (-113 + 76*eta) + (113904*Pi_p1by3*
128  pow((f*mc)/etap35,1.3333333333333333)*
129  (-81 + 4*eta)*chi)/pow(113 - 76*eta,2)))/
130  (32256.*pow(f,1.1666666666666667));
131 }
132 
133 /*
134  * Derivative of the amplitude with respect to \f$\eta\f$ (symm. mass ratio)
135  */
136 static REAL8 XLALSimInspiralTaylorF2RedSpinDerivAEta(
137  const REAL8 mc_msun, /**< chirp mass (M_sun) */
138  const REAL8 eta, /**< symmetric mass ratio */
139  const REAL8 chi, /**< reduced-spin parameter */
140  const REAL8 f /**< Fourier frequency (Hz) */
141 ) {
142  const REAL8 mc = mc_msun * LAL_MTSUN_SI;
143  const REAL8 etap2 = eta*eta;
144  const REAL8 etap3 = etap2*eta;
145  const REAL8 etap4 = etap3*eta;
146  const REAL8 etap5 = etap4*eta;
147  const REAL8 chip2 = chi*chi;
148  const REAL8 mcp2 = mc*mc;
149  const REAL8 eta_fac = (-113 + 76 * eta) * (-113 + 76 * eta) * (-113 + 76 * eta);
150  REAL8 eta_three_fifths = pow(eta, 0.6);
151  REAL8 f_mc_eta_one_third = cbrt((f*mc)/eta_three_fifths);
152 
153  return (pow(mc/eta_three_fifths,0.8333333333333334)*
154  (2235340800*f_mc_eta_one_third * f_mc_eta_one_third*
155  eta_three_fifths*eta_three_fifths*eta_fac*(-743 + 1386*eta) +
156  f*f*mcp2*Pi_p4by3*
157  eta_fac *
158  (257959397806029 - 36738273116160*LAL_GAMMA -
159  95047630643350*eta - 12126223216800*etap2 +
160  5541700903200*etap3 +
161  7823692800*Pi_p2*(-1920 + 451*eta) -
162  1940400*Pi_p4by3*
163  f_mc_eta_one_third*
164  (-5111593 - 2311072*eta + 64752*etap2)) +
165  369600*f*mc*Pi_p1by3*eta_three_fifths*
166  (12192768*LAL_PI*eta_fac +
167  35962920*Pi_p5by3*
168  f_mc_eta_one_third * f_mc_eta_one_third *
169  eta_fac -
170  28703808*eta_fac*chi -
171  71190*Pi_p2by3*
172  f_mc_eta_one_third * f_mc_eta_one_third *
173  (-6415515901 + 12944580756*eta -
174  10861276272*etap2 + 3401313728*etap3)*
175  chi + Pi_p1by3*
176  f_mc_eta_one_third *
177  (34636018030144*etap3 -
178  27532505500416*etap4 +
179  6407002215936*etap5 -
180  1442897*(-7266251 + 20575296*chip2) -
181  21696*etap2*(-4887203 + 102257316*chip2) +
182  76614*eta*(-320998087 + 907387488*chip2))) -
183  12246091038720*f*f*mcp2*Pi_p4by3*
184  eta_fac*log((64*f*mc*LAL_PI)/eta_three_fifths)))/
185  (1.5021490176e12*sqrt(30)*pow(f,1.1666666666666667)*pow(eta,1.7)*
186  eta_fac);
187 }
188 
189 /*
190  * Derivative of the amplitude with respect to the chirp mass
191  */
192 static REAL8 XLALSimInspiralTaylorF2RedSpinDerivAMChirp(
193  const REAL8 mc_msun, /**< chirp mass (M_sun) */
194  const REAL8 eta, /**< symmetric mass ratio */
195  const REAL8 chi, /**< reduced-spin parameter */
196  const REAL8 f /**< Fourier frequency (Hz) */
197 ) {
198  const REAL8 mc = mc_msun * LAL_MTSUN_SI;
199  const REAL8 etap2 = eta*eta;
200  const REAL8 etap3 = etap2*eta;
201  const REAL8 etap4 = etap3*eta;
202  const REAL8 chip2 = chi*chi;
203  const REAL8 mcp2 = mc*mc;
204  REAL8 eta_three_fifths = pow(eta, 0.6);
205 
206  return (sqrt(0.8333333333333334)*(5 -
207  (17*f*f*mcp2*Pi_p2*Pi_p2*(-320 + 451*eta))/
208  (96.*eta_three_fifths*eta_three_fifths) +
209  (3*Pi_p2by3*
210  pow((f*mc)/eta_three_fifths,0.6666666666666666)*
211  (743 + 924*eta))/224. +
212  (5*Pi_p3/Pi_p1by3*
213  pow((f*mc)/eta_three_fifths,1.6666666666666667)*
214  (-4757 + 4788*eta))/448. -
215  (19*pow(LAL_PI,3.3333333333333335)*
216  pow((f*mc)/eta_three_fifths,2.3333333333333335)*
217  (5111593 + 8088752*eta + 151088*etap2))/2.709504e6 +
218  (f*mc*Pi_p2*(-33047278387200*eta_three_fifths +
219  f*mc*(-1471974996766431 + 208183547658240*LAL_GAMMA +
220  3231619441873900*eta - 103072897342800*etap2 +
221  20935314523200*etap3)))/
222  (1.5021490176e12*eta_three_fifths*eta_three_fifths) +
223  (1243*f*mc*LAL_PI*chi)/(24.*eta_three_fifths) -
224  (565*Pi_p5by3*
225  pow((f*mc)/eta_three_fifths,1.6666666666666667)*
226  (502429 - 591368*eta + 1680*etap2)*chi)/
227  (5376.*(-113 + 76*eta)) +
228  (13*Pi_p4by3*
229  pow((f*mc)/eta_three_fifths,1.3333333333333333)*
230  (2490853280*etap2 - 112068617472*etap3 +
231  56201773824*etap4 +
232  1808*eta*(-1708561 + 7175952*chip2) -
233  12769*(-7266251 + 20575296*chip2)))/
234  (8.128512e6*pow(113 - 76*eta,2)) +
235  (14552*f*f*mcp2*Pi_p2*
236  log((64*f*mc*LAL_PI)/eta_three_fifths))/(315.*eta_three_fifths*eta_three_fifths)))/
237  (12.*pow(f,1.1666666666666667)*Pi_p2by3*
238  pow(mc/eta_three_fifths,0.16666666666666666)*pow(eta,0.1));
239 }
240 
241 /*
242  * Derivative of the phasae with respect to \f$\chi\f$ (reduced spin parameter)
243  */
244 static REAL8 XLALSimInspiralTaylorF2RedSpinDerivPsiChi(
245  const REAL8 mc_msun, /**< chirp mass (M_sun) */
246  const REAL8 eta, /**< symmetric mass ratio */
247  const REAL8 chi, /**< reduced-spin parameter */
248  const REAL8 f /**< Fourier frequency (Hz) */
249 ) {
250  const REAL8 mc = mc_msun * LAL_MTSUN_SI;
251  const REAL8 etap2 = eta*eta;
252  REAL8 eta_three_fifths = pow(eta, 0.6);
253 
254  return (113*((756*f*mc)/eta_three_fifths +
255  (320355*Pi_p1by3*
256  pow((f*mc)/eta_three_fifths,1.3333333333333333)*
257  (-81 + 4*eta)*chi)/pow(113 - 76*eta,2) -
258  (5*Pi_p2by3*
259  pow((f*mc)/eta_three_fifths,1.6666666666666667)*
260  (-146597 + 135856*eta + 17136*etap2)*
261  (1 + log(f / 40.)))/
262  (-113 + 76*eta)))/
263  (96768.*Pi_p2by3*
264  pow((f*mc)/eta_three_fifths,1.6666666666666667)*eta);
265 }
266 
267 /*
268  * Derivative of the phasae with respect to \f$\eta\f$ (symmetric mass ratio)
269  */
270 static REAL8 XLALSimInspiralTaylorF2RedSpinDerivPsiEta(
271  const REAL8 mc_msun, /**< chirp mass (M_sun) */
272  const REAL8 eta, /**< symmetric mass ratio */
273  const REAL8 chi, /**< reduced-spin parameter */
274  const REAL8 f /**< Fourier frequency (Hz) */
275 ) {
276  const REAL8 mc = mc_msun * LAL_MTSUN_SI;
277  const REAL8 etap2 = eta*eta;
278  const REAL8 etap3 = etap2*eta;
279  const REAL8 etap4 = etap3*eta;
280  const REAL8 etap5 = etap4*eta;
281  const REAL8 chip2 = chi*chi;
282  const REAL8 mcp2 = mc*mc;
283  REAL8 eta_three_fifths = pow(eta, 0.6);
284  const REAL8 eta_fac = -113 + 76 * eta;
285  const REAL8 fmc_fac = cbrt(f * mc / eta_three_fifths);
286 
287  return (-77616000*f*mc*LAL_PI*eta_fac*
288  (23187*LAL_PI*eta_fac*eta_fac -
289  113*(16565461 - 22282744*eta + 12261424*etap2)*chi)*
290  log(fmc_fac/
291  (2.*cbrt(5)*
292  cbrt(mc/eta_three_fifths))) +
293  (31046400*fmc_fac*fmc_fac*
294  eta_three_fifths * eta_three_fifths *eta_fac*eta_fac*eta_fac*(-743 + 1386*eta) -
295  f*f*mcp2*LAL_PI * Pi_p1by3 *
296  eta_fac * eta_fac * eta_fac *
297  (33984313019673 - 4592284139520*LAL_GAMMA -
298  12118079538950*eta - 413215941600*etap2 +
299  2083465692000*etap3 +
300  977961600*Pi_p2*(-3072 + 451*eta) +
301  323400*LAL_PI*Pi_p1by3*
302  fmc_fac*
303  (15419335 + 3633744*eta + 2132496*etap2)) -
304  18480*f*mc*Pi_p1by3*eta_three_fifths*
305  (-6096384*LAL_PI*eta_fac * eta_fac * eta_fac +
306  32461800*Pi_p2/Pi_p1by3*
307  fmc_fac*fmc_fac*
308  eta_fac * eta_fac * eta_fac +
309  14351904*eta_fac * eta_fac * eta_fac*chi -
310  158200*LAL_PI/Pi_p1by3*
311  fmc_fac*fmc_fac*
312  (-1871897093 + 3776925108*eta -
313  3079029456*etap2 + 931868224*etap3)*
314  chi - 5*Pi_p1by3*
315  fmc_fac*
316  (14990425815136*etap3 -
317  12186233587584*etap4 +
318  2866657264128*etap5 -
319  1442897*(-3058673 + 5143824*chip2) +
320  612912*eta*(-17749451 + 28355859*chip2) -
321  2712*etap2*(-202619251 + 204514632*chip2)))
322  + 1530761379840*f*f*mcp2*LAL_PI*Pi_p1by3*
323  eta_fac * eta_fac * eta_fac*log((64*f*mc*LAL_PI)/eta_three_fifths))/
324  (fmc_fac * fmc_fac *eta_three_fifths))/
325  (5.007163392e11*f*mc*LAL_PI*etap2*eta_fac * eta_fac * eta_fac);
326 }
327 
328 /*
329  * Derivative of the phase with respect to the chirp mass
330  */
331 static REAL8 XLALSimInspiralTaylorF2RedSpinDerivPsiMChirp(
332  const REAL8 mc_msun, /**< chirp mass (M_sun) */
333  const REAL8 eta, /**< symmetric mass ratio */
334  const REAL8 chi, /**< reduced-spin parameter */
335  const REAL8 f /**< Fourier frequency (Hz) */
336 ) {
337  const REAL8 mc = mc_msun * LAL_MTSUN_SI;
338  const REAL8 etap2 = eta*eta;
339  const REAL8 etap3 = etap2*eta;
340  const REAL8 etap4 = etap3*eta;
341  const REAL8 chip2 = chi*chi;
342  const REAL8 mcp2 = mc*mc;
343  const REAL8 mcp3 = mcp2*mc;
344  REAL8 eta_three_fifths = pow(eta, 0.6);
345 
346  return (cbrt((f*mc)/eta_three_fifths)*
347  (-93139200*pow(113 - 76*eta,2)*eta_three_fifths * eta_three_fifths *
348  (252 + Pi_p2by3*
349  pow((f*mc)/eta_three_fifths,0.6666666666666666)*
350  (743 + 924*eta)) +
351  f*f*mcp2*LAL_PI*LAL_PI*pow(113 - 76*eta,2)*
352  (10052469856691 - 1530761379840*LAL_GAMMA -
353  24236159077900*eta + 206607970800*etap2 -
354  462992376000*etap3 +
355  1955923200*Pi_p2*(-512 + 451*eta) -
356  184800*Pi_p4by3*
357  cbrt((f*mc)/eta_three_fifths)*
358  (-15419335 - 12718104*eta + 4975824*etap2)) +
359  9240*f*mc*LAL_PI*eta_three_fifths*
360  (16257024*LAL_PI*pow(113 - 76*eta,2) -
361  38271744*pow(113 - 76*eta,2)*chi -
362  5*cbrt(LAL_PI)*
363  cbrt((f*mc)/eta_three_fifths)*
364  (-20737091296*etap2 - 43167841920*etap3 +
365  25146116352*etap4 +
366  904*eta*(19183315 + 3587976*chip2) -
367  12769*(-3058673 + 5143824*chip2))) -
368  510253793280*f*f*mcp2*LAL_PI*LAL_PI*
369  pow(113 - 76*eta,2)*log((64*f*mc*LAL_PI)/eta_three_fifths)))/
370  (6.0085960704e11*f*f*mcp3*pow(LAL_PI,1.6666666666666667)*
371  pow(113 - 76*eta,2)*eta);
372 }
373 
374 /**
375  * @addtogroup LALSimInspiralTaylorF2ReducedSpin_c
376  * @{
377  */
378 
379 /*
380  * Compute the template-space metric of "reduced-spin" PN templates in
381  * Mchirp-eta-chi parameter space.
382  */
383 int XLALSimInspiralTaylorF2RedSpinMetricMChirpEtaChi(
384  REAL8 *gamma00, /**< template metric coeff. 00 in mChirp-eta-chi */
385  REAL8 *gamma01, /**< template metric coeff. 01/10 in mChirp-eta-chi */
386  REAL8 *gamma02, /**< template metric coeff. 02/20 in mChirp-eta-chi */
387  REAL8 *gamma11, /**< template metric coeff. 11 in mChirp-eta-chi */
388  REAL8 *gamma12, /**< template metric coeff. 12/21 in mChirp-eta-chi */
389  REAL8 *gamma22, /**< template metric coeff. 22 in mChirp-eta-chi */
390  const REAL8 mc, /**< chirp mass (in solar mass) */
391  const REAL8 eta, /**< symmetric mass ratio */
392  const REAL8 chi, /**< reduced-spin parameter */
393  const REAL8 fLow, /**< low-frequency cutoff (Hz) */
394  const REAL8FrequencySeries *Sh /**< PSD in strain per root Hertz */
395 ) {
396  REAL8Vector *A=NULL, *dAMc=NULL, *dAEta=NULL;
397  REAL8Vector *dAChi=NULL, *dAT0=NULL, *dAPhi0=NULL, *dPsiMc=NULL;
398  REAL8Vector *dPsiEta=NULL, *dPsiChi=NULL, *dPsiT0=NULL, *dPsiPhi0=NULL;
399 
400  const REAL8 df = Sh->deltaF;
401  REAL8 hSqr = 0.;
402  int s = 0;
403 
404  /* compute the Schwarzschild ISCO frequency */
405  REAL8 fCut = 1. / (LAL_PI*mc*pow(eta, -0.6) * sqrt(6 * 6 * 6) * LAL_MTSUN_SI);
406 
407  /* create a view of the PSD between fLow and fCut */
408  size_t nBins = (fCut - fLow) / df;
409  size_t k = nBins;
410  REAL8Vector Shdata = {nBins, Sh->data->data + (size_t) (fLow / df)}; /* copy the Vector, including its pointer to the actual data */
411  /* drop low-frequency samples */
412  Shdata.length = nBins; /* drop high-frequency samples */
413 
414  /* allocate memory for various vectors */
415  A = XLALCreateREAL8Vector(nBins);
416  dAMc = XLALCreateREAL8Vector(nBins);
417  dAEta = XLALCreateREAL8Vector(nBins);
418  dAChi = XLALCreateREAL8Vector(nBins);
419  dAT0 = XLALCreateREAL8Vector(nBins);
420  dAPhi0 = XLALCreateREAL8Vector(nBins);
421  dPsiMc = XLALCreateREAL8Vector(nBins);
422  dPsiEta = XLALCreateREAL8Vector(nBins);
423  dPsiChi = XLALCreateREAL8Vector(nBins);
424  dPsiT0 = XLALCreateREAL8Vector(nBins);
425  dPsiPhi0 = XLALCreateREAL8Vector(nBins);
426 
427  /* create a frequency vector from fLow to fCut with frequency resolution df */
428  for (;k--;) {
429  const REAL8 f = fLow + k * df;
430 
431  /* compute the amplitude of the frequency-domain waveform */
432  A->data[k] = XLALSimInspiralTaylorF2RedSpinAofF(mc, eta, chi, f);
433 
434  /* compute the square of the template norm */
435  hSqr += A->data[k] * A->data[k] / Shdata.data[k];
436 
437  /* compute the waveform deratives with respect to the parameters */
438  dAMc->data[k] = XLALSimInspiralTaylorF2RedSpinDerivAMChirp(mc, eta, chi, f);
439  dAEta->data[k] = XLALSimInspiralTaylorF2RedSpinDerivAEta(mc, eta, chi, f);
440  dAChi->data[k] = XLALSimInspiralTaylorF2RedSpinDerivAChi(mc, eta, chi, f);
441  dAT0->data[k] = 0.;
442  dAPhi0->data[k] = 0.;
443  dPsiMc->data[k] = XLALSimInspiralTaylorF2RedSpinDerivPsiMChirp(mc, eta, chi, f);
444  dPsiEta->data[k] = XLALSimInspiralTaylorF2RedSpinDerivPsiEta(mc, eta, chi, f);
445  dPsiChi->data[k] = XLALSimInspiralTaylorF2RedSpinDerivPsiChi(mc, eta, chi, f);
446  dPsiT0->data[k] = LAL_TWOPI * f;
447  dPsiPhi0->data[k] = 1.;
448  }
449  hSqr *= 4 * df;
450 
451  /* allocate memory, and initialize the Fisher matrix */
452  gsl_matrix * g = gsl_matrix_calloc (5, 5);
453 
454  /* compute the components of the Fisher matrix in coordinates mc, eta, chi, t0, phi0 */
455  gsl_matrix_set (g, 0,0, MetricCoeffs(A, dPsiMc, dPsiMc, dAMc, dAMc, &Shdata, hSqr, df));
456  gsl_matrix_set (g, 0,1, MetricCoeffs(A, dPsiMc, dPsiEta, dAMc, dAEta, &Shdata, hSqr, df));
457  gsl_matrix_set (g, 0,2, MetricCoeffs(A, dPsiMc, dPsiChi, dAMc, dAChi, &Shdata, hSqr, df));
458  gsl_matrix_set (g, 0,3, MetricCoeffs(A, dPsiMc, dPsiT0, dAMc, dAT0, &Shdata, hSqr, df));
459  gsl_matrix_set (g, 0,4, MetricCoeffs(A, dPsiMc, dPsiPhi0, dAMc, dAPhi0, &Shdata, hSqr, df));
460 
461  gsl_matrix_set (g, 1,0, gsl_matrix_get(g, 0,1));
462  gsl_matrix_set (g, 1,1, MetricCoeffs(A, dPsiEta, dPsiEta, dAEta, dAEta, &Shdata, hSqr, df));
463  gsl_matrix_set (g, 1,2, MetricCoeffs(A, dPsiEta, dPsiChi, dAEta, dAChi, &Shdata, hSqr, df));
464  gsl_matrix_set (g, 1,3, MetricCoeffs(A, dPsiEta, dPsiT0, dAEta, dAT0, &Shdata, hSqr, df));
465  gsl_matrix_set (g, 1,4, MetricCoeffs(A, dPsiEta, dPsiPhi0, dAEta, dAPhi0, &Shdata, hSqr, df));
466 
467  gsl_matrix_set (g, 2,0, gsl_matrix_get(g, 0,2));
468  gsl_matrix_set (g, 2,1, gsl_matrix_get(g, 1,2));
469  gsl_matrix_set (g, 2,2, MetricCoeffs(A, dPsiChi, dPsiChi, dAChi, dAChi, &Shdata, hSqr, df));
470  gsl_matrix_set (g, 2,3, MetricCoeffs(A, dPsiChi, dPsiT0, dAChi, dAT0, &Shdata, hSqr, df));
471  gsl_matrix_set (g, 2,4, MetricCoeffs(A, dPsiChi, dPsiPhi0, dAChi, dAPhi0, &Shdata, hSqr, df));
472 
473  gsl_matrix_set (g, 3,0, gsl_matrix_get(g, 0,3));
474  gsl_matrix_set (g, 3,1, gsl_matrix_get(g, 1,3));
475  gsl_matrix_set (g, 3,2, gsl_matrix_get(g, 2,3));
476  gsl_matrix_set (g, 3,3, MetricCoeffs(A, dPsiT0, dPsiT0, dAT0, dAT0, &Shdata, hSqr, df));
477  gsl_matrix_set (g, 3,4, MetricCoeffs(A, dPsiT0, dPsiPhi0, dAT0, dAPhi0, &Shdata, hSqr, df));
478 
479  gsl_matrix_set (g, 4,0, gsl_matrix_get(g, 0,4));
480  gsl_matrix_set (g, 4,1, gsl_matrix_get(g, 1,4));
481  gsl_matrix_set (g, 4,2, gsl_matrix_get(g, 2,4));
482  gsl_matrix_set (g, 4,3, gsl_matrix_get(g, 3,4));
483  gsl_matrix_set (g, 4,4, MetricCoeffs(A, dPsiPhi0, dPsiPhi0, dAPhi0, dAPhi0, &Shdata, hSqr, df));
484 
485  /* free the memory */
486  XLALDestroyREAL8Vector(A);
487  XLALDestroyREAL8Vector(dAMc);
488  XLALDestroyREAL8Vector(dAEta);
489  XLALDestroyREAL8Vector(dAChi);
490  XLALDestroyREAL8Vector(dAT0);
491  XLALDestroyREAL8Vector(dAPhi0);
492  XLALDestroyREAL8Vector(dPsiMc);
493  XLALDestroyREAL8Vector(dPsiEta);
494  XLALDestroyREAL8Vector(dPsiChi);
495  XLALDestroyREAL8Vector(dPsiT0);
496  XLALDestroyREAL8Vector(dPsiPhi0);
497 
498  {
499  /* Form submatrices g1, g2, g3, g4, defined as:
500  * g = [ g1 g2
501  * g4 g3 ] */
502  gsl_matrix_view g1v = gsl_matrix_submatrix (g, 0, 0, 3, 3);
503  gsl_matrix_view g2v = gsl_matrix_submatrix (g, 0, 3, 3, 2);
504  gsl_matrix_view g3v = gsl_matrix_submatrix (g, 3, 3, 2, 2);
505  gsl_matrix_view g4v = gsl_matrix_submatrix (g, 3, 0, 2, 3);
506 
507  gsl_matrix * g1 = gsl_matrix_calloc (3, 3);
508  gsl_matrix * g2 = gsl_matrix_calloc (3, 2);
509  gsl_matrix * g3 = gsl_matrix_calloc (2, 2);
510  gsl_matrix * g4 = gsl_matrix_calloc (2, 3);
511  gsl_matrix * g3invg4 = gsl_matrix_calloc (2, 3);
512  gsl_matrix * g2g3invg4 = gsl_matrix_calloc (3, 3);
513 
514  /* Project out the t0 and phi0 dimensions: gamma = g1 - g2 g3^{-1} g4 */
515  gsl_matrix * g3inv = gsl_matrix_calloc (2, 2);
516  gsl_permutation * p = gsl_permutation_calloc (2);
517 
518  gsl_matrix_memcpy (g1, &g1v.matrix);
519  gsl_matrix_memcpy (g2, &g2v.matrix);
520  gsl_matrix_memcpy (g3, &g3v.matrix);
521  gsl_matrix_memcpy (g4, &g4v.matrix);
522  gsl_matrix_free (g);
523 
524  gsl_linalg_LU_decomp (g3, p, &s);
525  gsl_linalg_LU_invert (g3, p, g3inv);
526  gsl_permutation_free (p);
527  gsl_matrix_free (g3);
528 
529  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, g3inv, g4, 0.0, g3invg4);
530  gsl_matrix_free (g4);
531  gsl_matrix_free (g3inv);
532  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, g2, g3invg4, 0.0, g2g3invg4);
533  gsl_matrix_free (g2);
534  gsl_matrix_free (g3invg4);
535 
536  gsl_matrix_sub (g1, g2g3invg4);
537  gsl_matrix_free (g2g3invg4);
538 
539  *gamma00 = gsl_matrix_get(g1, 0, 0);
540  *gamma01 = gsl_matrix_get(g1, 0, 1);
541  *gamma02 = gsl_matrix_get(g1, 0, 2);
542  *gamma11 = gsl_matrix_get(g1, 1, 1);
543  *gamma12 = gsl_matrix_get(g1, 1, 2);
544  *gamma22 = gsl_matrix_get(g1, 2, 2);
545  gsl_matrix_free (g1);
546  }
547 
548  return XLAL_SUCCESS;
549 }
550 
551 
552 /**
553  * Compute the template-space metric of "reduced-spin" PN templates in
554  * theta0, theta3, theta3s parameter space.
555  */
556 int XLALSimInspiralTaylorF2RedSpinComputeNoiseMoments(
557  REAL8Vector *momI_0, /**< noise moments: \f$momI_0(f) = \int_{f0}^f (f'/f0)^{(0-17)/3} df'\f$ */
558  REAL8Vector *momI_2, /**< noise moments: \f$momI_2(f) = \int_{f0}^f (f'/f0)^{(2-17)/3} df'\f$ */
559  REAL8Vector *momI_3, /**< noise moments: \f$momI_3(f) = \int_{f0}^f (f'/f0)^{(3-17)/3} df'\f$ */
560  REAL8Vector *momI_4, /**< noise moments: \f$momI_4(f) = \int_{f0}^f (f'/f0)^{(4-17)/3} df'\f$ */
561  REAL8Vector *momI_5, /**< noise moments: \f$momI_5(f) = \int_{f0}^f (f'/f0)^{(5-17)/3} df'\f$ */
562  REAL8Vector *momI_6, /**< noise moments: \f$momI_6(f) = \int_{f0}^f (f'/f0)^{(6-17)/3} df'\f$ */
563  REAL8Vector *momI_7, /**< noise moments: \f$momI_7(f) = \int_{f0}^f (f'/f0)^{(7-17)/3} df'\f$ */
564  REAL8Vector *momI_8, /**< noise moments: \f$momI_8(f) = \int_{f0}^f (f'/f0)^{(8-17)/3} df'\f$ */
565  REAL8Vector *momI_9, /**< noise moments: \f$momI_9(f) = \int_{f0}^f (f'/f0)^{(9-17)/3} df'\f$ */
566  REAL8Vector *momI_10, /**< noise moments: \f$momI_10(f) = \int_{f0}^f (f'/f0)^{(10-17)/3} df'\f$ */
567  REAL8Vector *momI_11, /**< noise moments: \f$momI_11(f) = \int_{f0}^f (f'/f0)^{(11-17)/3} df'\f$ */
568  REAL8Vector *momI_12, /**< noise moments: \f$momI_12(f) = \int_{f0}^f (f'/f0)^{(12-17)/3} df'\f$ */
569  REAL8Vector *momI_13, /**< noise moments: \f$momI_13(f) = \int_{f0}^f (f'/f0)^{(13-17)/3} df'\f$ */
570  REAL8Vector *momI_14, /**< noise moments: \f$momI_14(f) = \int_{f0}^f (f'/f0)^{(14-17)/3} df'\f$ */
571  REAL8Vector *momI_15, /**< noise moments: \f$momI_15(f) = \int_{f0}^f (f'/f0)^{(15-17)/3} df'\f$ */
572  REAL8Vector *momI_16, /**< noise moments: \f$momI_16(f) = \int_{f0}^f (f'/f0)^{(16-17)/3} df'\f$ */
573  REAL8Vector *momJ_5, /**< noise moments: \f$momJ_5(f) = \int_{f0}^f (f'/f0)^{(5-17)/3} log(f'/f0) df'\f$ */
574  REAL8Vector *momJ_6, /**< noise moments: \f$momJ_6(f) = \int_{f0}^f (f'/f0)^{(6-17)/3} log(f'/f0) df'\f$ */
575  REAL8Vector *momJ_7, /**< noise moments: \f$momJ_7(f) = \int_{f0}^f (f'/f0)^{(7-17)/3} log(f'/f0) df'\f$ */
576  REAL8Vector *momJ_8, /**< noise moments: \f$momJ_8(f) = \int_{f0}^f (f'/f0)^{(8-17)/3} log(f'/f0) df'\f$ */
577  REAL8Vector *momJ_9, /**< noise moments: \f$momJ_9(f) = \int_{f0}^f (f'/f0)^{(9-17)/3} log(f'/f0) df'\f$ */
578  REAL8Vector *momJ_10, /**< noise moments: \f$momJ_10(f) = \int_{f0}^f (f'/f0)^{(10-17)/3} log(f'/f0) df'\f$ */
579  REAL8Vector *momJ_11, /**< noise moments: \f$momJ_11(f) = \int_{f0}^f (f'/f0)^{(11-17)/3} log(f'/f0) df'\f$ */
580  REAL8Vector *momJ_12, /**< noise moments: \f$momJ_12(f) = \int_{f0}^f (f'/f0)^{(12-17)/3} log(f'/f0) df'\f$ */
581  REAL8Vector *momJ_13, /**< noise moments: \f$momJ_13(f) = \int_{f0}^f (f'/f0)^{(13-17)/3} log(f'/f0) df'\f$ */
582  REAL8Vector *momJ_14, /**< noise moments: \f$momJ_14(f) = \int_{f0}^f (f'/f0)^{(14-17)/3} log(f'/f0) df'\f$ */
583  REAL8Vector *momK_10, /**< noise moments: \f$momK_10(f) = \int_{f0}^f (f'/f0)^{(10-17)/3} log(f'/f0) log(f'/f0) df'\f$ */
584  REAL8Vector *momK_11, /**< noise moments: \f$momK_11(f) = \int_{f0}^f (f'/f0)^{(11-17)/3} log(f'/f0) log(f'/f0) df'\f$ */
585  REAL8Vector *momK_12, /**< noise moments: \f$momK_12(f) = \int_{f0}^f (f'/f0)^{(12-17)/3} log(f'/f0) log(f'/f0) df'\f$ */
586  REAL8Vector *Sh, /**< one sided PSD of the detector noise: Sh(f) for f = [fLow, fNyq] */
587  REAL8 fLow, /**< low frequency cutoff (Hz) */
588  REAL8 df) /**< frequency resolution of the psd vector (Hz) */
589 {
590  size_t i;
591  const size_t ilow = (size_t) (fLow / df); /* PSD starts at 0 Hz; others start at fLow */
592 
593  /* some minimum error checking */
594  if ( !momI_0 || !momI_2 || !momI_3 || !momI_4 || !momI_5 || !momI_6 || !momI_7 || !momI_8 || !momI_9 || !momI_10 || !momI_11 || !momI_12 || !momI_13 || !momI_14 || !momI_15 || !momI_16 || !momJ_5 || !momJ_6 || !momJ_7 || !momJ_8 || !momJ_9 || !momJ_10 || !momJ_11 || !momJ_12 || !momJ_13 || !momJ_14 || !momK_10 || !momK_11 || !momK_12) {
595  XLALPrintError("Moments not initialized");
596  XLAL_ERROR(XLAL_EFAULT);
597  }
598  if (momI_0->length > (Sh->length - ilow)) {
599  XLALPrintError("Sh not long enough to fill moment vectors");
600  XLAL_ERROR(XLAL_EDOM);
601  }
602 
603  momI_0->data[0] = 0.;
604  momI_2->data[0] = 0.;
605  momI_3->data[0] = 0.;
606  momI_4->data[0] = 0.;
607  momI_5->data[0] = 0.;
608  momI_6->data[0] = 0.;
609  momI_7->data[0] = 0.;
610  momI_8->data[0] = 0.;
611  momI_9->data[0] = 0.;
612  momI_10->data[0] = 0.;
613  momI_11->data[0] = 0.;
614  momI_12->data[0] = 0.;
615  momI_13->data[0] = 0.;
616  momI_14->data[0] = 0.;
617  momI_15->data[0] = 0.;
618  momI_16->data[0] = 0.;
619 
620  momJ_5->data[0] = 0.;
621  momJ_6->data[0] = 0.;
622  momJ_7->data[0] = 0.;
623  momJ_8->data[0] = 0.;
624  momJ_9->data[0] = 0.;
625  momJ_10->data[0] = 0.;
626  momJ_11->data[0] = 0.;
627  momJ_12->data[0] = 0.;
628  momJ_13->data[0] = 0.;
629  momJ_14->data[0] = 0.;
630 
631  momK_10->data[0] = 0.;
632  momK_11->data[0] = 0.;
633  momK_12->data[0] = 0.;
634 
635  const REAL8 fLowmSevenBythree = pow(fLow, -7./3.);
636  const REAL8 fLowFac = fLowmSevenBythree * df;
637  for (i=1; i < momI_0->length; i++) {
638  const REAL8 psdfac = Sh->data[i + ilow];
639  if (psdfac) {
640  const REAL8 fbyfLow = (i*df+fLow)/fLow;
641  const REAL8 logfbyfLow = log(fbyfLow);
642  const REAL8 logfbyfLowSq = logfbyfLow * logfbyfLow;
643  const REAL8 fbyfLowp2 = fbyfLow * fbyfLow;
644  const REAL8 fbyfLowp3 = fbyfLowp2 * fbyfLow;
645  const REAL8 fbyfLowp4 = fbyfLowp3 * fbyfLow;
646  const REAL8 fbyfLowp5 = fbyfLowp4 * fbyfLow;
647  const REAL8 cbrtfbyfLow = cbrt(fbyfLow);
648 
649  momI_0->data[i] = momI_0->data[i-1] + fLowFac/(fbyfLowp5*cbrtfbyfLow*cbrtfbyfLow)/psdfac; /* f^{-17/3} */
650  momI_2->data[i] = momI_2->data[i-1] + fLowFac/fbyfLowp5/psdfac; /* f^{-15/3} */
651  momI_3->data[i] = momI_3->data[i-1] + fLowFac/fbyfLowp5*cbrtfbyfLow/psdfac; /* f^{-14/3} */
652  momI_4->data[i] = momI_4->data[i-1] + fLowFac/fbyfLowp4/cbrtfbyfLow/psdfac; /* f^{-13/3} */
653  momI_5->data[i] = momI_5->data[i-1] + fLowFac/fbyfLowp4/psdfac; /* f^{-12/3} */
654  momI_6->data[i] = momI_6->data[i-1] + fLowFac/fbyfLowp4*cbrtfbyfLow/psdfac; /* f^{-11/3} */
655  momI_7->data[i] = momI_7->data[i-1] + fLowFac/fbyfLowp3/cbrtfbyfLow/psdfac; /* f^{-10/3} */
656  momI_8->data[i] = momI_8->data[i-1] + fLowFac/fbyfLowp3/psdfac; /* f^{-9/3} */
657  momI_9->data[i] = momI_9->data[i-1] + fLowFac/fbyfLowp3*cbrtfbyfLow/psdfac; /* f^{-8/3} */
658  momI_10->data[i] = momI_10->data[i-1] + fLowFac/fbyfLowp2/cbrtfbyfLow/psdfac; /* f^{-7/3} */
659  momI_11->data[i] = momI_11->data[i-1] + fLowFac/fbyfLowp2/psdfac; /* f^{-6/3} */
660  momI_12->data[i] = momI_12->data[i-1] + fLowFac/fbyfLowp2*cbrtfbyfLow/psdfac; /* f^{-5/3} */
661  momI_13->data[i] = momI_13->data[i-1] + fLowFac/fbyfLow/cbrtfbyfLow/psdfac; /* f^{-4/3} */
662  momI_14->data[i] = momI_14->data[i-1] + fLowFac/fbyfLow/psdfac; /* f^{-3/3} */
663  momI_15->data[i] = momI_15->data[i-1] + fLowFac/fbyfLow*cbrtfbyfLow/psdfac; /* f^{-2/3} */
664  momI_16->data[i] = momI_16->data[i-1] + fLowFac/cbrtfbyfLow/psdfac; /* f^{-1/3} */
665 
666  momJ_5->data[i] = momJ_5->data[i-1] + fLowFac*logfbyfLow/fbyfLowp4/psdfac; /* f^{-12/3} */
667  momJ_6->data[i] = momJ_6->data[i-1] + fLowFac*logfbyfLow/fbyfLowp4*cbrtfbyfLow/psdfac; /* f^{-11/3} */
668  momJ_7->data[i] = momJ_7->data[i-1] + fLowFac*logfbyfLow/fbyfLowp3/cbrtfbyfLow/psdfac; /* f^{-10/3} */
669  momJ_8->data[i] = momJ_8->data[i-1] + fLowFac*logfbyfLow/fbyfLowp3/psdfac; /* f^{-9/3} */
670  momJ_9->data[i] = momJ_9->data[i-1] + fLowFac*logfbyfLow/fbyfLowp3*cbrtfbyfLow/psdfac; /* f^{-8/3} */
671  momJ_10->data[i] = momJ_10->data[i-1] + fLowFac*logfbyfLow/fbyfLowp2/cbrtfbyfLow/psdfac; /* f^{-7/3} */
672  momJ_11->data[i] = momJ_11->data[i-1] + fLowFac*logfbyfLow/fbyfLowp2/psdfac; /* f^{-6/3} */
673  momJ_12->data[i] = momJ_12->data[i-1] + fLowFac*logfbyfLow/fbyfLowp2*cbrtfbyfLow/psdfac; /* f^{-5/3} */
674  momJ_13->data[i] = momJ_13->data[i-1] + fLowFac*logfbyfLow/fbyfLow/cbrtfbyfLow/psdfac; /* f^{-4/3} */
675  momJ_14->data[i] = momJ_14->data[i-1] + fLowFac*logfbyfLow/fbyfLow/psdfac; /* f^{-3/3} */
676 
677  momK_10->data[i] = momK_10->data[i-1] + fLowFac*logfbyfLowSq/fbyfLowp2/cbrtfbyfLow/psdfac; /* f^{-7/3} */
678  momK_11->data[i] = momK_11->data[i-1] + fLowFac*logfbyfLowSq/fbyfLowp2/psdfac; /* f^{-6/3} */
679  momK_12->data[i] = momK_12->data[i-1] + fLowFac*logfbyfLowSq/fbyfLowp2*cbrtfbyfLow/psdfac; /* f^{-5/3} */
680  }
681  }
682 
683  return XLAL_SUCCESS;
684 }
685 
686 
687 /**
688  * Compute the Fisher information matrix of "reduced-spin" PN templates in
689  * theta0, theta3, theta3s, t0, phi0 parameter space, for an SNR=1/sqrt(2) signal.
690  */
691 gsl_matrix *XLALSimInspiralTaylorF2RedSpinFisherMatrixChirpTimes(
692  const REAL8 theta0, /**< dimensionless parameter related to the chirp time by theta0 = 2 pi fLow tau0 */
693  const REAL8 theta3, /**< dimensionless parameter related to the chirp time by theta3 = -2 pi fLow tau3 */
694  const REAL8 theta3s, /**< dimensionless parameter related to the chirp time by theta3s = 2 pi fLow tau3s */
695  const REAL8 fLow, /**< low-frequency cutoff (Hz) */
696  const REAL8 df, /**< frequency resolution of the noise moment vectors (Hz) */
697  REAL8Vector *momI_0, /**< noise moments: \f$momI_0(f) = \int_{f0}^f (f'/f0)^{(0-17)/3} df'\f$ */
698  REAL8Vector *momI_2, /**< noise moments: \f$momI_2(f) = \int_{f0}^f (f'/f0)^{(2-17)/3} df'\f$ */
699  REAL8Vector *momI_3, /**< noise moments: \f$momI_3(f) = \int_{f0}^f (f'/f0)^{(3-17)/3} df'\f$ */
700  REAL8Vector *momI_4, /**< noise moments: \f$momI_4(f) = \int_{f0}^f (f'/f0)^{(4-17)/3} df'\f$ */
701  REAL8Vector *momI_5, /**< noise moments: \f$momI_5(f) = \int_{f0}^f (f'/f0)^{(5-17)/3} df'\f$ */
702  REAL8Vector *momI_6, /**< noise moments: \f$momI_6(f) = \int_{f0}^f (f'/f0)^{(6-17)/3} df'\f$ */
703  REAL8Vector *momI_7, /**< noise moments: \f$momI_7(f) = \int_{f0}^f (f'/f0)^{(7-17)/3} df'\f$ */
704  REAL8Vector *momI_8, /**< noise moments: \f$momI_8(f) = \int_{f0}^f (f'/f0)^{(8-17)/3} df'\f$ */
705  REAL8Vector *momI_9, /**< noise moments: \f$momI_9(f) = \int_{f0}^f (f'/f0)^{(9-17)/3} df'\f$ */
706  REAL8Vector *momI_10, /**< noise moments: \f$momI_10(f) = \int_{f0}^f (f'/f0)^{(10-17)/3} df'\f$ */
707  REAL8Vector *momI_11, /**< noise moments: \f$momI_11(f) = \int_{f0}^f (f'/f0)^{(11-17)/3} df'\f$ */
708  REAL8Vector *momI_12, /**< noise moments: \f$momI_12(f) = \int_{f0}^f (f'/f0)^{(12-17)/3} df'\f$ */
709  REAL8Vector *momI_13, /**< noise moments: \f$momI_13(f) = \int_{f0}^f (f'/f0)^{(13-17)/3} df'\f$ */
710  REAL8Vector *momI_14, /**< noise moments: \f$momI_14(f) = \int_{f0}^f (f'/f0)^{(14-17)/3} df'\f$ */
711  REAL8Vector *momI_15, /**< noise moments: \f$momI_15(f) = \int_{f0}^f (f'/f0)^{(15-17)/3} df'\f$ */
712  REAL8Vector *momI_16, /**< noise moments: \f$momI_16(f) = \int_{f0}^f (f'/f0)^{(16-17)/3} df'\f$ */
713  REAL8Vector *momJ_5, /**< noise moments: \f$momJ_5(f) = \int_{f0}^f (f'/f0)^{(5-17)/3} log(f'/f0) df'\f$ */
714  REAL8Vector *momJ_6, /**< noise moments: \f$momJ_6(f) = \int_{f0}^f (f'/f0)^{(6-17)/3} log(f'/f0) df'\f$ */
715  REAL8Vector *momJ_7, /**< noise moments: \f$momJ_7(f) = \int_{f0}^f (f'/f0)^{(7-17)/3} log(f'/f0) df'\f$ */
716  REAL8Vector *momJ_8, /**< noise moments: \f$momJ_8(f) = \int_{f0}^f (f'/f0)^{(8-17)/3} log(f'/f0) df'\f$ */
717  REAL8Vector *momJ_9, /**< noise moments: \f$momJ_9(f) = \int_{f0}^f (f'/f0)^{(9-17)/3} log(f'/f0) df'\f$ */
718  REAL8Vector *momJ_10, /**< noise moments: \f$momJ_10(f) = \int_{f0}^f (f'/f0)^{(10-17)/3} log(f'/f0) df'\f$ */
719  REAL8Vector *momJ_11, /**< noise moments: \f$momJ_11(f) = \int_{f0}^f (f'/f0)^{(11-17)/3} log(f'/f0) df'\f$ */
720  REAL8Vector *momJ_12, /**< noise moments: \f$momJ_12(f) = \int_{f0}^f (f'/f0)^{(12-17)/3} log(f'/f0) df'\f$ */
721  REAL8Vector *momJ_13, /**< noise moments: \f$momJ_13(f) = \int_{f0}^f (f'/f0)^{(13-17)/3} log(f'/f0) df'\f$ */
722  REAL8Vector *momJ_14, /**< noise moments: \f$momJ_14(f) = \int_{f0}^f (f'/f0)^{(14-17)/3} log(f'/f0) df'\f$ */
723  REAL8Vector *momK_10, /**< noise moments: \f$momK_14(f) = \int_{f0}^f (f'/f0)^{(14-17)/3} log(f'/f0) log(f'/f0) df'\f$ */
724  REAL8Vector *momK_11, /**< noise moments: \f$momK_15(f) = \int_{f0}^f (f'/f0)^{(15-17)/3} log(f'/f0) log(f'/f0) df'\f$ */
725  REAL8Vector *momK_12 /**< noise moments: \f$momK_16(f) = \int_{f0}^f (f'/f0)^{(16-17)/3} log(f'/f0) log(f'/f0) df'\f$ */
726 ) {
727 
728  if (theta0 <= 0) XLAL_ERROR_NULL(XLAL_EDOM);
729  if (theta3 <= 0) XLAL_ERROR_NULL(XLAL_EDOM);
730  if (fLow <= 0) XLAL_ERROR_NULL(XLAL_EDOM);
731  if (df <= 0) XLAL_ERROR_NULL(XLAL_EDOM);
732 
733  REAL8 theta0_p2 = theta0 * theta0;
734  REAL8 theta3_p2 = theta3 * theta3;
735  REAL8 theta3s_p2 = theta3s * theta3s;
736  REAL8 theta3_p3 = theta3_p2 * theta3;
737  REAL8 theta3_p4 = theta3_p3 * theta3;
738  REAL8 theta3_p5 = theta3_p4 * theta3;
739  REAL8 theta0_p1by3 = cbrt(theta0);
740  REAL8 theta3_p1by3 = cbrt(theta3);
741  REAL8 theta0_p2by3 = theta0 / theta0_p1by3;
742  REAL8 theta3_p2by3 = theta3 / theta3_p1by3;
743  REAL8 theta3_p5by3 = theta3_p2 / theta3_p1by3;
744  REAL8 theta0_p4by3 = theta0 * theta0_p1by3;
745  REAL8 theta3_p4by3 = theta3 * theta3_p1by3;
746  REAL8 theta0_p5by3 = theta0_p2 / theta0_p1by3;
747  REAL8 theta3_p8by3 = theta3_p3 / theta3_p1by3;
748  REAL8 theta0_p7by3 = theta0_p2 * theta0_p1by3;
749  REAL8 theta3_p10by3 = theta3_p3 * theta3_p1by3;
750  REAL8 theta3bytheta0_p2by3 = cbrt(theta3_p2/theta0_p2);
751 
752  /* harder to name combinations that come up a lot */
753  REAL8 theta03_fac = -152*2.15443469003188*Pi_p5by3*theta0_p2by3 + 565*theta3_p5by3;
754  REAL8 theta03_fac_p2 = theta03_fac * theta03_fac;
755  REAL8 theta03_fac_p3 = theta03_fac_p2 * theta03_fac;
756 
757  /* compute the frequency bin corresponding to the cutoff frequency
758  * (Schwarzschild test particle ISCO). Note that the first frequency
759  * bin correspond to a frequency of fLow */
760  REAL8 fISCO = (8*sqrt(0.6666666666666666)*fLow*LAL_PI*theta0)/(15.*theta3);
761  if (fISCO <= fLow) {
762  XLALPrintError("fISCO <= fLow");
763  XLAL_ERROR_NULL(XLAL_EDOM);
764  }
765  size_t iCut = (fISCO - fLow) / df;
766  if (iCut > momI_10->length - 1) {
767  XLALPrintWarning("moments do not cover fISCO (%g Hz); truncating to (%g Hz)", fISCO, (momI_10->length - 1) * df + fLow);
768  iCut = momI_10->length - 1;
769  }
770 
771  /* template norm */
772  REAL8 hSqr = 2.*momI_10->data[iCut];
773 
774  /*******************************************************************/
775  /* derivatives of the phasing coefficients w.r.t. parameters */
776  /*******************************************************************/
777  REAL8 dPhi0byDtheta0 = 0.6;
778 
779  REAL8 dPhi2byDtheta0 = (11088*LAL_PI + (743*tenbyPi_p2by3*theta3_p5by3)/theta0_p2by3)/(12096.*theta3);
780 
781  REAL8 dPhi4byDtheta0 = (5429*2.92401773821287*cbrt(LAL_PI/(2.*theta0_p2*theta3)))/16128. - (15293365*1.7099759466767)/(3.2514048e7*1.5874010519682*Pi_p4by3/theta3bytheta0_p2by3/theta3bytheta0_p2by3) + (617*Pi_p2)/(384.*theta3_p2) - (13596750*4.64158883361278*Pi_p7by3*theta3_p8by3*theta3s_p2)/(theta0_p2by3*theta03_fac_p3) + (91125*1.7099759466767*Pi_p2by3*theta3_p8by3*theta3s_p2)/(2.*1.5874010519682*theta0_p4by3*theta03_fac_p2) + (2052000*Pi_p4*theta3*theta3s_p2)/theta03_fac_p3;
782 
783  REAL8 dPhi5byDtheta0 = (5*1.7099759466767*theta3_p2by3*(-873377*2.15443469003188*theta3 + 2345552*2.15443469003188*theta3s + (344829849600*Pi_p10by3*theta0_p4by3*theta3s)/theta03_fac_p2))/(1.0934784e7*twoPi_p2by3*theta0_p5by3);
784 
785  REAL8 dPhi6byDtheta0 = (-888945361920*Pi_p5*theta0_p2 + 33057275328*4.64158883361278*Pi_p10by3*theta0_p4by3*theta3_p5by3 + 9694463631160*2.15443469003188*Pi_p5by3*theta0_p2by3*theta3_p10by3 - 352848545280*2.15443469003188*Pi_11by3*theta0_p2by3*theta3_p10by3 + 3*(-11072977443251 + 1530761379840*LAL_GAMMA)*theta3_p5 + 3004298035200*Pi_p2*theta3_p5 + 1530761379840*theta3_p5*log((10*theta3)/(LAL_PI*theta0)))/(9.61375371264e11*Pi_p2*theta0_p2*theta3_p3);
786 
787  REAL8 dPhi7byDtheta0 = (-5*theta3_p2by3*(12718104*4.64158883361278*Pi_p5by3*theta0_p2by3 + 77096675*2.15443469003188*theta3_p5by3))/(2.60112384e8*Pi_p4by3*theta0_p7by3);
788 
789  REAL8 dPhi2byDtheta3 = (-5544*LAL_PI*theta0 + 743*tenbyPi_p2by3*theta0_p1by3*theta3_p5by3)/(6048.*theta3_p2);
790 
791  REAL8 dPhi3byDtheta3 = -1.5;
792 
793  REAL8 dPhi4byDtheta3 = ((-52242624*Pi_p2*theta0_p4by3)/theta3_p3 - (2736216*4.64158883361278*Pi_p1by3*theta0_p2by3)/theta3_p4by3 + (15293365*2.15443469003188*theta3_p1by3)/Pi_p4by3 + (740710656000*2.15443469003188*Pi_p2by3*theta3_p5by3*(608*2.15443469003188*Pi_p5by3*theta0_p2by3 + 565*theta3_p5by3)*theta3s_p2)/theta03_fac_p3 - (7315660800*4.64158883361278*Pi_p7by3*theta0_p2by3*(456*2.15443469003188*Pi_p5by3*theta0_p2by3 + 3955*theta3_p5by3)*theta3s_p2)/theta03_fac_p3)/(1.6257024e7*theta0_p1by3);
794 
795  REAL8 dPhi5byDtheta3 = (5*(640*2.15443469003188*Pi_p10by3*theta0_p4by3*theta3_p5by3*(13950845*theta3 - 36141056*theta3s) - 4000*Pi_p5by3*theta0_p2by3*theta3_p10by3*(16594163*theta3 - 12773768*theta3s) + 2825*4.64158883361278*theta3_p5*(4366885*theta3 - 4691104*theta3s) + 1000194048*4.64158883361278*Pi_p5*theta0_p2*theta3s))/(387072.*Pi_p2by3*theta0_p2by3*pow(152*2.15443469003188*Pi_p5by3*theta0_p2by3*theta3 - 565*theta3_p8by3,2));
796 
797  REAL8 dPhi6byDtheta3 = (1333418042880*Pi_p5*theta0_p2 - 66114550656*4.64158883361278*Pi_p10by3*theta0_p4by3*theta3_p5by3 - 4847231815580*2.15443469003188*Pi_p5by3*theta0_p2by3*theta3_p10by3 + 176424272640*2.15443469003188*Pi_11by3*theta0_p2by3*theta3_p10by3 + 3*(11328104339891 - 1530761379840*LAL_GAMMA)*theta3_p5 - 3004298035200*Pi_p2*theta3_p5 - 1530761379840*theta3_p5*log((10*theta3)/(LAL_PI*theta0)))/(4.80687685632e11*Pi_p2*theta0*theta3_p4);
798 
799  REAL8 dPhi7byDtheta3 = (5*(17059968*Pi_p10by3 + (7267488*4.64158883361278*Pi_p5by3*theta3_p5by3)/theta0_p2by3 + (77096675*2.15443469003188*theta3_p10by3)/theta0_p4by3))/(1.48635648e8*Pi_p4by3*theta3_p2);
800 
801  REAL8 dPhi3byDtheta3s = 1.5;
802 
803  REAL8 dPhi4byDtheta3s = (675*theta3*(8*4.64158883361278*Pi_p7by3*theta0_p2by3 - 405*2.15443469003188*Pi_p2by3*theta3_p5by3)*theta3s)/(2.*theta0_p1by3*theta03_fac_p2);
804 
805  REAL8 dPhi5byDtheta3s = (5*(-2785343*tenbyPi_p2by3*theta3bytheta0_p2by3 + (816*LAL_PI*(-2373 + (3301450*theta3_p5by3)/theta03_fac))/theta3))/1.7313408e7;
806 
807  REAL8 dPhi8byDt0 = 2.*fLow*LAL_PI;
808 
809  REAL8 dPhi5byDphi0 = 1.;
810 
811  REAL8 dPhiL5byDtheta0 = (5*1.7099759466767*theta3_p2by3*(-873377*2.15443469003188*theta3 + 2345552*2.15443469003188*theta3s + (344829849600*Pi_p10by3*theta0_p4by3*theta3s)/theta03_fac_p2))/(1.0934784e7*twoPi_p2by3*theta0_p5by3);
812 
813  REAL8 dPhiL6byDtheta0 = (535*theta3_p2)/(336.*Pi_p2*theta0_p2);
814 
815  REAL8 dPhiL5byDtheta3 = (5*(640*2.15443469003188*Pi_p10by3*theta0_p4by3*theta3_p5by3*(13950845*theta3 - 36141056*theta3s) - 4000*Pi_p5by3*theta0_p2by3*theta3_p10by3*(16594163*theta3 - 12773768*theta3s) + 2825*4.64158883361278*theta3_p5*(4366885*theta3 - 4691104*theta3s) + 1000194048*4.64158883361278*Pi_p5*theta0_p2*theta3s))/(387072.*Pi_p2by3*theta0_p2by3*pow(152*2.15443469003188*Pi_p5by3*theta0_p2by3*theta3 - 565*theta3_p8by3,2));
816 
817  REAL8 dPhiL6byDtheta3 = (-535*theta3)/(168.*Pi_p2*theta0);
818 
819  REAL8 dPhiL5byDtheta3s = (5*(-2785343*tenbyPi_p2by3*theta3bytheta0_p2by3 + (816*LAL_PI*(-2373 + (3301450*theta3_p5by3)/theta03_fac))/theta3))/1.7313408e7;
820 
821  /*******************************************************************/
822  /* (theta0, theta0) component of the metric */
823  /*******************************************************************/
824 
825  REAL8 g_theta0theta0 = dPhi0byDtheta0 * dPhi0byDtheta0 * momI_0->data[iCut] +
826  dPhi0byDtheta0 * dPhi2byDtheta0 * momI_2->data[iCut] +
827  dPhi0byDtheta0 * dPhi4byDtheta0 * momI_4->data[iCut] +
828  dPhi0byDtheta0 * dPhi5byDtheta0 * momI_5->data[iCut] +
829  dPhi0byDtheta0 * dPhiL5byDtheta0 * momJ_5->data[iCut] +
830  dPhi0byDtheta0 * dPhi6byDtheta0 * momI_6->data[iCut] +
831  dPhi0byDtheta0 * dPhiL6byDtheta0 * momJ_6->data[iCut] +
832  dPhi0byDtheta0 * dPhi7byDtheta0 * momI_7->data[iCut] +
833  dPhi2byDtheta0 * dPhi0byDtheta0 * momI_2->data[iCut] +
834  dPhi2byDtheta0 * dPhi2byDtheta0 * momI_4->data[iCut] +
835  dPhi2byDtheta0 * dPhi4byDtheta0 * momI_6->data[iCut] +
836  dPhi2byDtheta0 * dPhi5byDtheta0 * momI_7->data[iCut] +
837  dPhi2byDtheta0 * dPhiL5byDtheta0 * momJ_7->data[iCut] +
838  dPhi2byDtheta0 * dPhi6byDtheta0 * momI_8->data[iCut] +
839  dPhi2byDtheta0 * dPhiL6byDtheta0 * momJ_8->data[iCut] +
840  dPhi2byDtheta0 * dPhi7byDtheta0 * momI_9->data[iCut] +
841  dPhi4byDtheta0 * dPhi0byDtheta0 * momI_4->data[iCut] +
842  dPhi4byDtheta0 * dPhi2byDtheta0 * momI_6->data[iCut] +
843  dPhi4byDtheta0 * dPhi4byDtheta0 * momI_8->data[iCut] +
844  dPhi4byDtheta0 * dPhi5byDtheta0 * momI_9->data[iCut] +
845  dPhi4byDtheta0 * dPhiL5byDtheta0 * momJ_9->data[iCut] +
846  dPhi4byDtheta0 * dPhi6byDtheta0 * momI_10->data[iCut] +
847  dPhi4byDtheta0 * dPhiL6byDtheta0 * momJ_10->data[iCut] +
848  dPhi4byDtheta0 * dPhi7byDtheta0 * momI_11->data[iCut] +
849  dPhi5byDtheta0 * dPhi0byDtheta0 * momI_5->data[iCut] +
850  dPhiL5byDtheta0 * dPhi0byDtheta0 * momJ_5->data[iCut] +
851  dPhi5byDtheta0 * dPhi2byDtheta0 * momI_7->data[iCut] +
852  dPhiL5byDtheta0 * dPhi2byDtheta0 * momJ_7->data[iCut] +
853  dPhi5byDtheta0 * dPhi4byDtheta0 * momI_9->data[iCut] +
854  dPhiL5byDtheta0 * dPhi4byDtheta0 * momJ_9->data[iCut] +
855  dPhi5byDtheta0 * dPhi5byDtheta0 * momI_10->data[iCut] +
856  dPhi5byDtheta0 * dPhiL5byDtheta0 * momJ_10->data[iCut] +
857  dPhiL5byDtheta0 * dPhi5byDtheta0 * momJ_10->data[iCut] +
858  dPhiL5byDtheta0 * dPhiL5byDtheta0 * momK_10->data[iCut] +
859  dPhi5byDtheta0 * dPhi6byDtheta0 * momI_11->data[iCut] +
860  dPhi5byDtheta0 * dPhiL6byDtheta0 * momJ_11->data[iCut] +
861  dPhiL5byDtheta0 * dPhi6byDtheta0 * momJ_11->data[iCut] +
862  dPhiL5byDtheta0 * dPhiL6byDtheta0 * momK_11->data[iCut] +
863  dPhi5byDtheta0 * dPhi7byDtheta0 * momI_12->data[iCut] +
864  dPhiL5byDtheta0 * dPhi7byDtheta0 * momJ_12->data[iCut] +
865  dPhi6byDtheta0 * dPhi0byDtheta0 * momI_6->data[iCut] +
866  dPhiL6byDtheta0 * dPhi0byDtheta0 * momJ_6->data[iCut] +
867  dPhi6byDtheta0 * dPhi2byDtheta0 * momI_8->data[iCut] +
868  dPhiL6byDtheta0 * dPhi2byDtheta0 * momJ_8->data[iCut] +
869  dPhi6byDtheta0 * dPhi4byDtheta0 * momI_10->data[iCut] +
870  dPhiL6byDtheta0 * dPhi4byDtheta0 * momJ_10->data[iCut] +
871  dPhi6byDtheta0 * dPhi5byDtheta0 * momI_11->data[iCut] +
872  dPhi6byDtheta0 * dPhiL5byDtheta0 * momJ_11->data[iCut] +
873  dPhiL6byDtheta0 * dPhi5byDtheta0 * momJ_11->data[iCut] +
874  dPhiL6byDtheta0 * dPhiL5byDtheta0 * momK_11->data[iCut] +
875  dPhi6byDtheta0 * dPhi6byDtheta0 * momI_12->data[iCut] +
876  dPhi6byDtheta0 * dPhiL6byDtheta0 * momJ_12->data[iCut] +
877  dPhiL6byDtheta0 * dPhi6byDtheta0 * momJ_12->data[iCut] +
878  dPhiL6byDtheta0 * dPhiL6byDtheta0 * momK_12->data[iCut] +
879  dPhi6byDtheta0 * dPhi7byDtheta0 * momI_13->data[iCut] +
880  dPhiL6byDtheta0 * dPhi7byDtheta0 * momJ_13->data[iCut] +
881  dPhi7byDtheta0 * dPhi0byDtheta0 * momI_7->data[iCut] +
882  dPhi7byDtheta0 * dPhi2byDtheta0 * momI_9->data[iCut] +
883  dPhi7byDtheta0 * dPhi4byDtheta0 * momI_11->data[iCut] +
884  dPhi7byDtheta0 * dPhi5byDtheta0 * momI_12->data[iCut] +
885  dPhi7byDtheta0 * dPhiL5byDtheta0 * momJ_12->data[iCut] +
886  dPhi7byDtheta0 * dPhi6byDtheta0 * momI_13->data[iCut] +
887  dPhi7byDtheta0 * dPhiL6byDtheta0 * momJ_13->data[iCut] +
888  dPhi7byDtheta0 * dPhi7byDtheta0 * momI_14->data[iCut] +
889  0. ;
890 
891  /*******************************************************************/
892  /* (theta0, theta3) component of the metric */
893  /*******************************************************************/
894 
895  REAL8 g_theta0theta3 = dPhi0byDtheta0 * dPhi2byDtheta3 * momI_2->data[iCut] +
896  dPhi0byDtheta0 * dPhi3byDtheta3 * momI_3->data[iCut] +
897  dPhi0byDtheta0 * dPhi4byDtheta3 * momI_4->data[iCut] +
898  dPhi0byDtheta0 * dPhi5byDtheta3 * momI_5->data[iCut] +
899  dPhi0byDtheta0 * dPhiL5byDtheta3 * momJ_5->data[iCut] +
900  dPhi0byDtheta0 * dPhi6byDtheta3 * momI_6->data[iCut] +
901  dPhi0byDtheta0 * dPhiL6byDtheta3 * momJ_6->data[iCut] +
902  dPhi0byDtheta0 * dPhi7byDtheta3 * momI_7->data[iCut] +
903  dPhi2byDtheta0 * dPhi2byDtheta3 * momI_4->data[iCut] +
904  dPhi2byDtheta0 * dPhi3byDtheta3 * momI_5->data[iCut] +
905  dPhi2byDtheta0 * dPhi4byDtheta3 * momI_6->data[iCut] +
906  dPhi2byDtheta0 * dPhi5byDtheta3 * momI_7->data[iCut] +
907  dPhi2byDtheta0 * dPhiL5byDtheta3 * momJ_7->data[iCut] +
908  dPhi2byDtheta0 * dPhi6byDtheta3 * momI_8->data[iCut] +
909  dPhi2byDtheta0 * dPhiL6byDtheta3 * momJ_8->data[iCut] +
910  dPhi2byDtheta0 * dPhi7byDtheta3 * momI_9->data[iCut] +
911  dPhi4byDtheta0 * dPhi2byDtheta3 * momI_6->data[iCut] +
912  dPhi4byDtheta0 * dPhi3byDtheta3 * momI_7->data[iCut] +
913  dPhi4byDtheta0 * dPhi4byDtheta3 * momI_8->data[iCut] +
914  dPhi4byDtheta0 * dPhi5byDtheta3 * momI_9->data[iCut] +
915  dPhi4byDtheta0 * dPhiL5byDtheta3 * momJ_9->data[iCut] +
916  dPhi4byDtheta0 * dPhi6byDtheta3 * momI_10->data[iCut] +
917  dPhi4byDtheta0 * dPhiL6byDtheta3 * momJ_10->data[iCut] +
918  dPhi4byDtheta0 * dPhi7byDtheta3 * momI_11->data[iCut] +
919  dPhi5byDtheta0 * dPhi2byDtheta3 * momI_7->data[iCut] +
920  dPhiL5byDtheta0 * dPhi2byDtheta3 * momJ_7->data[iCut] +
921  dPhi5byDtheta0 * dPhi3byDtheta3 * momI_8->data[iCut] +
922  dPhiL5byDtheta0 * dPhi3byDtheta3 * momJ_8->data[iCut] +
923  dPhi5byDtheta0 * dPhi4byDtheta3 * momI_9->data[iCut] +
924  dPhiL5byDtheta0 * dPhi4byDtheta3 * momJ_9->data[iCut] +
925  dPhi5byDtheta0 * dPhi5byDtheta3 * momI_10->data[iCut] +
926  dPhi5byDtheta0 * dPhiL5byDtheta3 * momJ_10->data[iCut] +
927  dPhiL5byDtheta0 * dPhi5byDtheta3 * momJ_10->data[iCut] +
928  dPhiL5byDtheta0 * dPhiL5byDtheta3 * momK_10->data[iCut] +
929  dPhi5byDtheta0 * dPhi6byDtheta3 * momI_11->data[iCut] +
930  dPhi5byDtheta0 * dPhiL6byDtheta3 * momJ_11->data[iCut] +
931  dPhiL5byDtheta0 * dPhi6byDtheta3 * momJ_11->data[iCut] +
932  dPhiL5byDtheta0 * dPhiL6byDtheta3 * momK_11->data[iCut] +
933  dPhi5byDtheta0 * dPhi7byDtheta3 * momI_12->data[iCut] +
934  dPhiL5byDtheta0 * dPhi7byDtheta3 * momJ_12->data[iCut] +
935  dPhi6byDtheta0 * dPhi2byDtheta3 * momI_8->data[iCut] +
936  dPhiL6byDtheta0 * dPhi2byDtheta3 * momJ_8->data[iCut] +
937  dPhi6byDtheta0 * dPhi3byDtheta3 * momI_9->data[iCut] +
938  dPhiL6byDtheta0 * dPhi3byDtheta3 * momJ_9->data[iCut] +
939  dPhi6byDtheta0 * dPhi4byDtheta3 * momI_10->data[iCut] +
940  dPhiL6byDtheta0 * dPhi4byDtheta3 * momJ_10->data[iCut] +
941  dPhi6byDtheta0 * dPhi5byDtheta3 * momI_11->data[iCut] +
942  dPhi6byDtheta0 * dPhiL5byDtheta3 * momJ_11->data[iCut] +
943  dPhiL6byDtheta0 * dPhi5byDtheta3 * momJ_11->data[iCut] +
944  dPhiL6byDtheta0 * dPhiL5byDtheta3 * momK_11->data[iCut] +
945  dPhi6byDtheta0 * dPhi6byDtheta3 * momI_12->data[iCut] +
946  dPhi6byDtheta0 * dPhiL6byDtheta3 * momJ_12->data[iCut] +
947  dPhiL6byDtheta0 * dPhi6byDtheta3 * momJ_12->data[iCut] +
948  dPhiL6byDtheta0 * dPhiL6byDtheta3 * momK_12->data[iCut] +
949  dPhi6byDtheta0 * dPhi7byDtheta3 * momI_13->data[iCut] +
950  dPhiL6byDtheta0 * dPhi7byDtheta3 * momJ_13->data[iCut] +
951  dPhi7byDtheta0 * dPhi2byDtheta3 * momI_9->data[iCut] +
952  dPhi7byDtheta0 * dPhi3byDtheta3 * momI_10->data[iCut] +
953  dPhi7byDtheta0 * dPhi4byDtheta3 * momI_11->data[iCut] +
954  dPhi7byDtheta0 * dPhi5byDtheta3 * momI_12->data[iCut] +
955  dPhi7byDtheta0 * dPhiL5byDtheta3 * momJ_12->data[iCut] +
956  dPhi7byDtheta0 * dPhi6byDtheta3 * momI_13->data[iCut] +
957  dPhi7byDtheta0 * dPhiL6byDtheta3 * momJ_13->data[iCut] +
958  dPhi7byDtheta0 * dPhi7byDtheta3 * momI_14->data[iCut] +
959  0. ;
960 
961  /*******************************************************************/
962  /* (theta0, theta3s) component of the metric */
963  /*******************************************************************/
964 
965  REAL8 g_theta0theta3s = dPhi0byDtheta0 * dPhi3byDtheta3s * momI_3->data[iCut] +
966  dPhi0byDtheta0 * dPhi4byDtheta3s * momI_4->data[iCut] +
967  dPhi0byDtheta0 * dPhi5byDtheta3s * momI_5->data[iCut] +
968  dPhi0byDtheta0 * dPhiL5byDtheta3s * momJ_5->data[iCut] +
969  dPhi2byDtheta0 * dPhi3byDtheta3s * momI_5->data[iCut] +
970  dPhi2byDtheta0 * dPhi4byDtheta3s * momI_6->data[iCut] +
971  dPhi2byDtheta0 * dPhi5byDtheta3s * momI_7->data[iCut] +
972  dPhi2byDtheta0 * dPhiL5byDtheta3s * momJ_7->data[iCut] +
973  dPhi4byDtheta0 * dPhi3byDtheta3s * momI_7->data[iCut] +
974  dPhi4byDtheta0 * dPhi4byDtheta3s * momI_8->data[iCut] +
975  dPhi4byDtheta0 * dPhi5byDtheta3s * momI_9->data[iCut] +
976  dPhi4byDtheta0 * dPhiL5byDtheta3s * momJ_9->data[iCut] +
977  dPhi5byDtheta0 * dPhi3byDtheta3s * momI_8->data[iCut] +
978  dPhiL5byDtheta0 * dPhi3byDtheta3s * momJ_8->data[iCut] +
979  dPhi5byDtheta0 * dPhi4byDtheta3s * momI_9->data[iCut] +
980  dPhiL5byDtheta0 * dPhi4byDtheta3s * momJ_9->data[iCut] +
981  dPhi5byDtheta0 * dPhi5byDtheta3s * momI_10->data[iCut] +
982  dPhi5byDtheta0 * dPhiL5byDtheta3s * momJ_10->data[iCut] +
983  dPhiL5byDtheta0 * dPhi5byDtheta3s * momJ_10->data[iCut] +
984  dPhiL5byDtheta0 * dPhiL5byDtheta3s * momK_10->data[iCut] +
985  dPhi6byDtheta0 * dPhi3byDtheta3s * momI_9->data[iCut] +
986  dPhiL6byDtheta0 * dPhi3byDtheta3s * momJ_9->data[iCut] +
987  dPhi6byDtheta0 * dPhi4byDtheta3s * momI_10->data[iCut] +
988  dPhiL6byDtheta0 * dPhi4byDtheta3s * momJ_10->data[iCut] +
989  dPhi6byDtheta0 * dPhi5byDtheta3s * momI_11->data[iCut] +
990  dPhi6byDtheta0 * dPhiL5byDtheta3s * momJ_11->data[iCut] +
991  dPhiL6byDtheta0 * dPhi5byDtheta3s * momJ_11->data[iCut] +
992  dPhiL6byDtheta0 * dPhiL5byDtheta3s * momK_11->data[iCut] +
993  dPhi7byDtheta0 * dPhi3byDtheta3s * momI_10->data[iCut] +
994  dPhi7byDtheta0 * dPhi4byDtheta3s * momI_11->data[iCut] +
995  dPhi7byDtheta0 * dPhi5byDtheta3s * momI_12->data[iCut] +
996  dPhi7byDtheta0 * dPhiL5byDtheta3s * momJ_12->data[iCut] +
997  0. ;
998 
999  /*******************************************************************/
1000  /* (theta0, t0) component of the metric */
1001  /*******************************************************************/
1002 
1003  REAL8 g_theta0t0 = dPhi0byDtheta0 * dPhi8byDt0 * momI_8->data[iCut] +
1004  dPhi2byDtheta0 * dPhi8byDt0 * momI_10->data[iCut] +
1005  dPhi4byDtheta0 * dPhi8byDt0 * momI_12->data[iCut] +
1006  dPhi5byDtheta0 * dPhi8byDt0 * momI_13->data[iCut] +
1007  dPhiL5byDtheta0 * dPhi8byDt0 * momJ_13->data[iCut] +
1008  dPhi6byDtheta0 * dPhi8byDt0 * momI_14->data[iCut] +
1009  dPhiL6byDtheta0 * dPhi8byDt0 * momJ_14->data[iCut] +
1010  dPhi7byDtheta0 * dPhi8byDt0 * momI_15->data[iCut] +
1011  0. ;
1012 
1013  /*******************************************************************/
1014  /* (theta0, phi0) component of the metric */
1015  /*******************************************************************/
1016 
1017  REAL8 g_theta0phi0 = dPhi0byDtheta0 * dPhi5byDphi0 * momI_5->data[iCut] +
1018  dPhi2byDtheta0 * dPhi5byDphi0 * momI_7->data[iCut] +
1019  dPhi4byDtheta0 * dPhi5byDphi0 * momI_9->data[iCut] +
1020  dPhi5byDtheta0 * dPhi5byDphi0 * momI_10->data[iCut] +
1021  dPhiL5byDtheta0 * dPhi5byDphi0 * momJ_10->data[iCut] +
1022  dPhi6byDtheta0 * dPhi5byDphi0 * momI_11->data[iCut] +
1023  dPhiL6byDtheta0 * dPhi5byDphi0 * momJ_11->data[iCut] +
1024  dPhi7byDtheta0 * dPhi5byDphi0 * momI_12->data[iCut] +
1025  0. ;
1026 
1027  /*******************************************************************/
1028  /* (theta3, theta3) component of the metric */
1029  /*******************************************************************/
1030 
1031  REAL8 g_theta3theta3 = dPhi2byDtheta3 * dPhi2byDtheta3 * momI_4->data[iCut] +
1032  dPhi2byDtheta3 * dPhi3byDtheta3 * momI_5->data[iCut] +
1033  dPhi2byDtheta3 * dPhi4byDtheta3 * momI_6->data[iCut] +
1034  dPhi2byDtheta3 * dPhi5byDtheta3 * momI_7->data[iCut] +
1035  dPhi2byDtheta3 * dPhiL5byDtheta3 * momJ_7->data[iCut] +
1036  dPhi2byDtheta3 * dPhi6byDtheta3 * momI_8->data[iCut] +
1037  dPhi2byDtheta3 * dPhiL6byDtheta3 * momJ_8->data[iCut] +
1038  dPhi2byDtheta3 * dPhi7byDtheta3 * momI_9->data[iCut] +
1039  dPhi3byDtheta3 * dPhi2byDtheta3 * momI_5->data[iCut] +
1040  dPhi3byDtheta3 * dPhi3byDtheta3 * momI_6->data[iCut] +
1041  dPhi3byDtheta3 * dPhi4byDtheta3 * momI_7->data[iCut] +
1042  dPhi3byDtheta3 * dPhi5byDtheta3 * momI_8->data[iCut] +
1043  dPhi3byDtheta3 * dPhiL5byDtheta3 * momJ_8->data[iCut] +
1044  dPhi3byDtheta3 * dPhi6byDtheta3 * momI_9->data[iCut] +
1045  dPhi3byDtheta3 * dPhiL6byDtheta3 * momJ_9->data[iCut] +
1046  dPhi3byDtheta3 * dPhi7byDtheta3 * momI_10->data[iCut] +
1047  dPhi4byDtheta3 * dPhi2byDtheta3 * momI_6->data[iCut] +
1048  dPhi4byDtheta3 * dPhi3byDtheta3 * momI_7->data[iCut] +
1049  dPhi4byDtheta3 * dPhi4byDtheta3 * momI_8->data[iCut] +
1050  dPhi4byDtheta3 * dPhi5byDtheta3 * momI_9->data[iCut] +
1051  dPhi4byDtheta3 * dPhiL5byDtheta3 * momJ_9->data[iCut] +
1052  dPhi4byDtheta3 * dPhi6byDtheta3 * momI_10->data[iCut] +
1053  dPhi4byDtheta3 * dPhiL6byDtheta3 * momJ_10->data[iCut] +
1054  dPhi4byDtheta3 * dPhi7byDtheta3 * momI_11->data[iCut] +
1055  dPhi5byDtheta3 * dPhi2byDtheta3 * momI_7->data[iCut] +
1056  dPhiL5byDtheta3 * dPhi2byDtheta3 * momJ_7->data[iCut] +
1057  dPhi5byDtheta3 * dPhi3byDtheta3 * momI_8->data[iCut] +
1058  dPhiL5byDtheta3 * dPhi3byDtheta3 * momJ_8->data[iCut] +
1059  dPhi5byDtheta3 * dPhi4byDtheta3 * momI_9->data[iCut] +
1060  dPhiL5byDtheta3 * dPhi4byDtheta3 * momJ_9->data[iCut] +
1061  dPhi5byDtheta3 * dPhi5byDtheta3 * momI_10->data[iCut] +
1062  dPhi5byDtheta3 * dPhiL5byDtheta3 * momJ_10->data[iCut] +
1063  dPhiL5byDtheta3 * dPhi5byDtheta3 * momJ_10->data[iCut] +
1064  dPhiL5byDtheta3 * dPhiL5byDtheta3 * momK_10->data[iCut] +
1065  dPhi5byDtheta3 * dPhi6byDtheta3 * momI_11->data[iCut] +
1066  dPhi5byDtheta3 * dPhiL6byDtheta3 * momJ_11->data[iCut] +
1067  dPhiL5byDtheta3 * dPhi6byDtheta3 * momJ_11->data[iCut] +
1068  dPhiL5byDtheta3 * dPhiL6byDtheta3 * momK_11->data[iCut] +
1069  dPhi5byDtheta3 * dPhi7byDtheta3 * momI_12->data[iCut] +
1070  dPhiL5byDtheta3 * dPhi7byDtheta3 * momJ_12->data[iCut] +
1071  dPhi6byDtheta3 * dPhi2byDtheta3 * momI_8->data[iCut] +
1072  dPhiL6byDtheta3 * dPhi2byDtheta3 * momJ_8->data[iCut] +
1073  dPhi6byDtheta3 * dPhi3byDtheta3 * momI_9->data[iCut] +
1074  dPhiL6byDtheta3 * dPhi3byDtheta3 * momJ_9->data[iCut] +
1075  dPhi6byDtheta3 * dPhi4byDtheta3 * momI_10->data[iCut] +
1076  dPhiL6byDtheta3 * dPhi4byDtheta3 * momJ_10->data[iCut] +
1077  dPhi6byDtheta3 * dPhi5byDtheta3 * momI_11->data[iCut] +
1078  dPhi6byDtheta3 * dPhiL5byDtheta3 * momJ_11->data[iCut] +
1079  dPhiL6byDtheta3 * dPhi5byDtheta3 * momJ_11->data[iCut] +
1080  dPhiL6byDtheta3 * dPhiL5byDtheta3 * momK_11->data[iCut] +
1081  dPhi6byDtheta3 * dPhi6byDtheta3 * momI_12->data[iCut] +
1082  dPhi6byDtheta3 * dPhiL6byDtheta3 * momJ_12->data[iCut] +
1083  dPhiL6byDtheta3 * dPhi6byDtheta3 * momJ_12->data[iCut] +
1084  dPhiL6byDtheta3 * dPhiL6byDtheta3 * momK_12->data[iCut] +
1085  dPhi6byDtheta3 * dPhi7byDtheta3 * momI_13->data[iCut] +
1086  dPhiL6byDtheta3 * dPhi7byDtheta3 * momJ_13->data[iCut] +
1087  dPhi7byDtheta3 * dPhi2byDtheta3 * momI_9->data[iCut] +
1088  dPhi7byDtheta3 * dPhi3byDtheta3 * momI_10->data[iCut] +
1089  dPhi7byDtheta3 * dPhi4byDtheta3 * momI_11->data[iCut] +
1090  dPhi7byDtheta3 * dPhi5byDtheta3 * momI_12->data[iCut] +
1091  dPhi7byDtheta3 * dPhiL5byDtheta3 * momJ_12->data[iCut] +
1092  dPhi7byDtheta3 * dPhi6byDtheta3 * momI_13->data[iCut] +
1093  dPhi7byDtheta3 * dPhiL6byDtheta3 * momJ_13->data[iCut] +
1094  dPhi7byDtheta3 * dPhi7byDtheta3 * momI_14->data[iCut] +
1095  0. ;
1096 
1097  /*******************************************************************/
1098  /* (theta3, theta3s) component of the metric */
1099  /*******************************************************************/
1100 
1101  REAL8 g_theta3theta3s = dPhi2byDtheta3 * dPhi3byDtheta3s * momI_5->data[iCut] +
1102  dPhi2byDtheta3 * dPhi4byDtheta3s * momI_6->data[iCut] +
1103  dPhi2byDtheta3 * dPhi5byDtheta3s * momI_7->data[iCut] +
1104  dPhi2byDtheta3 * dPhiL5byDtheta3s * momJ_7->data[iCut] +
1105  dPhi3byDtheta3 * dPhi3byDtheta3s * momI_6->data[iCut] +
1106  dPhi3byDtheta3 * dPhi4byDtheta3s * momI_7->data[iCut] +
1107  dPhi3byDtheta3 * dPhi5byDtheta3s * momI_8->data[iCut] +
1108  dPhi3byDtheta3 * dPhiL5byDtheta3s * momJ_8->data[iCut] +
1109  dPhi4byDtheta3 * dPhi3byDtheta3s * momI_7->data[iCut] +
1110  dPhi4byDtheta3 * dPhi4byDtheta3s * momI_8->data[iCut] +
1111  dPhi4byDtheta3 * dPhi5byDtheta3s * momI_9->data[iCut] +
1112  dPhi4byDtheta3 * dPhiL5byDtheta3s * momJ_9->data[iCut] +
1113  dPhi5byDtheta3 * dPhi3byDtheta3s * momI_8->data[iCut] +
1114  dPhiL5byDtheta3 * dPhi3byDtheta3s * momJ_8->data[iCut] +
1115  dPhi5byDtheta3 * dPhi4byDtheta3s * momI_9->data[iCut] +
1116  dPhiL5byDtheta3 * dPhi4byDtheta3s * momJ_9->data[iCut] +
1117  dPhi5byDtheta3 * dPhi5byDtheta3s * momI_10->data[iCut] +
1118  dPhi5byDtheta3 * dPhiL5byDtheta3s * momJ_10->data[iCut] +
1119  dPhiL5byDtheta3 * dPhi5byDtheta3s * momJ_10->data[iCut] +
1120  dPhiL5byDtheta3 * dPhiL5byDtheta3s * momK_10->data[iCut] +
1121  dPhi6byDtheta3 * dPhi3byDtheta3s * momI_9->data[iCut] +
1122  dPhiL6byDtheta3 * dPhi3byDtheta3s * momJ_9->data[iCut] +
1123  dPhi6byDtheta3 * dPhi4byDtheta3s * momI_10->data[iCut] +
1124  dPhiL6byDtheta3 * dPhi4byDtheta3s * momJ_10->data[iCut] +
1125  dPhi6byDtheta3 * dPhi5byDtheta3s * momI_11->data[iCut] +
1126  dPhi6byDtheta3 * dPhiL5byDtheta3s * momJ_11->data[iCut] +
1127  dPhiL6byDtheta3 * dPhi5byDtheta3s * momJ_11->data[iCut] +
1128  dPhiL6byDtheta3 * dPhiL5byDtheta3s * momK_11->data[iCut] +
1129  dPhi7byDtheta3 * dPhi3byDtheta3s * momI_10->data[iCut] +
1130  dPhi7byDtheta3 * dPhi4byDtheta3s * momI_11->data[iCut] +
1131  dPhi7byDtheta3 * dPhi5byDtheta3s * momI_12->data[iCut] +
1132  dPhi7byDtheta3 * dPhiL5byDtheta3s * momJ_12->data[iCut] +
1133  0. ;
1134 
1135  /*******************************************************************/
1136  /* (theta3, t0) component of the metric */
1137  /*******************************************************************/
1138 
1139  REAL8 g_theta3t0 = dPhi2byDtheta3 * dPhi8byDt0 * momI_10->data[iCut] +
1140  dPhi3byDtheta3 * dPhi8byDt0 * momI_11->data[iCut] +
1141  dPhi4byDtheta3 * dPhi8byDt0 * momI_12->data[iCut] +
1142  dPhi5byDtheta3 * dPhi8byDt0 * momI_13->data[iCut] +
1143  dPhiL5byDtheta3 * dPhi8byDt0 * momJ_13->data[iCut] +
1144  dPhi6byDtheta3 * dPhi8byDt0 * momI_14->data[iCut] +
1145  dPhiL6byDtheta3 * dPhi8byDt0 * momJ_14->data[iCut] +
1146  dPhi7byDtheta3 * dPhi8byDt0 * momI_15->data[iCut] +
1147  0. ;
1148 
1149  /*******************************************************************/
1150  /* (theta3, phi0) component of the metric */
1151  /*******************************************************************/
1152 
1153  REAL8 g_theta3phi0 = dPhi2byDtheta3 * dPhi5byDphi0 * momI_7->data[iCut] +
1154  dPhi3byDtheta3 * dPhi5byDphi0 * momI_8->data[iCut] +
1155  dPhi4byDtheta3 * dPhi5byDphi0 * momI_9->data[iCut] +
1156  dPhi5byDtheta3 * dPhi5byDphi0 * momI_10->data[iCut] +
1157  dPhiL5byDtheta3 * dPhi5byDphi0 * momJ_10->data[iCut] +
1158  dPhi6byDtheta3 * dPhi5byDphi0 * momI_11->data[iCut] +
1159  dPhiL6byDtheta3 * dPhi5byDphi0 * momJ_11->data[iCut] +
1160  dPhi7byDtheta3 * dPhi5byDphi0 * momI_12->data[iCut] +
1161  0. ;
1162 
1163  /*******************************************************************/
1164  /* (theta3s, theta3s) component of the metric */
1165  /*******************************************************************/
1166 
1167  REAL8 g_theta3stheta3s = dPhi3byDtheta3s * dPhi3byDtheta3s * momI_6->data[iCut] +
1168  dPhi3byDtheta3s * dPhi4byDtheta3s * momI_7->data[iCut] +
1169  dPhi3byDtheta3s * dPhi5byDtheta3s * momI_8->data[iCut] +
1170  dPhi3byDtheta3s * dPhiL5byDtheta3s * momJ_8->data[iCut] +
1171  dPhi4byDtheta3s * dPhi3byDtheta3s * momI_7->data[iCut] +
1172  dPhi4byDtheta3s * dPhi4byDtheta3s * momI_8->data[iCut] +
1173  dPhi4byDtheta3s * dPhi5byDtheta3s * momI_9->data[iCut] +
1174  dPhi4byDtheta3s * dPhiL5byDtheta3s * momJ_9->data[iCut] +
1175  dPhi5byDtheta3s * dPhi3byDtheta3s * momI_8->data[iCut] +
1176  dPhiL5byDtheta3s * dPhi3byDtheta3s * momJ_8->data[iCut] +
1177  dPhi5byDtheta3s * dPhi4byDtheta3s * momI_9->data[iCut] +
1178  dPhiL5byDtheta3s * dPhi4byDtheta3s * momJ_9->data[iCut] +
1179  dPhi5byDtheta3s * dPhi5byDtheta3s * momI_10->data[iCut] +
1180  dPhi5byDtheta3s * dPhiL5byDtheta3s * momJ_10->data[iCut] +
1181  dPhiL5byDtheta3s * dPhi5byDtheta3s * momJ_10->data[iCut] +
1182  dPhiL5byDtheta3s * dPhiL5byDtheta3s * momK_10->data[iCut] +
1183  0. ;
1184 
1185  /*******************************************************************/
1186  /* (theta3s, t0) component of the metric */
1187  /*******************************************************************/
1188 
1189  REAL8 g_theta3st0 = dPhi3byDtheta3s * dPhi8byDt0 * momI_11->data[iCut] +
1190  dPhi4byDtheta3s * dPhi8byDt0 * momI_12->data[iCut] +
1191  dPhi5byDtheta3s * dPhi8byDt0 * momI_13->data[iCut] +
1192  dPhiL5byDtheta3s * dPhi8byDt0 * momJ_13->data[iCut] +
1193  0. ;
1194 
1195  /*******************************************************************/
1196  /* (theta3s, phi0) component of the metric */
1197  /*******************************************************************/
1198 
1199  REAL8 g_theta3sphi0 = dPhi3byDtheta3s * dPhi5byDphi0 * momI_8->data[iCut] +
1200  dPhi4byDtheta3s * dPhi5byDphi0 * momI_9->data[iCut] +
1201  dPhi5byDtheta3s * dPhi5byDphi0 * momI_10->data[iCut] +
1202  dPhiL5byDtheta3s * dPhi5byDphi0 * momJ_10->data[iCut] +
1203  0. ;
1204 
1205  /*******************************************************************/
1206  /* (t0, t0) component of the metric */
1207  /*******************************************************************/
1208 
1209  REAL8 g_t0t0 = dPhi8byDt0 * dPhi8byDt0 * momI_16->data[iCut] +
1210  0. ;
1211 
1212  /*******************************************************************/
1213  /* (t0, phi0) component of the metric */
1214  /*******************************************************************/
1215 
1216  REAL8 g_t0phi0 = dPhi8byDt0 * dPhi5byDphi0 * momI_13->data[iCut] +
1217  0. ;
1218 
1219  /*******************************************************************/
1220  /* (phi0, phi0) component of the metric */
1221  /*******************************************************************/
1222 
1223  REAL8 g_phi0phi0 = dPhi5byDphi0 * dPhi5byDphi0 * momI_10->data[iCut] +
1224  0. ;
1225 
1226  ////////////////////////////////////////////////////////////////////////////////////////
1227  ////////////////////////////////////////////////////////////////////////////////////////
1228 
1229  /* allocate memory, and initialize the Fisher matrix */
1230  gsl_matrix * g = gsl_matrix_calloc (5, 5);
1231 
1232  /* compute the components of the Fisher matrix in coordinates theta0, theta3, theta3s, t0, phi0 */
1233  gsl_matrix_set (g, 0,0, g_theta0theta0/hSqr);
1234  gsl_matrix_set (g, 0,1, g_theta0theta3/hSqr);
1235  gsl_matrix_set (g, 0,2, g_theta0theta3s/hSqr);
1236  gsl_matrix_set (g, 0,3, g_theta0t0/hSqr);
1237  gsl_matrix_set (g, 0,4, g_theta0phi0/hSqr);
1238 
1239  gsl_matrix_set (g, 1,0, gsl_matrix_get(g, 0,1));
1240  gsl_matrix_set (g, 1,1, g_theta3theta3/hSqr);
1241  gsl_matrix_set (g, 1,2, g_theta3theta3s/hSqr);
1242  gsl_matrix_set (g, 1,3, g_theta3t0/hSqr);
1243  gsl_matrix_set (g, 1,4, g_theta3phi0/hSqr);
1244 
1245  gsl_matrix_set (g, 2,0, gsl_matrix_get(g, 0,2));
1246  gsl_matrix_set (g, 2,1, gsl_matrix_get(g, 1,2));
1247  gsl_matrix_set (g, 2,2, g_theta3stheta3s/hSqr);
1248  gsl_matrix_set (g, 2,3, g_theta3st0/hSqr);
1249  gsl_matrix_set (g, 2,4, g_theta3sphi0/hSqr);
1250 
1251  gsl_matrix_set (g, 3,0, gsl_matrix_get(g, 0,3));
1252  gsl_matrix_set (g, 3,1, gsl_matrix_get(g, 1,3));
1253  gsl_matrix_set (g, 3,2, gsl_matrix_get(g, 2,3));
1254  gsl_matrix_set (g, 3,3, g_t0t0/hSqr);
1255  gsl_matrix_set (g, 3,4, g_t0phi0/hSqr);
1256 
1257  gsl_matrix_set (g, 4,0, gsl_matrix_get(g, 0,4));
1258  gsl_matrix_set (g, 4,1, gsl_matrix_get(g, 1,4));
1259  gsl_matrix_set (g, 4,2, gsl_matrix_get(g, 2,4));
1260  gsl_matrix_set (g, 4,3, gsl_matrix_get(g, 3,4));
1261  gsl_matrix_set (g, 4,4, g_phi0phi0/hSqr);
1262 
1263  return g;
1264 }
1265 
1266 
1267 /**
1268  * Compute the template-space metric of "reduced-spin" PN templates in
1269  * theta0, theta3, theta3s parameter space.
1270  */
1271 int XLALSimInspiralTaylorF2RedSpinMetricChirpTimes(
1272  REAL8 *gamma00, /**< template metric coeff. 00 in theta0-theta3-theta3s*/
1273  REAL8 *gamma01, /**< template metric coeff. 01/10 in theta0-theta3-theta3s */
1274  REAL8 *gamma02, /**< template metric coeff. 02/20 in theta0-theta3-theta3s */
1275  REAL8 *gamma11, /**< template metric coeff. 11 in theta0-theta3-theta3s */
1276  REAL8 *gamma12, /**< template metric coeff. 12/21 in theta0-theta3-theta3s */
1277  REAL8 *gamma22, /**< template metric coeff. 22 in theta0-theta3-theta3s */
1278  const REAL8 theta0, /**< dimensionless parameter related to the chirp time by theta0 = 2 pi fLow tau0 */
1279  const REAL8 theta3, /**< dimensionless parameter related to the chirp time by theta3 = -2 pi fLow tau3 */
1280  const REAL8 theta3s, /**< dimensionless parameter related to the chirp time by theta3s = 2 pi fLow tau3s */
1281  const REAL8 fLow, /**< low-frequency cutoff (Hz) */
1282  const REAL8 df, /**< frequency resolution of the noise moment vectors (Hz) */
1283  REAL8Vector *momI_0, /**< noise moments: \f$momI_0(f) = \int_{f0}^f (f'/f0)^{(0-17)/3} df'\f$ */
1284  REAL8Vector *momI_2, /**< noise moments: \f$momI_2(f) = \int_{f0}^f (f'/f0)^{(2-17)/3} df'\f$ */
1285  REAL8Vector *momI_3, /**< noise moments: \f$momI_3(f) = \int_{f0}^f (f'/f0)^{(3-17)/3} df'\f$ */
1286  REAL8Vector *momI_4, /**< noise moments: \f$momI_4(f) = \int_{f0}^f (f'/f0)^{(4-17)/3} df'\f$ */
1287  REAL8Vector *momI_5, /**< noise moments: \f$momI_5(f) = \int_{f0}^f (f'/f0)^{(5-17)/3} df'\f$ */
1288  REAL8Vector *momI_6, /**< noise moments: \f$momI_6(f) = \int_{f0}^f (f'/f0)^{(6-17)/3} df'\f$ */
1289  REAL8Vector *momI_7, /**< noise moments: \f$momI_7(f) = \int_{f0}^f (f'/f0)^{(7-17)/3} df'\f$ */
1290  REAL8Vector *momI_8, /**< noise moments: \f$momI_8(f) = \int_{f0}^f (f'/f0)^{(8-17)/3} df'\f$ */
1291  REAL8Vector *momI_9, /**< noise moments: \f$momI_9(f) = \int_{f0}^f (f'/f0)^{(9-17)/3} df'\f$ */
1292  REAL8Vector *momI_10, /**< noise moments: \f$momI_10(f) = \int_{f0}^f (f'/f0)^{(10-17)/3} df'\f$ */
1293  REAL8Vector *momI_11, /**< noise moments: \f$momI_11(f) = \int_{f0}^f (f'/f0)^{(11-17)/3} df'\f$ */
1294  REAL8Vector *momI_12, /**< noise moments: \f$momI_12(f) = \int_{f0}^f (f'/f0)^{(12-17)/3} df'\f$ */
1295  REAL8Vector *momI_13, /**< noise moments: \f$momI_13(f) = \int_{f0}^f (f'/f0)^{(13-17)/3} df'\f$ */
1296  REAL8Vector *momI_14, /**< noise moments: \f$momI_14(f) = \int_{f0}^f (f'/f0)^{(14-17)/3} df'\f$ */
1297  REAL8Vector *momI_15, /**< noise moments: \f$momI_15(f) = \int_{f0}^f (f'/f0)^{(15-17)/3} df'\f$ */
1298  REAL8Vector *momI_16, /**< noise moments: \f$momI_16(f) = \int_{f0}^f (f'/f0)^{(16-17)/3} df'\f$ */
1299  REAL8Vector *momJ_5, /**< noise moments: \f$momJ_5(f) = \int_{f0}^f (f'/f0)^{(5-17)/3} log(f'/f0) df'\f$ */
1300  REAL8Vector *momJ_6, /**< noise moments: \f$momJ_6(f) = \int_{f0}^f (f'/f0)^{(6-17)/3} log(f'/f0) df'\f$ */
1301  REAL8Vector *momJ_7, /**< noise moments: \f$momJ_7(f) = \int_{f0}^f (f'/f0)^{(7-17)/3} log(f'/f0) df'\f$ */
1302  REAL8Vector *momJ_8, /**< noise moments: \f$momJ_8(f) = \int_{f0}^f (f'/f0)^{(8-17)/3} log(f'/f0) df'\f$ */
1303  REAL8Vector *momJ_9, /**< noise moments: \f$momJ_9(f) = \int_{f0}^f (f'/f0)^{(9-17)/3} log(f'/f0) df'\f$ */
1304  REAL8Vector *momJ_10, /**< noise moments: \f$momJ_10(f) = \int_{f0}^f (f'/f0)^{(10-17)/3} log(f'/f0) df'\f$ */
1305  REAL8Vector *momJ_11, /**< noise moments: \f$momJ_11(f) = \int_{f0}^f (f'/f0)^{(11-17)/3} log(f'/f0) df'\f$ */
1306  REAL8Vector *momJ_12, /**< noise moments: \f$momJ_12(f) = \int_{f0}^f (f'/f0)^{(12-17)/3} log(f'/f0) df'\f$ */
1307  REAL8Vector *momJ_13, /**< noise moments: \f$momJ_13(f) = \int_{f0}^f (f'/f0)^{(13-17)/3} log(f'/f0) df'\f$ */
1308  REAL8Vector *momJ_14, /**< noise moments: \f$momJ_14(f) = \int_{f0}^f (f'/f0)^{(14-17)/3} log(f'/f0) df'\f$ */
1309  REAL8Vector *momK_10, /**< noise moments: \f$momK_14(f) = \int_{f0}^f (f'/f0)^{(14-17)/3} log(f'/f0) log(f'/f0) df'\f$ */
1310  REAL8Vector *momK_11, /**< noise moments: \f$momK_15(f) = \int_{f0}^f (f'/f0)^{(15-17)/3} log(f'/f0) log(f'/f0) df'\f$ */
1311  REAL8Vector *momK_12 /**< noise moments: \f$momK_16(f) = \int_{f0}^f (f'/f0)^{(16-17)/3} log(f'/f0) log(f'/f0) df'\f$ */
1312 ) {
1313  gsl_matrix * g = XLALSimInspiralTaylorF2RedSpinFisherMatrixChirpTimes(
1314  theta0, theta3, theta3s, fLow, df, momI_0, momI_2, momI_3, momI_4,
1315  momI_5, momI_6, momI_7, momI_8, momI_9, momI_10, momI_11, momI_12,
1316  momI_13, momI_14, momI_15, momI_16, momJ_5, momJ_6, momJ_7, momJ_8,
1317  momJ_9, momJ_10, momJ_11, momJ_12, momJ_13, momJ_14, momK_10, momK_11,
1318  momK_12);
1319 
1320  int s = 0;
1321 
1322  if (!g)
1323  XLAL_ERROR(XLAL_FAILURE);
1324 
1325  /* Form submatrices g1, g2, g3, g4, defined as:
1326  * g = [ g1 g2
1327  * g4 g3 ] */
1328  gsl_matrix_view g1v = gsl_matrix_submatrix (g, 0, 0, 3, 3);
1329  gsl_matrix_view g2v = gsl_matrix_submatrix (g, 0, 3, 3, 2);
1330  gsl_matrix_view g3v = gsl_matrix_submatrix (g, 3, 3, 2, 2);
1331  gsl_matrix_view g4v = gsl_matrix_submatrix (g, 3, 0, 2, 3);
1332 
1333  gsl_matrix * g1 = gsl_matrix_calloc (3, 3);
1334  gsl_matrix * g2 = gsl_matrix_calloc (3, 2);
1335  gsl_matrix * g3 = gsl_matrix_calloc (2, 2);
1336  gsl_matrix * g4 = gsl_matrix_calloc (2, 3);
1337  gsl_matrix * g3invg4 = gsl_matrix_calloc (2, 3);
1338  gsl_matrix * g2g3invg4 = gsl_matrix_calloc (3, 3);
1339 
1340  /* Project out the t0 and phi0 dimensions: gamma = g1 - g2 g3^{-1} g4 */
1341  gsl_matrix * g3inv = gsl_matrix_calloc (2, 2);
1342  gsl_permutation * p = gsl_permutation_calloc (2);
1343 
1344  gsl_matrix_memcpy (g1, &g1v.matrix);
1345  gsl_matrix_memcpy (g2, &g2v.matrix);
1346  gsl_matrix_memcpy (g3, &g3v.matrix);
1347  gsl_matrix_memcpy (g4, &g4v.matrix);
1348  gsl_matrix_free (g);
1349 
1350  gsl_linalg_LU_decomp (g3, p, &s);
1351  gsl_linalg_LU_invert (g3, p, g3inv);
1352  gsl_permutation_free (p);
1353  gsl_matrix_free (g3);
1354 
1355  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, g3inv, g4, 0.0, g3invg4);
1356  gsl_matrix_free (g4);
1357  gsl_matrix_free (g3inv);
1358  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, g2, g3invg4, 0.0, g2g3invg4);
1359  gsl_matrix_free (g2);
1360  gsl_matrix_free (g3invg4);
1361 
1362  gsl_matrix_sub (g1, g2g3invg4);
1363  gsl_matrix_free (g2g3invg4);
1364 
1365  *gamma00 = gsl_matrix_get(g1, 0, 0);
1366  *gamma01 = gsl_matrix_get(g1, 0, 1);
1367  *gamma02 = gsl_matrix_get(g1, 0, 2);
1368  *gamma11 = gsl_matrix_get(g1, 1, 1);
1369  *gamma12 = gsl_matrix_get(g1, 1, 2);
1370  *gamma22 = gsl_matrix_get(g1, 2, 2);
1371  gsl_matrix_free (g1);
1372 
1373  return XLAL_SUCCESS;
1374 }
1375 
1376 
1377 /* compute theta0, theta3, theta3s from mc, eta, chi */
1378 void XLALSimInspiralTaylorF2RedSpinChirpTimesFromMchirpEtaChi(
1379  double *theta0, /**< dimensionless parameter related to the chirp time by theta0 = 2 pi fLow tau0 */
1380  double *theta3, /**< dimensionless parameter related to the chirp time by theta3 = -2 pi fLow tau3 */
1381  double *theta3s,/**< dimensionless parameter related to the chirp time by theta3s = 2 pi fLow tau3s */
1382  double mc, /**< chirp mass (M_sun) */
1383  double eta, /**< symmetric mass ratio */
1384  double chi, /**< reduced-spin parameter */
1385  double fLow) /**< low-frequency cutoff (Hz) */
1386 {
1387  *theta0 = 5./(128.*pow(LAL_PI*mc*LAL_MTSUN_SI*fLow,1.6666666666666667));
1388  *theta3 = cbrt(LAL_PI)/(4.*pow(mc*LAL_MTSUN_SI*fLow,0.6666666666666666)*pow(eta,0.6));
1389  *theta3s = (113.*chi)/(192.*pow(LAL_PI*mc*LAL_MTSUN_SI*fLow,0.6666666666666666)*pow(eta,0.6));
1390 }
1391 
1392 
1393 /* compute mc, eta, chi from theta0, theta3, theta3s */
1394 void XLALSimInspiralTaylorF2RedSpinMchirpEtaChiFromChirpTimes(
1395  double *mc, /**< chirp mass (M_sun) */
1396  double *eta, /**< symmetric mass ratio */
1397  double *chi, /**< reduced-spin parameter */
1398  double theta0, /**< dimensionless parameter related to the chirp time by theta0 = 2 pi fLow tau0 */
1399  double theta3, /**< dimensionless parameter related to the chirp time by theta3 = -2 pi fLow tau3 */
1400  double theta3s, /**< dimensionless parameter related to the chirp time by theta3s = 2 pi fLow tau3s */
1401  double fLow) /**< low-frequency cutoff (Hz) */
1402 {
1403  *mc = pow(2.*theta0*theta0*theta0/125., -1./5.)/(16.*LAL_PI*fLow)/LAL_MTSUN_SI;
1404  *eta = cbrt(16.*pow(LAL_PI,5.)*theta0*theta0/ (25.*pow(theta3,5.)));
1405  *chi = 48.*LAL_PI*theta3s/(113.*theta3);
1406 }
1407 
1408 /** @} */
g3
const double g3
Definition: LALSimIMRLackeyTidal2013.c:57
g2
const double g2
Definition: LALSimIMRLackeyTidal2013.c:56
g1
const double g1
Definition: LALSimIMRLackeyTidal2013.c:55
XLALSimInspiralTaylorF2RedSpinDerivAEta
static REAL8 XLALSimInspiralTaylorF2RedSpinDerivAEta(const REAL8 mc_msun, const REAL8 eta, const REAL8 chi, const REAL8 f)
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:136
XLALSimInspiralTaylorF2RedSpinDerivPsiEta
static REAL8 XLALSimInspiralTaylorF2RedSpinDerivPsiEta(const REAL8 mc_msun, const REAL8 eta, const REAL8 chi, const REAL8 f)
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:270
XLALSimInspiralTaylorF2RedSpinDerivPsiMChirp
static REAL8 XLALSimInspiralTaylorF2RedSpinDerivPsiMChirp(const REAL8 mc_msun, const REAL8 eta, const REAL8 chi, const REAL8 f)
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:331
MetricCoeffs
static REAL8 MetricCoeffs(REAL8Vector *A, REAL8Vector *dPsii, REAL8Vector *dPsij, REAL8Vector *dAi, REAL8Vector *dAj, REAL8Vector *Sh, REAL8 hSqr, REAL8 df)
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:50
XLALSimInspiralTaylorF2RedSpinDerivPsiChi
static REAL8 XLALSimInspiralTaylorF2RedSpinDerivPsiChi(const REAL8 mc_msun, const REAL8 eta, const REAL8 chi, const REAL8 f)
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:244
XLALSimInspiralTaylorF2RedSpinDerivAMChirp
static REAL8 XLALSimInspiralTaylorF2RedSpinDerivAMChirp(const REAL8 mc_msun, const REAL8 eta, const REAL8 chi, const REAL8 f)
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:192
XLALSimInspiralTaylorF2RedSpinAofF
static REAL8 XLALSimInspiralTaylorF2RedSpinAofF(const REAL8 mc_msun, const REAL8 eta, const REAL8 chi, const REAL8 f)
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:64
XLALSimInspiralTaylorF2RedSpinDerivAChi
static REAL8 XLALSimInspiralTaylorF2RedSpinDerivAChi(const REAL8 mc_msun, const REAL8 eta, const REAL8 chi, const REAL8 f)
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:111
g4
REAL8 g4
Definition: SEOBNRv2_opt_3derivs.h:420
s
int s
Definition: bh_qnmode.c:137
i
double i
Definition: bh_ringdown.c:118
LAL_PI
#define LAL_PI
LAL_TWOPI
#define LAL_TWOPI
LAL_MTSUN_SI
#define LAL_MTSUN_SI
LAL_GAMMA
#define LAL_GAMMA
REAL8
double REAL8
XLALSimInspiralTaylorF2RedSpinComputeNoiseMoments
int XLALSimInspiralTaylorF2RedSpinComputeNoiseMoments(REAL8Vector *momI_0, REAL8Vector *momI_2, REAL8Vector *momI_3, REAL8Vector *momI_4, REAL8Vector *momI_5, REAL8Vector *momI_6, REAL8Vector *momI_7, REAL8Vector *momI_8, REAL8Vector *momI_9, REAL8Vector *momI_10, REAL8Vector *momI_11, REAL8Vector *momI_12, REAL8Vector *momI_13, REAL8Vector *momI_14, REAL8Vector *momI_15, REAL8Vector *momI_16, REAL8Vector *momJ_5, REAL8Vector *momJ_6, REAL8Vector *momJ_7, REAL8Vector *momJ_8, REAL8Vector *momJ_9, REAL8Vector *momJ_10, REAL8Vector *momJ_11, REAL8Vector *momJ_12, REAL8Vector *momJ_13, REAL8Vector *momJ_14, REAL8Vector *momK_10, REAL8Vector *momK_11, REAL8Vector *momK_12, REAL8Vector *Sh, REAL8 fLow, REAL8 df)
Compute the template-space metric of "reduced-spin" PN templates in theta0, theta3,...
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:556
XLALSimInspiralTaylorF2RedSpinMchirpEtaChiFromChirpTimes
void XLALSimInspiralTaylorF2RedSpinMchirpEtaChiFromChirpTimes(double *mc, double *eta, double *chi, double theta0, double theta3, double theta3s, double fLow)
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:1394
XLALSimInspiralTaylorF2RedSpinMetricMChirpEtaChi
int XLALSimInspiralTaylorF2RedSpinMetricMChirpEtaChi(REAL8 *gamma00, REAL8 *gamma01, REAL8 *gamma02, REAL8 *gamma11, REAL8 *gamma12, REAL8 *gamma22, const REAL8 mc, const REAL8 eta, const REAL8 chi, const REAL8 fLow, const REAL8FrequencySeries *Sh)
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:383
XLALSimInspiralTaylorF2RedSpinFisherMatrixChirpTimes
gsl_matrix * XLALSimInspiralTaylorF2RedSpinFisherMatrixChirpTimes(const REAL8 theta0, const REAL8 theta3, const REAL8 theta3s, const REAL8 fLow, const REAL8 df, REAL8Vector *momI_0, REAL8Vector *momI_2, REAL8Vector *momI_3, REAL8Vector *momI_4, REAL8Vector *momI_5, REAL8Vector *momI_6, REAL8Vector *momI_7, REAL8Vector *momI_8, REAL8Vector *momI_9, REAL8Vector *momI_10, REAL8Vector *momI_11, REAL8Vector *momI_12, REAL8Vector *momI_13, REAL8Vector *momI_14, REAL8Vector *momI_15, REAL8Vector *momI_16, REAL8Vector *momJ_5, REAL8Vector *momJ_6, REAL8Vector *momJ_7, REAL8Vector *momJ_8, REAL8Vector *momJ_9, REAL8Vector *momJ_10, REAL8Vector *momJ_11, REAL8Vector *momJ_12, REAL8Vector *momJ_13, REAL8Vector *momJ_14, REAL8Vector *momK_10, REAL8Vector *momK_11, REAL8Vector *momK_12)
Compute the Fisher information matrix of "reduced-spin" PN templates in theta0, theta3,...
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:691
XLALSimInspiralTaylorF2RedSpinChirpTimesFromMchirpEtaChi
void XLALSimInspiralTaylorF2RedSpinChirpTimesFromMchirpEtaChi(double *theta0, double *theta3, double *theta3s, double mc, double eta, double chi, double fLow)
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:1378
XLALSimInspiralTaylorF2RedSpinMetricChirpTimes
int XLALSimInspiralTaylorF2RedSpinMetricChirpTimes(REAL8 *gamma00, REAL8 *gamma01, REAL8 *gamma02, REAL8 *gamma11, REAL8 *gamma12, REAL8 *gamma22, const REAL8 theta0, const REAL8 theta3, const REAL8 theta3s, const REAL8 fLow, const REAL8 df, REAL8Vector *momI_0, REAL8Vector *momI_2, REAL8Vector *momI_3, REAL8Vector *momI_4, REAL8Vector *momI_5, REAL8Vector *momI_6, REAL8Vector *momI_7, REAL8Vector *momI_8, REAL8Vector *momI_9, REAL8Vector *momI_10, REAL8Vector *momI_11, REAL8Vector *momI_12, REAL8Vector *momI_13, REAL8Vector *momI_14, REAL8Vector *momI_15, REAL8Vector *momI_16, REAL8Vector *momJ_5, REAL8Vector *momJ_6, REAL8Vector *momJ_7, REAL8Vector *momJ_8, REAL8Vector *momJ_9, REAL8Vector *momJ_10, REAL8Vector *momJ_11, REAL8Vector *momJ_12, REAL8Vector *momJ_13, REAL8Vector *momJ_14, REAL8Vector *momK_10, REAL8Vector *momK_11, REAL8Vector *momK_12)
Compute the template-space metric of "reduced-spin" PN templates in theta0, theta3,...
Definition: LALSimInspiralTaylorF2ReducedSpinMetric.c:1271
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLAL_ERROR_NULL
#define XLAL_ERROR_NULL(...)
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_SUCCESS
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFAULT
XLAL_EDOM
XLAL_EDOM
XLAL_FAILURE
XLAL_FAILURE
p
list p
REAL8FrequencySeries::data
REAL8Sequence * data
REAL8Vector::data
REAL8 * data