Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JQuadrature.hh
Go to the documentation of this file.
1#ifndef __JTOOLS__JQUADRATURE__
2#define __JTOOLS__JQUADRATURE__
3
4#include <cmath>
5
7#include "JTools/JElement.hh"
9
10
11/**
12 * \file
13 *
14 * Auxiliary classes for numerical integration.
15 * \author mdejong
16 */
17namespace JTOOLS {}
18namespace JPP { using namespace JTOOLS; }
19
20namespace JTOOLS {
21
22 /**
23 * Type definition of basic element for quadratures.
24 */
26
27
28 /**
29 * Type definition for numerical integration.
30 *
31 * \f[\displaystyle \int_{x_1}^{x_2} f(x) dx = \sum_{i=1}^{N} w_i f(x_i) \f]
32 *
33 * The abscissa and ordinate values of the collection can be used
34 * as abscissa and weight values of the summation to approximately
35 * determine the integral of the function.
36 */
38 public JCollection<JElement2D_t>
39 {
40 public:
41 /**
42 * Default constructor.
43 */
45 {}
46
47
48 /**
49 * General purpose constructor.
50 *
51 * The template argument should correspond to a function requiring two arguments.\n
52 * These two arguments should correspond to the lower and upper integration limit, respectively.\n
53 * The given function should return the value of the integral between the two integration limits.
54 *
55 * \param xmin minimal x
56 * \param xmax maximal x
57 * \param nx number of points
58 * \param integral integral
59 * \param eps precision
60 */
61 template<class JFunction_t>
62 JQuadrature(const double xmin,
63 const double xmax,
64 const int nx,
65 JFunction_t integral,
66 const double eps = 1.0e-4)
67 {
68 double Xmin = xmin;
69 double Xmax = xmax;
70
71 const double Vmin = integral(Xmin, Xmax) / (double) nx;
72
73 for (int i = 0; i != nx; ++i) {
74
75 for (double xmin = Xmin, xmax = Xmax; ; ) {
76
77 const double x = 0.5 * (xmin + xmax);
78 const double v = integral(Xmin, x);
79
80 if (fabs(Vmin - v) < eps * Vmin ||
81 fabs(xmax - xmin) < eps * (Xmax - Xmin)) {
82
83 const double __x = 0.5 * (Xmin + x);
84 const double __y = Vmin / integral(__x);
85
86 insert(JElement2D_t(__x,__y));
87
88 Xmin = x;
89 xmax = Xmax;
90
91 break;
92 }
93
94 if (v < Vmin)
95 xmin = x;
96 else
97 xmax = x;
98 }
99 }
100 }
101 };
102
103
104 /**
105 * Numerical integrator for \f$ W(x) = 1 \f$.
106 *
107 * Gauss-Legendre integration code is taken from reference:
108 * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
109 * Cambridge University Press.
110 */
112 public JQuadrature
113 {
114 public:
115 /**
116 * Constructor.
117 *
118 * \param n number of points
119 * \param eps precision
120 */
121 JGaussLegendre(const int n,
122 const double eps = 1.0e-12) :
124 {
125 resize(n);
126
127 const int M = (n + 1) / 2;
128
129 double p0, p1, p2, pp;
130
131 for (int i = 0; i < M; ++i) {
132
133 double z = cos(PI * (i+0.75) / (n+0.5));
134 double z1;
135
136 do {
137
138 p1 = 0.0;
139 p2 = 1.0;
140
141 // recurrence relation
142
143 for (int j = 0; j < n; ++j) {
144 p0 = p1;
145 p1 = p2;
146 p2 = ((2*j + 1) * z*p1 - j*p0) / (j+1);
147 }
148
149 pp = n * (z*p2 - p1) / (z*z - 1.0);
150
151 z1 = z;
152 z = z1 - p2/pp;
153
154 } while (fabs(z-z1) > eps);
155
156 const double y = 2.0 / ((1.0-z*z)*pp*pp);
157
158 at( i ) = JElement2D_t(-z,y);
159 at(n-i-1) = JElement2D_t(+z,y);
160 }
161 }
162 };
163
164
165 /**
166 * Numerical integrator for \f$ W(x) = x^{a} \, e^{-x} \f$.
167 *
168 * Gauss-Laguerre integration code is taken from reference:
169 * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
170 * Cambridge University Press.
171 */
173 public JQuadrature
174 {
175 public:
176 /**
177 * Constructor.
178 *
179 * \param n number of points
180 * \param alf power
181 * \param eps precision
182 */
183 JGaussLaguerre(const int n,
184 const double alf,
185 const double eps = 1.0e-12) :
187 {
188 const int number_of_iterations = 100;
189
190 double z1;
191 double p0, p1, p2, pp;
192
193 double z = (1.0 + alf) * (3.0 + 0.92*alf) / (1.0 + 2.4*n + 1.8*alf);
194
195 for (int i = 0; i < n; ++i) {
196
197 switch (i) {
198
199 case 0:
200 break;
201
202 case 1:
203 z += (15.0 + 6.25*alf) / (1.0 + 0.9*alf + 2.5*n);
204 break;
205
206 default:
207 const double ai = i - 1;
208 z += ((1.0+2.55*ai)/(1.9*ai) + (1.26*ai*alf)/(1.0+3.5*ai)) * (z - at(i-2).getX()) / (1.0 + 0.3*alf);
209 break;
210 }
211
212 for (int k = 0; k != number_of_iterations; ++k) {
213
214 p1 = 0.0;
215 p2 = 1.0;
216
217 // recurrence relation
218
219 for (int j = 0; j < n; ++j) {
220 p0 = p1;
221 p1 = p2;
222 p2 = ((2*j + 1 + alf - z) * p1 - (j + alf)*p0) / (j+1);
223 }
224
225 pp = (n*p2 - (n+alf)*p1) / z;
226
227 z1 = z;
228 z = z1 - p2/pp;
229
230 if (fabs(z-z1) < eps)
231 break;
232 }
233
234 const double y = -tgamma(alf+n) / tgamma((double) n) / (pp*n*p1);
235
236 insert(JElement2D_t(z,y));
237 }
238 }
239 };
240
241
242 /**
243 * Numerical integrator for \f$ W(x) = e^{-x^{2}} \f$.
244 *
245 * Gauss-Hermite integration code is taken from reference:
246 * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
247 * Cambridge University Press.
248 */
250 public JQuadrature
251 {
252 public:
253 /**
254 * Constructor.
255 *
256 * \param n number of points
257 * \param eps precision
258 */
259 JGaussHermite(const int n,
260 const double eps = 1.0e-12) :
262 {
263 resize(n);
264
265 const double pii = 1.0 / pow(PI,0.25);
266
267 const int number_of_iterations = 100;
268
269 const int M = (n + 1) / 2;
270
271 double p0, p1, p2, pp;
272 double z = 0.0;
273 double z1;
274
275 for (int i = 0; i < M; ++i) {
276
277 switch (i) {
278
279 case 0:
280 z = sqrt((double) (2*n+1)) - 1.85575 * pow((double) (2*n+1),-0.16667);
281 break;
282
283 case 1:
284 z -= 1.14 * pow((double) n,0.426) / z;
285 break;
286
287 case 2:
288 z = 1.86*z + 0.86*at( 0 ).getX();
289 break;
290
291 case 3:
292 z = 1.91*z + 0.91*at( 1 ).getX();
293 break;
294
295 default:
296 z = 2.00*z + 1.00*at(i-2).getX();
297 break;
298 }
299
300 for (int k = 0; k != number_of_iterations; ++k) {
301
302 p1 = 0.0;
303 p2 = pii;
304
305 // recurrence relation
306
307 for (int j = 0; j < n; ++j) {
308 p0 = p1;
309 p1 = p2;
310 p2 = z * sqrt(2.0/(double) (j+1)) * p1 - sqrt((double) j / (double) (j+1)) * p0;
311 }
312
313 pp = sqrt((double) (2*n)) * p1;
314
315 z1 = z;
316 z = z1 - p2/pp;
317
318 if (fabs(z-z1) < eps)
319 break;
320 }
321
322 const double y = 2.0 / (pp*pp);
323
324 at( i ) = JElement2D_t(-z,y);
325 at(n-i-1) = JElement2D_t(+z,y);
326 }
327 }
328 };
329
330
331 /**
332 * Numerical integrator for \f$ W(x) = (1 + g^{2} - 2gx)^{a} \f$, where \f$ g > 0 \f$.
333 *
334 * Henyey-Greenstein integration points and weights.
335 */
337 public JQuadrature
338 {
339 public:
340 /**
341 * Constructor.
342 *
343 * \param n number of points
344 * \param g angular dependence parameter
345 * \param a power
346 */
348 const double g,
349 const double a) :
351 {
352 const double b = -2*g * (a + 1.0);
353 const double ai = 1.0 / (a + 1.0);
354
355 const double ymin = pow(1.0 + g, 2*(a + 1.0)) / b;
356 const double ymax = pow(1.0 - g, 2*(a + 1.0)) / b;
357
358 const double dy = (ymax - ymin) / (n + 1);
359
360 for (double y = ymax - 0.5*dy; y > ymin; y -= dy) {
361
362 const double v = y*b;
363 const double w = pow(v, ai);
364 const double x = (1.0 + g*g - w) / (2*g);
365 const double dx = pow(v, -a*ai)*dy;
366
367 insert(JElement2D_t(x,dx));
368 }
369 }
370
371
372 /**
373 * Constructor.
374 *
375 * \param n number of points
376 * \param g angular dependence parameter
377 * \param a power
378 * \param xmin minimal value
379 * \param xmax maximal value
380 */
382 const double g,
383 const double a,
384 const double xmin,
385 const double xmax) :
387 {
388 const double b = -2*g * (a + 1.0);
389 const double ai = 1.0 / (a + 1.0);
390
391 const double ymin = pow(1.0 + g*g -2*g*xmin, a + 1.0) / b;
392 const double ymax = pow(1.0 + g*g -2*g*xmax, a + 1.0) / b;
393
394 const double dy = (ymax - ymin) / (n + 1);
395
396 for (double y = ymax - 0.5*dy; y > ymin; y -= dy) {
397
398 const double v = y*b;
399 const double w = pow(v, ai);
400 const double x = (1.0 + g*g - w) / (2*g);
401 const double dx = pow(v, -a*ai)*dy;
402
403 insert(JElement2D_t(x,dx));
404 }
405 }
406
407
408 /**
409 * Constructor for special case where a = -1.
410 *
411 * \param n number of points
412 * \param g angular dependence parameter
413 */
415 const double g) :
417 {
418 const double dy = 1.0 / (n + 1);
419 const double gi = log((1.0 + g*g) / (1.0 - g*g)) / (2.0*g);
420
421 for (double y = 1.0 - 0.5*dy; y > 0.0; y -= dy) {
422
423 const double v = -y*2.0*g*gi + log(1.0 + g*g);
424 const double w = exp(v);
425 const double x = (1.0 + g*g - w) / (2.0*g);
426 const double dx = w*gi*dy;
427
428 insert(JElement2D_t(x,dx));
429 }
430 }
431 };
432
433
434 /**
435 * Numerical integrator for \f$ W(x) = 1 + g \, x^{2} \f$, where \f$ g > 0 \f$.
436 *
437 * Rayleigh integration points and weights.
438 */
439 class JRayleigh :
440 public JQuadrature
441 {
442 public:
443 /**
444 * Constructor.
445 *
446 * \param n number of points
447 * \param g angular dependence parameter
448 */
449 JRayleigh(const int n,
450 const double g) :
452 {
453 const double dy = 1.0 / (n + 1);
454 const double gi = 3.0/g + 1.0;
455
456 // t^3 + 3pt + 2q = 0
457
458 const double p = 1.0/g;
459
460 for (double y = 0.5*dy; y < 1.0; y += dy) {
461
462 const double q = 0.5*gi - gi*y;
463
464 const double b = sqrt(q*q + p*p*p);
465 const double u = pow(-q + b, 1.0/3.0);
466 const double v = pow(+q + b, 1.0/3.0);
467
468 const double x = u - v;
469 const double dx = (u + v) / (3.0*b);
470
471 insert(JElement2D_t(x, dx*gi*dy));
472 }
473 }
474 };
475
476
477 /**
478 * Numerical integrator for \f$ W(x) = \left|x\right| / \sqrt{1 - x^{2}} \f$.
479 *
480 * Co-tangent integration points and weights.
481 */
482 class JCotangent :
483 public JQuadrature
484 {
485 public:
486 /**
487 * Constructor.
488 *
489 * \param n number of points
490 */
491 JCotangent(const int n) :
493 {
494 for (double ds = 1.0 / (n/2), sb = 0.5*ds; sb < 1.0; sb += ds) {
495
496 const double cb = sqrt((1.0 + sb)*(1.0 - sb));
497 const double dc = ds*sb/cb;
498
499 insert(JElement2D_t(+cb, dc));
500 insert(JElement2D_t(-cb, dc));
501 }
502 }
503 };
504
505
506 /**
507 * Numerical integrator for \f$ W(x) = \left|x\right| / \sqrt{1 - x^{2}} \f$ for \f$ x > 0 \f$ and \f$ W(x) = 1 \f$ for \f$ x \le 0 \f$.
508 *
509 * Bi-tangent integration points and weights.
510 */
511 class JBitangent :
512 public JQuadrature
513 {
514 public:
515 /**
516 * Constructor.
517 *
518 * \param n number of points
519 */
520 JBitangent(const int n) :
522 {
523 double sb, ds;
524 double cb = 0.0;
525 double dc = 0.0;
526
527 for (ds = 1.0 / (n/2), sb = 0.5*ds; sb < 1.0; sb += ds) {
528
529 cb = sqrt((1.0 + sb)*(1.0 - sb));
530 dc = ds*sb/cb;
531
532 insert(JElement2D_t(cb, dc));
533 }
534
535 for (dc = (cb + 1.0) / (n/2), cb -= 0.5*dc ; cb > -1.0; cb -= dc) {
536 insert(JElement2D_t(cb, dc));
537 }
538 }
539 };
540}
541
542#endif
General purpose class for a collection of sorted elements.
TPaveText * p1
The elements in a collection are sorted according to their abscissa values and a given distance opera...
Constants.
Numerical integrator for for and for .
JBitangent(const int n)
Constructor.
General purpose class for collection of elements, see: <a href="JTools.PDF";>Collection of elements....
Definition JSet.hh:22
void resize(typename container_type::size_type size)
virtual abscissa_type getX(int index) const override
pair_type insert(const value_type &element)
Numerical integrator for .
JCotangent(const int n)
Constructor.
Numerical integrator for .
JGaussHermite(const int n, const double eps=1.0e-12)
Constructor.
Numerical integrator for .
JGaussLaguerre(const int n, const double alf, const double eps=1.0e-12)
Constructor.
Numerical integrator for .
JGaussLegendre(const int n, const double eps=1.0e-12)
Constructor.
Numerical integrator for , where .
JHenyeyGreenstein(const int n, const double g, const double a)
Constructor.
JHenyeyGreenstein(const int n, const double g, const double a, const double xmin, const double xmax)
Constructor.
JHenyeyGreenstein(const int n, const double g)
Constructor for special case where a = -1.
Type definition for numerical integration.
JQuadrature()
Default constructor.
JQuadrature(const double xmin, const double xmax, const int nx, JFunction_t integral, const double eps=1.0e-4)
General purpose constructor.
Numerical integrator for , where .
JRayleigh(const int n, const double g)
Constructor.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary classes and methods for multi-dimensional interpolations and histograms.
JElement2D< double, double > JElement2D_t
Type definition of basic element for quadratures.
const int n
Definition JPolint.hh:791
int j
Definition JPolint.hh:801
2D Element.
Definition JPolint.hh:1131