function v = spval(sp,x)
% SPVAL	Evaluate spline.
%
%          v = spval(sp,x)
%
%  returns the vector  sp(x) .

% C. de Boor / latest change: May 23, 1989
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.


%  if necessary, sort  x 
xs=x(:)';lx=length(xs);tosort=0;
if (~isempty(find(diff(xs)<=0))),
   tosort=1;[xs,indx]=sort(xs);
end

%  take apart  spline
[t,a,n,k,d]=spbrk(sp);

%  and augment the knot sequence so that first and last knot each have
%  multiplicity  k .

index=find(diff(t)>0);addl=k-index(1);addr=index(length(index))-n;
if (addl>0|addr>0),
   t=[t(1)*ones(1,addl) t(:)' t(n+k)*ones(1,addr)];
   a=[zeros(d,addl) a zeros(d,addr)];
   n=n+addl+addr;
end

% for each data point, compute its knot interval
index=max([sorted(t(1:n),xs);k*ones(1,lx)]);

% set up and carry out in lockstep the first spline evaluation algorithm:

b=a(:,index);
if (k>1)
   ones(d,1)*index;dindex=ans(:);
   tx=ones(d*lx,1)*[2-k:k-1]+dindex*ones(1,2*(k-1));tx(:)=t(tx);
   ones(d,1)*xs;tx=tx-ans(:)*ones(1,2*(k-1));
   a=a(:);d*ones(d,1)*index+[1-d:0]'*ones(1,lx);dindex(:)=ans(:);
   b=d*ones(d*lx,1)*[1-k:0]+dindex*ones(1,k);b(:)=a(b);

   % the following loop is taken from  sprpp .

   for r=1:k-1;
      for i=1:k-r;
         b(:,i) =(tx(:,i+k-1).*b(:,i)-tx(:,i+r-1).*b(:,i+1))./...
                  (tx(:,i+k-1)-tx(:,i+r-1));
      end
   end

end

v=zeros(d,lx);v(:)=b(:,1);
if tosort>0,[junk,indx]=sort(indx); v=v(:,indx);end

% finally, zero out all values for points outside the knot interval
index=find(x<t(1)|x>t(n+1));
if (~isempty(index)),
   v(:,index)=zeros(d,length(index)); end
