function [Q,nfun]=da0glob(f,a,tolabs,tolrel,trace,varargin)
%DA0GLOB  Numerically evaluates an integral using adaptively
%   a Newton-Cotes 5/9 point rule. The routine is globally adaptive and the
%   fact that the routine makes use of two rules implies that it is doubly adaptive.
%   Except for some minor changes this routine is identical to the routine
%   COTEGLOB published by this author in 2002.
%
%   [Q,NFUN]=DA0GLOB('F',A) approximates the sum of the integrals of F(x) from
%   A(i) to A(i+1), i=1,...,length(A)-1, to a relative precision 10^-3. Note that length(A)
%   has to be at least 2. Q is the estimate of the integral and NFUN is the number 
%   of function values used by the code. '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,NFUN]=DA0GLOB('F',A,TOLABS) integrates to an absolute
%   error of TOLABS.
%
%   [Q,NFUN]=DA0GLOB('F',A,TOLABS,TOLREL) integrates to an absolute
%   error of max(TOLABS,EST*TOLREL), where EST is the current estimate of the
%   integral. The code also estimates the problems NOISE level
%   and it will not allow the effective absolute error to be below this
%   NOISE level: thus max([TOLABS,EST*TOLREL,NOISE])
%
%   [Q,NFUN]=DA0GLOB('F',A,TOLABS,TOLREL,TRACE) displays the left
%   and right end points of the current interval, the partial integral 
%   and the rule used (5 or 9).
%
%   [Q,NFUN]=DA0GLOB('F',A,TOLABS,TOLREL,TRACE,P1,P2,...) allows
%   coefficients P1, ... to be passed directly to the
%   function F:  G=F(X,P1,P2,...). To use default values
%   for TOLABS, TOLREL or TRACE, one may pass the empty matrix ([]).
%   The default values are: TOLABS=0, TOLREL=0.001 
%   (if both tolerances are empty matrices) and TRACE=0. 
%
%   The code also counts the number of function values and it may stop due
%   to the fact that the maximum number (10 000) has been reached.
%   Furtermore the code may stop processing an interval if the interval
%   size becomes too small giving a warning that this has happened. Finally
%   the code checks that the tolerance is not below the problem's noise
%   level. In such a case the stopping criterion is modified and a warning
%   given.

%   Modified by
%   Terje O. Espelid 15/11/03

%   References: Berntsen and Espelid, TOMS 1991
%               Espelid, in Espelid and Genz, 1992.
%               Espelid, Bit 2003
%               Espelid, Report 266 February 2004

global interval p2 w2 nw2

  if(nargin<3), tolabs=[];end
  if(nargin<4), tolrel=[]; end
  if(nargin<5), trace=[]; end
  if(isempty(tolrel)), tolrel=0; end
  if(isempty(tolabs)),tolabs=0;end
  if(tolabs==0)&(tolrel==0),tolrel=0.001;end
  if(isempty(trace)), trace=0; end
% 9 points Newton-Cotes rule in the interval [-1,1];
  p2=[-4,-3, -2,-1,0,1,2,3,4]/4;
% Weights:
  w2=[989,5888,-928,10496,-4540,10496,-928,5888,989]/14175;
% 8 Orthogonal Nullrules of degree  7,6,5,4,...,0 are given. All odd numbered are symmetric.
% All even numbered are antisymmetric. We given first the five elements.
nulw=[0.1101854769234527D-01,-0.8814838153876216D-01,0.3085193353856676D+00,-0.6170386707713351D+00, 0.7712983384641688D+00;
      0.4267465171184540D-01,-0.2560479102710724D+00,0.5974451239658356D+00,-0.5974451239658356D+00,0;
      0.1123675793894425D+00,-0.4775622124051307D+00,0.6180216866419335D+00,0.2809189484736109D-01,-0.5618378969472148D+00;
      0.2311270062743306D+00,-0.6355992672544093D+00 ,0.2311270062743311D+00,0.5200357641172426D+00,0;
      0.3911196434508170D+00,-0.5866794651762256D+00,-0.3073082912827879D+00,0.2514340565041003D+00,0.5028681130081833D+00;
      0.5561911416025476D+00,-0.2780955708012750D+00,-0.5164632029166498D+00,-0.3575514481730693D+00,0;
      0.6647755647017197D+00,0.1661938911754310D+00,-0.1899358756290700D+00,-0.4036137357117497D+00,-0.4748396890726778D+00;
      0.6455025992424890D+00,0.4841269494318643D+00,0.3227512996212468D+00,0.1613756498106190D+00,0];
