function spline = spmak(knots,coefs,sizec)
%SPMAK Put together a spline in B-form.
%
%   SPMAK(KNOTS,COEFS) puts together a spline from the knots and 
%   coefficients. COEFS should be d-by-n , with  d  the number of
%   components in each of the  n  spline coefficient supplied.  
%   The order of the spline is inferred as  k := length(KNOTS) - n .
%   Knot multiplicity is held to <= k , with the coefficients
%   corresponding to a B-spline with trivial support ignored.
%
%   SPMAK  will prompt for KNOTS and COEFS.
%
%   If KNOTS is a cell array of length  m , then COEFS is expected to
%   be an m- or (m+1)-dimensional array, and the corresponding m-variate
%   d-valued tensor product spline is returned, with  d  the first 
%   dimension of COEFS in case COEFS is (m+1)-dimensional, and 1 otherwise.
%   FEWER CHECKS ARE CARRIED OUT in this case.  
%
%   SPMAK(KNOTS,COEFS,SIXEC) uses SIZEC to specify the intended array
%   dimensions of COEFS, and may be needed for proper interpretation
%   of COEFS in case one or more of its trailing dimensions is a singleton
%   and thus COEFS appears to be of lower dimension.
%
%   For example, if the intent is to construct a 2-vector-valued bivariate
%   polynomial on the rectangle [-1 .. 1] x [0 .. 1], linear in the first
%   variable and constant in the second, say
%      coefs = zeros([2 2 1]); coefs(:,:,1) = [1 0;0 1];
%   then the straightforward
%      sp = spmak({[-1 -1 1 1],[0 1]},coefs);
%   will fail, producing a scalar-valued function of (illegal) order [2 0],
%   as will
%      sp = spmak({[-1 -1 1 1],[0 1]},coefs,size(coefs));
%   while proper use of that third argument, as in
%      sp = spmak({[-1 -1 1 1],[0 1]},coefs,[2 2 1]);
%   will succeed.
%
%   See also SPBRK, RSMAK, PPMAK, RPMAK, FNBRK.

%   Carl de Boor 9 dec 89
%   cb : February 4, 1991 (disallow the check for d>n)
%   cb : 20 oct 1994 (enforce knot multiplicity <= k )
%   cb :  9 may '95 (use .' instead of ')
%   cb : 15aug96 (check that basic interval is nontrivial)
%   cb : 18oct97 (change to structure, also cover multivariate spline)
%   cb : 07dec97 (reorder form for backward consistency, ensure knots is 1-row)
%   cb : 10may98 (standardize the help)
%   cb : 10apr00 (introduce optional third argument for multidimensional case)
%   cb : 23jun00 (improve the help)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc. 
%   $Revision: 1.10 $

if nargin==0;
   knots = input('Give the vector of knots  >');
   coefs = input('Give the array of B-spline coefficients  >');
end
if iscell(knots) % we are putting together a tensor-product spline

   m = length(knots); if nargin<3, sizec = size(coefs); end
   switch length(sizec)
      case m, sizec = [1, sizec]; coefs = reshape(coefs,sizec);
      case m+1
      otherwise
         error('COEFS must be a ([1+]length(KNOTS))-dimensional array.')
   end
   d = sizec(1); n = sizec(2:end);
   for j=m:-1:1
      k(j) = length(knots{j})-n(j);
      if k(j)<=0, error('There should be more knots than coefficients.'), end
      knots{j} = reshape(knots{j},1,k(j)+n(j));
   end

else            % we are putting together a univariate spline

   [d,n] = size(coefs);
   if isempty(coefs), error('The coefficient sequence is empty.'), end
   
   k = length(knots)-n;
   if k<=0, error('There should be more knots than coefficients.'), end
   diffk = diff(knots);
   if any(diffk<0)
      error('The knot sequence should be nondecreasing.')
   end
   if knots(1)==knots(end)
      error('The extreme knots should be different.')
   end
   
   index = find(knots(k+[1:n])-knots(1:n)>0); % throw out trivial B-splines
   if length(index)<n
      coefs = coefs(:,index); knots = knots([index n+[1:k]]); n = length(index);
   end
   knots = reshape(knots,1,n+k);
end

spline.form = 'B-';
spline.knots = knots;
spline.coefs = coefs;
spline.number = n;
spline.order = k;
spline.dim = d;
% spline = [11 d n coefs(:).' k knots(:).'];
