function p = spline_evaluation(c,xi,x) 
% SPLINE_EVALUATION: evaluation of a spline function
% p = spline_evaluation(c,xi,x) 
% values p(j) at x(j) of a spline function  
% with coefficients c and knots xi

% change to column vectors (if necessary)
c = c(:); xi = xi(:); x = x(:);

% degree
n = length(xi)-length(c)-1;

% add points and knots to avoid exceptional cases
c = [zeros(n+1,1); c; zeros(n+1,1)];
xi_max = (2*norm([xi; x],inf)+1)*ones(n+1,1);
xi = [-xi_max; xi; xi_max];

for j = 1:length(x)
   % interval: xi(kj) <= x(j) < xi(kj+1)
   kj = min(find(xi>x(j)))-1;
   % relevant control points
   pxj = c(kj-n:kj);
   % triangular scheme
   for m = 1:n 
      ind = [kj-n+m:kj];
      gamma = ((x(j) - xi(ind)) ./ (xi(ind+n+1-m) - xi(ind)));
      pxj = gamma.*pxj(2:end) + (1-gamma).*pxj(1:end-1);
   end 
   p(j) = pxj;
end
p = p(:); 
