[ldc,c,xi,lxi,k,x1,x2,pquad]=dppqad(ldc,c,xi,lxi,k,x1,x2,pquad);

function [ldc,c,xi,lxi,k,x1,x2,pquad]=dppqad(ldc,c,xi,lxi,k,x1,x2,pquad);
%***BEGIN PROLOGUE  DPPQAD
%***PURPOSE  Compute the integral on (X1,X2) of a K-th order B-spline
%            using the piecewise polynomial (PP) representation.
%***LIBRARY   SLATEC
%***CATEGORY  H2A2A1, E3, K6
%***TYPE      doubleprecision (PPQAD-S, DPPQAD-D)
%***KEYWORDS  B-SPLINE, DATA FITTING, INTERPOLATION, QUADRATURE, SPLINES
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%     Abstract    **** a doubleprecision routine ****
%         DPPQAD computes the integral on (X1,X2) of a K-th order
%         B-spline using the piecewise polynomial representation
%         (C,XI,LXI,K).  Here the Taylor expansion about the left
%         end point XI(J) of the J-th interval is integrated and
%         evaluated on subintervals of (X1,X2) which are formed by
%         included break points.  Integration outside (XI(1),XI(LXI+1))
%         is permitted.
%
%     Description of Arguments
%         Input      C,XI,X1,X2 are doubleprecision
%           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
%           X1,X2  - end points of quadrature interval, normally in
%                    XI(1) .LE. X .LE. XI(LXI+1)
%
%         Output     PQUAD is doubleprecision
%           PQUAD  - integral of the PP representation over (X1,X2)
%
%     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  DINTRV, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   800901  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890911  Removed unnecessary intrinsics.  (WRB)
%   890911  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)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  DPPQAD
%
persistent a aa bb dx flk i ii il il1 il2 ilo im left mf1 mf2 q s ss ta tb x ; 

if isempty(i), i=0; end;
if isempty(ii), ii=0; end;
if isempty(il), il=0; end;
if isempty(ilo), ilo=0; end;
if isempty(il1), il1=0; end;
if isempty(il2), il2=0; end;
if isempty(im), im=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(bb), bb=0; end;
if isempty(dx), dx=0; end;
if isempty(flk), flk=0; end;
if isempty(q), q=0; end;
if isempty(s), s=0; end;
if isempty(ss), ss=zeros(1,2); end;
if isempty(ta), ta=0; end;
if isempty(tb), tb=0; end;
if isempty(x), x=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  DPPQAD
pquad = 0.0d0;
if( k<1 )
%
%
xermsg('SLATEC','DPPQAD','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( lxi<1 ) ;
xermsg('SLATEC','DPPQAD','LXI DOES NOT SATISFY LXI.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','DPPQAD','LDC DOES NOT SATISFY LDC.GE.K',2,1);
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]=dintrv(xi,lxi,aa,ilo,il1,mf1);
[xi,lxi,bb,ilo,il2,mf2]=dintrv(xi,lxi,bb,ilo,il2,mf2);
q = 0.0d0;
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;
x = min(bb,tb);
for ii = 1 : 2;
ss(ii) = 0.0d0;
dx = x - xi(left);
if( dx~=0.0d0 )
s = c(k,left);
flk = k;
im = fix(k - 1);
il = fix(im);
for i = 1 : il;
s = s.*dx./flk + c(im,left);
im = fix(im - 1);
flk = flk - 1.0d0;
end; i = fix(il+1);
ss(ii) = s.*dx;
end;
x = a;
end; ii = fix(2+1);
q = q +(ss(1)-ss(2));
end; left = fix(il2+1);
if( x1>x2 )
q = -q;
end;
pquad = 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;
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 DPPSL