/* ----------------------------------- */
/* ADAPTIVE SIMPSON'S ALGORITHM

   Procedure for computing the definite
   integral of function f on a closed
   interval [a,b] using the adaptive
   Simpson procedure.                  */ 
/* ----------------------------------- */ 

double Simpson(double a, double b, double eps,
	       int level, int level_max)
{
  int i, j, k, kmax=1;
  double c, d, e, h, result;
  double one_simpson, two_simpson;
  double left_simpson, right_simpson;
 
  h = b-a;
  c = 0.5*(a+b);
  one_simpson = h*(f(a)+4.0*f(c)+f(b))/6.0;
  d = 0.5*(a+c);
  e = 0.5*(c+b);
  two_simpson = h*(f(a)+4.0*f(d)+2.0*f(c)+4.0*f(e)+f(b))/12.0;
  /* Check for level */
  if(level+1 >= level_max) {  
    result = two_simpson;
    printf("Maximum level reached\n"); 
  }
  else{
    /* Check for desired accuracy */
    if(fabs(two_simpson-one_simpson) < 15.0*eps)
      result = two_simpson + (two_simpson-one_simpson)/15.0;
    /* Divide further */
    else {
      left_simpson = Simpson(a,c,eps/2.0,level+1,level_max);   
      right_simpson = Simpson(c,b,eps/2.0,level+1,level_max);   
      result =  left_simpson + right_simpson;
    }
  }
  return(result);
}