GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GIntegral.hpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GIntegral.hpp - Integration class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2014 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.hpp
23  * @brief Integration class interface definition
24  * @author Juergen Knoedlseder
25  */
26 
27 #ifndef GINTEGRAL_HPP
28 #define GINTEGRAL_HPP
29 
30 /* __ Includes ___________________________________________________________ */
31 #include <string>
32 #include <vector>
33 #include "GBase.hpp"
34 #include "GFunction.hpp"
35 
36 
37 /***********************************************************************//**
38  * @class GIntegral
39  *
40  * @brief GIntegral class interface definition.
41  *
42  * This class allows to perform integration using various methods. The
43  * integrand is implemented by a derived class of GFunction.
44  ***************************************************************************/
45 class GIntegral : public GBase {
46 
47 public:
48 
49  // Constructors and destructors
50  explicit GIntegral(void);
51  explicit GIntegral(GFunction* kernel);
52  GIntegral(const GIntegral& integral);
53  virtual ~GIntegral(void);
54 
55  // Operators
56  GIntegral& operator=(const GIntegral& integral);
57 
58  // Methods
59  void clear(void);
60  GIntegral* clone(void) const;
61  std::string classname(void) const;
62  void max_iter(const int& iter);
63  const int& max_iter(void) const;
64  void fixed_iter(const int& iter);
65  const int& fixed_iter(void) const;
66  void eps(const double& eps);
67  const double& eps(void) const;
68  void silent(const bool& silent);
69  const bool& silent(void) const;
70  const int& iter(void) const;
71  const int& calls(void) const;
72  const bool& is_valid(void) const;
73  const std::string& message(void) const;
74  void kernel(GFunction* kernel);
75  const GFunction* kernel(void) const;
76  double romberg(std::vector<double> bounds,
77  const int& order = 5);
78  double romberg(const double& a, const double& b,
79  const int& order = 5);
80  double trapzd(const double& a, const double& b,
81  const int& n = 1, double result = 0.0);
82  double adaptive_simpson(const double& a, const double& b) const;
83  double gauss_kronrod(const double& a, const double& b) const;
84  std::string print(const GChatter& chatter = NORMAL) const;
85 
86 protected:
87  // Protected methods
88  void init_members(void);
89  void copy_members(const GIntegral& integral);
90  void free_members(void);
91  double polint(double* xa, double* ya, int n, double x, double *dy);
92  double adaptive_simpson_aux(const double& a, const double& b,
93  const double& eps, const double& S,
94  const double& fa, const double& fb,
95  const double& fc,
96  const int& bottom) const;
97  double rescale_error(double err,
98  const double& result_abs,
99  const double& result_asc) const;
100 
101  // Protected data area
102  GFunction* m_kernel; //!< Pointer to function kernel
103  double m_eps; //!< Requested relative integration precision
104  int m_max_iter; //!< Maximum number of iterations
105  int m_fix_iter; //!< Fixed number of iterations
106  bool m_silent; //!< Suppress integration warnings in console
107 
108  // Integrator results
109  mutable int m_iter; //!< Number of iterations used
110  mutable int m_calls; //!< Number of function calls used
111  mutable bool m_isvalid; //!< Integration result valid (true=yes)
112  mutable bool m_has_abserr; //!< Has absolute integration error
113  mutable bool m_has_relerr; //!< Has relative integration error
114  mutable double m_abserr; //!< Absolute integration error
115  mutable double m_relerr; //!< Absolute integration error
116  mutable std::string m_message; //!< Status message (if result is invalid)
117 
118 };
119 
120 
121 /***********************************************************************//**
122  * @brief Return class name
123  *
124  * @return String containing the class name ("GIntegral").
125  ***************************************************************************/
126 inline
127 std::string GIntegral::classname(void) const
128 {
129  return ("GIntegral");
130 }
131 
132 
133 /***********************************************************************//**
134  * @brief Return number of iterations
135  *
136  * @return Number of iterations.
137  ***************************************************************************/
138 inline
139 const int& GIntegral::iter(void) const
140 {
141  return m_iter;
142 }
143 
144 
145 /***********************************************************************//**
146  * @brief Set maximum number of iterations
147  *
148  * @param[in] iter Maximum number of iterations.
149  ***************************************************************************/
150 inline
151 void GIntegral::max_iter(const int& iter)
152 {
153  m_max_iter = iter;
154  return;
155 }
156 
157 
158 /***********************************************************************//**
159  * @brief Return maximum number of iterations
160  *
161  * @return Maximum number of iterations.
162  ***************************************************************************/
163 inline
164 const int& GIntegral::max_iter(void) const
165 {
166  return m_max_iter;
167 }
168 
169 
170 /***********************************************************************//**
171  * @brief Set fixed number of iterations
172  *
173  * @param[in] iter Fixed number of iterations.
174  *
175  * If the fixed number of iterations is set, the integration algorithm will
176  * always performed the given number of iterations, irrespectively of the
177  * precision that is reached. This feature is relevant for computing
178  * numerical derivates from numerically integrated functions.
179  ***************************************************************************/
180 inline
181 void GIntegral::fixed_iter(const int& iter)
182 {
183  m_fix_iter = iter;
184  return;
185 }
186 
187 
188 /***********************************************************************//**
189  * @brief Return fixed number of iterations
190  *
191  * @return Fixed number of iterations.
192  ***************************************************************************/
193 inline
194 const int& GIntegral::fixed_iter(void) const
195 {
196  return m_fix_iter;
197 }
198 
199 
200 /***********************************************************************//**
201  * @brief Set relative precision
202  *
203  * @param[in] eps Relative precision.
204  ***************************************************************************/
205 inline
206 void GIntegral::eps(const double& eps)
207 {
208  m_eps = eps;
209  return;
210 }
211 
212 
213 /***********************************************************************//**
214  * @brief Get relative precision
215  *
216  * @return Relative precision.
217  ***************************************************************************/
218 inline
219 const double& GIntegral::eps(void) const
220 {
221  return m_eps;
222 }
223 
224 
225 /***********************************************************************//**
226  * @brief Get number of function calls
227  *
228  * @return Number of function calls.
229  ***************************************************************************/
230 inline
231 const int& GIntegral::calls(void) const
232 {
233  return m_calls;
234 }
235 
236 
237 /***********************************************************************//**
238  * @brief Set silence flag
239  *
240  * @param[in] silent Silence flag.
241  ***************************************************************************/
242 inline
243 void GIntegral::silent(const bool& silent)
244 {
245  m_silent = silent;
246  return;
247 }
248 
249 
250 /***********************************************************************//**
251  * @brief Get silence flag
252  *
253  * @return True is class is silent, false otherwise.
254  ***************************************************************************/
255 inline
256 const bool& GIntegral::silent(void) const
257 {
258  return m_silent;
259 }
260 
261 
262 /***********************************************************************//**
263  * @brief Set kernel
264  *
265  * @param[in] kernel Kernel.
266  *
267  * Sets the kernel for which the integral should be determined.
268  ***************************************************************************/
269 inline
271 {
272  m_kernel = kernel;
273  return;
274 }
275 
276 
277 /***********************************************************************//**
278  * @brief Get kernel
279  *
280  * @return Kernel.
281  ***************************************************************************/
282 inline
283 const GFunction* GIntegral::kernel(void) const
284 {
285  return m_kernel;
286 }
287 
288 
289 /***********************************************************************//**
290  * @brief Signal if integration result is valid
291  *
292  * @return True is integration result is valid.
293  ***************************************************************************/
294 inline
295 const bool& GIntegral::is_valid(void) const
296 {
297  return m_isvalid;
298 }
299 
300 
301 /***********************************************************************//**
302  * @brief Return integration status message
303  *
304  * @return Integration status message.
305  ***************************************************************************/
306 inline
307 const std::string& GIntegral::message(void) const
308 {
309  return m_message;
310 }
311 
312 #endif /* GINTEGRAL_HPP */
const double & eps(void) const
Get relative precision.
Definition: GIntegral.hpp:219
const int & fixed_iter(void) const
Return fixed number of iterations.
Definition: GIntegral.hpp:194
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
int m_iter
Number of iterations used.
Definition: GIntegral.hpp:109
const bool & silent(void) const
Get silence flag.
Definition: GIntegral.hpp:256
void init_members(void)
Initialise class members.
Definition: GIntegral.cpp:987
Definition of interface for all GammaLib classes.
double polint(double *xa, double *ya, int n, double x, double *dy)
Perform Polynomial interpolation.
Definition: GIntegral.cpp:1066
double rescale_error(double err, const double &result_abs, const double &result_asc) const
Rescale errors for Gauss-Kronrod integration.
Definition: GIntegral.cpp:1208
std::string classname(void) const
Return class name.
Definition: GIntegral.hpp:127
int m_max_iter
Maximum number of iterations.
Definition: GIntegral.hpp:104
const std::string & message(void) const
Return integration status message.
Definition: GIntegral.hpp:307
GIntegral class interface definition.
Definition: GIntegral.hpp:45
double gauss_kronrod(const double &a, const double &b) const
Gauss-Kronrod integration.
Definition: GIntegral.cpp:734
int m_fix_iter
Fixed number of iterations.
Definition: GIntegral.hpp:105
GIntegral & operator=(const GIntegral &integral)
Assignment operator.
Definition: GIntegral.cpp:304
void clear(void)
Clear integral.
Definition: GIntegral.cpp:334
const int & iter(void) const
Return number of iterations.
Definition: GIntegral.hpp:139
bool m_isvalid
Integration result valid (true=yes)
Definition: GIntegral.hpp:111
GIntegral(void)
Void constructor.
Definition: GIntegral.cpp:230
Single parameter function abstract base class definition.
Interface class for all GammaLib classes.
Definition: GBase.hpp:52
double trapzd(const double &a, const double &b, const int &n=1, double result=0.0)
Perform Trapezoidal integration.
Definition: GIntegral.cpp:564
bool m_has_abserr
Has absolute integration error.
Definition: GIntegral.hpp:112
std::string m_message
Status message (if result is invalid)
Definition: GIntegral.hpp:116
GChatter
Definition: GTypemaps.hpp:33
bool m_has_relerr
Has relative integration error.
Definition: GIntegral.hpp:113
const int & calls(void) const
Get number of function calls.
Definition: GIntegral.hpp:231
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&#39;s method.
Definition: GIntegral.cpp:1143
const GFunction * kernel(void) const
Get kernel.
Definition: GIntegral.hpp:283
double m_relerr
Absolute integration error.
Definition: GIntegral.hpp:115
void copy_members(const GIntegral &integral)
Copy class members.
Definition: GIntegral.cpp:1016
double m_eps
Requested relative integration precision.
Definition: GIntegral.hpp:103
bool m_silent
Suppress integration warnings in console.
Definition: GIntegral.hpp:106
std::string print(const GChatter &chatter=NORMAL) const
Print integral information.
Definition: GIntegral.cpp:918
virtual ~GIntegral(void)
Destructor.
Definition: GIntegral.cpp:282
Single parameter function abstract base class.
Definition: GFunction.hpp:44
const int & max_iter(void) const
Return maximum number of iterations.
Definition: GIntegral.hpp:164
GFunction * m_kernel
Pointer to function kernel.
Definition: GIntegral.hpp:102
int m_calls
Number of function calls used.
Definition: GIntegral.hpp:110
const bool & is_valid(void) const
Signal if integration result is valid.
Definition: GIntegral.hpp:295
double m_abserr
Absolute integration error.
Definition: GIntegral.hpp:114
void free_members(void)
Delete class members.
Definition: GIntegral.cpp:1043
GIntegral * clone(void) const
Clone integral.
Definition: GIntegral.cpp:352
double adaptive_simpson(const double &a, const double &b) const
Adaptive Simpson&#39;s integration.
Definition: GIntegral.cpp:676