function spnew = sprfn(sp,addknots)
%SPRFN Insert additional knots into B-form of a spline.
%
%   SPNEW = SPRFN(SP,ADDKNOTS)  insert the specified ADDKNOTS into the B-form 
%   of the spline contained in SP.
%
%   No knot multiplicity will be increased beyond the order of the spline.
%
%   If SP describes an m-variate spline, then ADDKNOTS is expected to be
%   a cell array with m entries, any of which may be empty if no refinement
%   in the corresponding knot sequence is wanted.
%
%   See also FNRFN, PPRFN.

%   Carl de Boor 10 dec 95; 14 feb 96; 26aug96(error when d>1); 21sep96
%   07dec96 (check against ADDKNOTS being empty)
%   05oct97 (replace use of ANS, to help compilation)
%   01nov97 (also handle tensor products)
%   07dec97 (permit ADDKNOTS to be 1-column matrix)
%   10may98 (standardize the help)
%   06may99 (correct failure to reshape coefs properly)
%   06sep99 (permit ADDKNOTS to be a matrix)
%   09apr00 (deal with fact that matlab drops trailing singleton dims)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.12 $

if iscell(addknots)    % we must be dealing with a multivariate spline

   [t,a,n,k,d] = spbrk(sp);

   if ~iscell(t)
    error('If ADDKNOTS is a cell-array, then SP must be a multivariate spline.')
   end
   
   m = length(n);
   coefs = a; sizec = [d,n];
   for i=m:-1:1   % carry out coordinatewise knot refinement
      if ~isempty(addknots{i})
         spi = sprfn1(spmak(t{i},reshape(coefs,prod(sizec(1:m)),sizec(m+1))),...
                                     addknots{i});
         t{i} = spi.knots; sizec(m+1) = spi.number; 
         coefs = reshape(spi.coefs, sizec);
      end
      if m>1
         coefs = permute(coefs,[1,m+1,2:m]); sizec(2:m+1) = sizec([m+1,2:m]);
      end
 
   end
   % At this point, COEFS contains the tensor-product B-spline coefficients;
   % also, the various knot sequences will have been updated. 
   % It remains to return information:
   spnew = spmak(t, coefs,sizec);

else             % univariate spline refinement
   spnew = sprfn1(sp,addknots);
end

function spnew = sprfn1(sp,addknots)
%SPRFN1 Insert additional knots into B-form of a univariate spline.

if isempty(addknots), spnew = sp; return, end
addknots = sort(addknots(:).'); ladd = length(addknots); 

[t,a,n,k,d] = spbrk(sp);

% retain only distinct points, but record their input multiplicity
index = [1 find(diff(addknots)>0)+1];
inmults = diff([index ladd+1]); sortedadds = addknots;
addknots = addknots(index); ladd = length(index);

% compute the current multiplicity of ADDKNOTS in T.
indexr = sorted(t, addknots);
temp = n+k - sorted(-t,-addknots); mults = indexr -  temp(ladd:-1:1);
% ... then reduce INMULTS to make certain that output has knots of
%     multiplicity at most  k .
excess = subplus(inmults+mults-k);
if any(excess>0)
   inmults = inmults - excess;
   index = find(inmults>0);
   if isempty(index)
      fprintf('WARNING (from SPRFN): ')
      fprintf('All additional knots occur already to full multiplicity.\n')
      spnew = sp; return
   end
   fprintf('WARNING (from SPRFN): ')
   fprintf('Insertion of some knot(s) would have caused excess multiplicity.\n')
   if length(index)<ladd
      addknots = addknots(index); mults = mults(index);
      inmults = inmults(index); ladd = length(addknots);
   end
end

% if the endknot multiplicity is to be increased and/or there are knots
% outside the current basic interval, do it now
lamin=1;
index = find(addknots<t(1));
if ~isempty(index) % there are knots to be put to the left
   totals = cumsum(inmults(index));
   t = [sortedadds(1:totals) t];
   a = [zeros(d,totals) a]; n = n + totals; indexr = indexr + totals;
   lamin = 1+length(index);
elseif addknots(1)==t(1)
   t = t([ones(1,inmults(1)) 1:(n+k)]);
   a = [zeros(d,inmults(1)) a]; n = n+inmults(1); indexr = indexr + inmults(1);
   lamin = 2;
end
lamax = ladd;
index = find(addknots>t(n+k));
if ~isempty(index) % there are knots to be put to the right
   totals = cumsum(inmults(index));
   t = [t sortedadds(length(sortedadds)+((1-totals):0))];
   a = [a zeros(d,totals)]; n = n + totals;
   lamax = ladd-length(index);
elseif addknots(ladd)==t(n+k)
   t = t([1:(n+k) (n+k)*ones(1,inmults(ladd))]);
   a = [a zeros(d,inmults(ladd))]; n = n+inmults(ladd);
   lamax = ladd-1;
end

% Increase endknot multiplicity to  k , to avoid difficulties.
% (Taken from SPVAL)
index = find(diff(t)>0); addl = k-index(1); addr = index(end)-n;
if ( addl>0 | addr>0 )
   t = t([ones(1,addl) 1:(n+k) (n+k)*ones(1,addr)]);
   a = [zeros(d,addl) a zeros(d,addr)];
   n = n+addl+addr; indexr = indexr+addl;
end

% Ready for knot insertion, one at a time.
for la=lamin:lamax
   for mm=1:inmults(la)
      newa = a(:,[1:indexr(la)-k+1 indexr(la)-k+1:n]);
      newt = [t((1:indexr(la))) addknots(la) t((indexr(la)+1:n+k))];

      for j=(indexr(la)-k+2):(indexr(la)-mults(la))
         newa(:,j) = ...
             (a(:,j-1)*(t(j+k-1)-addknots(la)) + a(:,j)*(addknots(la)-t(j)))/...
             (          t(j+k-1)                                     -t(j) );
      end
      t = newt; a = newa; n = n+1; indexr = indexr+1; mults(la) = mults(la)+1;
   end
end

if addl>0|addr>0 % remove again those additional end knots and coefficients
   a(:,[1:addl n+((1-addr):0)]) = [];
   t([1:addl n+k+((1-addr):0)]) = []; n = n - addl-addr;
end

spnew = spmak(t,a);
