GammaLib 2.1.0.dev
Loading...
Searching...
No Matches
GIntegral.cpp
Go to the documentation of this file.
1/***************************************************************************
2 * GIntegral.cpp - Integration class *
3 * ----------------------------------------------------------------------- *
4 * copyright (C) 2010-2023 by Juergen Knoedlseder *
5 * ----------------------------------------------------------------------- *
6 * *
7 * This program is free software: you can redistribute it and/or modify *
8 * it under the terms of the GNU General Public License as published by *
9 * the Free Software Foundation, either version 3 of the License, or *
10 * (at your option) any later version. *
11 * *
12 * This program is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have received a copy of the GNU General Public License *
18 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19 * *
20 ***************************************************************************/
21/**
22 * @file GIntegral.cpp
23 * @brief Integration class implementation
24 * @author Juergen Knoedlseder
25 */
26
27/* __ Includes ___________________________________________________________ */
28#include <cmath> // std::abs()
29#include <vector>
30#include <algorithm> // std::sort
31#include <limits> // numeric_limits
32#include "GIntegral.hpp"
33#include "GException.hpp"
34#include "GTools.hpp"
35#include "GFunction.hpp"
36
37/* __ Method name definitions ____________________________________________ */
38#define G_ROMBERG "GIntegral::romberg(double&, double&, int&)"
39#define G_TRAPZD "GIntegral::trapzd(double&, double&, int&, double)"
40#define G_POLINT "GIntegral::polint(double*, double*, int, double, double*)"
41
42/* __ Macros _____________________________________________________________ */
43
44/* __ Coding definitions _________________________________________________ */
45
46/* __ Debug definitions __________________________________________________ */
47
48/* __ Constants __________________________________________________________ */
49namespace gammalib {
50
51 // Gauss-Kronrod abscissae, common to the 10-, 21-, 43- and 87-point rule
52 const double gkx1[5] = {
53 0.973906528517171720077964012084452,
54 0.865063366688984510732096688423493,
55 0.679409568299024406234327365114874,
56 0.433395394129247190799265943165784,
57 0.148874338981631210884826001129720
58 };
59
60 // Gauss-Kronrod weights of the 10-point rule
61 const double gkw10[5] = {
62 0.066671344308688137593568809893332,
63 0.149451349150580593145776339657697,
64 0.219086362515982043995534934228163,
65 0.269266719309996355091226921569469,
66 0.295524224714752870173892994651338
67 };
68
69 // Gauss-Kronrod abscissae, common to the 21-, 43- and 87-point rule
70 const double gkx2[5] = {
71 0.995657163025808080735527280689003,
72 0.930157491355708226001207180059508,
73 0.780817726586416897063717578345042,
74 0.562757134668604683339000099272694,
75 0.294392862701460198131126603103866
76 };
77
78 // Gauss-Kronrod weights of the 21-point rule for abscissae gkx1
79 const double gkw21a[5] = {
80 0.032558162307964727478818972459390,
81 0.075039674810919952767043140916190,
82 0.109387158802297641899210590325805,
83 0.134709217311473325928054001771707,
84 0.147739104901338491374841515972068
85 };
86
87 // Gauss-Kronrod weights of the 21-point rule for abscissae gkx2
88 const double gkw21b[6] = {
89 0.011694638867371874278064396062192,
90 0.054755896574351996031381300244580,
91 0.093125454583697605535065465083366,
92 0.123491976262065851077958109831074,
93 0.142775938577060080797094273138717,
94 0.149445554002916905664936468389821
95 };
96
97 // Gauss-Kronrod abscissae, common to the 43- and 87-point rule
98 const double gkx3[11] = {
99 0.999333360901932081394099323919911,
100 0.987433402908088869795961478381209,
101 0.954807934814266299257919200290473,
102 0.900148695748328293625099494069092,
103 0.825198314983114150847066732588520,
104 0.732148388989304982612354848755461,
105 0.622847970537725238641159120344323,
106 0.499479574071056499952214885499755,
107 0.364901661346580768043989548502644,
108 0.222254919776601296498260928066212,
109 0.074650617461383322043914435796506
110 };
111
112 // Gauss-Kronrod weights of the 43-point rule for abscissae gkx1, gkx3
113 const double gkw43a[10] = {
114 0.016296734289666564924281974617663,
115 0.037522876120869501461613795898115,
116 0.054694902058255442147212685465005,
117 0.067355414609478086075553166302174,
118 0.073870199632393953432140695251367,
119 0.005768556059769796184184327908655,
120 0.027371890593248842081276069289151,
121 0.046560826910428830743339154433824,
122 0.061744995201442564496240336030883,
123 0.071387267268693397768559114425516
124 };
125
126 // Gauss-Kronrod weights of the 43-point formula for abscissae gkx3
127 const double gkw43b[12] = {
128 0.001844477640212414100389106552965,
129 0.010798689585891651740465406741293,
130 0.021895363867795428102523123075149,
131 0.032597463975345689443882222526137,
132 0.042163137935191811847627924327955,
133 0.050741939600184577780189020092084,
134 0.058379395542619248375475369330206,
135 0.064746404951445885544689259517511,
136 0.069566197912356484528633315038405,
137 0.072824441471833208150939535192842,
138 0.074507751014175118273571813842889,
139 0.074722147517403005594425168280423
140 };
141
142 // Gauss-Kronrod abscissae, of the 87-point rule
143 const double gkx4[22] = {
144 0.999902977262729234490529830591582,
145 0.997989895986678745427496322365960,
146 0.992175497860687222808523352251425,
147 0.981358163572712773571916941623894,
148 0.965057623858384619128284110607926,
149 0.943167613133670596816416634507426,
150 0.915806414685507209591826430720050,
151 0.883221657771316501372117548744163,
152 0.845710748462415666605902011504855,
153 0.803557658035230982788739474980964,
154 0.757005730685495558328942793432020,
155 0.706273209787321819824094274740840,
156 0.651589466501177922534422205016736,
157 0.593223374057961088875273770349144,
158 0.531493605970831932285268948562671,
159 0.466763623042022844871966781659270,
160 0.399424847859218804732101665817923,
161 0.329874877106188288265053371824597,
162 0.258503559202161551802280975429025,
163 0.185695396568346652015917141167606,
164 0.111842213179907468172398359241362,
165 0.037352123394619870814998165437704
166 };
167
168 // Gauss-Kronrod weights of the 87-point rule for abscissae gkx1, gkx2, gkx3
169 const double gkw87a[21] = {
170 0.008148377384149172900002878448190,
171 0.018761438201562822243935059003794,
172 0.027347451050052286161582829741283,
173 0.033677707311637930046581056957588,
174 0.036935099820427907614589586742499,
175 0.002884872430211530501334156248695,
176 0.013685946022712701888950035273128,
177 0.023280413502888311123409291030404,
178 0.030872497611713358675466394126442,
179 0.035693633639418770719351355457044,
180 0.000915283345202241360843392549948,
181 0.005399280219300471367738743391053,
182 0.010947679601118931134327826856808,
183 0.016298731696787335262665703223280,
184 0.021081568889203835112433060188190,
185 0.025370969769253827243467999831710,
186 0.029189697756475752501446154084920,
187 0.032373202467202789685788194889595,
188 0.034783098950365142750781997949596,
189 0.036412220731351787562801163687577,
190 0.037253875503047708539592001191226
191 };
192
193 // Gauss-Kronrod weights of the 87-point formula for abscissae gkx4
194 const double gkw87b[23] = {
195 0.000274145563762072350016527092881,
196 0.001807124155057942948341311753254,
197 0.004096869282759164864458070683480,
198 0.006758290051847378699816577897424,
199 0.009549957672201646536053581325377,
200 0.012329447652244853694626639963780,
201 0.015010447346388952376697286041943,
202 0.017548967986243191099665352925900,
203 0.019938037786440888202278192730714,
204 0.022194935961012286796332102959499,
205 0.024339147126000805470360647041454,
206 0.026374505414839207241503786552615,
207 0.028286910788771200659968002987960,
208 0.030052581128092695322521110347341,
209 0.031646751371439929404586051078883,
210 0.033050413419978503290785944862689,
211 0.034255099704226061787082821046821,
212 0.035262412660156681033782717998428,
213 0.036076989622888701185500318003895,
214 0.036698604498456094498018047441094,
215 0.037120549269832576114119958413599,
216 0.037334228751935040321235449094698,
217 0.037361073762679023410321241766599
218 };
219
220} // end gammalib namespace
221
222
223/*==========================================================================
224 = =
225 = Constructors/destructors =
226 = =
227 ==========================================================================*/
228
229/***********************************************************************//**
230 * @brief Void constructor
231 ***************************************************************************/
233{
234 // Initialise members
235 init_members();
236
237 // Return
238 return;
239}
240
241
242/***********************************************************************//**
243 * @brief Function kernel constructor
244 *
245 * @param[in] kernel Pointer to function kernel.
246 *
247 * The function kernel constructor assigns the function kernel pointer in
248 * constructing the object.
249 ***************************************************************************/
251{
252 // Initialise members
253 init_members();
254
255 // Set function kernel
257
258 // Return
259 return;
260}
261
262
263/***********************************************************************//**
264 * @brief Copy constructor
265 *
266 * @param[in] integral Integral.
267 ***************************************************************************/
269{
270 // Initialise members
271 init_members();
272
273 // Copy members
274 copy_members(integral);
275
276 // Return
277 return;
278}
279
280
281/***********************************************************************//**
282 * @brief Destructor
283 ***************************************************************************/
285{
286 // Free members
287 free_members();
288
289 // Return
290 return;
291}
292
293
294/*==========================================================================
295 = =
296 = Operators =
297 = =
298 ==========================================================================*/
299
300/***********************************************************************//**
301 * @brief Assignment operator
302 *
303 * @param[in] integral Integral.
304 * @return Integral.
305 ***************************************************************************/
307{
308 // Execute only if object is not identical
309 if (this != &integral) {
310
311 // Free members
312 free_members();
313
314 // Initialise integral
315 init_members();
316
317 // Copy members
318 copy_members(integral);
319
320 } // endif: object was not identical
321
322 // Return
323 return *this;
324}
325
326
327/*==========================================================================
328 = =
329 = Public methods =
330 = =
331 ==========================================================================*/
332
333/***********************************************************************//**
334 * @brief Clear integral
335 ***************************************************************************/
337{
338 // Free members
339 free_members();
340
341 // Initialise private members
342 init_members();
343
344 // Return
345 return;
346}
347
348
349/***********************************************************************//**
350 * @brief Clone integral
351 *
352 * @return Pointer to deep copy of integral.
353 ***************************************************************************/
355{
356 return new GIntegral(*this);
357}
358
359
360/***********************************************************************//**
361 * @brief Perform Romberg integration
362 *
363 * @param[in] bounds Integration boundaries.
364 * @param[in] order Integration order (default: 5)
365 *
366 * @exception GException::invalid_argument
367 * Integration order incompatible with number of iterations.
368 *
369 * Returns the integral of the integrand, computed over a number of
370 * intervals [a0,a1], [a1,a2], ... that are given as an unordered vector
371 * by the @p bounds argument.
372 *
373 * Integration is performed by Romberg's method of order 2*order, where
374 *
375 * order=1 is equivalent to the trapezoidal rule,
376 * order=2 is equivalent to Simpson's rule, and
377 * order=3 is equivalent to Boole's rule.
378 *
379 * The number of iterations is limited by m_max_iter. m_eps specifies the
380 * requested fractional accuracy. By default it is set to 1e-6.
381 ***************************************************************************/
382double GIntegral::romberg(std::vector<double> bounds, const int& order)
383{
384 // Sort integration boundaries in ascending order
385 std::sort(bounds.begin(), bounds.end());
386
387 // Initialise integral
388 double value = 0.0;
389
390 // Initialise integration status information
391 int calls = 0;
392
393 // Add integral of all intervals
394 for (int i = 0; i < bounds.size()-1; ++i) {
395 value += romberg(bounds[i], bounds[i+1], order);
396 calls += m_calls;
397 }
398
399 // Set integration status information
400 m_calls = calls;
401
402 // Return value
403 return value;
404}
405
406
407/***********************************************************************//**
408 * @brief Perform Romberg integration
409 *
410 * @param[in] a Left integration boundary.
411 * @param[in] b Right integration boundary.
412 * @param[in] order Integration order (default: 5)
413 *
414 * @exception GException::invalid_argument
415 * Integration order incompatible with number of iterations.
416 *
417 * Returns the integral of the integrand from a to b. Integration is
418 * performed by Romberg's method of order 2*order, where
419 *
420 * order=1 is equivalent to the trapezoidal rule,
421 * order=2 is equivalent to Simpson's rule, and
422 * order=3 is equivalent to Boole's rule.
423 *
424 * The number of iterations is limited by m_max_iter. m_eps specifies the
425 * requested fractional accuracy. By default it is set to 1e-6.
426 ***************************************************************************/
427double GIntegral::romberg(const double& a, const double& b, const int& order)
428{
429 // Initialise result and status
430 double result = 0.0;
431
432 // Initialise integration status information
433 m_isvalid = true;
434 m_calls = 0;
435 m_has_abserr = false;
436 m_has_relerr = false;
437
438 // Continue only if integration range is valid
439 if (b > a) {
440
441 // Initialise variables
442 bool converged = false;
443 double dss = 0.0;
444
445 // Determine (maximum) number of iterations
447
448 // Check whether maximum number of iterations is compliant with
449 // order
450 if (order > max_iter) {
451 std::string msg = "Requested integration order "+
452 gammalib::str(order)+" is larger than the "
453 "maximum number of iterations "+
454 gammalib::str(max_iter)+". Either reduce the "
455 "integration order or increase the (maximum) "
456 "number of iterations.";
458 }
459
460 // Allocate temporal storage
461 double* s = new double[max_iter+2];
462 double* h = new double[max_iter+2];
463
464 // Initialise step size
465 h[1] = 1.0;
466 s[0] = 0.0;
467
468 // Iterative loop
469 for (m_iter = 1; m_iter <= max_iter; ++m_iter) {
470
471 // Integration using Trapezoid rule
472 s[m_iter] = trapzd(a, b, m_iter, s[m_iter-1]);
473
474 // Compile option: Check for NaN/Inf
475 #if defined(G_NAN_CHECK)
476 if (is_notanumber(s[m_iter]) || is_infinite(s[m_iter])) {
477 m_message = "*** ERROR: GIntegral::romberg"
478 "(a="+gammalib::str(a)+", b="+gammalib::str(b)+""
479 ", k="+gammalib::str(k)+"): NaN/Inf encountered"
480 " (s["+gammalib::str(m_iter)+"]="
481 ""+gammalib::str(s[m_iter])+")";
482 std::cout << m_message << std::endl;
483 m_isvalid = false;
484 }
485 #endif
486
487 // Starting from iteration order on, use polynomial interpolation
488 if (m_iter >= order) {
489
490 // Compute result using polynom interpolation
491 result = polint(&h[m_iter-order], &s[m_iter-order],
492 order, 0.0, &dss);
493
494 // If a fixed number of iterations has been requested then
495 // check whether we reached the final one; otherwise check
496 // whether we reached the requested precision.
497 if (m_fix_iter > 0) {
498 if (m_iter == max_iter) {
499 converged = true;
500 }
501 }
502 else if (std::abs(dss) <= m_eps * std::abs(result)) {
503 converged = true;
504 }
505 if (converged) {
506 m_has_abserr = true;
507 m_abserr = std::abs(dss);
508 if (std::abs(result) > 0) {
509 m_has_relerr = true;
510 m_relerr = m_abserr / std::abs(result);
511 }
512 break;
513 }
514
515 } // endif: polynomial interpolation performed
516
517 // Reduce step size
518 h[m_iter+1]= 0.25 * h[m_iter];
519
520 } // endfor: iterative loop
521
522 // Free temporal storage
523 delete [] s;
524 delete [] h;
525
526 // Set status and optionally dump warning
527 if (!converged) {
528 m_isvalid = false;
529 m_message = "Integration uncertainty "+
530 gammalib::str(std::abs(dss))+
531 " exceeds absolute tolerance of "+
532 gammalib::str(m_eps * std::abs(result))+
533 " after "+gammalib::str(m_iter)+
534 " iterations. Result "+
535 gammalib::str(result)+
536 " is inaccurate.";
537 if (!m_silent) {
538 std::string origin = "GIntegral::romberg("+
539 gammalib::str(a)+", "+
540 gammalib::str(b)+", "+
541 gammalib::str(order)+")";
543 }
544 }
545
546 } // endif: integration range was valid
547
548 // Return result
549 return result;
550}
551
552
553/***********************************************************************//**
554 * @brief Perform Trapezoidal integration
555 *
556 * @param[in] a Left integration boundary.
557 * @param[in] b Right integration boundary.
558 * @param[in] n Number of steps.
559 * @param[in] result Result from a previous trapezoidal integration step.
560 * @return Integration results.
561 *
562 * @exception GException::invalid_value
563 * Function kernel not set.
564 *
565 * The method performs a trapezoidal integration of the function kernel for
566 * the integration interval [a,b].
567 *
568 * If @p n = 1 the integral is computed using
569 *
570 * \f[
571 * \int_a^b f(x) \, dx = \frac{1}{2} (b-a) (f(a) + f(b))
572 * \f]
573 *
574 * For @p n > 1 the integral is computed using
575 *
576 * \f[
577 * \int_a^b f(x) \, dx = \frac{1}{2} \left[{\tt result} +
578 * \frac{b-a}{2^{n-1}}
579 * \sum_{i=0}^{2^{n-1}-1} f \left( a + (0.5+i) \frac{b-a}{2^{n-1}} \right) \right]
580 *
581 * \f]
582 *
583 * where \f${\tt result}\f$ is the integration result from a previous call
584 * to the method with @p n = @p n - 1.
585 ***************************************************************************/
586double GIntegral::trapzd(const double& a, const double& b, const int& n,
587 double result)
588{
589 // Throw an exception if the instance has no kernel
590 if (m_kernel == NULL) {
591 std::string msg = "Function kernel not set. Please set a function "
592 "kernel before calling the method.";
594 }
595
596 // Handle case of identical boundaries
597 if (a == b) {
598 result = 0.0;
599 }
600
601 // ... otherwise use trapeziodal rule
602 else {
603
604 // Case A: Only a single step is requested
605 if (n == 1) {
606
607 // Evaluate integrand at boundaries
608 double y_a = m_kernel->eval(a);
609 double y_b = m_kernel->eval(b);
610 m_calls += 2;
611
612 // Compute result
613 result = 0.5*(b-a)*(y_a + y_b);
614
615 } // endif: only a single step was requested
616
617 // Case B: More than a single step is requested
618 else {
619
620 // Compute step level 2^(n-1)
621 int it = 1;
622 for (int j = 1; j < n-1; ++j) {
623 it <<= 1;
624 }
625
626 // Verify that step level is valid
627 if (it == 0) {
628 m_isvalid = false;
629 m_message = "Invalid step level "+gammalib::str(it)+
630 " encountered for"
631 " a="+gammalib::str(a)+
632 ", b="+gammalib::str(b)+
633 ", n="+gammalib::str(n)+
634 ", result="+gammalib::str(result)+
635 ". Looks like n is too large.";
636 if (!m_silent) {
638 }
639 }
640
641 // Set step size
642 double tnm = double(it);
643 double del = (b-a)/tnm;
644
645 // Verify that step is >0
646 if (del == 0) {
647 m_isvalid = false;
648 m_message = "Invalid step size "+gammalib::str(del)+
649 " encountered for"
650 " a="+gammalib::str(a)+
651 ", b="+gammalib::str(b)+
652 ", n="+gammalib::str(n)+
653 ", result="+gammalib::str(result)+
654 ". Step is too small to make sense.";
655 if (!m_silent) {
657 }
658 }
659
660 // Sum up values
661 double x = a + 0.5*del;
662 double sum = 0.0;
663 for (int j = 0; j < it; ++j, x+=del) {
664
665 // Evaluate integrand
666 double y = m_kernel->eval(x);
667 m_calls++;
668
669 // Add integrand
670 sum += y;
671
672 } // endfor: looped over steps
673
674 // Set result
675 result = 0.5*(result + (b-a)*sum/tnm);
676 }
677
678 } // endelse: trapeziodal rule was applied
679
680 // Return result
681 return result;
682}
683
684
685/***********************************************************************//**
686 * @brief Adaptive Simpson's integration
687 *
688 * @param[in] a Left integration boundary.
689 * @param[in] b Right integration boundary.
690 *
691 * Integrates the function using an adaptive Simpson's rule. The initial
692 * interval [a,b] is split into two sub-intervals [a,c] and [c,b] for which
693 * the integral is computed using
694 *
695 * \f[
696 * \frac{b-a}{6} f(a) + 4f(c) + f(b)
697 * \f]
698 *
699 * where \f$c=(a+b)/2\f$ is the mid-point of interval [a,b]. Each
700 * sub-interval is then recursively divided into sub-interval and the process
701 * is repeated. Dividing of sub-intervals is stopped when the difference
702 * between subsequent intervals falls below the relative tolerance specified
703 * by eps(). The maximum recursion depth is set by the max_iter() method.
704 *
705 * I almost do not dare to confess, but the code has been taken from
706 * http://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method
707 * It's really pretty simple ...
708 ***************************************************************************/
709double GIntegral::adaptive_simpson(const double& a, const double& b) const
710{
711 // Initialise integration status information
712 m_isvalid = true;
713 m_calls = 0;
715 m_has_abserr = false;
716 m_has_relerr = false;
717
718 // Compute mid-point c
719 double c = 0.5*(a + b); //!< Mid-point of interval [a,b]
720 double h = b - a; //!< Length of interval [a,b]
721
722 // Evaluate function at boundaries and mid-point c
723 double fa = m_kernel->eval(a);
724 double fb = m_kernel->eval(b);
725 double fc = m_kernel->eval(c);
726 m_calls += 3;
727
728 // Compute integral using Simpson's rule
729 double S = (h/6.0) * (fa + 4.0*fc + fb);
730
731 // Initialise absolute precision
732 double epsilon = (std::abs(S) > 0) ? m_eps * S : m_eps;
733
734 // Call recursive auxiliary function
735 double value = adaptive_simpson_aux(a, b, epsilon, S, fa, fb, fc, m_max_iter);
736
737 // Deduce the number of iterations from the iteration counter
739
740 // If result is not valid, set and output status message
741 if (!m_isvalid) {
742 m_message = "Integration uncertainty exceeds relative tolerance "
743 "of "+gammalib::str(m_eps)+" and absolute tolerance of "
744 ""+gammalib::str(epsilon)+" after "+gammalib::str(m_iter)+
745 " iterations and "+gammalib::str(m_calls)+" function "
746 "calls. Result "+gammalib::str(value)+" inaccurate.";
747 if (!m_silent) {
748 std::string origin = "GIntegral::adaptive_simpson("+
749 gammalib::str(a)+", "+
750 gammalib::str(b)+")";
752 }
753 }
754
755 // Return result
756 return value;
757}
758
759
760/***********************************************************************//**
761 * @brief Adaptive Gauss-Lobatto-Kronrod integration
762 *
763 * @param[in] a Left integration boundary.
764 * @param[in] b Right integration boundary.
765 * @return Integral of kernel over interval [a,b].
766 *
767 * Integrates the function using an adaptive Gauss-Lobatto method with a
768 * Kronrod extension.
769 ***************************************************************************/
770double GIntegral::adaptive_gauss_kronrod(const double& a, const double& b) const
771{
772 // Set constants
773 const double alpha = std::sqrt(2.0/3.0);
774 const double beta = 1.0/std::sqrt(5.0);
775 const double x1 = 0.942882415695480;
776 const double x2 = 0.641853342345781;
777 const double x3 = 0.236383199662150;
778 const double x[12] = {0,-x1,-alpha,-x2,-beta,-x3,0.0,x3,beta,x2,alpha,x1};
779 const double eps = std::numeric_limits<double>::epsilon();
780
781 // Set tolerance, assuring that the tolerance is not smaller than the
782 // machine precision
783 double tol = (m_eps < 10.0*eps) ? 10.0 * eps : m_eps;
784
785 // Initialise integration status information
786 m_isvalid = true;
787 m_iter = 0;
788 m_calls = 0;
789 m_has_abserr = false;
790 m_has_relerr = false;
791 m_terminate = true;
792
793 // Allocate result storage array
794 double y[13];
795
796 // Compute midpoint and step size
797 double m = 0.5 * (a+b);
798 double h = 0.5 * (b-a);
799
800 // Evaluate function at end points
801 double fa = m_kernel->eval(a);
802 double fb = m_kernel->eval(b);
803 m_calls += 2;
804
805 // Store end points
806 y[0] = fa;
807 y[12] = fb;
808
809 // Evaluate function at intermediate points
810 for (int i = 1; i < 12; ++i) {
811 y[i] = m_kernel->eval(m + x[i]*h);
812 m_calls++;
813 }
814
815 // 4-point Gauss-Lobatto formula
816 double i2 = (h/6.0)*(y[0]+y[12]+5.0*(y[4]+y[8]));
817
818 // 7-point Konrod extension
819 double i1 = (h/1470.0) * ( 77.0 * (y[0]+y[12])+
820 432.0 * (y[2]+y[10])+
821 625.0 * (y[4]+y[8])+
822 672.0 * y[6]);
823
824 // 13-point Konrod extension
825 double is = h*(0.0158271919734802 * (y[0]+y[12])+
826 0.0942738402188500 * (y[1]+y[11])+
827 0.155071987336585 * (y[2]+y[10])+
828 0.188821573960182 * (y[3]+y[9])+
829 0.199773405226859 * (y[4]+y[8])+
830 0.224926465333340 * (y[5]+y[7])+
831 0.242611071901408 * y[6]);
832
833 // Estimate errors
834 double erri1 = std::abs(i1-is);
835 double erri2 = std::abs(i2-is);
836
837 // Check on convergence
838 double r = (erri2 != 0.0) ? erri1/erri2 : 1.0;
839 double toler = (r > 0.0 && r < 1.0) ? tol/r : tol;
840 if (is == 0.0) {
841 is = b-a;
842 }
843
844 // Call helper
845 is = std::abs(is);
846 double value = adaptive_gauss_kronrod_aux(a, b, fa, fb, is, toler);
847
848 // Return result
849 return value;
850}
851
852
853/***********************************************************************//**
854 * @brief Gauss-Kronrod integration
855 *
856 * @param[in] a Left integration boundary.
857 * @param[in] b Right integration boundary.
858 *
859 ***************************************************************************/
860double GIntegral::gauss_kronrod(const double& a, const double& b) const
861{
862 // Initialise integration status information
863 m_isvalid = true;
864 m_iter = 0;
865 m_calls = 0;
866 m_has_abserr = false;
867 m_has_relerr = false;
868
869 // Initialise integration result
870 double result = 0.0;
871
872 // Allocate some arrays
873 double fv1[5];
874 double fv2[5];
875 double fv3[5];
876 double fv4[5];
877 double savfun[21];
878
879 // Main code loop (so that we can exit using break)
880 do {
881
882 // Tolerance check
883 if (m_eps < 1.12e-14) {
884 m_isvalid = false;
885 m_message = "Requested relative tolerance of "+gammalib::str(m_eps)+
886 " cannot be acheived. Please relax the integration "
887 "precision.";
888 if (!m_silent) {
889 std::string origin = "GIntegral::gauss_kronrod("+
890 gammalib::str(a)+", "+
891 gammalib::str(b)+")";
893 }
894 }
895
896 // Compute function at mid-point
897 double h = 0.5 * (b - a);
898 double abs_h = std::abs(h);
899 double c = 0.5 * (b + a);
900 double f_c = m_kernel->eval(c);
901 m_calls++;
902
903 // Compute the integral using the 10- and 21-point formulae
904 m_iter++;
905 double res10 = 0;
906 double res21 = gammalib::gkw21b[5] * f_c;
907 double resabs = gammalib::gkw21b[5] * std::abs(f_c);
908 for (int k = 0; k < 5; ++k) {
909 double x = h * gammalib::gkx1[k];
910 double fval1 = m_kernel->eval(c+x);
911 double fval2 = m_kernel->eval(c-x);
912 double fval = fval1 + fval2;
913 m_calls += 2;
914 res10 += gammalib::gkw10[k] * fval;
915 res21 += gammalib::gkw21a[k] * fval;
916 resabs += gammalib::gkw21a[k] * (std::abs(fval1) + std::abs(fval2));
917 savfun[k] = fval;
918 fv1[k] = fval1;
919 fv2[k] = fval2;
920 }
921 for (int k = 0; k < 5; ++k) {
922 double x = h * gammalib::gkx2[k];
923 double fval1 = m_kernel->eval(c+x);
924 double fval2 = m_kernel->eval(c-x);
925 double fval = fval1 + fval2;
926 m_calls += 2;
927 res21 += gammalib::gkw21b[k] * fval;
928 resabs += gammalib::gkw21b[k] * (std::abs(fval1) + std::abs(fval2));
929 savfun[k+5] = fval;
930 fv3[k] = fval1;
931 fv4[k] = fval2;
932 }
933 resabs *= abs_h;
934 double mean = 0.5 * res21;
935 double resasc = gammalib::gkw21b[5] * std::abs(f_c - mean);
936 for (int k = 0; k < 5; ++k) {
937 resasc += (gammalib::gkw21a[k] *
938 (std::abs(fv1[k] - mean) + std::abs(fv2[k] - mean)) +
940 (std::abs(fv3[k] - mean) + std::abs(fv4[k] - mean)));
941 }
942 resasc *= abs_h ;
943 result = res21 * h;
944 double error = rescale_error((res21 - res10) * h, resabs, resasc);
945
946 // Test for convergence */
947 if (error < m_eps * std::abs(result)) {
948 m_has_abserr = true;
949 m_abserr = error;
950 if (std::abs(result) > 0) {
951 m_has_relerr = true;
952 m_relerr = error / std::abs(result);
953 }
954 break;
955 }
956
957 // Compute the integral using the 43-point formula
958 m_iter++;
959 double res43 = gammalib::gkw43b[11] * f_c;
960 for (int k = 0; k < 10; ++k) {
961 res43 += savfun[k] * gammalib::gkw43a[k];
962 }
963 for (int k = 0; k < 11; ++k) {
964 double x = h * gammalib::gkx3[k];
965 double fval = (m_kernel->eval(c+x) +
966 m_kernel->eval(c-x));
967 m_calls += 2;
968 res43 += fval * gammalib::gkw43b[k];
969 savfun[k+10] = fval;
970 }
971 result = res43 * h;
972
973 // Test for convergence */
974 error = rescale_error((res43 - res21) * h, resabs, resasc);
975 if (error < m_eps * std::abs(result)) {
976 m_has_abserr = true;
977 m_abserr = error;
978 if (std::abs(result) > 0) {
979 m_has_relerr = true;
980 m_relerr = error / std::abs(result);
981 }
982 break;
983 }
984
985 // Compute the integral using the 87-point formula
986 m_iter++;
987 double res87 = gammalib::gkw87b[22] * f_c;
988 for (int k = 0; k < 21; ++k) {
989 res87 += savfun[k] * gammalib::gkw87a[k];
990 }
991 for (int k = 0; k < 22; ++k) {
992 double x = h * gammalib::gkx4[k];
993 res87 += gammalib::gkw87b[k] *
994 (m_kernel->eval(c+x) +
995 m_kernel->eval(c-x));
996 m_calls += 2;
997 }
998 result = res87 * h ;
999
1000 // Test for convergence */
1001 error = rescale_error ((res87 - res43) * h, resabs, resasc);
1002 if (error < m_eps * std::abs(result)) {
1003 m_has_abserr = true;
1004 m_abserr = error;
1005 if (std::abs(result) > 0) {
1006 m_has_relerr = true;
1007 m_relerr = error / std::abs(result);
1008 }
1009 break;
1010 }
1011
1012 // Failed to converge
1013 m_isvalid = false;
1014 m_message = "Integration uncertainty "+gammalib::str(error)+" exceeds "
1015 "absolute tolerance of "+
1016 gammalib::str(m_eps * std::abs(result))+" after "+
1017 gammalib::str(m_iter)+" iterations and "+
1018 gammalib::str(m_calls)+" function calls. Result "+
1019 gammalib::str(result)+" inaccurate.";
1020 if (!m_silent) {
1021 std::string origin = "GIntegral::gauss_kronrod("+
1022 gammalib::str(a)+", "+
1023 gammalib::str(b)+")";
1025 }
1026
1027 } while (false); // end of main loop
1028
1029 // Return result
1030 return result;
1031}
1032
1033
1034/***********************************************************************//**
1035 * @brief Print integral information
1036 *
1037 * @param[in] chatter Chattiness.
1038 * @return String containing integral information.
1039 ***************************************************************************/
1040std::string GIntegral::print(const GChatter& chatter) const
1041{
1042 // Initialise result string
1043 std::string result;
1044
1045 // Continue only if chatter is not silent
1046 if (chatter != SILENT) {
1047
1048 // Append header
1049 result.append("=== GIntegral ===");
1050
1051 // Append information
1052 result.append("\n"+gammalib::parformat("Relative precision"));
1053 result.append(gammalib::str(eps()));
1054 if (m_has_abserr) {
1055 result.append("\n"+gammalib::parformat("Absolute error"));
1056 result.append(gammalib::str(m_abserr));
1057 }
1058 if (m_has_relerr) {
1059 result.append("\n"+gammalib::parformat("Relative error"));
1060 result.append(gammalib::str(m_relerr));
1061 }
1062 result.append("\n"+gammalib::parformat("Function calls"));
1063 result.append(gammalib::str(calls()));
1064 result.append("\n"+gammalib::parformat("Iterations"));
1065 result.append(gammalib::str(iter()));
1066 if (m_fix_iter > 0) {
1067 result.append(" (fixed: ");
1068 result.append(gammalib::str(fixed_iter()));
1069 result.append(")");
1070 }
1071 else {
1072 result.append(" (maximum: ");
1073 result.append(gammalib::str(max_iter()));
1074 result.append(")");
1075 }
1076
1077 // Append status information
1078 result.append("\n"+gammalib::parformat("Status"));
1079 if (is_valid()) {
1080 result.append("Result accurate.");
1081 }
1082 else {
1083 result.append(message());
1084 }
1085 if (silent()) {
1086 result.append("\n"+gammalib::parformat("Warnings")+"suppressed");
1087 }
1088 else {
1089 result.append("\n"+gammalib::parformat("Warnings"));
1090 result.append("in standard output");
1091 }
1092
1093 } // endif: chatter was not silent
1094
1095 // Return result
1096 return result;
1097}
1098
1099
1100/*==========================================================================
1101 = =
1102 = Protected methods =
1103 = =
1104 ==========================================================================*/
1105
1106/***********************************************************************//**
1107 * @brief Initialise class members
1108 ***************************************************************************/
1110{
1111 // Initialise members
1112 m_kernel = NULL;
1113 m_eps = 1.0e-6;
1114 m_max_iter = 20;
1115 m_fix_iter = 0;
1116 m_message.clear();
1117 m_silent = false;
1118
1119 // Initialise results
1120 m_iter = 0;
1121 m_calls = 0;
1122 m_isvalid = true;
1123 m_has_abserr = false;
1124 m_has_relerr = false;
1125 m_abserr = 0.0;
1126 m_relerr = 0.0;
1127
1128 // Return
1129 return;
1130}
1131
1132
1133/***********************************************************************//**
1134 * @brief Copy class members
1135 *
1136 * @param[in] integral Integral.
1137 ***************************************************************************/
1139{
1140 // Copy attributes
1141 m_kernel = integral.m_kernel;
1142 m_eps = integral.m_eps;
1143 m_max_iter = integral.m_max_iter;
1144 m_fix_iter = integral.m_fix_iter;
1145 m_message = integral.m_message;
1146 m_silent = integral.m_silent;
1147
1148 // Copy results
1149 m_iter = integral.m_iter;
1150 m_calls = integral.m_calls;
1151 m_isvalid = integral.m_isvalid;
1152 m_has_abserr = integral.m_has_abserr;
1153 m_has_relerr = integral.m_has_relerr;
1154 m_abserr = integral.m_abserr;
1155 m_relerr = integral.m_relerr;
1156
1157 // Return
1158 return;
1159}
1160
1161
1162/***********************************************************************//**
1163 * @brief Delete class members
1164 ***************************************************************************/
1166{
1167 // Return
1168 return;
1169}
1170
1171
1172/***********************************************************************//**
1173 * @brief Perform Polynomial interpolation
1174 *
1175 * @param[in] xa Pointer to array of X values.
1176 * @param[in] ya Pointer to array of Y values.
1177 * @param[in] n Number of elements in arrays.
1178 * @param[in] x X value at which interpolations should be performed.
1179 * @param[out] dy Error estimate for interpolated values.
1180 *
1181 * Given arrays xa[1,..,n] and ya[1,..,n], and given a value x, this
1182 * method returns a value y, and an error estimate dy. If P(x) is the
1183 * polynomial of degree n-1, then the returned value y=P(x).
1184 *
1185 * @todo Implement exceptions instead of screen dump.
1186 * @todo Use std::vector for xa and ya and start at 0
1187 ***************************************************************************/
1188double GIntegral::polint(double* xa, double* ya, int n, double x, double* dy)
1189{
1190 // Initialise result
1191 double y = 0.0;
1192
1193 // Allocate temporary memory
1194 std::vector<double> c(n, 0.0);
1195 std::vector<double> d(n, 0.0);
1196
1197 // Compute initial distance to first node
1198 double dif = std::abs(x-xa[1]);
1199
1200 // Find index ns of the closest table entry
1201 int ns = 0;
1202 for (int i = 0; i < n; ++i) {
1203 double dift = std::abs(x-xa[i+1]);
1204 if (dift < dif) {
1205 ns = i;
1206 dif = dift;
1207 }
1208 c[i] = ya[i+1];
1209 d[i] = ya[i+1];
1210 }
1211
1212 // Get initial approximation to y
1213 y = ya[ns+1];
1214 ns--;
1215
1216 // Loop over each column of the tableau
1217 for (int m = 1; m < n; ++m) {
1218
1219 // Update current c's and d's
1220 for (int i = 0; i < n-m; ++i) {
1221 double ho = xa[i+1] - x;
1222 double hp = xa[i+m+1] - x;
1223 double w = c[i+1] - d[i];
1224 double den = ho - hp;
1225 if (den == 0.0) {
1226 m_isvalid = false;
1227 m_message = "Invalid step size "+gammalib::str(den)+
1228 " encountered. Two values in xa array are"
1229 " identical.";
1230 if (!m_silent) {
1232 }
1233 }
1234 den = w/den;
1235 d[i] = hp*den;
1236 c[i] = ho*den;
1237 }
1238
1239 // Compute y correction
1240 *dy = (2*(ns+1) < (n-m)) ? c[ns+1] : d[ns--];
1241
1242 // Update y
1243 y += *dy;
1244
1245 } // endfor: looped over columns of tableau
1246
1247 // Return
1248 return y;
1249}
1250
1251
1252/***********************************************************************//**
1253 * @brief Auxiliary function for adaptive Simpson's method
1254 *
1255 * @param[in] a Left integration boundary.
1256 * @param[in] b Right integration boundary.
1257 * @param[in] eps Precision.
1258 * @param[in] S Integral of last computation.
1259 * @param[in] fa Function value at left integration boundary.
1260 * @param[in] fb Function value at right integration boundary.
1261 * @param[in] fc Function value at mid-point of interval [a,b]
1262 * @param[in] bottom Iteration counter (stop when 0)
1263 *
1264 * Implements a recursive auxiliary method for the adative_simpson()
1265 * integrator.
1266 ***************************************************************************/
1267double GIntegral::adaptive_simpson_aux(const double& a, const double& b,
1268 const double& eps, const double& S,
1269 const double& fa, const double& fb,
1270 const double& fc,
1271 const int& bottom) const
1272{
1273 // Store the iteration counter
1274 if (bottom < m_iter) {
1275 m_iter = bottom;
1276 }
1277
1278 // Compute mid-point c bet
1279 double c = 0.5*(a + b); //!< Mid-point of interval [a,b]
1280 double h = b - a; //!< Length of interval [a,b]
1281 double d = 0.5*(a + c); //!< Mid-point of interval [a,c]
1282 double e = 0.5*(c + b); //!< Mid-point of interval [c,b]
1283
1284 // Evaluate function at mid-points d and e
1285 double fd = m_kernel->eval(d);
1286 double fe = m_kernel->eval(e);
1287 m_calls += 2;
1288
1289 // Compute integral using Simpson's rule for the left and right interval
1290 double h12 = h / 12.0;
1291 double Sleft = h12 * (fa + 4.0*fd + fc);
1292 double Sright = h12 * (fc + 4.0*fe + fb);
1293 double S2 = Sleft + Sright;
1294
1295 // Allocate result
1296 double value;
1297
1298 // If converged then compute the result ...
1299 if (std::abs(S2 - S) <= 15.0 * eps) {
1300// if (std::abs(S2 - S) <= 15.0 * m_eps * std::abs(S2)) {
1301 value = S2 + (S2 - S)/15.0;
1302 }
1303
1304 // ... else if the maximum recursion depth was reached then compute the
1305 // result and signal result invalidity
1306 else if (bottom <= 0) {
1307 value = S2 + (S2 - S)/15.0;
1308 m_isvalid = false;
1309 }
1310
1311 // ... otherwise call this method recursively
1312 else {
1313 value = adaptive_simpson_aux(a, c, 0.5*eps, Sleft, fa, fc, fd, bottom-1) +
1314 adaptive_simpson_aux(c, b, 0.5*eps, Sright, fc, fb, fe, bottom-1);
1315 }
1316
1317 // Return result
1318 return value;
1319}
1320
1321
1322/***********************************************************************//**
1323 * @brief Adaptive Gauss-Lobatto-Kronrod integration helper
1324 *
1325 * @param[in] a Left integration boundary.
1326 * @param[in] b Right integration boundary.
1327 * @param[in] fa Function value at the left integration boundary.
1328 * @param[in] fb Function value at the right integration boundary.
1329 * @param[in] is 13-point Kronrod extension of previous step.
1330 * @param[in] toler Tolerance.
1331 * @return Integral of kernel over interval [a,b].
1332 *
1333 * Integrates the function using an adaptive Gauss-Lobatto method with a
1334 * Kronrod extension.
1335 ***************************************************************************/
1337 const double& b,
1338 const double& fa,
1339 const double& fb,
1340 const double& is,
1341 const double& toler) const
1342{
1343 // Set constants
1344 const double alpha = std::sqrt(2.0/3.0);
1345 const double beta = 1.0/std::sqrt(5.0);
1346
1347 // Initialise value
1348 double value = 0.0;
1349
1350 // Compute midpoint and step size
1351 double m = 0.5*(a+b);
1352 double h = 0.5*(b-a);
1353
1354 // Compute intermediate points
1355 double mll = m - alpha*h;
1356 double ml = m - beta*h;
1357 double mr = m + beta*h;
1358 double mrr = m + alpha*h;
1359
1360 // Evaluate function at intermediate points
1361 double fmll = m_kernel->eval(mll);
1362 double fml = m_kernel->eval(ml);
1363 double fm = m_kernel->eval(m);
1364 double fmr = m_kernel->eval(mr);
1365 double fmrr = m_kernel->eval(mrr);
1366 m_calls += 5;
1367
1368 // 4-point Gauss-Lobatto formula
1369 double i2 = h/6.0*(fa+fb+5.0*(fml+fmr));
1370
1371 // 7-point Konrod extension
1372 double i1 = h/1470.0*(77.0*(fa+fb)+432.0*(fmll+fmrr)+625.0*(fml+fmr)+672.0*fm);
1373
1374 // Check for convergence
1375 if (std::abs(i1-i2) <= toler*is || mll <= a || b <= mrr) {
1376 if ((mll <= a || b <= mrr) && m_terminate) {
1377 m_isvalid = false;
1378 m_terminate = false;
1379 }
1380 value = i1;
1381 }
1382
1383 // ... otherwise subdivide interval
1384 else {
1385 value = adaptive_gauss_kronrod_aux(a, mll, fa, fmll, is, toler) +
1386 adaptive_gauss_kronrod_aux(mll, ml, fmll, fml, is, toler) +
1387 adaptive_gauss_kronrod_aux(ml, m, fml, fm, is, toler) +
1388 adaptive_gauss_kronrod_aux(m, mr, fm, fmr, is, toler) +
1389 adaptive_gauss_kronrod_aux(mr, mrr, fmr, fmrr, is, toler) +
1390 adaptive_gauss_kronrod_aux(mrr, b, fmrr, fb, is, toler);
1391 }
1392
1393 // Return result
1394 return value;
1395}
1396
1397
1398/***********************************************************************//**
1399 * @brief Rescale errors for Gauss-Kronrod integration
1400 *
1401 * @param[in] err Error estimate.
1402 * @param[in] result_abs ???.
1403 * @param[in] result_asc ???.
1404 * @return Rescaled error estimate.
1405 ***************************************************************************/
1407 const double& result_abs,
1408 const double& result_asc) const
1409{
1410 // Take absolute value of error
1411 err = std::abs(err);
1412
1413 // ...
1414 if (result_asc != 0.0 && err != 0.0) {
1415 double scale = std::pow((200.0 * err / result_asc), 1.5);
1416 if (scale < 1.0) {
1417 err = result_asc * scale ;
1418 }
1419 else {
1420 err = result_asc ;
1421 }
1422 }
1423 if (result_abs > 2.2250738585072014e-308 / (50.0 * 2.2204460492503131e-16)) {
1424 double min_err = 50.0 * 2.2204460492503131e-16 * result_abs;
1425 if (min_err > err) {
1426 err = min_err ;
1427 }
1428 }
1429
1430 // Return error
1431 return err;
1432}
Exception handler interface definition.
Single parameter function abstract base class definition.
#define G_ROMBERG
Definition GIntegral.cpp:38
#define G_TRAPZD
Definition GIntegral.cpp:39
#define G_POLINT
Definition GIntegral.cpp:40
Integration class interface definition.
Gammalib tools definition.
GChatter
Definition GTypemaps.hpp:33
@ SILENT
Definition GTypemaps.hpp:34
double sum(const GVector &vector)
Computes vector sum.
Definition GVector.cpp:1012
Single parameter function abstract base class.
Definition GFunction.hpp:44
virtual double eval(const double &x)=0
GIntegral class interface definition.
Definition GIntegral.hpp:46
const bool & silent(void) const
Get silence flag.
const std::string & message(void) const
Return integration status message.
std::string m_message
Status message (if result is invalid)
int m_calls
Number of function calls used.
const double & eps(void) const
Get relative precision.
double adaptive_simpson_aux(const double &a, const double &b, const double &eps, const double &S, const double &fa, const double &fb, const double &fc, const int &bottom) const
Auxiliary function for adaptive Simpson's method.
const int & iter(void) const
Return number of iterations.
double adaptive_gauss_kronrod_aux(const double &a, const double &b, const double &fa, const double &fb, const double &is, const double &toler) const
Adaptive Gauss-Lobatto-Kronrod integration helper.
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
double m_relerr
Absolute integration error.
void copy_members(const GIntegral &integral)
Copy class members.
double m_eps
Requested relative integration precision.
const int & fixed_iter(void) const
Return fixed number of iterations.
const int & calls(void) const
Get number of function calls.
bool m_has_abserr
Has absolute integration error.
bool m_has_relerr
Has relative integration error.
GIntegral * clone(void) const
Clone integral.
double trapzd(const double &a, const double &b, const int &n=1, double result=0.0)
Perform Trapezoidal integration.
void clear(void)
Clear integral.
bool m_terminate
Signals termination of subdivision.
bool m_silent
Suppress integration warnings in console.
const GFunction * kernel(void) const
Get function kernel.
int m_iter
Number of iterations used.
double gauss_kronrod(const double &a, const double &b) const
Gauss-Kronrod integration.
virtual ~GIntegral(void)
Destructor.
int m_max_iter
Maximum number of iterations.
void init_members(void)
Initialise class members.
double adaptive_simpson(const double &a, const double &b) const
Adaptive Simpson's integration.
GFunction * m_kernel
Pointer to function kernel.
std::string print(const GChatter &chatter=NORMAL) const
Print integral information.
GIntegral(void)
Void constructor.
void free_members(void)
Delete class members.
GIntegral & operator=(const GIntegral &integral)
Assignment operator.
double adaptive_gauss_kronrod(const double &a, const double &b) const
Adaptive Gauss-Lobatto-Kronrod integration.
const int & max_iter(void) const
Return maximum number of iterations.
double rescale_error(double err, const double &result_abs, const double &result_asc) const
Rescale errors for Gauss-Kronrod integration.
const bool & is_valid(void) const
Signal if integration result is valid.
double m_abserr
Absolute integration error.
bool m_isvalid
Integration result valid (true=yes)
double polint(double *xa, double *ya, int n, double x, double *dy)
Perform Polynomial interpolation.
int m_fix_iter
Fixed number of iterations.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition GTools.cpp:1162
const double gkw43b[12]
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition GTools.cpp:508
const double gkw43a[10]
const double gkw87b[23]
const double gkx2[5]
Definition GIntegral.cpp:70
const double gkw21b[6]
Definition GIntegral.cpp:88
const double gkw10[5]
Definition GIntegral.cpp:61
const double gkw87a[21]
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition GTools.cpp:1405
const double gkx1[5]
Definition GIntegral.cpp:52
const double gkx3[11]
Definition GIntegral.cpp:98
const double gkw21a[5]
Definition GIntegral.cpp:79
const double gkx4[22]