function [answer,errest,flag] = Adapt(f,a,b,abserr,relerr)
%  Estimate the definite integral of f(x) from a to b using an
%  adaptive quadrature scheme based on Gauss-Kronrod (3,7) formulas.  
%  REQUIRES Add.m and Quad.m.
%
% Input parameters:
%    f     = name of the function defining f(x).
%   a, b   = end points of integration interval.
%   abserr = absolute error tolerance desired.
%   relerr = relative error tolerance desired.
%
% Output parameters:
%   answer = computed estimate of the integral.
%   errest = estimate of the absolute error in answer.
%   flag   = 0  for normal return;
%          = 1  insufficient storage in queue (120);
%          = 2  too many function evaluations (3577);

maxq = 120;  
queue = zeros(maxq,4);
maxfe = 3577;

%  Test input data.
if relerr < 10*eps
  error('relerr < 10*eps is not allowed in Adapt.')
end
if abserr < 0
  error('abserr <= 0 is not allowed in Adapt.')
end

%  Initialization.
length = 0;
top = 1;
bottom = 1;
nfe = 0;

%  Form an initial approximation answer to the integral over [a,b].
%  If it is not sufficiently accurate, initialize the queue and
%  begin the main loop.

[q,e,nfe] = Quad(f,a,b,nfe);
answer = q;
errest = e;
if abs(errest) > max(abserr,relerr * abs(answer))
  [queue,length,bottom] = Add(queue,q,e,a,b,length,bottom);
end

%  Main loop; if queue is empty then return, else subdivide the
%  top entry.
while length > 0

  %  Remove the top entries from the queue.
  q =     queue(top,1);
  e =     queue(top,2);
  alpha = queue(top,3);
  beta =  queue(top,4);
  length = length - 1;
  if top < maxq
    top = top + 1;
  else
    top = 1;
  end

  h = 0.5*(beta-alpha);
  [ql,el,nfe] = Quad(f,alpha,alpha+h,nfe);
  [qr,er,nfe] = Quad(f,alpha+h,beta,nfe);

  %  Update answer and the error estimate.
  answer = answer + ((ql + qr) - q);
  errest = errest + ((el + er) - e);

  %  Test for failures.
  if length >= maxq - 1
    flag = 1;
    return;
  end
  if nfe >= maxfe
    flag = 2;
    return;
  end

  %  Test for convergence.
  tol = max(abserr, relerr * abs(answer));
  if abs(errest) <= tol
    flag = 0;
    return;
  end

  %  Add new subintervals to queue if errors are too big.
  tol = tol * h / (b - a);
  if abs(el) > tol
    [queue,length,bottom] = Add(queue,ql,el,alpha,alpha+h,length,bottom);
  end
  if abs(er) > tol
    [queue,length,bottom] = Add(queue,qr,er,alpha+h,beta,length,bottom);
  end
end
flag = 0;