%Barycentric Lagrange interpolation.
%VecX, VecY: interpolation data. x0: observation point.
function value1=barycentric(VecX,VecY,x0)
%input checking
n=length(VecX);
if length(VecY) ~=n; display('Interpolation data inconsistent'); return;end
% calculating the weights
V=bsxfun(@minus,VecX,VecX.');
V(1:n+1:end) = 1;
VV = exp(sum(log(abs(V))));
w = 1./(prod(sign(V)).*VV);
% calculating the interpolant
num=sum(w.*VecY./(x0-VecX));
den=sum(w./(x0-VecX));
value1=num/den;
end