[fun,a,b,err,ansmlv,ierr]=gaus8(fun,a,b,err,ansmlv,ierr);

function [fun,a,b,err,ansmlv,ierr]=gaus8(fun,a,b,err,ansmlv,ierr);
%***BEGIN PROLOGUE  GAUS8
%***PURPOSE  Integrate a real function of one variable over a finite
%            interval using an adaptive 8-point Legendre-Gauss
%            algorithm.  Intended primarily for high accuracy
%            integration or integration of smooth functions.
%***LIBRARY   SLATEC
%***CATEGORY  H2A1A1
%***TYPE      SINGLE PRECISION (GAUS8-S, DGAUS8-D)
%***KEYWORDS  ADAPTIVE QUADRATURE, AUTOMATIC INTEGRATOR,
%             GAUSS QUADRATURE, NUMERICAL INTEGRATION
%***AUTHOR  Jones, R. E., (SNLA)
%***DESCRIPTION
%
%     Abstract
%        GAUS8 integrates real functions of one variable over finite
%        intervals using an adaptive 8-point Legendre-Gauss algorithm.
%        GAUS8 is intended primarily for high accuracy integration
%        or integration of smooth functions.
%
%     Description of Arguments
%
%        Input--
%        FUN - name of external function to be integrated.  This name
%              must be in an EXTERNAL statement in the calling program.
%              FUN must be a REAL function of one REAL argument.  The
%              value of the argument to FUN is the variable of
%              integration which ranges from A to B.
%        A   - lower limit of integration
%        B   - upper limit of integration (may be less than A)
%        ERR - is a requested pseudorelative error tolerance.  Normally
%              pick a value of ABS(ERR) so that STOL .LT. ABS(ERR) .LE.
%              1.0E-3 where STOL is the single precision unit roundoff
%              R1MACH(4).  ANS will normally have no more error than
%              ABS(ERR) times the integral of the absolute value of
%              FUN(X).  Usually, smaller values for ERR yield more
%              accuracy and require more function evaluations.
%
%              A negative value for ERR causes an estimate of the
%              absolute error in ANS to be returned in ERR.  Note that
%              ERR must be a variable (not a constant) in this case.
%              Note also that the user must reset the value of ERR
%              before making any more calls that use the variable ERR.
%
%        Output--
%        ERR - will be an estimate of the absolute error in ANS if the
%              input value of ERR was negative.  (ERR is unchanged if
%              the input value of ERR was non-negative.)  The estimated
%              error is solely for information to the user and should
%              not be used as a correction to the computed integral.
%        ANS - computed value of integral
%        IERR- a status code
%            --Normal codes
%               1 ANS most likely meets requested error tolerance,
%                 or A=B.
%              -1 A and B are too nearly equal to allow normal
%                 integration.  ANS is set to zero.
%            --Abnormal code
%               2 ANS probably does not meet requested error tolerance.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  I1MACH, R1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   810223  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%***end PROLOGUE  GAUS8
persistent aa ae anib area c ce ee ef eps est firstCall gl glr gr gt h hh k kml kmx l lmn lmx lr mxl nbits nib nlmn nlmx sq2 tol vl vr w1 w2 w3 w4 x x1 x2 x3 x4 ; if isempty(firstCall),firstCall=1;end; 

