LALSimulation  5.4.0.1-fe68b98
LALSimInspiralEOBPANumericalFuncs.c
Go to the documentation of this file.
1 #include <lal/LALSimIMR.h>
2 #include <lal/LALSimInspiral.h>
3 #include <lal/Date.h>
4 #include <lal/TimeSeries.h>
5 #include <lal/Units.h>
6 #include <lal/LALAdaptiveRungeKuttaIntegrator.h>
7 #include <lal/SphericalHarmonics.h>
8 #include <lal/LALSimSphHarmMode.h>
9 #include <LALSimInspiralWaveformFlags.h>
10 #include <lal/LALDict.h>
11 
12 #include <gsl/gsl_sf_gamma.h>
13 #include <gsl/gsl_matrix.h>
14 #include <gsl/gsl_poly.h>
15 #include <gsl/gsl_spline.h>
16 #include <gsl/gsl_roots.h>
17 #include <gsl/gsl_deriv.h>
18 
19 #include "LALSimInspiralEOBPostAdiabatic.h"
20 
21 /**
22  * Function which calculates the 4-order first finite difference
23  * derivative of a numerical function
24  */
25 int
26 XLALFDDerivative1Order4(
27  REAL8Vector *XVec,
28  /**<< An array of X values */
29  REAL8Vector *YVec,
30  /**<< An array of Y values */
31  REAL8Vector *derivativeVec
32  /**<< OUTPUT, the derivative dY/dX */
33 )
34 {
35  REAL8 fourthOrderCoeffs[5][5] = {
36  {-25./12., 4., -3., 4./3., -1./4.},
37  {-1./4., -5./6., 3./2., -1./2., 1./12.},
38  {1./12., -2./3., 0, 2./3., -1./12.},
39  {-1./12., 1./2., -3./2., 5./6., 1./4.},
40  {1./4., -4./3., 3., -4., 25./12.}
41  };
42 
43  UINT4 vecLength;
44  vecLength = XVec->length;
45 
46  REAL8 h;
47  h = fabs(XVec->data[0] - XVec->data[1]);
48 
49  UINT4 i;
50  UINT4 j;
51 
52  for (i = 0; i <= vecLength-1; i++)
53  {
54  if (i == 0)
55  {
56  for (j = 0; j <= 4; j++)
57  {
58  derivativeVec->data[i] += (
59  fourthOrderCoeffs[0][j] * YVec->data[j]
60  );
61  }
62  }
63  else if (i == 1)
64  {
65  for (j = 0; j <= 4; j++)
66  {
67  derivativeVec->data[i] += (
68  fourthOrderCoeffs[1][j] * YVec->data[j]
69  );
70  }
71  }
72  else if (i == vecLength-2)
73  {
74  for (j = 0; j <= 4; j++)
75  {
76  derivativeVec->data[i] += (
77  fourthOrderCoeffs[3][j] * YVec->data[i+j-3]
78  );
79  }
80  }
81  else if (i == vecLength-1)
82  {
83  for (j = 0; j <= 4; j++)
84  {
85  derivativeVec->data[i] += (
86  fourthOrderCoeffs[4][j] * YVec->data[i+j-4]
87  );
88  }
89  }
90  else
91  {
92  for (j = 0; j <= 4; j++)
93  {
94  derivativeVec->data[i] += (
95  fourthOrderCoeffs[2][j] * YVec->data[i+j-2]
96  );
97  }
98  }
99 
100  derivativeVec->data[i] /= h;
101  }
102 
103  return XLAL_SUCCESS;
104 }
105 
106 /**
107  * Function which calculates the 2-order first finite difference
108  * derivative of a numerical function
109  */
110 int
111 XLALFDDerivative1Order2(
112  REAL8Vector *XVec,
113  /**<< An array of X values */
114  REAL8Vector *YVec,
115  /**<< An array of Y values */
116  REAL8Vector *derivativeVec
117  /**<< OUTPUT, the derivative dY/dX */
118 )
119 {
120  REAL8 secondOrderCoeffs[3][3] = {
121  {-3./2., 2, -1./2.},
122  {-1./2., 0, 1./2.},
123  {1./2., -2., 3./2.}
124  };
125 
126  UINT4 vecLength;
127  vecLength = XVec->length;
128 
129 
130  REAL8 h;
131  h = fabs(XVec->data[0] - XVec->data[1]);
132 
133  UINT4 i;
134  UINT4 j;
135 
136  for (i = 0; i <= vecLength-1; i++)
137  {
138  if (i == 0)
139  {
140  for (j = 0; j <= 2; j++)
141  {
142  derivativeVec->data[i] += (
143  secondOrderCoeffs[0][j] * YVec->data[j]
144  );
145  }
146  }
147  else if (i == vecLength-1)
148  {
149  for (j = 0; j <= 2; j++)
150  {
151  derivativeVec->data[i] += (
152  secondOrderCoeffs[2][j] * YVec->data[i+j-2]
153  );
154  }
155  }
156  else
157  {
158  for (j = 0; j <= 2; j++)
159  {
160  derivativeVec->data[i] += (
161  secondOrderCoeffs[1][j] * YVec->data[i+j-1]
162  );
163  }
164  }
165 
166  derivativeVec->data[i] /= h;
167  }
168 
169  return XLAL_SUCCESS;
170 }
171 
172 /**
173  * Function which calculates the 6-order first finite difference
174  * derivative of a numerical function
175  */
176 int
177 XLALFDDerivative1Order6(
178  REAL8Vector *XVec,
179  /**<< An array of X values */
180  REAL8Vector *YVec,
181  /**<< An array of Y values */
182  REAL8Vector *derivativeVec
183  /**<< OUTPUT, the derivative dY/dX */
184 )
185 {
186  REAL8 sixthOrderCoeffs[7][7] = {
187  {-49./20., 6., -15./2., 20./3., -15./4., 6./5., -1./6.},
188  {-1./6., -77./60., 5./2., -5./3., 5./6., -1./4., 1./30.},
189  {1./30., -2./5., -7./12., 4./3., -1./2., 2./15., -1./60.},
190  {-1./60., 3./20., -3./4., 0, 3./4., -3./20., 1./60.},
191  {1./60., -2./15., 1./2., -4./3., 7./12., 2./5., -1./30.},
192  {-1./30., 1./4., -5./6., 5./3., -5./2., 77./60., 1./6.},
193  {1./6., -6./5., 15./4., -20./3., 15./2., -6., 49./20.}
194  };
195 
196  UINT4 vecLength;
197  vecLength = XVec->length;
198 
199 
200  REAL8 h;
201  h = fabs(XVec->data[0] - XVec->data[1]);
202 
203  UINT4 i;
204  UINT4 j;
205 
206  for (i = 0; i <= vecLength-1; i++)
207  {
208  if (i == 0)
209  {
210  for (j = 0; j <= 6; j++)
211  {
212  derivativeVec->data[i] += (
213  sixthOrderCoeffs[0][j] * YVec->data[j]
214  );
215  }
216  }
217  else if (i == 1)
218  {
219  for (j = 0; j <= 6; j++)
220  {
221  derivativeVec->data[i] += (
222  sixthOrderCoeffs[1][j] * YVec->data[j]
223  );
224  }
225  }
226  else if (i == 2)
227  {
228  for (j = 0; j <= 6; j++)
229  {
230  derivativeVec->data[i] += (
231  sixthOrderCoeffs[2][j] * YVec->data[j]
232  );
233  }
234  }
235  else if (i == vecLength-3)
236  {
237  for (j = 0; j <= 6; j++)
238  {
239  derivativeVec->data[i] += (
240  sixthOrderCoeffs[4][j] * YVec->data[i+j-4]
241  );
242  }
243  }
244  else if (i == vecLength-2)
245  {
246  for (j = 0; j <= 6; j++)
247  {
248  derivativeVec->data[i] += (
249  sixthOrderCoeffs[5][j] * YVec->data[i+j-5]
250  );
251  }
252  }
253  else if (i == vecLength-1)
254  {
255  for (j = 0; j <= 6; j++)
256  {
257  derivativeVec->data[i] += (
258  sixthOrderCoeffs[6][j] * YVec->data[i+j-6]
259  );
260  }
261  }
262  else
263  {
264  for (j = 0; j <= 6; j++)
265  {
266  derivativeVec->data[i] += (
267  sixthOrderCoeffs[3][j] * YVec->data[i+j-3]
268  );
269  }
270  }
271 
272  derivativeVec->data[i] /= h;
273  }
274 
275  return XLAL_SUCCESS;
276 }
277 
278 /**
279  * Function which calculates the 8-order first finite difference
280  * derivative of a numerical function
281  */
282 int
283 XLALFDDerivative1Order8(
284  REAL8Vector *XVec,
285  /**<< An array of X values */
286  REAL8Vector *YVec,
287  /**<< An array of Y values */
288  REAL8Vector *derivativeVec
289  /**<< OUTPUT, the derivative dY/dX */
290 )
291 {
292  REAL8 eightOrderCoeffs[9][9] = {
293  {-761./280., 8., -14., 56./3., -35./2., 56./5., -14./3., 8./7., -1./8.},
294  {-1./8., -223./140., 7./2., -7./2., 35./12., -7./4., 7./10., -1./6., 1./56.},
295  {1./56., -2./7., -19./20., 2., -5./4., 2./3., -1./4., 2./35., -1./168.},
296  {-1./168., 1./14., -1./2., -9./20., 5./4., -1./2., 1./6., -1./28., 1./280.},
297  {1./280., -4./105., 1./5., -4./5., 0, 4./5., -1./5., 4./105., -1./280.},
298  {-1./280., 1./28., -1./6., 1./2., -5./4., 9./20., 1./2., -1./14., 1./168.},
299  {1./168., -2./35., 1./4., -2./3., 5./4., -2., 19./20., 2./7., -1./56.},
300  {-1./56., 1./6., -7./10., 7./4., -35./12., 7./2., -7./2., 223./140., 1./8.},
301  {1./8., -8./7., 14./3., -56./5., 35./2., -56./3., 14., -8., 761./280.}
302  };
303 
304  UINT4 vecLength;
305  vecLength = XVec->length;
306 
307 
308  REAL8 h;
309  h = fabs(XVec->data[0] - XVec->data[1]);
310 
311  UINT4 i;
312  UINT4 j;
313 
314  for (i = 0; i <= vecLength-1; i++)
315  {
316  if (i == 0)
317  {
318  for (j = 0; j <= 8; j++)
319  {
320  derivativeVec->data[i] += (
321  eightOrderCoeffs[0][j] * YVec->data[j]
322  );
323  }
324  }
325  else if (i == 1)
326  {
327  for (j = 0; j <= 8; j++)
328  {
329  derivativeVec->data[i] += (
330  eightOrderCoeffs[1][j] * YVec->data[j]
331  );
332  }
333  }
334  else if (i == 2)
335  {
336  for (j = 0; j <= 8; j++)
337  {
338  derivativeVec->data[i] += (
339  eightOrderCoeffs[2][j] * YVec->data[j]
340  );
341  }
342  }
343  else if (i == 3)
344  {
345  for (j = 0; j <= 8; j++)
346  {
347  derivativeVec->data[i] += (
348  eightOrderCoeffs[3][j] * YVec->data[j]
349  );
350  }
351  }
352  else if (i == vecLength-4)
353  {
354  for (j = 0; j <= 8; j++)
355  {
356  derivativeVec->data[i] += (
357  eightOrderCoeffs[5][j] * YVec->data[i+j-5]
358  );
359  }
360  }
361  else if (i == vecLength-3)
362  {
363  for (j = 0; j <= 8; j++)
364  {
365  derivativeVec->data[i] += (
366  eightOrderCoeffs[6][j] * YVec->data[i+j-6]
367  );
368  }
369  }
370  else if (i == vecLength-2)
371  {
372  for (j = 0; j <= 8; j++)
373  {
374  derivativeVec->data[i] += (
375  eightOrderCoeffs[7][j] * YVec->data[i+j-7]
376  );
377  }
378  }
379  else if (i == vecLength-1)
380  {
381  for (j = 0; j <= 8; j++)
382  {
383  derivativeVec->data[i] += (
384  eightOrderCoeffs[8][j] * YVec->data[i+j-8]
385  );
386  }
387  }
388  else
389  {
390  for (j = 0; j <= 8; j++)
391  {
392  derivativeVec->data[i] += (
393  eightOrderCoeffs[4][j] * YVec->data[i+j-4]
394  );
395  }
396  }
397 
398  derivativeVec->data[i] /= h;
399  }
400 
401  return XLAL_SUCCESS;
402 }
403 
404 /**
405  * Function which calculates the 3-order cumulative derivative of a
406  * numerical function
407  */
408 int
409 XLALCumulativeIntegral3(
410  REAL8Vector *XVec,
411  /**<< An array of X values */
412  REAL8Vector *YVec,
413  /**<< An array of Y values */
414  REAL8Vector *integralVec
415  /**<< OUTPUT, the integral of Y dX */
416 )
417 {
418  UINT4 vecLength;
419  vecLength = XVec->length;
420 
421  REAL8Vector *XVecExt = XLALCreateREAL8Vector(vecLength+2);
422  REAL8Vector *YVecExt = XLALCreateREAL8Vector(vecLength+2);
423 
424  memset(XVecExt->data, 0, XVecExt->length * sizeof(REAL8));
425  memset(YVecExt->data, 0, YVecExt->length * sizeof(REAL8));
426 
427  REAL8 *X0, *X1, *X2, *X3;
428  REAL8 *Y0, *Y1, *Y2, *Y3;
429 
430  REAL8 a, b, c, d, e, h, g, z;
431  const REAL8 oo12 = 0.08333333333333333;
432 
433  UINT4 i;
434 
435  for (i=1; i < vecLength+1; i++)
436  {
437  XVecExt->data[i] = XVec->data[i-1];
438  YVecExt->data[i] = YVec->data[i-1];
439  }
440 
441  XVecExt->data[0] = XVec->data[3];
442  XVecExt->data[vecLength+1] = XVec->data[vecLength-4];
443 
444  YVecExt->data[0] = YVec->data[3];
445  YVecExt->data[vecLength+1] = YVec->data[vecLength-4];
446 
447  X0 = &XVecExt->data[0];
448  X1 = &XVecExt->data[1];
449  X2 = &XVecExt->data[2];
450  X3 = &XVecExt->data[3];
451 
452  Y0 = &YVecExt->data[0];
453  Y1 = &YVecExt->data[1];
454  Y2 = &YVecExt->data[2];
455  Y3 = &YVecExt->data[3];
456 
457  for (i=0; i < vecLength-1; i++)
458  {
459  a = X1[i] - X0[i];
460  b = X2[i] - X1[i];
461  c = X3[i] - X2[i];
462  d = Y1[i] - Y0[i];
463  e = Y2[i] - Y1[i];
464  h = Y3[i] - Y2[i];
465  g = 0.5 * (Y1[i]+Y2[i]);
466  z = (
467  b * g
468  + oo12 * b * b * (
469  c * b * (2 * c + b) * (c + b) * d
470  - a * c * (c - a) * (2 * c + 2 * a + 3 * b) * e
471  - a * b * (2 * a + b) * (a + b) * h
472  ) / (a * c * (a + b) * (c + b) * (c + a + b))
473  );
474  integralVec->data[i+1] = integralVec->data[i] + z;
475  }
476  XLALDestroyREAL8Vector(XVecExt);
477  XLALDestroyREAL8Vector(YVecExt);
478 
479  return XLAL_SUCCESS;
480 }
XLALFDDerivative1Order6
int XLALFDDerivative1Order6(REAL8Vector *XVec, REAL8Vector *YVec, REAL8Vector *derivativeVec)
Function which calculates the 6-order first finite difference derivative of a numerical function.
Definition: LALSimInspiralEOBPANumericalFuncs.c:177
XLALFDDerivative1Order2
int XLALFDDerivative1Order2(REAL8Vector *XVec, REAL8Vector *YVec, REAL8Vector *derivativeVec)
Function which calculates the 2-order first finite difference derivative of a numerical function.
Definition: LALSimInspiralEOBPANumericalFuncs.c:111
XLALFDDerivative1Order4
int XLALFDDerivative1Order4(REAL8Vector *XVec, REAL8Vector *YVec, REAL8Vector *derivativeVec)
Function which calculates the 4-order first finite difference derivative of a numerical function.
Definition: LALSimInspiralEOBPANumericalFuncs.c:26
XLALCumulativeIntegral3
int XLALCumulativeIntegral3(REAL8Vector *XVec, REAL8Vector *YVec, REAL8Vector *integralVec)
Function which calculates the 3-order cumulative derivative of a numerical function.
Definition: LALSimInspiralEOBPANumericalFuncs.c:409
XLALFDDerivative1Order8
int XLALFDDerivative1Order8(REAL8Vector *XVec, REAL8Vector *YVec, REAL8Vector *derivativeVec)
Function which calculates the 8-order first finite difference derivative of a numerical function.
Definition: LALSimInspiralEOBPANumericalFuncs.c:283
LALSimInspiralEOBPostAdiabatic.h
d
static vector d(const double L_norm, const double J_norm, const vector roots)
Internal function that returns the coefficients "d_0", "d_2" and "d_4" from 1703.03967 corresponding ...
Definition: LALSimInspiralFDPrecAngles_internals.c:604
c
#define c
Definition: LALSimInspiralHGimri.c:35
i
double i
Definition: bh_ringdown.c:118
e
double e
Definition: bh_ringdown.c:117
REAL8
double REAL8
UINT4
uint32_t UINT4
a
static const INT4 a
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLAL_SUCCESS
XLAL_SUCCESS
REAL8Vector::data
REAL8 * data