function w = quasi_interpolation(xi,t)
% QUASI_INTERPOLATION: weights for quasi-interpolation
% w = quasi_interpolation(xi,t)
% functional sum_{i=1}^{n+1} w(k,i)f(t(k,i)) for knots xi

% convert to column vector (if necessary)
xi = xi(:);

% degree 
n = length(xi)-size(t,1)-1; 

% weights via Marsden's identity
if n==0
   w = ones(size(t));
else
   for k=1:size(t,1)
      tk = repmat(t(k,:),n+1,1);
      A = (tk'-tk).^n;
      x = repmat(xi(k+1:k+n),1,n+1)-tk(1:n,:);
      if n>1
         x = prod(x);
      end 
      w(k,:) = x / A;
   end
end

