function [Result, Output] = quadg (fun,a,b,opts,varargin)
%QUADG  Numerical quadrature using open methods.
%   Q = QUADG(fun,a,b) approximates the integral of fun
%   from a to b using some type of Gaussian quadrature.
%
%   The default type is a Gauss-Kronrod, 7(15) order.
%
%   Q = QUADG(fun,a,b,options) allows integrations option to
%   be specified.  
%
%   Q = QUADG(fun,a,b,options,p1,p2,...) will pass the extra 
%   parameters p1, p2, ... to fun.
%
%   Options is a structure specifying any of the following options:
%
%       Field name              Parameter           Possibilites and {Default}
%
%       options.tol             Tolerance           [{1e-6}, any fractional number]
%       options.diplay          Display on/off      [{'off'}, 'on']
%       options.nodes           nodes               [vector, {Gauss-Kronrod nodes}]
%       options.weights         weights             [matrix, {Gauss-Kronrod weights}]
%
%   Using options=[] will use the default options.
%
%   Tolerance is based on the relative error.  With double-precision  
%   arithmetic, feasible values for this parameter range 
%   from 1e-15 to 1.  The default is 1e-6.  
%
%   Setting the display on shows the integrand samples and a histogram
%   of their distribution, as well as a table reporting the progress of 
%   the integration as it proceeds.  Note that the display will slow down 
%   the performance of the quadrature.
%
%   The nodes and weights parameters are discussed in the 
%   "Advanced Features" section in the source code just after
%   the help entry.  Type 'edit quadg' to see this note.  There's 
%   no need to use this option unless you like playing with other 
%   people's source code and know what you're doing (and can 
%   figure out what I was doing).
%
%   Example:
%
%      f = inline('cos(x).*sin(x.*x)');
%      options.display = 'on';
%      quadg(f,-2*pi,2*pi,options)
%    
%      Increasing the range of this quadrature will cause most recursive 
%      to slow dramatically and run out of stack space.
%
%   See the source code for additional comments.  
%
%   See also QUADL, QUAD, DBLQUAD, INLINE, @.

%   ========================================================
%   ADVANCED FEATURES -- DEFINING YOUR OWN WEIGHTS AND NODES
%   --------------------------------------------------------
%   You can also provide your own nodes and weights to this
%   function, allowing you to use this code as a harness for
%   any integration method you like.  Your weights and nodes 
%   should be defined for an integration over [-1,+1]
%
%   This is an advanced feature.  Unless you're really interested
%   in using your own quadrature, there's no need to worry about this.
%
%   The QUADG function has two alternate algorithms which you can 
%   access in this way.  For example:
%
%       [options.nodes options.weights] = quadg('gausskronrod')
%
%   returns values for a Gauss-Kronrod 7(15) order quadrature,
%
%       [options.nodes options.weights] = quadg('gausslegendre',p)
%
%   returns values for a p-th order Gauss-Legendre quadrature.  Then:
%
%       quadg(f,a,b,options)
%
%   will integrate f using the algorithm specified.
%
%   options.nodes and options.weights should be defined as a one
%   column vector and a two row matrix, respectively, so that:
%
%       2*options.weights*feval(fun,a)
%
%   returns a 2 by 1 vector of the lower and higher order estimates 
%   for an integration over [-1,+1].  For example:
%
%       >> [opts.nodes opts.weights] = quadg('gausslegendre',2);
%
%       >> opts.nodes
%        ans =
%           -0.5774
%            0.5774
%           -0.7746
%                 0
%            0.7746
%
%       >> opts.weights
%        ans =
%            0.5000    0.5000         0         0         0
%                 0         0    0.2778    0.4444    0.2778
% 
%       >> 2*opts.weights*cos(opts.nodes)
%        ans =
%           1.6758
%           1.6830
%     
%   These are the 2-nd and 3-rd order estimates for the integral of cosine
%   taken over the range [-1,+1]

