function [b,knots] = knotinsertion(b,knots,bart)
n = size(b,2)-1;
k = length(knots)-n-1;
h = find(knots > bart,1)-1; % bart is in [u(h),u(h+1))
if (isempty(h))
  h = n+1;
end
newb = [b(:,1:h),b(:,h:n+1)];
for i = h-k+2:h
  val = knots(i+k-1) == knots(i);
  a = (1-val)*(bart-knots(i))/(knots(i+k-1)-knots(i)+val);
  newb(:,i) = (1-a)*b(:,i-1)+a*b(:,i);
end
b = newb;
knots = [knots(1:h),bart,knots(h+1:end)];
%!demo
%! b = rand(2,4);
%! t = linspace(0,1,51);
%! k = 3;
%! [S,knots] = BScurve(b,k,t);
%! plot(b(1,:),b(2,:),'o',S(1,:),S(2,:))
%! for i = 1:5
%!   bart = rand;
%!   [b,knots] = knotinsertion(b,knots,bart);
%! end
%! [S1,knots] = BScurve(b,k,t,knots);
%! hold on
%! plot(b(1,:),b(2,:),'x',S1(1,:),S1(2,:),'+')
%! legend('original control points','original B-spline curve',...
%! 'new control points (5 knot insertions)','new B-spline curve')