/* aquad.c adaptive quadrature  numerical integration       */
/*         for a desired accuracy on irregular functions.   */
/* define the function to be integrated as:                 */
/*  double f(double x)                                      */
/*  {                                                       */
/*    // compute value                                      */
/*    return value;                                         */
/*  }                                                       */
/*                                                          */
/* integrate from xmin to xmax                              */
/* approximate desired absolute error in all intervals, eps */
/* accuracy absolutely not guaranteed                       */

#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))

double aquad(double f(double x), double xmin, double xmax, double eps)
{
  double area, temp, part, h;
  double err;
  int nmax = 2000;
  double sbin[2000];
  double ebin[2000];
  double abin[2000];
  int    fbin[2000];
  int i, j, k, n, nn, done;
  int kmax = 20;

  n=32; /* initial number of bins */
  h = (xmax-xmin)/(double)n;
  for(i=0; i<n; i++)
  {
    sbin[i] = xmin+i*h;
    ebin[i] = xmin+(i+1)*h;
    fbin[i] = 0;
  }
  k = 0;
  done = 0;
  nn = n; /* next available bin */
  while(!done)
  {
    done = 1;
    k++;
    if(k>=kmax) break; /* quit if more than kmax subdivisions */
    area = 0.0;
    for(i=0; i<n; i++)
    {
      if(fbin[i]==1) /* this interval finished */
      {
        area = area + abin[i]; /* accumulate total area each pass */
        continue;
      }
      temp = f((sbin[i]+ebin[i])/2.0); /* two integration methods */
      part = f((3.0*sbin[i]+ebin[i])/4.0);
      part = part + f((sbin[i]+3.0*ebin[i])/4.0);
      abin[i] = (part+2.0*temp)*(ebin[i]-sbin[i])/4.0;
      area = area + abin[i];
      err = abs(temp-part/2.0);
      if(err*1.414 < eps) /* heuristic  */
      {
        fbin[i] = 1; /* this interval finished */
      }
      else
      {
        done = 0; /* must keep dividing */
        if(nn>=nmax) /* out of space, quit */
        {
          done = 1;
          for(j=i+1; j<n; j++) area = area + abin[j];
          break; /* result not correct */
	}
        else /* divide interval into two halves */
        {
          sbin[nn] = (sbin[i]+ebin[i])/2.0;
          ebin[nn] = ebin[i];
          fbin[nn] = 0;
          ebin[i] = sbin[nn];         
          nn++;
	}
      }
    } /* end for i */
    n = nn;
  } /* end while */  
  return area;
} /* end aquad.c */