Jpp - the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JQuadrature.hh
Go to the documentation of this file.
1 #ifndef __JTOOLS__JQUADRATURE__
2 #define __JTOOLS__JQUADRATURE__
3 
4 #include <cmath>
5 
6 #include "JTools/JConstants.hh"
7 #include "JTools/JElement.hh"
8 #include "JTools/JCollection.hh"
9 
10 
11 /**
12  * \file
13  *
14  * Auxiliary classes for numerical integration.
15  * \author mdejong
16  */
17 namespace JTOOLS {}
18 namespace JPP { using namespace JTOOLS; }
19 
20 namespace 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  */
37  class JQuadrature :
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.
52  * These two arguments should correspond to the lower and upper integration limit, respectively.
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 <tt>W(x) = 1</tt>.
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  */
111  class JGaussLegendre :
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) :
123  JQuadrature()
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 <tt>W(x) = x^a e^(-x)</tt>.
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  */
172  class JGaussLaguerre :
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) :
186  JQuadrature()
187  {
188  const int number_of_iterations = 100;
189 
190  double z, z1;
191  double p0, p1, p2, pp;
192 
193  for (int i = 0; i < n; ++i) {
194 
195  switch (i) {
196 
197  case 0:
198  z = (1.0 + alf) * (3.0 + 0.92*alf) / (1.0 + 2.4*n + 1.8*alf);
199  break;
200 
201  case 1:
202  z += (15.0 + 6.25*alf) / (1.0 + 0.9*alf + 2.5*n);
203  break;
204 
205  default:
206  const double ai = i - 1;
207  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);
208  break;
209  }
210 
211 
212  int k;
213 
214  for (k = 0; k != number_of_iterations; ++k) {
215 
216  p1 = 0.0;
217  p2 = 1.0;
218 
219  // recurrence relation
220 
221  for (int j = 0; j < n; ++j) {
222  p0 = p1;
223  p1 = p2;
224  p2 = ((2*j + 1 + alf - z) * p1 - (j + alf)*p0) / (j+1);
225  }
226 
227  pp = (n*p2 - (n+alf)*p1) / z;
228 
229  z1 = z;
230  z = z1 - p2/pp;
231 
232  if (fabs(z-z1) < eps)
233  break;
234  }
235 
236  const double y = -tgamma(alf+n) / tgamma((double) n) / (pp*n*p1);
237 
238  insert(JElement2D_t(z,y));
239  }
240  }
241  };
242 
243 
244  /**
245  * Numerical integrator for <tt>W(x) = e^-(x^2)</tt>.
246  *
247  * Gauss-Hermite integration code is taken from reference:
248  * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
249  * Cambridge University Press.
250  */
251  class JGaussHermite :
252  public JQuadrature
253  {
254  public:
255  /**
256  * Constructor.
257  *
258  * \param n number of points
259  * \param eps precision
260  */
261  JGaussHermite(const int n,
262  const double eps = 1.0e-12) :
263  JQuadrature()
264  {
265  resize(n);
266 
267  const double pii = 1.0 / pow(PI,0.25);
268 
269  const int number_of_iterations = 100;
270 
271  const int M = (n + 1) / 2;
272 
273  double p0, p1, p2, pp;
274  double z = 0.0;
275  double z1;
276 
277  for (int i = 0; i < M; ++i) {
278 
279  switch (i) {
280 
281  case 0:
282  z = sqrt((double) (2*n+1)) - 1.85575 * pow((double) (2*n+1),-0.16667);
283  break;
284 
285  case 1:
286  z -= 1.14 * pow((double) n,0.426) / z;
287  break;
288 
289  case 2:
290  z = 1.86*z + 0.86*at( 0 ).getX();
291  break;
292 
293  case 3:
294  z = 1.91*z + 0.91*at( 1 ).getX();
295  break;
296 
297  default:
298  z = 2.00*z + 1.00*at(i-2).getX();
299  break;
300  }
301 
302  for (int k = 0; k != number_of_iterations; ++k) {
303 
304  p1 = 0.0;
305  p2 = pii;
306 
307  // recurrence relation
308 
309  for (int j = 0; j < n; ++j) {
310  p0 = p1;
311  p1 = p2;
312  p2 = z * sqrt(2.0/(double) (j+1)) * p1 - sqrt((double) j / (double) (j+1)) * p0;
313  }
314 
315  pp = sqrt((double) (2*n)) * p1;
316 
317  z1 = z;
318  z = z1 - p2/pp;
319 
320  if (fabs(z-z1) < eps)
321  break;
322  }
323 
324  const double y = 2.0 / (pp*pp);
325 
326  at( i ) = JElement2D_t(-z,y);
327  at(n-i-1) = JElement2D_t(+z,y);
328  }
329  }
330  };
331 
332 
333  /**
334  * Numerical integrator for <tt>W(x) = (1 + g^2 - 2gx)^a</tt>.
335  * For this, <tt>g > 0</tt>.
336  *
337  * Henyey-Greenstein integration points and weights.
338  */
340  public JQuadrature
341  {
342  public:
343  /**
344  * Constructor.
345  *
346  * \param n number of points
347  * \param g angular dependence parameter
348  * \param a power
349  */
350  JHenyeyGreenstein(const int n,
351  const double g,
352  const double a) :
353  JQuadrature()
354  {
355  const double b = -2*g * (a + 1.0);
356  const double ai = 1.0 / (a + 1.0);
357 
358  const double ymin = pow(1.0 + g, 2*(a + 1.0)) / b;
359  const double ymax = pow(1.0 - g, 2*(a + 1.0)) / b;
360 
361  const double dy = (ymax - ymin) / (n + 1);
362 
363  for (double y = ymax - 0.5*dy; y > ymin; y -= dy) {
364 
365  const double v = y*b;
366  const double w = pow(v, ai);
367  const double x = (1.0 + g*g - w) / (2*g);
368  const double dx = pow(v, -a*ai)*dy;
369 
370  insert(JElement2D_t(x,dx));
371  }
372  }
373 
374 
375  /**
376  * Constructor.
377  *
378  * \param n number of points
379  * \param g angular dependence parameter
380  * \param a power
381  * \param xmin minimal value
382  * \param xmax maximal value
383  */
384  JHenyeyGreenstein(const int n,
385  const double g,
386  const double a,
387  const double xmin,
388  const double xmax) :
389  JQuadrature()
390  {
391  const double b = -2*g * (a + 1.0);
392  const double ai = 1.0 / (a + 1.0);
393 
394  const double ymin = pow(1.0 + g*g -2*g*xmin, a + 1.0) / b;
395  const double ymax = pow(1.0 + g*g -2*g*xmax, a + 1.0) / b;
396 
397  const double dy = (ymax - ymin) / (n + 1);
398 
399  for (double y = ymax - 0.5*dy; y > ymin; y -= dy) {
400 
401  const double v = y*b;
402  const double w = pow(v, ai);
403  const double x = (1.0 + g*g - w) / (2*g);
404  const double dx = pow(v, -a*ai)*dy;
405 
406  insert(JElement2D_t(x,dx));
407  }
408  }
409 
410 
411  /**
412  * Constructor for special case where a = -1.
413  *
414  * \param n number of points
415  * \param g angular dependence parameter
416  */
417  JHenyeyGreenstein(const int n,
418  const double g) :
419  JQuadrature()
420  {
421  const double dy = 1.0 / (n + 1);
422  const double gi = log((1.0 + g*g) / (1.0 - g*g)) / (2.0*g);
423 
424  for (double y = 1.0 - 0.5*dy; y > 0.0; y -= dy) {
425 
426  const double v = -y*2.0*g*gi + log(1.0 + g*g);
427  const double w = exp(v);
428  const double x = (1.0 + g*g - w) / (2.0*g);
429  const double dx = w*gi*dy;
430 
431  insert(JElement2D_t(x,dx));
432  }
433  }
434  };
435 
436 
437  /**
438  * Numerical integrator for <tt>W(x) = (1 + g*x*x)</tt>.
439  * For this, <tt>g > 0</tt>.
440  *
441  * Rayleigh integration points and weights.
442  */
443  class JRayleigh :
444  public JQuadrature
445  {
446  public:
447  /**
448  * Constructor.
449  *
450  * \param n number of points
451  * \param g angular dependence parameter
452  */
453  JRayleigh(const int n,
454  const double g) :
455  JQuadrature()
456  {
457  const double dy = 1.0 / (n + 1);
458  const double gi = 3.0/g + 1.0;
459 
460  // t^3 + 3pt + 2q = 0
461 
462  const double p = 1.0/g;
463 
464  for (double y = 0.5*dy; y < 1.0; y += dy) {
465 
466  const double q = 0.5*gi - gi*y;
467 
468  const double b = sqrt(q*q + p*p*p);
469  const double u = pow(-q + b, 1.0/3.0);
470  const double v = pow(+q + b, 1.0/3.0);
471 
472  const double x = u - v;
473  const double dx = (u + v) / (3.0*b);
474 
475  insert(JElement2D_t(x, dx*gi*dy));
476  }
477  }
478  };
479 
480 
481  /**
482  * Numerical integrator for <tt>W(x) = |x| / sqrt(1 - x*x)</tt>
483  *
484  * Co-tangent integration points and weights.
485  */
486  class JCotangent :
487  public JQuadrature
488  {
489  public:
490  /**
491  * Constructor.
492  *
493  * \param n number of points
494  */
495  JCotangent(const int n) :
496  JQuadrature()
497  {
498  for (double ds = 1.0 / (n/2), sb = 0.5*ds; sb < 1.0; sb += ds) {
499 
500  const double cb = sqrt((1.0 + sb)*(1.0 - sb));
501  const double dc = ds*sb/cb;
502 
503  insert(JElement2D_t(+cb, dc));
504  insert(JElement2D_t(-cb, dc));
505  }
506  }
507  };
508 
509 
510  /**
511  * Numerical integrator for <tt>W(x) = |x| / sqrt(1 - x*x), x > 0 </tt> and <tt>W(x) = 1, x <= 0</tt>.
512  *
513  * Bi-tangent integration points and weights.
514  */
515  class JBitangent :
516  public JQuadrature
517  {
518  public:
519  /**
520  * Constructor.
521  *
522  * \param n number of points
523  */
524  JBitangent(const int n) :
525  JQuadrature()
526  {
527  double sb, ds;
528  double cb = 0.0;
529  double dc = 0.0;
530 
531  for (ds = 1.0 / (n/2), sb = 0.5*ds; sb < 1.0; sb += ds) {
532 
533  cb = sqrt((1.0 + sb)*(1.0 - sb));
534  dc = ds*sb/cb;
535 
536  insert(JElement2D_t(cb, dc));
537  }
538 
539  for (dc = (cb + 1.0) / (n/2), cb -= 0.5*dc ; cb > -1.0; cb -= dc) {
540  insert(JElement2D_t(cb, dc));
541  }
542  }
543  };
544 }
545 
546 #endif
data_type w[N+1][M+1]
Definition: JPolint.hh:741
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
Numerical integrator for W(x) = 1.
Definition: JQuadrature.hh:111
TPaveText * p1
do $JPP JMEstimator M
Definition: JMEstimator.sh:37
General purpose class for collection of elements, see: &lt;a href=&quot;JTools.PDF&quot;;&gt;Collection of elements...
Definition: JCollection.hh:73
The elements in a collection are sorted according to their abscissa values and a given distance opera...
JGaussHermite(const int n, const double eps=1.0e-12)
Constructor.
Definition: JQuadrature.hh:261
Numerical integrator for W(x) = e^-(x^2).
Definition: JQuadrature.hh:251
JGaussLegendre(const int n, const double eps=1.0e-12)
Constructor.
Definition: JQuadrature.hh:121
Numerical integrator for W(x) = |x| / sqrt(1 - x*x)
Definition: JQuadrature.hh:486
JGaussLaguerre(const int n, const double alf, const double eps=1.0e-12)
Constructor.
Definition: JQuadrature.hh:183
JQuadrature()
Default constructor.
Definition: JQuadrature.hh:44
void resize(typename container_type::size_type size)
Resize collection.
Definition: JCollection.hh:789
Numerical integrator for W(x) = (1 + g*x*x).
Definition: JQuadrature.hh:443
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` typeset -Z 4 STRING JOpera1D -f hydrophone.root
JHenyeyGreenstein(const int n, const double g, const double a, const double xmin, const double xmax)
Constructor.
Definition: JQuadrature.hh:384
JRayleigh(const int n, const double g)
Constructor.
Definition: JQuadrature.hh:453
Numerical integrator for W(x) = |x| / sqrt(1 - x*x), x &gt; 0 and W(x) = 1, x &lt;= 0. ...
Definition: JQuadrature.hh:515
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:98
JQuadrature(const double xmin, const double xmax, const int nx, JFunction_t integral, const double eps=1.0e-4)
General purpose constructor.
Definition: JQuadrature.hh:62
JBitangent(const int n)
Constructor.
Definition: JQuadrature.hh:524
static const double PI
Mathematical constants.
General purpose class for a collection of sorted elements.
Type definition for numerical integration.
Definition: JQuadrature.hh:37
Numerical integrator for W(x) = x^a e^(-x).
Definition: JQuadrature.hh:172
JHenyeyGreenstein(const int n, const double g)
Constructor for special case where a = -1.
Definition: JQuadrature.hh:417
p2
Definition: module-Z:fit.sh:72
Constants.
Numerical integrator for W(x) = (1 + g^2 - 2gx)^a.
Definition: JQuadrature.hh:339
then JCalibrateToT a
Definition: JTuneHV.sh:116
alias put_queue eval echo n
Definition: qlib.csh:19
2D Element.
Definition: JElement.hh:46
JTOOLS::JElement2D< double, double > JElement2D_t
Definition: JLED.hh:27
int j
Definition: JPolint.hh:666
JCotangent(const int n)
Constructor.
Definition: JQuadrature.hh:495
pair_type insert(const value_type &element)
Insert element.
Definition: JCollection.hh:318
data_type v[N+1][M+1]
Definition: JPolint.hh:740
double u[N+1]
Definition: JPolint.hh:739
JHenyeyGreenstein(const int n, const double g, const double a)
Constructor.
Definition: JQuadrature.hh:350