[t,bcoef,n,k,x1,x2,bquad,work]=bsqad(t,bcoef,n,k,x1,x2,bquad,work);

function [t,bcoef,n,k,x1,x2,bquad,work]=bsqad(t,bcoef,n,k,x1,x2,bquad,work);
%***BEGIN PROLOGUE  BSQAD
%***PURPOSE  Compute the integral of a K-th order B-spline using the
%            B-representation.
%***LIBRARY   SLATEC
%***CATEGORY  H2A2A1, E3, K6
%***TYPE      SINGLE PRECISION (BSQAD-S, DBSQAD-D)
%***KEYWORDS  INTEGRAL OF B-SPLINES, QUADRATURE
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%     Abstract
%         BSQAD computes the integral on (X1,X2) of a K-th order
%         B-spline using the B-representation (T,BCOEF,N,K).  Orders
%         K as high as 20 are permitted by applying a 2, 6, or 10
%         point Gauss formula on subintervals of (X1,X2) which are
%         formed by included (distinct) knots.
%
%         If orders K greater than 20 are needed, use BFQAD with
%         F(X) = 1.
%
%     Description of Arguments
%         Input
%           T      - knot array of length N+K
%           BCOEF  - B-spline coefficient array of length N
%           N      - length of coefficient array
%           K      - order of B-spline, 1 .LE. K .LE. 20
%           X1,X2  - end points of quadrature interval in
%                    T(K) .LE. X .LE. T(N+1)
%
%         Output
%           BQUAD  - integral of the B-spline over (X1,X2)
%           WORK   - work vector of length 3*K
%
%     Error Conditions
%         Improper input is a fatal error
%
%***REFERENCES  D. E. Amos, Quadrature subroutines for splines and
%                 B-splines, Report SAND79-1825, Sandia Laboratories,
%                 December 1979.
%***ROUTINES CALLED  BVALU, INTRV, 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  BSQAD
%
persistent a aa b bb bma bpa c1 firstCall gpts gwts gx i il1 il2 ilo inbv jf left m mf mflag np1 npk q summlv ta tb y1 y2 ; if isempty(firstCall),firstCall=1;end; 

if isempty(i), i=0; end;
if isempty(il1), il1=0; end;
if isempty(il2), il2=0; end;
if isempty(ilo), ilo=0; end;
if isempty(inbv), inbv=0; end;
if isempty(jf), jf=0; end;
if isempty(left), left=0; end;
if isempty(m), m=0; end;
if isempty(mf), mf=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(b), b=0; end;
if isempty(bb), bb=0; end;
if isempty(bma), bma=0; end;
if isempty(bpa), bpa=0; end;
if isempty(c1), c1=0; end;
if isempty(gpts), gpts=zeros(1,9); end;
if isempty(gwts), gwts=zeros(1,9); end;
if isempty(gx), gx=0; end;
if isempty(q), q=0; end;
if isempty(summlv), summlv=zeros(1,5); end;
if isempty(ta), ta=0; end;
if isempty(tb), tb=0; end;
if isempty(y1), y1=0; end;
if isempty(y2), y2=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,[]);
%
if firstCall,   gpts(1) =[5.77350269189625764e-01];  end;
if firstCall,  gpts(2) =[2.38619186083196909e-01];  end;
if firstCall,  gpts(3) =[6.61209386466264514e-01];  end;
if firstCall,  gpts(4) =[9.32469514203152028e-01];  end;
if firstCall,  gpts(5) =[1.48874338981631211e-01];  end;
if firstCall,  gpts(6) =[4.33395394129247191e-01];  end;
if firstCall, gpts(7) =[6.79409568299024406e-01];  end;
if firstCall,  gpts(8) =[8.65063366688984511e-01];  end;
if firstCall,  gpts(9)=[9.73906528517171720e-01];  end;
if firstCall,   gwts(1) =[1.00000000000000000e+00];  end;
if firstCall,  gwts(2) =[4.67913934572691047e-01];  end;
if firstCall,  gwts(3) =[3.60761573048138608e-01];  end;
if firstCall,  gwts(4) =[1.71324492379170345e-01];  end;
if firstCall,  gwts(5) =[2.95524224714752870e-01];  end;
if firstCall,  gwts(6) =[2.69266719309996355e-01];  end;
if firstCall, gwts(7) =[2.19086362515982044e-01];  end;
if firstCall,  gwts(8) =[1.49451349150580593e-01];  end;
if firstCall,  gwts(9)=[6.66713443086881376e-02];  end;
firstCall=0;
%
%***FIRST EXECUTABLE STATEMENT  BSQAD
bquad = 0.0e0;
if( k<1 || k>20 )
xermsg('SLATEC','BSQAD','K DOES NOT SATISFY 1.LE.K.LE.20',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','BSQAD','N DOES NOT SATISFY N.GE.K',2,1);
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);
%     SELECTION OF 2, 6, OR 10 POINT GAUSS FORMULA
jf = 0;
mf = 1;
if( k>4 )
jf = 1;
mf = 3;
if( k>12 )
jf = 4;
mf = 5;
end;
end;
%
for i = 1 : mf;
summlv(i) = 0.0e0;
end; i = fix(mf+1);
ilo = 1;
inbv = 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;
for left = il1 : il2;
ta = t(left);
tb = t(left+1);
if( ta~=tb )
a = max(aa,ta);
b = min(bb,tb);
bma = 0.5e0.*(b-a);
bpa = 0.5e0.*(b+a);
for m = 1 : mf;
c1 = bma.*gpts(jf+m);
gx = -c1 + bpa;
[y2 ,t,bcoef,n,k,dumvar6,gx,inbv,work]=bvalu(t,bcoef,n,k,0,gx,inbv,work);
gx = c1 + bpa;
[y1 ,t,bcoef,n,k,dumvar6,gx,inbv,work]=bvalu(t,bcoef,n,k,0,gx,inbv,work);
summlv(m) = summlv(m) +(y1+y2).*bma;
end; m = fix(mf+1);
end;
end; left = fix(il2+1);
q = 0.0e0;
for m = 1 : mf;
q = q + gwts(jf+m).*summlv(m);
end; m = fix(mf+1);
if( x1>x2 )
q = -q;
end;
bquad = 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','BSQAD','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;
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 BSRH