[t,a,n,k,ldc,c,xi,lxi,work]=bsppp(t,a,n,k,ldc,c,xi,lxi,work);

function [t,a,n,k,ldc,c,xi,lxi,work]=bsppp(t,a,n,k,ldc,c,xi,lxi,work);
%***BEGIN PROLOGUE  BSPPP
%***PURPOSE  Convert the B-representation of a B-spline to the piecewise
%            polynomial (PP) form.
%***LIBRARY   SLATEC
%***CATEGORY  E3, K6
%***TYPE      SINGLE PRECISION (BSPPP-S, DBSPPP-D)
%***KEYWORDS  B-SPLINE, PIECEWISE POLYNOMIAL
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%     Written by Carl de Boor and modified by D. E. Amos
%
%     Abstract
%         BSPPP is the BSPLPP routine of the reference.
%
%         BSPPP converts the B-representation (T,A,N,K) to the
%         piecewise polynomial (PP) form (C,XI,LXI,K) for use with
%         PPVAL.  Here XI(*), the break point array of length LXI, is
%         the knot array T(*) with multiplicities removed.  The columns
%         of the matrix C(I,J) contain the right Taylor derivatives
%         for the polynomial expansion about XI(J) for the intervals
%         XI(J) .LE. X .LE. XI(J+1), I=1,K, J=1,LXI.  Function PPVAL
%         makes this evaluation at a specified point X in
%         XI(1) .LE. X .LE. XI(LXI(1) .LE. X .LE. XI+1)
%
%     Description of Arguments
%         Input
%          T       - knot vector of length N+K
%          A       - B-spline coefficient vector of length N
%          N       - number of B-spline coefficients
%                    N = sum of knot multiplicities-K
%          K       - order of the B-spline, K .GE. 1
%          LDC     - leading dimension of C, LDC .GE. K
%
%         Output
%          C       - matrix of dimension at least (K,LXI) containing
%                    right derivatives at break points
%          XI      - XI break point vector of length LXI+1
%          LXI     - number of break points, LXI .LE. N-K+1
%          WORK    - work vector of length K*(N+3)
%
%     Error Conditions
%         Improper input is a fatal error
%
%***REFERENCES  Carl de Boor, Package for calculating with B-splines,
%                 SIAM Journal on Numerical Analysis 14, 3 (June 1977),
%                 pp. 441-472.
%***ROUTINES CALLED  BSPDR, BSPEV, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   800901  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   890831  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  BSPPP
%
persistent ileft inev nk ; 

if isempty(ileft), ileft=0; end;
if isempty(inev), inev=0; end;
if isempty(nk), nk=0; end;
%     DIMENSION T(N+K),XI(LXI+1),C(LDC,*)
%     HERE, * = THE FINAL VALUE OF THE OUTPUT PARAMETER LXI.
t_shape=size(t);t=reshape(t,1,[]);
a_shape=size(a);a=reshape(a,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
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  BSPPP
if( k<1 )
xermsg('SLATEC','BSPPP','K DOES NOT SATISFY K.GE.1',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
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( n<k ) ;
xermsg('SLATEC','BSPPP','N DOES NOT SATISFY N.GE.K',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
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','BSPPP','LDC DOES NOT SATISFY LDC.GE.K',2,1);
else;
k_orig=k;    [t,a,n,k,dumvar5,work]=bspdr(t,a,n,k,k,work);    k(dumvar5~=k_orig)=dumvar5(dumvar5~=k_orig);
lxi = 0;
xi(1) = t(k);
inev = 1;
nk = fix(n.*k + 1);
for ileft = k : n;
if( t(ileft+1)~=t(ileft) )
lxi = fix(lxi + 1);
xi(lxi+1) = t(ileft+1);
k_orig=k;    [t,dumvar2,n,k,dumvar5,xi(lxi),inev,c(sub2ind(size(c),1,lxi):end),dumvar9]=bspev(t,work(sub2ind(size(work),max(1,1)):end),n,k,k,xi(lxi),inev,c(sub2ind(size(c),1,lxi):end),work(sub2ind(size(work),max(nk,1)):end));    k(dumvar5~=k_orig)=dumvar5(dumvar5~=k_orig);   dumvar2i=find((work(sub2ind(size(work),max(1,1)):end))~=(dumvar2));dumvar9i=find((work(sub2ind(size(work),max(nk,1)):end))~=(dumvar9));   work(1-1+dumvar2i)=dumvar2(dumvar2i); work(nk-1+dumvar9i)=dumvar9(dumvar9i); 
end;
end; ileft = fix(n+1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
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;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
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 BSPVD