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

function [fun,a,b,err,ansmlv,ierr,k]=dqnc79(fun,a,b,err,ansmlv,ierr,k);
%***BEGIN PROLOGUE  DQNC79
%***PURPOSE  Integrate a function using a 7-point adaptive Newton-Cotes
%            quadrature rule.
%***LIBRARY   SLATEC
%***CATEGORY  H2A1A1
%***TYPE      doubleprecision (QNC79-S, DQNC79-D)
%***KEYWORDS  ADAPTIVE QUADRATURE, INTEGRATION, NEWTON-COTES
%***AUTHOR  Kahaner, D. K., (NBS)
%           Jones, R. E., (SNLA)
%***DESCRIPTION
%
%     Abstract  *** a doubleprecision routine ***
%       DQNC79 is a general purpose program for evaluation of
%       one dimensional integrals of user defined functions.
%       DQNC79 will pick its own points for evaluation of the
%       integrand and these will vary from problem to problem.
%       Thus, DQNC79 is not designed to integrate over data sets.
%       Moderately smooth integrands will be integrated efficiently
%       and reliably.  For problems with strong singularities,
%       oscillations etc., the user may wish to use more sophis-
%       ticated routines such as those in QUADPACK.  One measure
%       of the reliability of DQNC79 is the output parameter K,
%       giving the number of integrand evaluations that were needed.
%
%     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 your calling
%              program.  You must write a Fortran function to evaluate
%              FUN.  This should be of the form
%                    doubleprecision function FUN (X)
%              C
%              C     X can vary from A to B
%              C     FUN(X) should be finite for all X on interval.
%              C
%                    FUN = ...
%                    RETURN
%                    end
%       A    - lower limit of integration
%       B    - upper limit of integration (may be less than A)
%       ERR  - is a requested error tolerance.  Normally, pick a value
%              0 .LT. ERR .LT. 1.0D-8.
%
%     --Output--
%       ANS  - computed value of the integral.  Hopefully, ANS is
%              accurate to within ERR * integral of ABS(FUN(X)).
%       IERR - a status code
%            - Normal codes
%               1  ANS most likely meets requested error tolerance.
%              -1  A equals B, or 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.
%       K    - the number of function evaluations actually used to do
%              the integration.  A value of K .GT. 1000 indicates a
%              difficult problem; other programs may be more efficient.
%              DQNC79 will gracefully give up if K exceeds 2000.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D1MACH, I1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   790601  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)
%   920218  Code redone to parallel QNC79.  (WRB)
%   930120  Increase array size 80->99, and KMX 2000->5000 for SUN -r8
%           wordlength.  (RWC)
%***end PROLOGUE  DQNC79
%     .. Scalar Arguments ..
persistent aa ae area bank blocal c ce ee ef eps f f1 f2 f3 f4 f5 f6 f7 first firstCall gt hh i kml kmx l lmn lmx lr nbits nib nlmn nlmx q13 q7 q7l q7r sq2 test tol vl vr w1 w2 w3 w4 ; if isempty(firstCall),firstCall=1;end; 

