% Copyright (C) 1998 Walter Gautschi
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2, or (at your option) any
% later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
% FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
% for more details.
%
% You should have received a copy of the GNU General Public License
% along with Octave; see the file COPYING.  If not, write to the Free
% Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

%QUADL  Numerically evaluate integral using adaptive
%   Lobatto rule.
%
%   Q=QUADL('F',A,B) approximates the integral of
%   F(X) from A to B to machine precision.  'F' is a
%   string containing the name of the function. The
%   function F must return a vector of output values if
%   given a vector of input values.
%
%   Q=QUADL('F',A,B,TOL) integrates to a relative
%   error of TOL.
%
%   Q=QUADL('F',A,B,TOL,TRACE) displays the left
%   end point of the current interval, the interval
%   length, and the partial integral.
%
%   Q=QUADL('F',A,B,TOL,TRACE,P1,P2,...) allows
%   coefficients P1, ... to be passed directly to the
%   function F:  G=F(X,P1,P2,...). To use default values
%   for TOL or TRACE, one may pass the empty matrix ([]).
%
% W. Gander and W. Gautschi, Adaptive Quadrature - Revisited,
% BIT Vol. 40, No. 1, March 2000, pp. 84--101.
% http://www.inf.ethz.ch/personal/gander/

%   Author: Walter Gautschi
%   Date: 08/03/98
%   Reference: Gander, Computermathematik, Birkhaeuser, 1992.

% 2003-08-05 Shai Ayal
%   * permission from author to release as GPL
% 2004-02-10 Paul Kienzle
%   * renamed to quadl for compatibility
%   * replace global variable terminate2 with local function need_warning
%   * add paper ref to docs

function Q=quadl(f,a,b,tol,trace,varargin)
  need_warning(1);
  if(nargin<4), tol=[]; end;
  if(nargin<5), trace=[]; end;
  if(isempty(tol)), tol=eps; end;
  if(isempty(trace)), trace=0; end;
  if tol < eps
    tol = eps;
  end
  m=(a+b)/2; h=(b-a)/2;
  alpha=sqrt(2/3); beta=1/sqrt(5);
  x1=.942882415695480; x2=.641853342345781;
  x3=.236383199662150;
  x=[a,m-x1*h,m-alpha*h,m-x2*h,m-beta*h,m-x3*h,m,m+x3*h,...
    m+beta*h,m+x2*h,m+alpha*h,m+x1*h,b];
  y=feval(f,x,varargin{:});
  fa=y(1); fb=y(13);
  i2=(h/6)*(y(1)+y(13)+5*(y(5)+y(9)));
  i1=(h/1470)*(77*(y(1)+y(13))+432*(y(3)+y(11))+ ...
     625*(y(5)+y(9))+672*y(7));
  is=h*(.0158271919734802*(y(1)+y(13))+.0942738402188500 ...
     *(y(2)+y(12))+.155071987336585*(y(3)+y(11))+ ...
     .188821573960182*(y(4)+y(10))+.199773405226859 ...
     *(y(5)+y(9))+.224926465333340*(y(6)+y(8))...
     +.242611071901408*y(7));
  s=sign(is); if(s==0), s=1; end;
  erri1=abs(i1-is);
  erri2=abs(i2-is);
  R=1; if(erri2~=0), R=erri1/erri2; end;
  if(R>0 & R<1), tol=tol/R; end;
  is=s*abs(is)*tol/eps;
  if(is==0), is=b-a, end;
  Q=adaptlobstp(f,a,b,fa,fb,is,trace,varargin{:});



function Q=adaptlobstp(f,a,b,fa,fb,is,trace,varargin)
%ADAPTLOBSTP  Recursive function used by ADAPTLOB.
%
%   Q = ADAPTLOBSTP('F',A,B,FA,FB,IS,TRACE) tries to
%   approximate the integral of F(X) from A to B to
%   an appropriate relative error. The argument 'F' is
%   a string containing the name of f.  The remaining
%   arguments are generated by ADAPTLOB or by recursion.
%
%   See also ADAPTLOB.

%   Walter Gautschi, 08/03/98

  h=(b-a)/2; m=(a+b)/2;
  alpha=sqrt(2/3); beta=1/sqrt(5);
  mll=m-alpha*h; ml=m-beta*h; mr=m+beta*h; mrr=m+alpha*h;
  x=[mll,ml,m,mr,mrr];
  y=feval(f,x,varargin{:});
  fmll=y(1); fml=y(2); fm=y(3); fmr=y(4); fmrr=y(5);
  i2=(h/6)*(fa+fb+5*(fml+fmr));
  i1=(h/1470)*(77*(fa+fb)+432*(fmll+fmrr)+625*(fml+fmr) ...
     +672*fm);
  if(is+(i1-i2)==is) | (mll<=a) | (b<=mrr),
   if ((m <= a) | (b<=m)) & need_warning;
       warning(['Interval contains no more machine number. ',...
                'Required tolerance may not be met.']);
       need_warning(0);
    end;
    Q=i1;
    if(trace), disp([a b-a Q]), end;
  else
    Q=adaptlobstp(f,a,mll,fa,fmll,is,trace,varargin{:})+...
      adaptlobstp(f,mll,ml,fmll,fml,is,trace,varargin{:})+...
      adaptlobstp(f,ml,m,fml,fm,is,trace,varargin{:})+...
      adaptlobstp(f,m,mr,fm,fmr,is,trace,varargin{:})+...
      adaptlobstp(f,mr,mrr,fmr,fmrr,is,trace,varargin{:})+...
      adaptlobstp(f,mrr,b,fmrr,fb,is,trace,varargin{:});
  end

function r = need_warning(v)
  persistent w = [];
  if nargin == 0, r = w;
  else w=v; end