Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
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>
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
20
21/**
22 * Function which calculates the 4-order first finite difference
23 * derivative of a numerical function
24 */
25int
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 */
110int
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 */
176int
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 */
282int
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 */
408int
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}
int XLALFDDerivative1Order6(REAL8Vector *XVec, REAL8Vector *YVec, REAL8Vector *derivativeVec)
Function which calculates the 6-order first finite difference derivative of a numerical function.
int XLALFDDerivative1Order2(REAL8Vector *XVec, REAL8Vector *YVec, REAL8Vector *derivativeVec)
Function which calculates the 2-order first finite difference derivative of a numerical function.
int XLALFDDerivative1Order4(REAL8Vector *XVec, REAL8Vector *YVec, REAL8Vector *derivativeVec)
Function which calculates the 4-order first finite difference derivative of a numerical function.
int XLALCumulativeIntegral3(REAL8Vector *XVec, REAL8Vector *YVec, REAL8Vector *integralVec)
Function which calculates the 3-order cumulative derivative of a numerical function.
int XLALFDDerivative1Order8(REAL8Vector *XVec, REAL8Vector *YVec, REAL8Vector *derivativeVec)
Function which calculates the 8-order first finite difference derivative of a numerical function.
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 ...
#define c
double i
Definition: bh_ringdown.c:118
double e
Definition: bh_ringdown.c:117
double REAL8
uint32_t UINT4
static const INT4 a
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLAL_SUCCESS
REAL8 * data