[f,w,p1,p2,p3,p4,kp,a,b,result,abserr,resabs,resasc]=dqk15w(f,w,p1,p2,p3,p4,kp,a,b,result,abserr,resabs,resasc);

function [f,w,p1,p2,p3,p4,kp,a,b,result,abserr,resabs,resasc]=dqk15w(f,w,p1,p2,p3,p4,kp,a,b,result,abserr,resabs,resasc);
%***BEGIN PROLOGUE  DQK15W
%***PURPOSE  To compute I = Integral of F*W over (A,B), with error
%                           estimate
%                       J = Integral of ABS(F*W) over (A,B)
%***LIBRARY   SLATEC (QUADPACK)
%***CATEGORY  H2A2A2
%***TYPE      doubleprecision (QK15W-S, DQK15W-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 driver program.
%
%              W      - doubleprecision
%                       function subprogram defining the integrand
%                       WEIGHT function W(X). The actual name for W
%                       needs to be declared E X T E R N A L in the
%                       calling program.
%
%              P1, P2, P3, P4 - doubleprecision
%                       Parameters in the WEIGHT function
%
%              KP     - Integer
%                       Key for indicating the type of WEIGHT function
%
%              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 equal or exceed ABS(I-RESULT)
%
%              RESABS - doubleprecision
%                       Approximation to the integral of ABS(F)
%
%              RESASC - doubleprecision
%                       Approximation to the integral of ABS(F-I/(B-A))
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D1MACH
%***REVISION HISTORY  (YYMMDD)
%   810101  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  DQK15W
%
persistent absc absc1 absc2 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(absc1), absc1=0; end;
if isempty(absc2), absc2=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 GAUSS-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 GAUSS-KRONROD RULE
%
%           WG     - WEIGHTS OF THE 7-POINT GAUSS RULE
%
if firstCall,   xgk(1) =[0.9914553711208126d+00];  end;
if firstCall,  xgk(2) =[0.9491079123427585d+00];  end;
if firstCall,  xgk(3) =[0.8648644233597691d+00];  end;
if firstCall,  xgk(4) =[0.7415311855993944d+00];  end;
if firstCall,  xgk(5) =[0.5860872354676911d+00];  end;
if firstCall,  xgk(6) =[0.4058451513773972d+00];  end;
if firstCall,  xgk(7)=[0.2077849550078985d+00];  end;
if firstCall,  xgk(8)=[0.0000000000000000d+00];  end;
%
if firstCall,   wgk(1) =[0.2293532201052922d-01];  end;
if firstCall,  wgk(2) =[0.6309209262997855d-01];  end;
if firstCall,  wgk(3) =[0.1047900103222502d+00];  end;
if firstCall,  wgk(4) =[0.1406532597155259d+00];  end;
if firstCall,  wgk(5) =[0.1690047266392679d+00];  end;
if firstCall,  wgk(6) =[0.1903505780647854d+00];  end;
if firstCall,  wgk(7)=[0.2044329400752989d+00];  end;
if firstCall,  wgk(8)=[0.2094821410847278d+00];  end;
%
if firstCall,   wg(1) =[0.1294849661688697d+00];  end;
if firstCall,  wg(2) =[0.2797053914892767d+00];  end;
if firstCall,  wg(3) =[0.3818300505051889d+00];  end;
if firstCall,  wg(4)=[0.4179591836734694d+00];  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*W 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  DQK15W
[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 ERROR.
%
fc = f(centr).*w(centr,p1,p2,p3,p4,kp);
resg = wg(4).*fc;
resk = wgk(8).*fc;
resabs = abs(resk);
for j = 1 : 3;
jtw = fix(j.*2);
absc = hlgth.*xgk(jtw);
absc1 = centr - absc;
absc2 = centr + absc;
fval1 = f(absc1).*w(absc1,p1,p2,p3,p4,kp);
fval2 = f(absc2).*w(absc2,p1,p2,p3,p4,kp);
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);
absc1 = centr - absc;
absc2 = centr + absc;
fval1 = f(absc1).*w(absc1,p1,p2,p3,p4,kp);
fval2 = f(absc2).*w(absc2,p1,p2,p3,p4,kp);
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 DQK21