function [yi, w] = fh(x, y, xi, d)
%------------------------------------------------------------------
% This function computes the F-H's interpolant for equidistant nodes x.
%
% INPUTS:
% (x,y) - two vectors: the first one contains the
%         nodes and the second one contains the function value 
%         at the nodes
% xi - vector of points, where we want to approximate function
% d - polynomial degree of F-H's interpolant (0 <= d <= n-1)
%
% OUTPUTS:
% yi - vector of function approximated values in xi
% w - vector of F-H's interpolant weights
%----------------------------------------------------------------
% Authors: 
%      Stefano De Marchi and Cinzia Bandiziol
%      Department of Pure and Applied Math.
%      University of Padova
%   
%           
% Version: May 2015		
%-----------------------------------------------------------------		

n=length(x);
m=length(xi);


% Compute matrix Z
Z=zeros(n,n-d);
vect_diagonals=ones(n-d,1);
t=(2:2:n-d);
vect_diagonals(t)=-1;
index_diagonals=(-d:0);
mat_diagonals=repmat(vect_diagonals,1,size(index_diagonals,2));
Z=spdiags(mat_diagonals,index_diagonals,Z);



A2=repmat(x',n,1);%nxn matrix

% Get distances between nodes and interpolation points
P1=A2-A2';

% Find all of the elements where the interpolation point is on a node
[a,b]=find(P1==0);
if(size(a,1)>n)
fprintf('<< Wrong input: points in x must be distinct >>\n');
      yi(i) = NaN;
return;
end

% Compute matrix G
G=ones(n-d,n);
P2=eye(n);
P=P1+P2;
Q=ones(n,n)./P;
for k=1:n-d
	G(k,:)=prod(Q(k:k+d,:),1);
end

% Matrix of weights
w=sum(Z'.*G,1)';


% Compute matrix A, Aij=xi(i)-x(j)
A1=repmat(xi,1,n);% mxn matrix
A3=repmat(x',m,1);% mxn matrix
A=A1-A3;
B=ones(m,n)./A;

% Evaluate denominator of F-H's interpolant at xi
den=B*w;

% Evaluate numerator of F-H's interpolant at xi
num=B*(w.*y);

% Evaluate F-H's interpolant at xi
yi=num./den;

[a,b]=find(A==0);

% Replace 0s with the given exact values. 
yi(a)=y(b);
