[f,a,b,result,abserr,resabs,resasc]=dqk15(f,a,b,result,abserr,resabs,resasc);

function [f,a,b,result,abserr,resabs,resasc]=dqk15(f,a,b,result,abserr,resabs,resasc);
%***BEGIN PROLOGUE  DQK15
%***PURPOSE  To compute I = Integral of F over (A,B), with error
%                           estimate
%                       J = integral of ABS(F) over (A,B)
%***LIBRARY   SLATEC (QUADPACK)
%***CATEGORY  H2A1A2
%***TYPE      doubleprecision (QK15-S, DQK15-D)
%***KEYWORDS  15-POINT GAUSS-KRONROD RULES, QUADPACK, QUADRATURE
%***AUTHOR  Piessens, Robert
%             Applied Mathematics and Programming Division
%             K. U. Leuven
%           de Doncker, Elise
%             Applied Mathematics and Programming Division
%             K. U. Leuven
%***DESCRIPTION
%
%           Integration rules
%           Standard fortran subroutine
%           doubleprecision version
%
%           PARAMETERS
%            ON ENTRY
%              F      - doubleprecision
%                       function subprogram defining the integrand
%                       function F(X). The actual name for F needs to be
%                       Declared E X T E R N A L in the calling program.
%
%              A      - doubleprecision
%                       Lower limit of integration
%
%              B      - doubleprecision
%                       Upper limit of integration
%
%            ON RETURN
%              RESULT - doubleprecision
%                       Approximation to the integral I
%                       Result is computed by applying the 15-POINT
%                       KRONROD RULE (RESK) obtained by optimal addition
%                       of abscissae to the 7-POINT GAUSS RULE(RESG).
%
%              ABSERR - doubleprecision
%                       Estimate of the modulus of the absolute error,
%                       which should not exceed ABS(I-RESULT)
%
%              RESABS - doubleprecision
%                       Approximation to the integral J
%
%              RESASC - doubleprecision
%                       Approximation to the integral of ABS(F-I/(B-A))
%                       over (A,B)
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D1MACH
%***REVISION HISTORY  (YYMMDD)
%   800101  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)
%***end PROLOGUE  DQK15
%
persistent absc centr dhlgth epmach fc firstCall fsum fv1 fv2 fval1 fval2 hlgth j jtw jtwm1 resg resk reskh uflow wg wgk xgk ; if isempty(firstCall),firstCall=1;end; 