if isempty(gt), gt=zeros(1,2); end;
%     .. function Arguments ..
%     .. Local Scalars ..
if isempty(ae), ae=0; end;
if isempty(area), area=0; end;
if isempty(bank), bank=0; end;
if isempty(blocal), blocal=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(q13), q13=0; end;
if isempty(q7), q7=0; end;
if isempty(q7l), q7l=0; end;
if isempty(sq2), sq2=0; end;
if isempty(test), test=0; end;
if isempty(tol), tol=0; 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(i), i=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(nbits), nbits=0; end;
if isempty(nib), nib=0; end;
if isempty(nlmn), nlmn=0; end;
if isempty(nlmx), nlmx=0; end;
if isempty(first), first=false; end;
%     .. Local Arrays ..
if isempty(aa), aa=zeros(1,99); end;
if isempty(f), f=zeros(1,13); end;
if isempty(f1), f1=zeros(1,99); end;
if isempty(f2), f2=zeros(1,99); end;
if isempty(f3), f3=zeros(1,99); end;
if isempty(f4), f4=zeros(1,99); end;
if isempty(f5), f5=zeros(1,99); end;
if isempty(f6), f6=zeros(1,99); end;
if isempty(f7), f7=zeros(1,99); end;
if isempty(hh), hh=zeros(1,99); end;
if isempty(q7r), q7r=zeros(1,99); end;
if isempty(vl), vl=zeros(1,99); end;
if isempty(lr), lr=zeros(1,99); end;
%     .. External Functions ..
%     .. External Subroutines ..
%     .. Intrinsic Functions ..
%     .. Save statement ..
%     .. Data statements ..
if firstCall,   kml=[7];  end;
if firstCall,  kmx=[5000];  end;
if firstCall,  nlmn=[2];  end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DQNC79
gt(:)=0;
if( first )
w1 = 41.0d0./140.0d0;
w2 = 216.0d0./140.0d0;
w3 = 27.0d0./140.0d0;
w4 = 272.0d0./140.0d0;
nbits = fix(d1mach(5).*i1mach(14)./0.30102000d0);
nlmx = fix(min(99,fix((nbits.*4)./5)));
sq2 = sqrt(2.0d0);
end;
first = false;
ansmlv = 0.0d0;
ierr = 1;
ce = 0.0d0;
while (1);
if( a~=b )
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 )
break;
end;
nib = fix(0.5d0 - log(c)./log(2.0d0));
lmx = fix(min(nlmx,nbits-nib-4));
if( lmx<2 )
break;
end;
lmn = fix(min(lmn,lmx));
end;
end;
end;
tol = max(abs(err),2.0d0.^(5-nbits));
if( err==0.0d0 )
tol = sqrt(d1mach(4));
end;
eps = tol;
hh(1) =(b-a)./12.0d0;
aa(1) = a;
lr(1) = 1;
for i = 1 : 2: 11 ;
f(i) = fun(a+(i-1).*hh(1));
end; i = fix(11 +1);
blocal = b;
f(13) = fun(blocal);
k = 7;
l = 1;
area = 0.0d0;
q7 = 0.0d0;
ef = 256.0d0./255.0d0;
bank = 0.0d0;
%
%     Compute refined estimates, estimate the error, etc.
%
while( true );
gt(:)=0;
for i = 2 : 2: 12 ;
f(i) = fun(aa(l)+(i-1).*hh(l));
end; i = fix(12 +1);
k = fix(k + 6);
%
%     Compute left and right half estimates
%
q7l = hh(l).*((w1.*(f(1)+f(7))+w2.*(f(2)+f(6)))+(w3.*(f(3)+f(5))+w4.*f(4)));
q7r(l) = hh(l).*((w1.*(f(7)+f(13))+w2.*(f(8)+f(12)))+(w3.*(f(9)+f(11))+w4.*f(10)));
%
%     Update estimate of integral of absolute value
%
area = area +(abs(q7l)+abs(q7r(l))-abs(q7));
%
%     Do not bother to test convergence before minimum refinement level
%
if( l>=lmn )
%
%     Estimate the error in new value for whole interval, Q13
%
q13 = q7l + q7r(l);
ee = abs(q7-q13).*ef;
%
%     Compute nominal allowed error
%
ae = eps.*area;
%
%     Borrow from bank account, but not too much
%
test = min(ae+0.8d0.*bank,10.0d0.*ae);
%
%     Don't ask for excessive accuracy
%
test = max([test,tol.*abs(q13),0.00003d0.*tol.*area]);
%
%     Now, did this interval pass or not?
%
if( ee<=test )
%
%     On good intervals accumulate the theoretical estimate
%
ce = ce +(q7-q13)./255.0d0;
else;
%
%     Consider the left half of next deeper level
%
if( k>kmx )
lmx = fix(min(kml,lmx));
end;
if( l<lmx )
gt(2)=1;
end;
%
%     Have hit maximum refinement level -- penalize the cumulative error
%
if(gt(2)==0)
ce = ce +(q7-q13);
end;
end;
%
%     Update the bank account.  Don't go into debt.
%
if(gt(2)==0)
bank = bank +(ae-ee);
if( bank<0.0d0 )
bank = 0.0d0;
end;
%
%     Did we just finish a left half or a right half?
%
if( lr(l)<=0 )
%
%     Proceed to right half at this level
%
vl(l) = q13;
else;
%
%     Left and right halves are done, so go back up a level
%
vr = q13;
while( l>1 );
if( l<=17 )
ef = ef.*sq2;
end;
eps = eps.*2.0d0;
l = fix(l - 1);
if( lr(l)<=0 )
vl(l) = vl(l+1) + vr;
gt(1)=1;
break;
else;
vr = vl(l+1) + vr;
end;
end;
if(gt(1)==0)
break;
end;
end;
gt(1)=0;
q7 = q7r(l-1);
lr(l) = 1;
aa(l) = aa(l) + 12.0d0.*hh(l);
f(1) = f1(l);
f(3) = f2(l);
f(5) = f3(l);
f(7) = f4(l);
f(9) = f5(l);
f(11) = f6(l);
f(13) = f7(l);
continue;
end;
gt(2)=0;
end;
l = fix(l + 1);
eps = eps.*0.5d0;
if( l<=17 )
ef = ef./sq2;
end;
hh(l) = hh(l-1).*0.5d0;
lr(l) = -1;
aa(l) = aa(l-1);
q7 = q7l;
f1(l) = f(7);
f2(l) = f(8);
f3(l) = f(9);
f4(l) = f(10);
f5(l) = f(11);
f6(l) = f(12);
f7(l) = f(13);
f(13) = f(7);
f(11) = f(6);
f(9) = f(5);
f(7) = f(4);
f(5) = f(3);
f(3) = f(2);
end;
%
%     Exit
%
ansmlv = vr;
if( abs(ce)>2.0d0.*tol.*area )
ierr = 2;
xermsg('SLATEC','DQNC79','ANS is probably insufficiently accurate.',2,1);
end;
return;
end;
break;
end;
ierr = -1;
xermsg('SLATEC','DQNC79',['A and B are too nearly equal to allow normal integration. $$','ANS is set to zero and IERR to -1.'],-1,-1);
end %subroutine dqnc79
%DECK DQNG