Hex  1.0
Hydrogen-electron collision solver
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
gausskronrod.h
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
2  * *
3  * / / / / __ \ \ / / *
4  * / /__ / / / _ \ \ \/ / *
5  * / ___ / | |/_/ / /\ \ *
6  * / / / / \_\ / / \ \ *
7  * *
8  * Jakub Benda (c) 2014 *
9  * Charles University in Prague *
10  * *
11 \* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
12 
13 #ifndef HEX_GAUSSKRONROD
14 #define HEX_GAUSSKRONROD
15 
16 #include <cmath>
17 
18 #include <gsl/gsl_integration.h>
19 
20 #include "misc.h"
21 #include "special.h"
22 
35 template <typename Functor> class GaussKronrod : public RadialFunction<double>
36 {
37  private:
38 
40  Functor Integrand;
41 
43  double Result;
44 
46  double AbsErr;
47 
49  bool Ok;
50 
52  double EpsAbs;
53 
55  double EpsRel;
56 
58  size_t Limit;
59 
61  gsl_integration_workspace * Workspace;
62 
64  std::string Status;
65 
66  public:
67 
68  // constructor
69  GaussKronrod (Functor f, size_t limit = 1000)
70  : Integrand(f), Result(special::constant::Nan), AbsErr(special::constant::Nan),
71  Ok(false), EpsAbs(0.), EpsRel(1e-5), Limit(limit),
72  Workspace(gsl_integration_workspace_alloc(Limit)) {}
73 
74  // destructor
76  {
77  gsl_integration_workspace_free(Workspace);
78  }
79 
80  // evaluate the integrand
81  inline double operator() (double x) const
82  {
83  return Integrand(x);
84  }
85 
92  static double eval (double x, void * ptr_radf)
93  {
94  RadialFunction<double> const & f = (* static_cast<RadialFunction<double> const *>(ptr_radf));
95  return f(x);
96  }
97 
121  bool integrate (double a, double b)
122  {
123  // reset
124  Result = AbsErr = special::constant::Nan;
125  Ok = true;
126 
127  // skip empty intervals
128  if (a == b)
129  {
130  Ok = true;
131  Result = 0;
132  Status = "";
133  return Ok;
134  }
135 
136  //
137  if (std::isnan(a) or std::isnan(b))
138  {
139  Ok = false;
140  Result = special::constant::Nan;
141  Status = "Some of the bounds is not finite.";
142  return false;
143  }
144 
145  // check order of bounds
146  if (a > b)
147  {
148  integrate(b, a);
149  Result = -Result;
150  return Ok;
151  }
152 
153  // setup integrand
154  gsl_function F;
155  F.function = &GaussKronrod::eval;
156  F.params = this;
157 
158  // setup integrator
159  int err = GSL_SUCCESS, err1 = GSL_SUCCESS, err2 = GSL_SUCCESS;
160  double Result1, Result2, AbsErr1, AbsErr2;
161 
162  // use correct integrator
163  if (std::isfinite(a) and std::isfinite(b)) /* -∞ < a < b < +∞ */
164  {
165  err = gsl_integration_qag
166  (
167  &F, a, b,
168  EpsAbs, EpsRel, Limit,
169  GSL_INTEG_GAUSS51,
170  Workspace,
171  &Result, &AbsErr
172  );
173  Status = (err != GSL_SUCCESS ? gsl_strerror(err) : "");
174  }
175  else if (std::isfinite(a) and not std::isfinite(b)) /* -∞ < a < b = +∞ */
176  {
177  err = gsl_integration_qagiu
178  (
179  &F, a,
180  EpsAbs, EpsRel, Limit,
181  Workspace,
182  &Result, &AbsErr
183  );
184  Status = (err != GSL_SUCCESS ? gsl_strerror(err) : "");
185  }
186  else if (not std::isfinite(a) and std::isfinite(b)) /* -∞ = a < b < +∞ */
187  {
188  err = gsl_integration_qagil
189  (
190  &F, b,
191  EpsAbs, EpsRel, Limit,
192  Workspace,
193  &Result, &AbsErr
194  );
195  Status = (err != GSL_SUCCESS ? gsl_strerror(err) : "");
196  }
197  else /* -∞ = a < b = +∞ */
198  {
199  err1 = gsl_integration_qagiu
200  (
201  &F, a,
202  EpsAbs, EpsRel, Limit,
203  Workspace,
204  &Result1, &AbsErr1
205  );
206 
207  err2 = gsl_integration_qagil
208  (
209  &F, b,
210  EpsAbs, EpsRel, Limit,
211  Workspace,
212  &Result2, &AbsErr2
213  );
214 
215  Result = Result1 + Result2;
216  AbsErr = AbsErr1 + AbsErr2;
217 
218  Status = (err1 != GSL_SUCCESS or err2 != GSL_SUCCESS ? std::string(gsl_strerror(err1)) + std::string(" & ") + std::string(gsl_strerror(err2)) : std::string());
219  }
220 
221  Ok = (err == GSL_SUCCESS and err1 == GSL_SUCCESS and err2 == GSL_SUCCESS);
222  return Ok;
223  }
224 
225  //
226  // getters & setters
227  //
228 
229  bool ok() const { return Ok; }
230  double result() const { return Result; }
231  double error() const { return AbsErr; }
232  std::string const & status() const { return Status; }
233 
234  void setEpsAbs (double epsabs) { EpsAbs = epsabs; }
235  void setEpsRel (double epsrel) { EpsRel = epsrel; }
236 
237  double epsabs () const { return EpsAbs; }
238  double epsrel () const { return EpsRel; }
239 };
240 
241 #endif
~GaussKronrod()
Definition: gausskronrod.h:75
double epsrel() const
Definition: gausskronrod.h:238
static double eval(double x, void *ptr_radf)
Evaluates radial function supplied in a pointer.
Definition: gausskronrod.h:92
double epsabs() const
Definition: gausskronrod.h:237
double F(double k, int l, double r, double sigma)
Definition: hydrogen.cpp:163
double operator()(double x) const
Evaluate the function.
Definition: gausskronrod.h:81
const double e
Definition: special.h:44
double result() const
Definition: gausskronrod.h:230
double error() const
Definition: gausskronrod.h:231
bool ok() const
Definition: gausskronrod.h:229
bool integrate(double a, double b)
Compute the integral.
Definition: gausskronrod.h:121
Numerical integrator.
Definition: gausskronrod.h:35
GaussKronrod(Functor f, size_t limit=1000)
Definition: gausskronrod.h:69
std::string const & status() const
Definition: gausskronrod.h:232
void setEpsAbs(double epsabs)
Definition: gausskronrod.h:234
Base class for radial functions.
Definition: special.h:536
const double Nan
Definition: special.h:68
void setEpsRel(double epsrel)
Definition: gausskronrod.h:235