if isempty(absc), absc=0; end;
if isempty(centr), centr=0; end;
if isempty(dhlgth), dhlgth=0; end;
if isempty(epmach), epmach=0; end;
if isempty(fc), fc=0; end;
if isempty(fsum), fsum=0; end;
if isempty(fval1), fval1=0; end;
if isempty(fval2), fval2=0; end;
if isempty(fv1), fv1=zeros(1,7); end;
if isempty(fv2), fv2=zeros(1,7); end;
if isempty(hlgth), hlgth=0; end;
if isempty(resg), resg=0; end;
if isempty(resk), resk=0; end;
if isempty(reskh), reskh=0; end;
if isempty(uflow), uflow=0; end;
if isempty(wg), wg=zeros(1,4); end;
if isempty(wgk), wgk=zeros(1,8); end;
if isempty(xgk), xgk=zeros(1,8); end;
if isempty(j), j=0; end;
if isempty(jtw), jtw=0; end;
if isempty(jtwm1), jtwm1=0; end;
%
%
%           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
%           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
%           CORRESPONDING WEIGHTS ARE GIVEN.
%
%           XGK    - ABSCISSAE OF THE 15-POINT KRONROD RULE
%                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 7-POINT
%                    GAUSS RULE
%                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
%                    ADDED TO THE 7-POINT GAUSS RULE
%
%           WGK    - WEIGHTS OF THE 15-POINT KRONROD RULE
%
%           WG     - WEIGHTS OF THE 7-POINT GAUSS RULE
%
%
% GAUSS QUADRATURE WEIGHTS AND KRONROD QUADRATURE ABSCISSAE AND WEIGHTS
% AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON,
% BELL LABS, NOV. 1981.
%
if firstCall,   wg(1)=[0.129484966168869693270611432679082d0];  end;
if firstCall,   wg(2)=[0.279705391489276667901467771423780d0];  end;
if firstCall,   wg(3)=[0.381830050505118944950369775488975d0];  end;
if firstCall,   wg(4)=[0.417959183673469387755102040816327d0];  end;
%
if firstCall,   xgk(1)=[0.991455371120812639206854697526329d0];  end;
if firstCall,   xgk(2)=[0.949107912342758524526189684047851d0];  end;
if firstCall,   xgk(3)=[0.864864423359769072789712788640926d0];  end;
if firstCall,   xgk(4)=[0.741531185599394439863864773280788d0];  end;
if firstCall,   xgk(5)=[0.586087235467691130294144838258730d0];  end;
if firstCall,   xgk(6)=[0.405845151377397166906606412076961d0];  end;
if firstCall,   xgk(7)=[0.207784955007898467600689403773245d0];  end;
if firstCall,   xgk(8)=[0.000000000000000000000000000000000d0];  end;
%
if firstCall,   wgk(1)=[0.022935322010529224963732008058970d0];  end;
if firstCall,   wgk(2)=[0.063092092629978553290700663189204d0];  end;
if firstCall,   wgk(3)=[0.104790010322250183839876322541518d0];  end;
if firstCall,   wgk(4)=[0.140653259715525918745189590510238d0];  end;
if firstCall,   wgk(5)=[0.169004726639267902826583426598550d0];  end;
if firstCall,   wgk(6)=[0.190350578064785409913256402421014d0];  end;
if firstCall,   wgk(7)=[0.204432940075298892414161999234649d0];  end;
if firstCall,   wgk(8)=[0.209482141084727828012999174891714d0];  end;
firstCall=0;
%
%
%           LIST OF MAJOR VARIABLES
%           -----------------------
%
%           CENTR  - MID POINT OF THE INTERVAL
%           HLGTH  - HALF-LENGTH OF THE INTERVAL
%           ABSC   - ABSCISSA
%           FVAL*  - FUNCTION VALUE
%           RESG   - RESULT OF THE 7-POINT GAUSS FORMULA
%           RESK   - RESULT OF THE 15-POINT KRONROD FORMULA
%           RESKH  - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B),
%                    I.E. TO I/(B-A)
%
%           MACHINE DEPENDENT CONSTANTS
%           ---------------------------
%
%           EPMACH IS THE LARGEST RELATIVE SPACING.
%           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
%
%***FIRST EXECUTABLE STATEMENT  DQK15
[epmach ]=d1mach(4);
[uflow ]=d1mach(1);
%
centr = 0.5d+00.*(a+b);
hlgth = 0.5d+00.*(b-a);
dhlgth = abs(hlgth);
%
%           COMPUTE THE 15-POINT KRONROD APPROXIMATION TO
%           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
%
fc = f(centr);
resg = fc.*wg(4);
resk = fc.*wgk(8);
resabs = abs(resk);
for j = 1 : 3;
jtw = fix(j.*2);
absc = hlgth.*xgk(jtw);
fval1 = f(centr-absc);
fval2 = f(centr+absc);
fv1(jtw) = fval1;
fv2(jtw) = fval2;
fsum = fval1 + fval2;
resg = resg + wg(j).*fsum;
resk = resk + wgk(jtw).*fsum;
resabs = resabs + wgk(jtw).*(abs(fval1)+abs(fval2));
end; j = fix(3+1);
for j = 1 : 4;
jtwm1 = fix(j.*2 - 1);
absc = hlgth.*xgk(jtwm1);
fval1 = f(centr-absc);
fval2 = f(centr+absc);
fv1(jtwm1) = fval1;
fv2(jtwm1) = fval2;
fsum = fval1 + fval2;
resk = resk + wgk(jtwm1).*fsum;
resabs = resabs + wgk(jtwm1).*(abs(fval1)+abs(fval2));
end; j = fix(4+1);
reskh = resk.*0.5d+00;
resasc = wgk(8).*abs(fc-reskh);
for j = 1 : 7;
resasc = resasc + wgk(j).*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh));
end; j = fix(7+1);
result = resk.*hlgth;
resabs = resabs.*dhlgth;
resasc = resasc.*dhlgth;
abserr = abs((resk-resg).*hlgth);
if( resasc~=0.0d+00 && abserr~=0.0d+00 )
abserr = resasc.*min(0.1d+01,(0.2d+03.*abserr./resasc).^1.5d+00);
end;
if( resabs>uflow./(0.5d+02.*epmach) )
abserr = max((epmach.*0.5d+02).*resabs,abserr);
end;
end
%DECK DQK15I