23#include <lal/Integrate.h>
24#include <lal/Interpolate.h>
25#include <lal/LALDatatypes.h>
26#include <lal/XLALError.h>
29#define UNUSED __attribute__ ((unused))
54 INT4 n = 1 << (refinement - 1);
55 REAL8 dx = (xmax - xmin)/n;
58 for (i = 0; i < n; ++i)
67 integral += (xmax - xmin)*sum/n;
74 y_0 =
f (xmin, params);
77 y_1 =
f (xmax, params);
80 integral = (xmax - xmin)*(y_1 + y_0)/2;
171 if (!((b > 0 &&
a > 0) || (b < 0 &&
a < 0)))
192 REAL8 dx1 = (xmax - xmin)/(3*n);
196 for (i = 0; i < n; ++i)
201 y =
f (ChangeOfVariables (
x,
a, b, &jac), params);
207 y =
f (ChangeOfVariables (
x,
a, b, &jac), params);
213 integral += (xmax - xmin)*sum/n;
218 REAL8 x = (xmax + xmin)/2;
221 y =
f (ChangeOfVariables (
x,
a, b, &jac), params);
224 integral = (xmax - xmin)*
y*jac;
241 const REAL8 epsilon = 1e-15;
242 enum { MaxSteps = 20 };
247 REAL8 integral[MaxSteps + 1];
248 REAL8 stepSize[MaxSteps + 1];
271 refineFactor = 1.0/9.0;
279 for (refinement = 0; refinement < MaxSteps; ++refinement)
281 temp = Algorithm (temp,
f, params, xmin, xmax, type, refinement);
285 integral[refinement] = temp;
287 if (refinement >= Order)
297 if (fabs(dy) < epsilon*fabs(
y))
303 integral[refinement + 1] = temp;
304 stepSize[refinement + 1] = refineFactor*stepSize[refinement];
static REAL8 XLALREAL8Trapezoid(REAL8 integral, REAL8(*f)(REAL8 x, void *params), void *params, REAL8 xmin, REAL8 xmax, IntegralType type, int refinement)
static REAL8 DEqualsAPlusXSq(REAL8 x, REAL8 a, REAL8 UNUSED b, REAL8 *jac)
static REAL8 DEqualsInvX(REAL8 x, REAL8 UNUSED a, REAL8 UNUSED b, REAL8 *jac)
static INT4 ThreePow(INT4 n)
static REAL8 DEqualsMinusLogX(REAL8 x, REAL8 UNUSED a, REAL8 UNUSED b, REAL8 *jac)
static REAL8 XLALREAL8Midpoint(REAL8 integral, REAL8(*f)(REAL8 x, void *params), void *params, REAL8 a, REAL8 b, IntegralType type, int refinement)
static REAL8 DEqualsBMinusXSq(REAL8 x, REAL8 UNUSED a, REAL8 b, REAL8 *jac)
static REAL8 DEqualsX(REAL8 x, REAL8 UNUSED a, REAL8 UNUSED b, REAL8 *jac)
static double f(double theta, double y, double xi)
IntegralType
Type of integral.
REAL8 XLALREAL8RombergIntegrate(REAL8(*f)(REAL8 x, void *params), void *params, REAL8 xmin, REAL8 xmax, IntegralType type)
@ InfiniteDomainExp
integrate infinite domain with exponential falloff
@ ClosedInterval
evaluate integral on a closed interval
@ SingularLowerLimit
integrate an inv sqrt singularity at lower limit
@ SingularUpperLimit
integrate an inv sqrt singularity at upper limit
@ InfiniteDomainPow
integrate infinite domain with power-law falloff
@ OpenInterval
evaluate integral on an open interval
REAL8 XLALREAL8PolynomialInterpolation(REAL8 *yout, REAL8 xtarget, REAL8 *y, REAL8 *x, UINT4 n)
double REAL8
Double precision real floating-point number (8 bytes).
int32_t INT4
Four-byte signed integer.
#define XLAL_ERROR_REAL8(...)
Macro to invoke a failure from a XLAL routine returning a REAL8.
#define xlalErrno
Modifiable lvalue containing the XLAL error number.
@ XLAL_EFAULT
Invalid pointer.
@ XLAL_EFUNC
Internal function call failed bit: "or" this with existing error number.
@ XLAL_EMAXITER
Exceeded maximum number of iterations.
@ XLAL_EDOM
Input domain error.
@ XLAL_EINVAL
Invalid argument.