function I = adaptiveQuadG(fcn,a,b,tol)
% Gander-Gautschi adaptive quadrature with robust
% machine independent termination criterion, &
% avoids arbitrary limits on depth of recursion.
% Using the 24- & 48-point Gauss-Legendre quadrature
% See: Algorithm 8.1 on p.358 Heath
% Usage: I = adaptiveQuadG('sin',0,pi/2,eps)
% K. Ming Leung, 03/15/03

I1 = gauss24pts(fcn, a, b);
I2 = gauss48pts(fcn, a, b);
m = a + (b-a)/2;
if m <= a | m >= b
    error('Tolerance may not be met');
    I = I2;
end
if abs(I2-I1) < tol % works better in Matlab
    I = I2;
else
    I = adaptiveQuadG(fcn,a,m,tol) + adaptiveQuadG(fcn,m,b,tol);
end