%   =============================================
%   MISCELLANEOUS COMMENTS
%   ---------------------------------------------
%   IMO, the iterative imlementation here is interesting to see; it's what 
%   makes this code fast.  
%
%   One way to test the accuracy of this code is to use the Symbolic Toolbox
%   to compute analytical results for some representative test integrands, 
%   and compare those to the results this code produces.
%
%   Pros of this code:
%       - fast and accurate (esp the default Gauss-Kronrod algorithm)
%       - nice graphical display
%       - can use your own nodes and weights 
%       - interesting iterative implementation (contrasts with the
%         typical recursive implementations often used)
%
%   Cons:
%       - though they works well in practice, the stopping criteria 
%         are simple (see QUADL for better ideas)
%
%   I have not compared this to the quadrature functions in MATLAB, but the more tools
%   you have in your quadrature collection the better off you'll be.  
%
%   If someone wants to clean up this code, please feel free to do so and
%   let me know; I'll gladly accept any useful modifications and incorporate them
%   into this file.  It's sloppy in some places, and there are other places where I'm
%   using KRON and SPDIAGS when I shouldn't.
%
%   Another more sophisticated enhancement would be to add code calculating 
%   the Kronrod extension nodes for any order quadrature.  I'm not sure if 
%   these extended nodes can be derived in a straightforward way, however. 
%
%   This function hasn't been extensively tested, and I haven't put in any "tweaks"
%   to account for bad inputs or any type of error checking.  
%
%   Please note that this file is *not* officially supported by The MathWorks. 
%
%   Nabeel Azar
%   nabeel@mathworks.com
%   nabeel@ieee.org
%
%   $Revision: 0.1 $  $Date: 1998 $

%______________________________________________________________
%  Input checking:

%   Are they calling it to just obtain nodes and weights?
%   If so, calculate them and then return
if isa(fun,'char')
    switch(lower(fun))
    case 'gausslegendre'
        if nargin>1; order = a; else; order = 7; end
        [Result Output] = gausslegendre(order);
        return
    case 'gausskronrod'
        [Result Output] = gausskronrod(7);
        return
    otherwise
        %   Do nothing and keep going
        msg = sprintf(['Assuming your first input is the name of\n',...
                'the function you want to integrate.  If you are\n',...
                'using MATLAB 6.0 or later, you should use\n',...
                'function handles to specify the function.']);
        warning(msg)
    end
end

%   If we're here, we're performing a quadrature.
%   We need at least 3 inputs
msg = nargchk(3,Inf,nargin);
if msg
    error(msg)
end

switch nargin
case 3; 
    %  Only 3 inputs, no options, no params
    Options = optionscheck([]); 
otherwise
    %   Options given, 
    Options = optionscheck(opts);
end
ExtraParameters = varargin;

%  End input checking
%--------------------------------------------------------------

%   Dynamically set MaxIter so no interval can 
%   be narrower than EPS.  The algorithm
%   typically stops well below MaxIter,
%   and throws a warning if it hits it.  
MaxIter = floor(log2(abs(b-a)/eps));
MaxIter = min(MaxIter,1024);

%   Set loop parameters and orient a and b as columns
Iter = 0;
Result = 0;
Tol = abs(Options.tol);
OutputOn = (nargout>=2);
a = a(:); b = b(:);

%  If verbosity is on, setup the figure window
if strcmp(Options.display,'on')
   %  Using verbose display also stored the information
   %  found in the second output argument.
   OutputOn = 1;
   %  Predefine storage for the data sets
   AllNodes = [];  AllIntegrand = [];  
   verbosity = 1;
   [Fig Axis Listbox HistogramAxis Spaces] = initializefigure(Options,fun,a,b);
   %  Start runtime counter
   t0 = cputime;
elseif OutputOn==1
   %  Predefine storage for the data sets
   AllNodes = [];  AllIntegrand = [];  
   verbosity = 0;
   t0 = cputime;
else
   verbosity = 0;
end


%______________________________________________________________
%  3)  Get the integration nodes and weights
%      for the first level and the second level
%      of checks.
X = Options.nodes;
W = Options.weights;
%--------------------------------------------------------------

%______________________________________________________________
%  4)  Begin the WHILE loop over which the integrand
%      is estimated.
%
%      This is where the real work takes place.  The code in this
%      loop isn't too difficult to follow.  Everything
%      up to the "if OutputOn" part is the actual quadrature;
%      the rest of the loop is just for options.
while ~isempty(a) & Iter<=MaxIter
   %  Increment iteration
   Iter = Iter + 1;

   %  Obtain nodes for the integration
   Nodes = scalex(X,a,b);
               
   %  Evaluate the integrand
   Integrand = feval(fun,Nodes,ExtraParameters{:});

   %  Scale integrand and compute the estimates.
   Estimates = W*Integrand*spdiags((b-a).',0,length(a),length(a));

   %  Find where the estimates are close and 
   %  remove the regions that had good estimates.
   %  This tolerance checking is simple, but works ok in practice.
   GoodRegions = abs((Estimates(1,:) - Estimates(2,:))./Estimates(2,:)) <=Tol;
   a(GoodRegions) = []; b(GoodRegions) = [];

   %  Update the integral to include the good regions.
   Result = Result + Estimates(2,:)*GoodRegions(:);

   %  Halve the remaining intervals
   %  Keep the a and b as row vectors.
   Middle_Terms = (a+b)/2;
   a = [a Middle_Terms].'; a = a(:).';
   b = [Middle_Terms b].'; b = b(:).';