function Q=adaptlob(f,a,b,tol,trace,varargin)
%ADAPTLOB  Numerically evaluate integral using adaptive
%   Lobatto rule.
%
%   Q=ADAPTLOB('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=ADAPTLOB('F',A,B,TOL) integrates to a relative
%   error of TOL.
%
%   Q=ADAPTLOB('F',A,B,TOL,TRACE) displays the left
%   end point of the current interval, the interval
%   length, and the partial integral.
%
%   Q=ADAPTLOB('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 ([]).
%
%   See also ADAPTLOBSTP.

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

global termination2

  termination2 = 0;
  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

global termination2

  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)) & (termination2==0);
       warning(['Interval contains no more machine number. ',...
                'Required tolerance may not be met.']);
       termination2 =1;
    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;