function [o1,o2,o3,o4,o5] = sdiff(i1,i2,i3,i4,i5,i6);
% function [a,u] = sdiff(a,u[,n]); (function or curve)
%
% function [a,u,v] = sdiff(a,u,v,n); (tensor product function)
%
% function [a1,a2,a3,u,v] = sdiff(a1,a2,a3,u,v,n); (tensor product surface)
%
% differentiation
%
% input:  a,a1,a2,a3 ... coefficient matrix
%         u,v        ... knot vector
%         n          ... order of differentiation, optional [default n = 1]
%
% output: a,a1,a2,a3 ... new coefficient matrix
%         u,v        ... new knot vector

% written by A. Bernzott and K. Holthaus in 1993
%                                                                            
% changed by Juergen Koch, June 1994

if (nargin == 2) i3 = 1; end

if nargin == 2 | nargin == 3   % spline function or spline curve

  for i = 1:i3

    [m,la] = size(i1); lu = length(i2); d = lu - la - 1;

    if (d == 0) disp('sdiff: n > d !'); end

    % find indices j with 2 <= j <= la and u_j < u_{j+d} 
    J = [];
    for j = 2:la
      if i2(j) < i2(j+d) J = [J j]; end
    end

    k = J(length(J)); 			% k = J_last
    i1 = d*(i1(:,J) - i1(:,J-1))./(ones(m,1)*(i2(J+d) - i2(J)));
    i2 = i2([J k+1:k+d]);
  end

  o1 = i1; o2 = i2;
  
elseif nargin == 4              % tensor product spline function

  [o1,o2] = sdiff(i1',i2,i4(1));
  [o1,o3] = sdiff(o1',i3,i4(2));
  
elseif nargin == 6              % tensor product spline surface
    
  [o1,o4,o5] = sdiff(i1,i4,i5,i6);
  [o2,o4,o5] = sdiff(i2,i4,i5,i6);
  [o3,o4,o5] = sdiff(i3,i4,i5,i6);
    
else
  disp('sdiff: nargin <> 2,3,4,6 !');
end