% Then we use the symmetry
  nw2=[nulw(1:8,1:5),nulw(1:8,4:-1:1)];
% and finally the antisymmetry
  nw2(2:2:8,6:9)=-nw2(2:2:8,6:9);
% Thus: p2 is the points, w2 is the weights and nw2 is the nullrules associated with the 9 points Newton-Cotes rule.
% The coordinates the four new points in a five point approximation
  Z=[-3,-1,1,3]/4;
% Define the number of intervals given
  nint=length(a)-1;
  if nint < 1, warning('No interval is given.');Q=0;nfun=0;return,end
% Compute half of the interval lengths and centers of each interval
  h=(a(2:nint+1)-a(1:nint))/2;c=(a(2:nint+1)+a(1:nint))/2;
% Initialize: the function count, the total error, too small intervals
% (logical: interval) and the maximum number of function values allowed.
  nfun=0;total=0;interval=0;nmax=10000;singular=0;
  
  for i=1:nint  
% Compute y and start the count of the number of function values
     x=c(i)+h(i)*p2; y=feval(f,x,varargin{:});nfun=nfun+9;
% Check if some of the function values are infinite. In case thre are: put
% the value(s) to zero.
     IN=find(isinf(y)); y(IN)=0;singular=max(1-isempty(IN),singular);
% Apply the degree 9 Cotes rule on the interval and
% initialize the datastructure by storing: error estimate, integral, absolute value, 
% function values, left and right point of the interval and
% signalling by a flag that a nine point rule has been used.
     [error(i),is(i),isa(i)]=cotes9(a(i),a(i+1),y,32);
     fu(i,1:9)=y;lp(i)=a(i);rp(i)=a(i+1);flag(i)=1;
  end
% Compute global values
  Q=sum(is);isabs=sum(isa);total=sum(error);noise=50*eps*isabs;

% global error test and check the number of function evaluations.

  while (total> max([abs(Q)*tolrel,tolabs,noise]))&(nfun<nmax-3)

%  Find the interval the greatest error estimate.
      [err,i]=max(error);
%  Prepare to process this interval: either by subdivision or by using a 9
%  point rule
      aa=lp(i);bb=rp(i);cc=(aa+bb)/2;total=total-err;Q=Q-is(i);isabs=isabs-isa(i);
      if flag(i)==1
%  On this interval a 9 point rule has been used. Subdivide it and use all nine function values
%  to compute: first the left part  of the integral
         z(1:5)=fu(i,1:5);[err,inte,inteabs]=cotes5(aa,cc,z(1:5));zz(1:5)=fu(i,5:9);
         error(i)=err;is(i)=inte;isa(i)=inteabs;total=total+err;Q=Q+inte;isabs=isabs+inteabs;
         flag(i)=0;fu(i,1:2:9)=z(1:5);fu(i,2:2:8)=0;lp(i)=aa;rp(i)=cc;
%  and then the right part of the integral.
         nint=nint+1;i=nint;[err,inte,inteabs]=cotes5(cc,bb,zz(1:5));
         error(i)=err;is(i)=inte;isa(i)=inteabs;total=total+err;Q=Q+inte;isabs=isabs+inteabs;
         flag(i)=0;fu(i,1:2:9)=zz(1:5);fu(i,2:2:8)=0;lp(i)=cc;rp(i)=bb;
      else
%  On this interval a 5 point rule has been used only. Compute f in 4 new
%  points and apply the 9 point rule
         hh=(bb-aa)/2;x4=cc+hh*Z;y4=feval(f,x4,varargin{:});nfun=nfun+4;
         IN=find(isinf(y)); y(IN)=0;singular=max(1-isempty(IN),singular);
         fu(i,2:2:8)=y4(1:4);y(1:9)=fu(i,1:9);[err,inte,inteabs]=cotes9(aa,bb,y(1:9),32);
         error(i)=err;is(i)=inte;isa(i)=inteabs;flag(i)=1;
         total=total+err;Q=Q+inte;isabs=isabs+inteabs;
      end
