Go to the documentation of this file.00001 #include "vnl_adaptsimpson_integral.h"
00002 #include <vcl_iostream.h>
00003 #include <vcl_cmath.h>
00004
00005 double vnl_adaptsimpson_integral::int_fnct_(double* x)
00006 {
00007 return pfnct_->f_(*x);
00008 }
00009
00010 double vnl_adaptsimpson_integral::integral(vnl_integrant_fnct* f, double a,
00011 double b, double acury)
00012 {
00013
00014 pfnct_ = f;
00015
00016 return adaptivesimpson(&vnl_adaptsimpson_integral::int_fnct_, a, b, acury, 0, depth_);
00017 }
00018
00019 double vnl_adaptsimpson_integral::adaptivesimpson(double(*f)(double*),
00020 double a, double b, double eps,
00021 int level, int level_max)
00022 {
00023 double c, d, e, h, result;
00024 double one_simpson, two_simpson;
00025 double left_simpson, right_simpson;
00026
00027 h = b-a;
00028 c = 0.5*(a+b);
00029 one_simpson = h*(f(&a)+4.0*f(&c)+f(&b))/6.0;
00030 d = 0.5*(a+c);
00031 e = 0.5*(c+b);
00032 two_simpson = h*(f(&a)+4.0*f(&d)+2.0*f(&c)+4.0*f(&e)+f(&b))/12.0;
00033
00034 if (level+1 >= level_max) {
00035 result = two_simpson;
00036 vcl_cerr<< "Maximum level reached\n";
00037 }
00038 else {
00039
00040 if (vcl_fabs(two_simpson-one_simpson) < 15.0*eps)
00041 result = two_simpson + (two_simpson-one_simpson)/15.0;
00042
00043 else {
00044 left_simpson = adaptivesimpson(f,a,c,eps/2.0,level+1,level_max);
00045 right_simpson = adaptivesimpson(f,c,b,eps/2.0,level+1,level_max);
00046 result = left_simpson + right_simpson;
00047 }
00048 }
00049 return result;
00050 }
00051