inte_cauchy_cern.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file inte_cauchy_cern.h
24  \brief File defining \ref o2scl::inte_cauchy_cern
25 */
26 #ifndef O2SCL_CERN_CAUCHY_H
27 #define O2SCL_CERN_CAUCHY_H
28 
29 #include <o2scl/inte.h>
30 #include <o2scl/inte_gauss_cern.h>
31 
32 #ifndef DOXYGEN_NO_O2NS
33 namespace o2scl {
34 #endif
35 
36  /** \brief Cauchy principal value integration (CERNLIB)
37 
38  The location of the singularity must be specified before-hand in
39  inte_cauchy_cern::s, and the singularity must not be at one of the
40  endpoints. Note that when integrating a function of the form \f$
41  \frac{f(x)}{(x-s)} \f$, the denominator \f$ (x-s) \f$ must be
42  specified in the argument \c func to integ(). This is different
43  from how the \ref inte_qawc_gsl operates.
44 
45  The method from \ref Longman58 is used for the decomposition of
46  the integral, and the resulting integrals are computed using
47  a user-specified base integration object.
48 
49  The uncertainty in the integral is not calculated, and is always
50  given as zero. The default base integration object is of type
51  \ref inte_gauss_cern. This is the CERNLIB default, but can be
52  modified by calling set_inte(). If the singularity is outside
53  the region of integration, then the result from the base
54  integration object is returned without calling the error
55  handler.
56 
57  Possible errors for integ() and integ_err():
58  - exc_einval - Singularity is on an endpoint
59  - exc_efailed - Could not reach requested accuracy
60  (occurs only if inte::err_nonconv is true)
61 
62  \note Currently \o2 supports only types \c double and
63  \c long \c double for the floating point type \c fp_t .
64 
65  This function is based on the CERNLIB routines RCAUCH and
66  DCAUCH which are documented at
67  http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d104/top.html
68  */
69  template<class func_t, class fp_t=double,
70  class weights_t=inte_gauss_coeffs_double>
71  class inte_cauchy_cern : public inte<func_t,fp_t> {
72 
73  public:
74 
75  const fp_t *w, *x;
76  weights_t wgts;
77 
79  it=&def_inte;
80  w=&(wgts.w[0]);
81  x=&(wgts.x[0]);
82  }
83 
84  /// The singularity (must be set before calling integ() or integ_err())
85  fp_t s;
86 
87  /** \brief Set the base integration object to use (default is \ref
88  inte_cauchy_cern::def_inte of type \ref inte_gauss_cern)
89  */
91  it=&i;
92  return 0;
93  }
94 
95  /** \brief Integrate function \c func from \c a to \c b
96  */
97  virtual int integ_err(func_t &func, fp_t a, fp_t b,
98  fp_t &res, fp_t &err) {
99  fp_t y1, y2, y3, y4;
100  size_t itx=0;
101 
102  err=0.0;
103 
104  if (s==a || s==b) {
105  O2SCL_CONV2_RET("Singularity on boundary in ",
106  "inte_cauchy_cern::integ_err().",
107  exc_einval,this->err_nonconv);
108  } else if ((s<a && s<b) || (s>a && s>b)) {
109  return it->integ_err(func,a,b,res,err);
110  }
111  fp_t h, b0;
112  if (2.0*s<a+b) {
113  h=it->integ(func,2.0*s-a,b);
114  b0=s-a;
115  } else {
116  h=it->integ(func,a,2.0*s-b);
117  b0=b-s;
118  }
119  fp_t c=0.005/b0;
120  fp_t bb=0.0;
121  bool loop1=true;
122  while(loop1==true) {
123  itx++;
124  fp_t s8, s16;
125  fp_t aa=bb;
126  bb=b0;
127  bool loop2=true;
128  while(loop2==true) {
129  fp_t c1=(bb+aa)/2.0;
130  fp_t c2=(bb-aa)/2.0;
131  fp_t c3=s+c1;
132  fp_t c4=s-c1;
133  s8=0.0;
134  s16=0.0;
135  for(int i=0;i<4;i++) {
136  fp_t u=c2*x[i];
137  y1=func(c3+u);
138  y2=func(c4-u);
139  y3=func(c3-u);
140  y4=func(c4+u);
141  s8+=w[i]*((y1+y2)+(y3+y4));
142  }
143  s8*=c2;
144  for(int i=4;i<12;i++) {
145  fp_t u=c2*x[i];
146  y1=func(c3+u);
147  y2=func(c4-u);
148  y3=func(c3-u);
149  y4=func(c4+u);
150  s16+=w[i]*((y1+y2)+(y3+y4));
151  }
152  s16*=c2;
153 
154  if (std::abs(s16-s8)<=this->tol_rel*(1.0+std::abs(s16))) {
155  loop2=false;
156  } else {
157  bb=c1;
158  if ((1.0+std::abs(c*c2))==1.0) {
159  loop2=false;
160  } else {
161  this->last_iter=itx;
162  O2SCL_CONV2_RET("Could not reach required accuracy in cern_",
163  "cauchy::integ()",exc_efailed,this->err_nonconv);
164  }
165  }
166  }
167  h+=s16;
168  if (bb==b0) loop1=false;
169  }
170  this->last_iter=itx;
171  res=h;
172  return o2scl::success;
173  }
174 
175  /// Default integration object
177 
178  protected:
179 
180 #ifndef DOXYGEN_INTERNAL
181 
182  /// The base integration object
184 
185 #endif
186 
187  };
188 
189 #ifndef DOXYGEN_NO_O2NS
190 }
191 #endif
192 
193 #endif
o2scl::inte_cauchy_cern
Cauchy principal value integration (CERNLIB)
Definition: inte_cauchy_cern.h:71
O2SCL_CONV2_RET
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:303
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::success
@ success
Success.
Definition: err_hnd.h:47
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::inte_cauchy_cern::s
fp_t s
The singularity (must be set before calling integ() or integ_err())
Definition: inte_cauchy_cern.h:85
o2scl::inte< func_t, double >::last_iter
size_t last_iter
The most recent number of iterations taken.
Definition: inte.h:70
o2scl::inte< func_t, double >::tol_rel
double tol_rel
The maximum relative uncertainty in the value of the integral (default )
Definition: inte.h:75
o2scl::inte
Base integration class [abstract base].
Definition: inte.h:51
o2scl::inte_cauchy_cern::it
inte< func_t, fp_t > * it
The base integration object.
Definition: inte_cauchy_cern.h:183
o2scl::inte_gauss_cern
Gaussian quadrature (CERNLIB)
Definition: inte_gauss_cern.h:293
o2scl::inte_cauchy_cern::def_inte
inte_gauss_cern< func_t, fp_t > def_inte
Default integration object.
Definition: inte_cauchy_cern.h:176
o2scl::inte< func_t, double >::err_nonconv
bool err_nonconv
If true, call the error handler if the routine does not converge or reach the desired tolerance (defa...
Definition: inte.h:88
o2scl::inte_cauchy_cern::set_inte
int set_inte(inte< func_t, fp_t > &i)
Set the base integration object to use (default is inte_cauchy_cern::def_inte of type inte_gauss_cern...
Definition: inte_cauchy_cern.h:90
o2scl::inte_cauchy_cern::integ_err
virtual int integ_err(func_t &func, fp_t a, fp_t b, fp_t &res, fp_t &err)
Integrate function func from a to b.
Definition: inte_cauchy_cern.h:97

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).