function [ q, n ] = trap_adapt ( f, a, b, tol, h_min, n_max )

%% trap_adapt applies adaptive trapezoidal quadrature.
%
%  Discussion:
%
%    Use a sequence of trapezoid rules to estimate the integral of
%    f(x) from a to b, with an estimated absolute error of no more
%    than tol.
%
%  Input:
%
%    function f(x): the function to be integrated.
%
%    real a, b: the left and right endpoints.
%
%    real tol: the error tolerance.
%
%    real h_min: the smallest allowable stepsize.  The computation
%    fails if h falls below this limit.
%
%    integer n_max: the maximum number of subintervals.  The computation
%    fails if the number of subintervals exceeds this limit.
%
%  Output:
%
%    real q: the estimate for the integral.
%
%    integer n: the number of subintervals used.
%

%
%  Start with what seems like a small suggested step.
%
  h = ( b - a ) / 100.0;

  n = 0;
  q = 0.0;
  x0 = a;
%
%  Estimate the integral from x0 to x0+h.
%
  while ( x0 < b )
%
%  Using a stepsize of h, estimate the integral and the error.
%  If the error is small enough, accept the estimate, and advance x0.
%  Otherwise, decrease h and try again.
%
    while ( true )

      if ( n_max <= n )
        error ( 'trap_adapt - Fatal error!  N too large!' )
        return
      end
%
%  If the step would go past b, cut it back.
%
      if ( b < x0 + h )
        h = ( b - x0 );
        x1 = x0 + h / 2.0;
        x2 = b;
      else
        x1 = x0 + h / 2.0;
        x2 = x0 + h;
      end
%
%  Compute integral and error estimates using 1 and 2 trapezoids.
%
      q1 = h * (       f(x0)         +       f(x2) ) / 2.0;
      q2 = h * ( 0.5 * f(x0) + f(x1) + 0.5 * f(x2) ) / 2.0 ;
      e1 = abs ( 4.0 * ( q2 - q1 ) / 3.0 );
      e2 = abs (       ( q2 - q1 ) / 3.0 );
%
%  Decide if h can be increased, or is about right, or needs to be reduced.
%
      if ( 8.0 * e2 <= tol * h / ( b - a ) )
        h = h * 2.0;
        q = q + q2;
        x0 = x2;
        n = n + 1;
        fprintf ( 'At x0 = %g, try bigger h = %g\n', x0, h );
        break
      elseif ( e2 <= tol * h / ( b - a ) )
        h = h;
        q = q + q2;
        x0 = x2;
        n = n + 1;
        break;
      else
        h = h / 2.0;
        fprintf ( 'At x0 =  %g, try smaller h = %g\n', x0, h );
        if ( h < h_min )
          error ( 'h too small!' );
        end
      end

    end

  end

  return
end

