function [Q,cnt] = quad(funfcn,a,b,tol,trace,varargin)
%QUAD   Numerically evaluate integral, low order method.
%   Q = QUAD('F',A,B) approximates the integral of F(X) from A to B to
%   within a relative error of 1e-3 using an adaptive recursive
%   Simpson's rule.  'F' is a string containing the name of the
%   function.  Function F must return a vector of output values if given
%   a vector of input values.  Q = Inf is returned if an excessive
%   recursion level is reached, indicating a possibly singular integral.
%
%   Q = QUAD('F',A,B,TOL) integrates to a relative error of TOL.  Use
%   a two element tolerance, TOL = [rel_tol abs_tol], to specify a
%   combination of relative and absolute error.
%
%   Q = QUAD('F',A,B,TOL,TRACE) integrates to a relative error of TOL and
%   for non-zero TRACE traces the function evaluations with a point plot
%   of the integrand. 
%
%   Q = QUAD('F',A,B,TOL,TRACE,P1,P2,...) allows parameters P1, P2, ...
%   to be passed directly to function F:   G = F(X,P1,P2,...).
%   To use default values for TOL or TRACE, you may pass in the empty
%   matrix ([]).
%
%   See also QUAD8, DBLQUAD.

%   C.B. Moler, 3-22-87.
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.16 $  $Date: 1997/12/02 15:35:23 $

% [Q,cnt] = quad(F,a,b,tol) also returns a function evaluation count.

if nargin < 4, tol = [1.e-3 0]; trace = 0; end
if nargin < 5, trace = 0; end
c = (a + b)/2;
if isempty(tol), tol = [1.e-3 0]; end
if length(tol)==1, tol = [tol 0]; end
if isempty(trace), trace = 0; end

% Top level initialization
x = [a b c a:(b-a)/10:b];
y = feval(funfcn,x,varargin{:});
fa = y(1);
fb = y(2);
fc = y(3);
if trace
    lims = [min(x) max(x) min(y) max(y)];
    if ~isreal(y)
        lims(3) = min(min(real(y)),min(imag(y)));
        lims(4) = max(max(real(y)),max(imag(y)));
    end
    ind = find(~isfinite(lims));
    if ~isempty(ind)
        [mind,nind] = size(ind);
        lims(ind) = 1.e30*(-ones(mind,nind) .^ rem(ind,2));
    end
    axis(lims);
    plot([c b],[fc fb],'.')
    hold on
    if ~isreal(y)
        plot([c b],imag([fc fb]),'+')
    end
    drawnow
end
lev = 1;

% Adaptive, recursive Simpson's quadrature
if ~isreal(y)
    Q0 = 1e30;
else
    Q0 = inf;
end
recur_lev_excess = 0;  
[Q,cnt,recur_lev_excess] = ...
    quadstp(funfcn,a,b,tol,lev,fa,fc,fb,Q0,trace,recur_lev_excess,varargin{:});
cnt = cnt + 3;
if trace
    hold off
    axis('auto');
end
 if (recur_lev_excess > 1)
   warning(sprintf('Recursion level limit reached %d times.', ...
         recur_lev_excess ))
   end
   
%------------------------------------------------------------
function [Q,cnt,recur_lev_excess] = quadstp(FunFcn,a,b,tol,lev,fa,fc,fb,Q0,trace,recur_lev_excess,varargin)
%QUADSTP Recursive function used by QUAD.
%   [Q,cnt] = quadstp(F,a,b,tol,lev,fa,fc,fb,Q0) tries to
%   approximate the integral of f(x) from a to b to within a
%   relative error of tol.  F is a string containing the name
%   of f.  The remaining arguments are generated by quad or by
%   the recursion.  lev is the recursion level.
%   fa = f(a). fc = f((a+b)/2). fb = f(b).
%   Q0 is an approximate value of the integral.
%   See also QUAD and QUAD8.

%   C.B. Moler, 3-22-87.

LEVMAX = 10;
if lev > LEVMAX
      % Print warning the first time the recursion level is exceeded.
      if ~recur_lev_excess
        warning('Recursion level limit reached in quad. Singularity likely.')
        recur_lev_excess = 1;
      else
        recur_lev_excess = recur_lev_excess + 1;
      end
    Q = Q0;
    cnt = 0;
    c = (a + b)/2;
    if trace
        yc = feval(FunFcn,c,varargin{:});
        plot(c,yc,'.');    
        if ~isreal(yc)
            plot(c,imag(yc),'+');
        end
        drawnow
    end
else
    % Evaluate function at midpoints of left and right half intervals.
    h = b - a;
    c = (a + b)/2;
    x = [a+h/4, b-h/4];
    f = feval(FunFcn,x,varargin{:});
    if trace
        plot(x,f,'.');
        if ~isreal(f)
            plot(x,imag(f),'+')
        end
        drawnow
    end
    cnt = 2;
    
    % Simpson's rule for half intervals.
    Q1 = h*(fa + 4*f(1) + fc)/12;
    Q2 = h*(fc + 4*f(2) + fb)/12;
    Q = Q1 + Q2;
    
    % Recursively refine approximations.
    if abs(Q - Q0) > tol(1)*abs(Q) + tol(2);
      
        [Q1,cnt1,recur_lev_excess] = quadstp(FunFcn,a,c,tol/2,lev+1,fa,f(1),fc,Q1,trace,recur_lev_excess,varargin{:});
        [Q2,cnt2,recur_lev_excess] = quadstp(FunFcn,c,b,tol/2,lev+1,fc,f(2),fb,Q2,trace,recur_lev_excess,varargin{:});
        Q = Q1 + Q2;
        cnt = cnt + cnt1 + cnt2;
    end
end
