function y=LagrangeInterp(xd,yd,x)
% Solmu 2/2019, Heikki Apiola
% Suoritetaan Lagrangen interpolaatio datalla (xd,yd)
% ja lasketaan polynomin arvot vektorilla x.
% Esim: 
% xd=-3.5:3.5;yd=sin(xd);
% x=linspace(-4.5,4.5,200); % Laskentapisteet
% y=LagrangeInterp(xd,yd,x);
% plot(xd,yd,'o',x,y);grid on
% Esimerkin tapauksessa interpolointi antaa varsin
% hyvän approksimaation; sin tällä välillä muistuttaa
% paritonta astetta olevaa polynomia.
N=length(xd);
y=zeros(size(x));
for i=1:N
    y=y+yd(i)*L(i,xd,x);
end    
% Alifunktio:
function  y  = L(k,xdata,x)
% xdata-vektoriin liittyvä Lagrangen kantafunktio L_k
% laskettuna vektorilla x.
n=length(xdata);
y=ones(size(x));
for j=[1:k-1 k+1:n]
    y=y.*((x-xdata(j))./(xdata(k)-xdata(j)));
end