% Redefine the global noise level
      noise=50*eps*isabs;
  end
  if (nfun>=nmax-3) 
       warning(['Stopping: maximum number of f-evaluations. ',...
                'Required tolerance may not be met.']);
  end
  if (interval==1)
         warning(['Interval too small. ',...
                'Required tolerance may not be met.']);
  end
  if max(abs(Q)*tolrel,tolabs) < noise
       warning(['Stopping: the tolerance is below the noise level of the problem. ',...
                'Required tolerance may not be met.']);
  end
  if singular==1
      warning(['Singularity probably detected. Required tolerance may not be met.'])
  end
  if trace
     for i=1:nint
         disp([lp(i),rp(i),is(i),5+flag(i)*4])
     end
  end
  Q=sum(is);
  
function [err,Q,Qabs]=cotes9(a,b,y,D)
%
%   [err,Q,Qabs] = COTES9(A,B,Y,D) applies a 9 point Newton-Cotes rule to 
%   approximate the integral of F(X) from A to B.
%   The argument Y contains the 9 function values to be used in the
%   approximation and D is a constant to be used in the error estimate.
%
%   See also COTES5
%   Modified
%   T O Espelid, 11/11/2003

global interval p2 w2 nw2
  h=(b-a)/2;c=(b+a)/2;x=c+h*p2;Q=(h*w2)*y';
  Qabs=(h*abs(w2))*abs(y');
% define the local noise level
  noise=50*Qabs*eps;
%  Compute the error estimates
  e=(h*nw2)*y';e2=e.^2;E2= (e2(1:2:7)+e2(2:2:8));E=sqrt(E2);
  Emin=min(E(2:4));
  if Emin==0,
     rmax=2;
  else
     r=E(1:3)./E(2:4);
     if sum(isinf(r))>0
        rmax=2;
     else
        rmax=max(r);
     end
  end
  if rmax > 1
     err=D*max(E(1:4));
  elseif 0.25 < rmax
         err=D*rmax*E(1);
      else
         err=D*4*rmax*rmax*E(1);
  end
  
% If the highest degree null rules are on the local noise level, then put the error to zero.
  if (E(1)<noise) & (E(2) < noise), err=0; end;
%check if interval has become too small: unable to distingwish between the
%endpoints and the two points close to the endpoints.
  if ((x(2) <= a) | (b<=x(8)))
       interval =1;
%In order to avoid handling this interval again: put the error to zero.
       err=0;
   end

function [err,Q,Qabs]=cotes5(a,b,y)
%
%
%   [err,Q,Qabs] = COTES5(A,B,Y) applies a 5 point Newton-Cotes rule to 
%   approximate the integral of F(X) from A to B. 
%   The argument Y contains the 5 function values to be used in the
%   approximation.
%
%   See also COTES9
%   Modified
%   T O Espelid, 11/11/2003

% The 5 points Cotes rule in the interval [-1,1];
  p1=[-1,-.5, 0,.5,1];
% Cotes weights: it should sum to 2.
  w1=[7,32,12,32,7]/45;
% Nullrules: 4 nullrules of degrees 3,2,1 and 0 for the 5 point Cotes rule.
  nw1=[0.1271031188518578,-0.5084124754074312,0.7626187131111468,-0.5084124754074312,0.1271031188518578;
         0.3362832433427012,-0.6725664866854024,0,0.6725664866854024,-0.3362832433427012;
         0.5684224278099782,-0.2842112139049893,-0.5684224278099774,-0.2842112139049893,0.5684224278099782;
         0.6725664866854023,0.3362832433427017,0,-0.3362832433427017,-0.6725664866854023];
  h=(b-a)/2;c=(b+a)/2;x=c+h*p1;Q=(h*w1)*y';Qabs=(h*w1)*abs(y');
% define the local noise level 
  noise=50*Qabs*eps;E=abs((h*nw1)*y');
  Emin=min(E(2:4));
  if Emin==0,
     rmax=2;
  else
     r=E(1:3)./E(2:4);
     if sum(isinf(r))>0
        rmax=2;
     else
        rmax=max(r);
     end
  end
  C=32;
  if rmax > 1
     err=C*max(E(1:4));
  elseif 0.5 < rmax
         err=C*rmax*E(2);
      else
         err=8*C*(rmax^4)*E(2);
  end
% If the highest degree null rules are on the local noise level, then put the error to zero.
  if (E(1)<noise) & (E(2) < noise), err=0; end;
