Loading [MathJax]/extensions/TeX/AMSsymbols.js
LAL 7.7.0.1-3a66518
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
Integrate.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2007 Jolien Creighton, Drew Keppel
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/* ---------- see Integrate.h for doxygen documentation ---------- */
21#include <math.h>
22#include <string.h>
23#include <lal/Integrate.h>
24#include <lal/Interpolate.h>
25#include <lal/LALDatatypes.h>
26#include <lal/XLALError.h>
27
28#ifdef __GNUC__
29#define UNUSED __attribute__ ((unused))
30#else
31#define UNUSED
32#endif
33
34static REAL8
36 REAL8 integral,
37 REAL8 (*f)(REAL8 x, void *params),
38 void *params,
39 REAL8 xmin,
40 REAL8 xmax,
41 IntegralType type,
42 int refinement
43 )
44{
45 switch (type)
46 {
47 default:
48 break;
49 }
50
51 if (refinement)
52 {
53 INT4 i;
54 INT4 n = 1 << (refinement - 1);
55 REAL8 dx = (xmax - xmin)/n;
56 REAL8 x = xmin + dx/2;
57 REAL8 sum = 0;
58 for (i = 0; i < n; ++i)
59 {
60 REAL8 y;
61 y = f (x, params);
62 if (xlalErrno)
64 sum += y;
65 x += dx;
66 }
67 integral += (xmax - xmin)*sum/n;
68 integral /= 2;
69 }
70 else
71 {
72 REAL8 y_0;
73 REAL8 y_1;
74 y_0 = f (xmin, params);
75 if (xlalErrno)
77 y_1 = f (xmax, params);
78 if (xlalErrno)
80 integral = (xmax - xmin)*(y_1 + y_0)/2;
81 }
82
83 return integral;
84}
85
86
87static INT4
89{
90 INT4 p = 1;
91 while (n-- > 0)
92 {
93 p *= 3;
94 }
95 return p;
96}
97
98
99static REAL8
100DEqualsX (REAL8 x, REAL8 UNUSED a, REAL8 UNUSED b, REAL8 *jac)
101{
102 *jac = 1;
103 return x;
104}
105
106static REAL8
107DEqualsInvX (REAL8 x, REAL8 UNUSED a, REAL8 UNUSED b, REAL8 *jac)
108{
109 REAL8 invx = 1/x;
110 *jac = invx*invx;
111 return invx;
112}
113
114static REAL8
116{
117 *jac = 2*x;
118 return a + x*x;
119}
120
121static REAL8
123{
124 *jac = 2*x;
125 return b - x*x;
126}
127
128static REAL8
129DEqualsMinusLogX (REAL8 x, REAL8 UNUSED a, REAL8 UNUSED b, REAL8 *jac)
130{
131 *jac = 1/x;
132 return -log(x);
133}
134
135static REAL8
137 REAL8 integral,
138 REAL8 (*f)(REAL8 x, void *params),
139 void *params,
140 REAL8 a,
141 REAL8 b,
142 IntegralType type,
143 int refinement
144 )
145{
146 REAL8 (*ChangeOfVariables)(REAL8, REAL8, REAL8, REAL8 *);
147 REAL8 xmax;
148 REAL8 xmin;
149
150 switch (type)
151 {
152 case OpenInterval:
153 ChangeOfVariables = DEqualsX;
154 xmax = b;
155 xmin = a;
156 break;
157
159 ChangeOfVariables = DEqualsAPlusXSq;
160 xmax = sqrt(b - a);
161 xmin = 0;
162 break;
163
165 ChangeOfVariables = DEqualsBMinusXSq;
166 xmax = sqrt(b - a);
167 xmin = 0;
168 break;
169
171 if (!((b > 0 && a > 0) || (b < 0 && a < 0)))
173 ChangeOfVariables = DEqualsInvX;
174 xmax = 1/a;
175 xmin = 1/b;
176 break;
177
179 ChangeOfVariables = DEqualsMinusLogX;
180 xmax = exp(-a);
181 xmin = 0;
182 break;
183
184 default: /* unrecognized type */
186 }
187
188 if (refinement)
189 {
190 INT4 i;
191 INT4 n = ThreePow (refinement - 1);
192 REAL8 dx1 = (xmax - xmin)/(3*n);
193 REAL8 dx2 = 2*dx1;
194 REAL8 x = xmin + dx1/2;
195 REAL8 sum = 0;
196 for (i = 0; i < n; ++i)
197 {
198 REAL8 y;
199 REAL8 jac;
200
201 y = f (ChangeOfVariables (x, a, b, &jac), params);
202 if (xlalErrno)
204 sum += y*jac;
205 x += dx2;
206
207 y = f (ChangeOfVariables (x, a, b, &jac), params);
208 if (xlalErrno)
210 sum += y*jac;
211 x += dx1;
212 }
213 integral += (xmax - xmin)*sum/n;
214 integral /= 3;
215 }
216 else
217 {
218 REAL8 x = (xmax + xmin)/2;
219 REAL8 y;
220 REAL8 jac;
221 y = f (ChangeOfVariables (x, a, b, &jac), params);
222 if (xlalErrno)
224 integral = (xmax - xmin)*y*jac;
225 }
226
227 return integral;
228}
229
230
231
232REAL8
234 REAL8 (*f)(REAL8 x, void *params),
235 void *params,
236 REAL8 xmin,
237 REAL8 xmax,
238 IntegralType type
239 )
240{
241 const REAL8 epsilon = 1e-15;
242 enum { MaxSteps = 20 };
243 enum { Order = 4 };
244
245 REAL8 (*Algorithm)(REAL8, REAL8 (*)(REAL8, void *), void *, REAL8, REAL8, IntegralType, int);
246
247 REAL8 integral[MaxSteps + 1];
248 REAL8 stepSize[MaxSteps + 1];
249 REAL8 refineFactor;
250 REAL8 temp = 0;
251 int refinement;
252
253 if (f == NULL)
255 if (xmax <= xmin)
257
258 switch (type)
259 {
260 case ClosedInterval:
261 Algorithm = XLALREAL8Trapezoid;
262 refineFactor = 0.25;
263 break;
264
265 case OpenInterval:
270 Algorithm = XLALREAL8Midpoint;
271 refineFactor = 1.0/9.0;
272 break;
273
274 default: /* unrecognized type */
276 }
277
278 stepSize[0] = 1;
279 for (refinement = 0; refinement < MaxSteps; ++refinement)
280 {
281 temp = Algorithm (temp, f, params, xmin, xmax, type, refinement);
282 if (xlalErrno)
284
285 integral[refinement] = temp;
286
287 if (refinement >= Order)
288 {
289 REAL8 dy;
290 REAL8 y;
291
292 /* extrapolate to continuum limit (stepSize = 0) */
293 dy = XLALREAL8PolynomialInterpolation (&y, 0, integral + refinement - Order, stepSize + refinement - Order, Order + 1);
294 if (xlalErrno)
296
297 if (fabs(dy) < epsilon*fabs(y))
298 {
299 return y;
300 }
301 }
302
303 integral[refinement + 1] = temp;
304 stepSize[refinement + 1] = refineFactor*stepSize[refinement];
305 }
306
308}
static REAL8 XLALREAL8Trapezoid(REAL8 integral, REAL8(*f)(REAL8 x, void *params), void *params, REAL8 xmin, REAL8 xmax, IntegralType type, int refinement)
Definition: Integrate.c:35
static REAL8 DEqualsAPlusXSq(REAL8 x, REAL8 a, REAL8 UNUSED b, REAL8 *jac)
Definition: Integrate.c:115
static REAL8 DEqualsInvX(REAL8 x, REAL8 UNUSED a, REAL8 UNUSED b, REAL8 *jac)
Definition: Integrate.c:107
static INT4 ThreePow(INT4 n)
Definition: Integrate.c:88
static REAL8 DEqualsMinusLogX(REAL8 x, REAL8 UNUSED a, REAL8 UNUSED b, REAL8 *jac)
Definition: Integrate.c:129
static REAL8 XLALREAL8Midpoint(REAL8 integral, REAL8(*f)(REAL8 x, void *params), void *params, REAL8 a, REAL8 b, IntegralType type, int refinement)
Definition: Integrate.c:136
static REAL8 DEqualsBMinusXSq(REAL8 x, REAL8 UNUSED a, REAL8 b, REAL8 *jac)
Definition: Integrate.c:122
static REAL8 DEqualsX(REAL8 x, REAL8 UNUSED a, REAL8 UNUSED b, REAL8 *jac)
Definition: Integrate.c:100
static double f(double theta, double y, double xi)
Definition: XLALMarcumQ.c:258
IntegralType
Type of integral.
Definition: Integrate.h:136
REAL8 XLALREAL8RombergIntegrate(REAL8(*f)(REAL8 x, void *params), void *params, REAL8 xmin, REAL8 xmax, IntegralType type)
Definition: Integrate.c:233
@ InfiniteDomainExp
integrate infinite domain with exponential falloff
Definition: Integrate.h:142
@ ClosedInterval
evaluate integral on a closed interval
Definition: Integrate.h:137
@ SingularLowerLimit
integrate an inv sqrt singularity at lower limit
Definition: Integrate.h:139
@ SingularUpperLimit
integrate an inv sqrt singularity at upper limit
Definition: Integrate.h:140
@ InfiniteDomainPow
integrate infinite domain with power-law falloff
Definition: Integrate.h:141
@ OpenInterval
evaluate integral on an open interval
Definition: Integrate.h:138
REAL8 XLALREAL8PolynomialInterpolation(REAL8 *yout, REAL8 xtarget, REAL8 *y, REAL8 *x, UINT4 n)
Definition: Interpolate.c:140
double REAL8
Double precision real floating-point number (8 bytes).
int32_t INT4
Four-byte signed integer.
static const INT4 a
Definition: Random.c:79
#define XLAL_ERROR_REAL8(...)
Macro to invoke a failure from a XLAL routine returning a REAL8.
Definition: XLALError.h:752
#define xlalErrno
Modifiable lvalue containing the XLAL error number.
Definition: XLALError.h:571
@ XLAL_EFAULT
Invalid pointer.
Definition: XLALError.h:408
@ XLAL_EFUNC
Internal function call failed bit: "or" this with existing error number.
Definition: XLALError.h:462
@ XLAL_EMAXITER
Exceeded maximum number of iterations.
Definition: XLALError.h:455
@ XLAL_EDOM
Input domain error.
Definition: XLALError.h:410
@ XLAL_EINVAL
Invalid argument.
Definition: XLALError.h:409