function [c,tau] = knot_insertion(c,tau,t) 
% KNOT_INSERTION: knot insertion for splines
% [c,tau] = knot_insertion(c,tau,t) 
% change control points c(k,:) and knots tau
% by adding knots t(j) in interior knot intervals with
% multiplicity at most equal to the degree (=> n>0)

% convert to column vectors (if necessary)
tau = tau(:); t = t(:);

% number of control points, dimension, and degree 
[m,dim] = size(c); n = length(tau) - m - 1;

% add knots successively 
for tj = t'
   % relevant indices
   ind = find(tau(2:end-n-1)<tj & tau(n+2:end-1)>tj)+1;
   if length(ind)>0
      % convex combinations
      gamma = repmat((tj-tau(ind))./(tau(ind+n)-tau(ind)),1,dim);
     c = [c(1:min(ind-1),:); 
          (1-gamma).*c(ind-1,:)+gamma.*c(ind,:); 
          c(max(ind):end,:)];
     tau = sort([tau; tj]);
   else 
      display(['knot ' num2str(tj) ' not admissible'])
   end
end

