[f,ldc,c,xi,lxi,k,id,x1,x2,tol,quad,ierr]=pfqad(f,ldc,c,xi,lxi,k,id,x1,x2,tol,quad,ierr);

function [f,ldc,c,xi,lxi,k,id,x1,x2,tol,quad,ierr]=pfqad(f,ldc,c,xi,lxi,k,id,x1,x2,tol,quad,ierr);
%***BEGIN PROLOGUE  PFQAD
%***PURPOSE  Compute the integral on (X1,X2) of a product of a function
%            F and the ID-th derivative of a B-spline,
%            (PP-representation).
%***LIBRARY   SLATEC
%***CATEGORY  H2A2A1, E3, K6
%***TYPE      SINGLE PRECISION (PFQAD-S, DPFQAD-D)
%***KEYWORDS  B-SPLINE, DATA FITTING, INTERPOLATION, QUADRATURE, SPLINES
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%     Abstract
%         PFQAD computes the integral on (X1,X2) of a product of a
%         function F and the ID-th derivative of a B-spline, using the
%         PP-representation (C,XI,LXI,K). (X1,X2) is normally a sub-
%         interval of XI(1) .LE. X .LE. XI(LXI+1).  An integration rou-
%         tine, PPGQ8(a modification of GAUS8), integrates the product
%         on sub-intervals of (X1,X2) formed by the included break
%         points.  Integration outside of (XI(1),XI(LXI+1)) is permitted
%         provided F is defined.
%
%     Description of Arguments
%         Input
%           F      - external function of one argument for the
%                    integrand PF(X)=F(X)*PPVAL(LDC,C,XI,LXI,K,ID,X,
%                    INPPV)
%           LDC    - leading dimension of matrix C, LDC .GE. K
%           C(I,J) - right Taylor derivatives at XI(J), I=1,K , J=1,LXI
%           XI(*)  - break point array of length LXI+1
%           LXI    - number of polynomial pieces
%           K      - order of B-spline, K .GE. 1
%           ID     - order of the spline derivative, 0 .LE. ID .LE. K-1
%                    ID=0 gives the spline function
%           X1,X2  - end points of quadrature interval, normally in
%                    XI(1) .LE. X .LE. XI(LXI+1)
%           TOL    - desired accuracy for the quadrature, suggest
%                    10.*STOL .LT. TOL .LE. 0.1 where STOL is the single
%                    precision unit roundoff for the machine = R1MACH(4)
%
%         Output
%           QUAD   - integral of PF(X) on (X1,X2)
%           IERR   - a status code
%                    IERR=1 normal return
%                         2 some quadrature does not meet the
%                           requested tolerance
%
%     Error Conditions
%         TOL not greater than the single precision unit roundoff or
%         less than 0.1 is a fatal error.
%         Some quadrature does not meet the requested tolerance.
%
%***REFERENCES  D. E. Amos, Quadrature subroutines for splines and
%                 B-splines, Report SAND79-1825, Sandia Laboratories,
%                 December 1979.
%***ROUTINES CALLED  INTRV, PPGQ8, R1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   800901  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)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  PFQAD
%
persistent a aa ansmlv b bb iflg il1 il2 ilo inppv left mf1 mf2 q ta tb wtol ; 

if isempty(iflg), iflg=0; end;
if isempty(ilo), ilo=0; end;
if isempty(il1), il1=0; end;
if isempty(il2), il2=0; end;
if isempty(inppv), inppv=0; end;
if isempty(left), left=0; end;
if isempty(mf1), mf1=0; end;
if isempty(mf2), mf2=0; end;
if isempty(a), a=0; end;
if isempty(aa), aa=0; end;
if isempty(ansmlv), ansmlv=0; end;
if isempty(b), b=0; end;
if isempty(bb), bb=0; end;
if isempty(q), q=0; end;
if isempty(ta), ta=0; end;
if isempty(tb), tb=0; end;
if isempty(wtol), wtol=0; end;
xi_shape=size(xi);xi=reshape(xi,1,[]);
c_shape=size(c);c=reshape([c(:).',zeros(1,ceil(numel(c)./prod([ldc])).*prod([ldc])-numel(c))],abs(ldc),[]);
%
%***FIRST EXECUTABLE STATEMENT  PFQAD
ierr = 1;
quad = 0.0e0;
if( k<1 )
xermsg('SLATEC','PFQAD','K DOES NOT SATISFY K.GE.1',2,1);
xi_shape=zeros(xi_shape);xi_shape(:)=xi(1:numel(xi_shape));xi=xi_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
elseif( ldc<k ) ;
xermsg('SLATEC','PFQAD','LDC DOES NOT SATISFY LDC.GE.K',2,1);
xi_shape=zeros(xi_shape);xi_shape(:)=xi(1:numel(xi_shape));xi=xi_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
elseif( id<0 || id>=k ) ;
xermsg('SLATEC','PFQAD','ID DOES NOT SATISFY 0.LE.ID.LT.K',2,1);
xi_shape=zeros(xi_shape);xi_shape(:)=xi(1:numel(xi_shape));xi=xi_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
elseif( lxi<1 ) ;
xermsg('SLATEC','PFQAD','LXI DOES NOT SATISFY LXI.GE.1',2,1);
else;
[wtol ]=r1mach(4);
if( tol<wtol || tol>0.1e0 )
%
xermsg('SLATEC','PFQAD',['TOL IS LESS THAN THE SINGLE PRECISION TOLERANCE OR ','GREATER THAN 0.1'],2,1);
xi_shape=zeros(xi_shape);xi_shape(:)=xi(1:numel(xi_shape));xi=xi_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
else;
aa = min(x1,x2);
bb = max(x1,x2);
if( aa==bb )
xi_shape=zeros(xi_shape);xi_shape(:)=xi(1:numel(xi_shape));xi=xi_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
end;
ilo = 1;
[xi,lxi,aa,ilo,il1,mf1]=intrv(xi,lxi,aa,ilo,il1,mf1);
[xi,lxi,bb,ilo,il2,mf2]=intrv(xi,lxi,bb,ilo,il2,mf2);
q = 0.0e0;
inppv = 1;
for left = il1 : il2;
ta = xi(left);
a = max(aa,ta);
if( left==1 )
a = aa;
end;
tb = bb;
if( left<lxi )
tb = xi(left+1);
end;
b = min(bb,tb);
[f,ldc,c,xi,lxi,k,id,a,b,inppv,tol,ansmlv,iflg]=ppgq8(f,ldc,c,xi,lxi,k,id,a,b,inppv,tol,ansmlv,iflg);
if( iflg>1 )
ierr = 2;
end;
q = q + ansmlv;
end; left = fix(il2+1);
if( x1>x2 )
q = -q;
end;
quad = q;
xi_shape=zeros(xi_shape);xi_shape(:)=xi(1:numel(xi_shape));xi=xi_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
end;
end;
xi_shape=zeros(xi_shape);xi_shape(:)=xi(1:numel(xi_shape));xi=xi_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
end
%DECK PGSF