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

function [fun,a,b,err,ansmlv,ierr]=dgaus8(fun,a,b,err,ansmlv,ierr);
%***BEGIN PROLOGUE  DGAUS8
%***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      doubleprecision (GAUS8-S, DGAUS8-D)
%***KEYWORDS  ADAPTIVE QUADRATURE, AUTOMATIC INTEGRATOR,
%             GAUSS QUADRATURE, NUMERICAL INTEGRATION
%***AUTHOR  Jones, R. E., (SNLA)
%***DESCRIPTION
%
%     Abstract  *** a doubleprecision routine ***
%        DGAUS8 integrates real functions of one variable over finite
%        intervals using an adaptive 8-point Legendre-Gauss algorithm.
%        DGAUS8 is intended primarily for high accuracy integration
%        or integration of smooth functions.
%
%        The maximum number of significant digits obtainable in ANS
%        is the smaller of 18 and the number of digits carried in
%        doubleprecision arithmetic.
%
%     Description of Arguments
%
%        Input--* FUN, A, B, ERR are doubleprecision *
%        FUN - name of external function to be integrated.  This name
%              must be in an EXTERNAL statement in the calling program.
%              FUN must be a doubleprecision function of one DOUBLE
%              PRECISION 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 DTOL .LT. ABS(ERR) .LE.
%              1.0D-3 where DTOL is the larger of 1.0D-18 and the
%              doubleprecision unit roundoff D1MACH(4).  ANS will
%              normally have no more error than ABS(ERR) times the
%              integral of the absolute value of FUN(X).  Usually,
%              smaller values of 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,ANS are doubleprecision *
%        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  D1MACH, I1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   810223  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890911  Removed unnecessary intrinsics.  (WRB)
%   890911  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  DGAUS8
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,60); 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,60); 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,60); end;
if isempty(hh), hh=zeros(1,60); end;
if isempty(sq2), sq2=0; end;
if isempty(tol), tol=0; end;
if isempty(vl), vl=zeros(1,60); 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))));doubleprecision :: g8 ;
if firstCall,   x1 =[1.83434642495649805d-01];  end;
if firstCall,  x2 =[5.25532409916328986d-01];  end;
if firstCall,  x3 =[7.96666477413626740d-01];  end;
if firstCall,  x4=[9.60289856497536232d-01];  end;
if firstCall,   w1 =[3.62683783378361983d-01];  end;
if firstCall,  w2 =[3.13706645877887287d-01];  end;
if firstCall,  w3 =[2.22381034453374471d-01];  end;
if firstCall,  w4=[1.01228536290376259d-01];  end;
if firstCall,   sq2=[1.41421356d0];  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  DGAUS8
%
%     Initialize
%
[k ]=i1mach(14);
anib = d1mach(5).*k./0.30102000d0;
nbits = fix(anib);
nlmx = fix(min(60,fix((nbits.*5)./8)));
ansmlv = 0.0d0;
ierr = 1;
ce = 0.0d0;
if( a==b )
if( err<0.0d0 )
err = ce;
end;
else;
lmx = fix(nlmx);
lmn = fix(nlmn);
if( b~=0.0d0 )
if( (abs(1.0d0).*sign(b)).*a>0.0d0 )
c = abs(1.0d0-a./b);
if( c<=0.1d0 )
if( c<=0.0d0 )
if( err<0.0d0 )
err = ce;
end;
return;
else;
anib = 0.5d0 - log(c)./0.69314718d0;
nib = fix(anib);
lmx = fix(min(nlmx,nbits-nib-7));
if( lmx<1 )
ierr = -1;
xermsg('SLATEC','DGAUS8',['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.0d0 )
err = ce;
end;
return;
else;
lmn = fix(min(lmn,lmx));
end;
end;
end;
end;
end;
tol = max(abs(err),2.0d0.^(5-nbits))./2.0d0;
if( err==0.0d0 )
tol = sqrt(d1mach(4));
end;
eps = tol;
hh(1) =(b-a)./4.0d0;
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.0d0.*hh(l),2.0d0.*hh(l));
k = 8;
area = abs(est);
ef = 0.5d0;
mxl = 0;
%
%     Compute refined estimates, estimate the error, etc.
%
while( true );
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.0d0.*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.5d0;
ef = ef./sq2;
hh(l) = hh(l-1).*0.5d0;
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;
gt=0;
while( l>1 );
l = fix(l - 1);
eps = eps.*2.0d0;
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.0d0.*hh(l);
end;
%
%     Exit
%
ansmlv = vr;
if((mxl~=0) &&(abs(ce)>2.0d0.*tol.*area) )
ierr = 2;
xermsg('SLATEC','DGAUS8','ANS is probably insufficiently accurate.',3,1);
end;
if( err<0.0d0 )
err = ce;
end;
end;
end %subroutine dgaus8
%DECK DGBCO