1 #include <iostream>
2 #include <cmath>
3 #include "vnl_adaptsimpson_integral.h"
4
int_fnct_(double * x)5 double vnl_adaptsimpson_integral::int_fnct_(double* x)
6 {
7 return pfnct_->f_(*x);
8 }
9
integral(vnl_integrant_fnct * f,double a,double b,double acury)10 double vnl_adaptsimpson_integral::integral(vnl_integrant_fnct* f, double a,
11 double b, double acury)
12 {
13 //set the function
14 pfnct_ = f;
15
16 return adaptivesimpson(&vnl_adaptsimpson_integral::int_fnct_, a, b, acury, 0, depth_);
17 }
18
adaptivesimpson(double (* f)(double *),double a,double b,double eps,int level,int level_max)19 double vnl_adaptsimpson_integral::adaptivesimpson(double(*f)(double*),
20 double a, double b, double eps,
21 int level, int level_max)
22 {
23 double c, d, e, h, result;
24 double one_simpson, two_simpson;
25 double left_simpson, right_simpson;
26
27 h = b-a;
28 c = 0.5*(a+b);
29 one_simpson = h*(f(&a)+4.0*f(&c)+f(&b))/6.0;
30 d = 0.5*(a+c);
31 e = 0.5*(c+b);
32 two_simpson = h*(f(&a)+4.0*f(&d)+2.0*f(&c)+4.0*f(&e)+f(&b))/12.0;
33 /* Check for level */
34 if (level+1 >= level_max) {
35 result = two_simpson;
36 std::cerr<< "Maximum level reached\n";
37 }
38 else {
39 /* Check for desired accuracy */
40 if (std::fabs(two_simpson-one_simpson) < 15.0*eps)
41 result = two_simpson + (two_simpson-one_simpson)/15.0;
42 /* Divide further */
43 else {
44 left_simpson = adaptivesimpson(f,a,c,eps/2.0,level+1,level_max);
45 right_simpson = adaptivesimpson(f,c,b,eps/2.0,level+1,level_max);
46 result = left_simpson + right_simpson;
47 }
48 }
49 return result;
50 }
51