if isempty(k), k=0; end;
if isempty(kml), kml=0; end;
if isempty(kmx), kmx=0; end;
if isempty(l), l=0; end;
if isempty(lmn), lmn=0; end;
if isempty(lmx), lmx=0; end;
if isempty(lr), lr=zeros(1,30); end;
if isempty(mxl), mxl=0; end;
if isempty(nbits), nbits=0; end;
if isempty(nib), nib=0; end;
if isempty(nlmn), nlmn=0; end;
if isempty(nlmx), nlmx=0; end;
if isempty(gt), gt=0; end;
if isempty(aa), aa=zeros(1,30); end;
if isempty(ae), ae=0; end;
if isempty(anib), anib=0; end;
if isempty(area), area=0; end;
if isempty(c), c=0; end;
if isempty(ce), ce=0; end;
if isempty(ee), ee=0; end;
if isempty(ef), ef=0; end;
if isempty(eps), eps=0; end;
if isempty(est), est=0; end;
if isempty(gl), gl=0; end;
if isempty(glr), glr=0; end;
if isempty(gr), gr=zeros(1,30); end;
if isempty(hh), hh=zeros(1,30); end;
if isempty(sq2), sq2=0; end;
if isempty(tol), tol=0; end;
if isempty(vl), vl=zeros(1,30); end;
if isempty(vr), vr=0; end;
if isempty(w1), w1=0; end;
if isempty(w2), w2=0; end;
if isempty(w3), w3=0; end;
if isempty(w4), w4=0; end;
if isempty(x1), x1=0; end;
if isempty(x2), x2=0; end;
if isempty(x3), x3=0; end;
if isempty(x4), x4=0; end;
if isempty(x), x=0; end;
if isempty(h), h=0; end;
% g8= @(x,h)  h*((w1*(fun(x-x1*h)+fun(x+x1*h))+w2*(fun(x-x2*h)+fun(x+x2*h)))+(w3*(fun(x-x3*h)+fun(x+x3*h))+w4*(fun(x-x4*h)+fun(x+x4*h))));real :: g8 ;
if firstCall,   x1 =[1.83434642495649805e-01];  end;
if firstCall,  x2 =[5.25532409916328986e-01];  end;
if firstCall,  x3 =[7.96666477413626740e-01];  end;
if firstCall,  x4=[9.60289856497536232e-01];  end;
if firstCall,   w1 =[3.62683783378361983e-01];  end;
if firstCall,  w2 =[3.13706645877887287e-01];  end;
if firstCall,  w3 =[2.22381034453374471e-01];  end;
if firstCall,  w4=[1.01228536290376259e-01];  end;
if firstCall,   sq2=[1.41421356e0];  end;
if firstCall,   nlmn=[1];  end;
if firstCall,  kmx=[5000];  end;
if firstCall,  kml=[6];  end;
firstCall=0;
g8= @(x,h)  h.*((w1.*(fun(x-x1.*h)+fun(x+x1.*h))+w2.*(fun(x-x2.*h)+fun(x+x2.*h)))+(w3.*(fun(x-x3.*h)+fun(x+x3.*h))+w4.*(fun(x-x4.*h)+fun(x+x4.*h))));
%***FIRST EXECUTABLE STATEMENT  GAUS8
%
%     Initialize
%
[k ]=i1mach(11);
anib = r1mach(5).*k./0.30102000e0;
nbits = fix(anib);
nlmx = fix(min(30,fix((nbits.*5)./8)));
ansmlv = 0.0e0;
ierr = 1;
ce = 0.0e0;
if( a==b )
if( err<0.0e0 )
err = ce;
end;
else;
lmx = fix(nlmx);
lmn = fix(nlmn);
if( b~=0.0e0 )
if( (abs(1.0e0).*sign(b)).*a>0.0e0 )
c = abs(1.0e0-a./b);
if( c<=0.1e0 )
if( c<=0.0e0 )
if( err<0.0e0 )
err = ce;
end;
return;
else;
anib = 0.5e0 - log(c)./0.69314718e0;
nib = fix(anib);
lmx = fix(min(nlmx,nbits-nib-7));
if( lmx<1 )
ierr = -1;
xermsg('SLATEC','GAUS8',['A and B are too nearly equal to allow normal integration. $$','ANS is set to zero and IERR to -1.'],1,-1);
if( err<0.0e0 )
err = ce;
end;
return;
else;
lmn = fix(min(lmn,lmx));
end;
end;
end;
end;
end;
tol = max(abs(err),2.0e0.^(5-nbits))./2.0e0;
if( err==0.0e0 )
tol = sqrt(r1mach(4));
end;
eps = tol;
hh(1) =(b-a)./4.0e0;
aa(1) = a;
lr(1) = 1;
l = 1;
g8= @(x,h)  h.*((w1.*(fun(x-x1.*h)+fun(x+x1.*h))+w2.*(fun(x-x2.*h)+fun(x+x2.*h)))+(w3.*(fun(x-x3.*h)+fun(x+x3.*h))+w4.*(fun(x-x4.*h)+fun(x+x4.*h))));
est = g8(aa(l)+2.0e0.*hh(l),2.0e0.*hh(l));
k = 8;
area = abs(est);
ef = 0.5e0;
mxl = 0;
%
%     Compute refined estimates, estimate the error, etc.
%
while( true );
gt=0;
g8= @(x,h)  h.*((w1.*(fun(x-x1.*h)+fun(x+x1.*h))+w2.*(fun(x-x2.*h)+fun(x+x2.*h)))+(w3.*(fun(x-x3.*h)+fun(x+x3.*h))+w4.*(fun(x-x4.*h)+fun(x+x4.*h))));
gl = g8(aa(l)+hh(l),hh(l));
g8= @(x,h)  h.*((w1.*(fun(x-x1.*h)+fun(x+x1.*h))+w2.*(fun(x-x2.*h)+fun(x+x2.*h)))+(w3.*(fun(x-x3.*h)+fun(x+x3.*h))+w4.*(fun(x-x4.*h)+fun(x+x4.*h))));
gr(l) = g8(aa(l)+3.0e0.*hh(l),hh(l));
k = fix(k + 16);
area = area +(abs(gl)+abs(gr(l))-abs(est));
%     IF (L .LT. LMN) GO TO 11
glr = gl + gr(l);
ee = abs(est-glr).*ef;
ae = max(eps.*area,tol.*abs(glr));
if( ee>ae )
%
%     Consider the left half of this level
%
if( k>kmx )
lmx = fix(kml);
end;
if( l>=lmx )
mxl = 1;
else;
l = fix(l + 1);
eps = eps.*0.5e0;
ef = ef./sq2;
hh(l) = hh(l-1).*0.5e0;
lr(l) = -1;
aa(l) = aa(l-1);
est = gl;
continue;
end;
end;
ce = ce +(est-glr);
if( lr(l)<=0 )
%
%     Proceed to right half at this level
%
vl(l) = glr;
else;
%
%     Return one level
%
vr = glr;
while( l>1 );
l = fix(l - 1);
eps = eps.*2.0e0;
ef = ef.*sq2;
if( lr(l)<=0 )
vl(l) = vl(l+1) + vr;
gt=1;
break;
else;
vr = vl(l+1) + vr;
end;
end;
if(gt==0)
break;
end;
end;
est = gr(l-1);
lr(l) = 1;
aa(l) = aa(l) + 4.0e0.*hh(l);
end;
%
%     Exit
%
ansmlv = vr;
if((mxl~=0) &&(abs(ce)>2.0e0.*tol.*area) )
ierr = 2;
xermsg('SLATEC','GAUS8','ANS is probably insufficiently accurate.',3,1);
end;
if( err<0.0e0 )
err = ce;
end;
end;
end %subroutine gaus8
%DECK GENBUN