[f,t,bcoef,n,k,id,x1,x2,tol,quad,ierr,work]=bfqad(f,t,bcoef,n,k,id,x1,x2,tol,quad,ierr,work);

function [f,t,bcoef,n,k,id,x1,x2,tol,quad,ierr,work]=bfqad(f,t,bcoef,n,k,id,x1,x2,tol,quad,ierr,work);
persistent a aa ansmlv b bb iflg il1 il2 ilo inbv left mflag np1 npk q ta tb wtol ; 

if isempty(inbv), inbv=0; end;
%***BEGIN PROLOGUE  BFQAD
%***PURPOSE  Compute the integral of a product of a function and a
%            derivative of a B-spline.
%***LIBRARY   SLATEC
%***CATEGORY  H2A2A1, E3, K6
%***TYPE      SINGLE PRECISION (BFQAD-S, DBFQAD-D)
%***KEYWORDS  INTEGRAL OF B-SPLINE, QUADRATURE
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%     Abstract
%         BFQAD computes the integral on (X1,X2) of a product of a
%         function F and the ID-th derivative of a K-th order B-spline,
%         using the B-representation (T,BCOEF,N,K).  (X1,X2) must be
%         a subinterval of T(K) .LE. X <= T(N+1).  An integration
%         routine BSGQ8 (a modification
%         of GAUS8), integrates the product on sub-
%         intervals of (X1,X2) formed by included (distinct) knots.
%
%     Description of Arguments
%         Input
%           F      - external function of one argument for the
%                    integrand BF(X)=F(X)*BVALU(T,BCOEF,N,K,ID,X,INBV,
%                    WORK)
%           T      - knot array of length N+K
%           BCOEF  - coefficient array of length N
%           N      - length of coefficient array
%           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 in
%                    T(K) .LE. X .LE. T(N+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 BF(X) on (X1,X2)
%           IERR   - a status code
%                    IERR=1  normal return
%                         2  some quadrature on (X1,X2) does not meet
%                            the requested tolerance.
%           WORK   - work vector of length 3*K
%
%     Error Conditions
%         X1 or X2 not in T(K) .LE. X .LE. T(N+1) is a fatal error.
%         TOL not greater than the single precision unit roundoff or
%         less than 0.1 is a fatal error.
%         Some quadrature fails to meet the requested tolerance.
%
%***REFERENCES  D. E. Amos, Quadrature subroutines for splines and
%                 B-splines, Report SAND79-1825, Sandia Laboratories,
%                 December 1979.
%***ROUTINES CALLED  BSGQ8, INTRV, 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  BFQAD
%
%
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(left), left=0; end;
if isempty(mflag), mflag=0; end;
if isempty(npk), npk=0; end;
if isempty(np1), np1=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;
t_shape=size(t);t=reshape(t,1,[]);
bcoef_shape=size(bcoef);bcoef=reshape(bcoef,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
%***FIRST EXECUTABLE STATEMENT  BFQAD
ierr = 1;
quad = 0.0e0;
if( k<1 )
xermsg('SLATEC','BFQAD','K DOES NOT SATISFY K.GE.1',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
elseif( n<k ) ;
xermsg('SLATEC','BFQAD','N DOES NOT SATISFY N.GE.K',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
elseif( id<0 || id>=k ) ;
xermsg('SLATEC','BFQAD','ID DOES NOT SATISFY 0 .LE. ID .LT. K',2,1);
else;
[wtol ]=r1mach(4);
if( tol<wtol || tol>0.1e0 )
xermsg('SLATEC','BFQAD',['TOL IS LESS THAN THE SINGLE PRECISION TOLERANCE OR ','GREATER THAN 0.1'],2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
else;
aa = min(x1,x2);
bb = max(x1,x2);
if( aa>=t(k) )
np1 = fix(n + 1);
if( bb<=t(np1) )
if( aa==bb )
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
npk = fix(n + k);
%
ilo = 1;
[t,npk,aa,ilo,il1,mflag]=intrv(t,npk,aa,ilo,il1,mflag);
[t,npk,bb,ilo,il2,mflag]=intrv(t,npk,bb,ilo,il2,mflag);
if( il2>=np1 )
il2 = fix(n);
end;
inbv = 1;
q = 0.0e0;
for left = il1 : il2;
ta = t(left);
tb = t(left+1);
if( ta~=tb )
a = max(aa,ta);
b = min(bb,tb);
[f,t,bcoef,n,k,id,a,b,inbv,tol,ansmlv,iflg,work]=bsgq8(f,t,bcoef,n,k,id,a,b,inbv,tol,ansmlv,iflg,work);
if( iflg>1 )
ierr = 2;
end;
q = q + ansmlv;
end;
end; left = fix(il2+1);
if( x1>x2 )
q = -q;
end;
quad = q;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
end;
%
%
xermsg('SLATEC','BFQAD','X1 OR X2 OR BOTH DO NOT SATISFY T(K).LE.X.LE.T(N+1)',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
end;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end
%DECK BIE