function pp=bspline(t,window)
% BSPLINE	Generate the B-spline for a given knot sequence.
%
%                 [pp=] bspline(t[,window])
%
%  plots the B-spline with knot sequence  t  as well as the polynomial
%  pieces of which it is composed, in the window specified.
%  Also returns its pp form if requested.

% C. de Boor / latest change: May 22, 1989
% C. de Boor / latest change: 11 May 1992 (correct misprint)
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.

%  Generate the spline description

k=length(t)-1;
if (k>1)
   adds=ones(1,k-1);
   tt=[adds*t(1) t(:)' adds*t(k+1)];j=k+1;
   a=[adds*0 1 adds*0]; 

   %  From this, generate the pp description

   inter=find( diff(tt)>0 );l=length(inter);
   tx=ones(l,1)*[2-k:k-1]+inter'*ones(1,2*(k-1));tx(:)=tt(tx);
   tx=tx-tt(inter)'*ones(1,2*(k-1));
   b=ones(l,1)*[1-k:0]+inter'*ones(1,k);b(:)=a(b);
   c=sprpp(tx,b);x=[tt(inter) tt(2*k)];
else
  l=1;x=t;c=1;
end

if (nargout>0), pp=ppmak(x,c,1); end

%  Now generate a mesh ...

step=100;xx=x(1)+[-10:step+10]*(x(l+1)-x(1))/step;nstep=length(xx);

%  and generate and plot the polynomial pieces

if (nargin>1), subplot(window); else, clf; end
xxx=[xx(1) xx(nstep)]; yyy=[-1 2];axis([xxx yyy]);hold on;
for j=1:k+1; plot([t(j) t(j)],yyy);end
bspl=zeros(1,step+21);find(xx>=x(1));jh=ans(1);
jsmax=3;js=jsmax;
for j=1:l;js=js+1;if (js>jsmax), js=1;end
   jl=jh; find(xx>=x(j+1));jh=ans(1);
   pval=polyval(c(j,:),xx-x(j));bspl(jl:jh)=pval(jl:jh);
   if (js==1),     plot(xx,pval,':');
   elseif (js==2), plot(xx,pval,'-.');
   else,           plot(xx,pval,'--');

   end
end
plot(xx,bspl);hold off;
