function newtonuv_polynom(zobraz_vse,x,f)
%
% Sestroji Newtonuv interpolacni polynom 
% pro funkci zadanou tabulkou [x;f(x)]
%
% pro zobraz_vse=1 zobrazi vsechny dilci polynomy
% nejsou-li zadane vektory x,f => zadaji se mysi 

clc;
presah=10;
pocet_deleni=1000;
clear texty;
barva=['r','m','g','c','k'];
styl=['- ';'--';'-.';': '];

figure(1);
clf;
hold on;
zoom on;
set(figure(1),'Units','centimeters');
set(figure(1),'Position',[10 3 18 15]);


if ~exist('zobraz_vse')
  zobraz_vse=0;
end
if ~exist('x') || ~exist('f')
  [x,f]=zadej_body;
  clc;
else
  xmin=min(x);
  xmax=max(x);
  fmin=min(f);
  fmax=max(f);
  axis([xmin-presah/100*(xmax-xmin) xmax+presah/100*(xmax-xmin) ...
	fmin-presah/100*(fmax-fmin) fmax+presah/100*(fmax-fmin)]); 
  for i=1:max(size(x))
    plot(x(i),f(i),'ro');
  end  
end  

[sortx,~]=sort(x);
%f=f(indexy);
n=max(size(x));
deleni=(sortx(n)-sortx(1))/pocet_deleni;
zoom on;

disp('----------------------------------------------------------------');
fprintf('  Newtonuv interpolacni polynom pro funkci danou tabulkou\n');
disp(' ');
fprintf(' |    x(k)    |    f(k)    |\n');
fprintf(' ===========================\n');
for i=1:n
  fprintf(' | %9.4f  | %9.4f  |\n',x(i),f(i));
end  

A=zeros(n);
A(:,1)=f';
for k=2:n
  for i=k:n
    A(i,k)=(A(i,k-1)-A(i-1,k-1))/(x(i)-x(i-k+1));
  end
end

[xk,fk]=polynom(x,f,deleni,presah/2);

xmin=min(xk);
xmax=max(xk);
fmin=min(fk);
fmax=max(fk);
if zobraz_vse==1
  for i=1:n
  pomf=f(1:i);
  pomx=x(1:i);
  [~,fk]=polynom(pomx,pomf,deleni,presah/2);
  if fmin>min(fk)
    fmin=min(fk);
  end
  if fmax<max(fk)
    fmax=max(fk);
  end
 end
end
axis([xmin-presah/100*(xmax-xmin) xmax+presah/100*(xmax-xmin) ...
      fmin-presah/100*(fmax-fmin) fmax+presah/100*(fmax-fmin)]); 
disp(' ');
disp('Stiskni klavesu');
disp(' ');
pause;

if zobraz_vse==1
 set(figure(1),'Position',[7 3 28 15]);
 for i=1:n
  clf;
  hold on;
  axis([xmin-presah/100*(xmax-xmin) xmax+presah/100*(xmax-xmin) ...
        fmin-presah/100*(fmax-fmin) fmax+presah/100*(fmax-fmin)]); 
  for k=1:min(i,n-1)
   pomf=f(1:k);
   pomx=x(1:k);
   [xk,fk]=polynom(pomx,pomf,deleni,presah/2);
   color=strcat(barva(mod(k-1,5)+1),styl(mod(floor((k-1)/5),5)+1,:));
   if k==1
    plot([xk-(xmax-xmin)/20,xk+(xmax-xmin)/20 ],[fk,fk],color);
   else
    plot(xk,fk,color);
   end  
   popis=sprintf('polynom stupne %d            ', k-1);
   texty(k,:)=popis(1:25); 
  end
  if i==n
   popis=sprintf('Vysledny polynom            ');
   texty(i,:)=popis(1:25); 
   [xk,fk,koeficienty]=polynom(x,f,deleni,presah/2);  
   plot(xk,fk,'b-');
  end  
  
  legend(texty,'Location','bestoutside');
  
  bbb=plot(x,f,'ro');
  set(get(get(bbb,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
  pr=predpis_newton(x,f,diag(A),i);
  disp(pr);
  if i~=n
    pause;
  end
 end 
else
  [xk,fk,koeficienty]=polynom(x,f,deleni,presah/2);  
  plot(xk,fk,'b-');
  bbb=plot(x,f,'ro'); 
  set(get(get(bbb,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end

n=max(size(koeficienty));

text='N(x) = ';
if koeficienty(1)<0 
  text=[text,'- '];
end  
for i=n-1:-1:0
  c1=num2str(abs(koeficienty(n-i)));
  c2=num2str(i);
  if i>1
    if koeficienty(n-i+1)>=0
      text=[text, c1,' x^', c2, ' + '];
    else       
      text=[text, c1,' x^', c2, ' - '];
    end   
  else
    if i==1
      if koeficienty(n)>=0
        text=[text, c1,' x',' + '];
      else       
        text=[text, c1,' x',' - '];
      end   
    else  
      text=[text, c1];    
    end
  end 
end
disp(' ');
disp(text);
disp(' ');

% koeficienty = polyfit(x,f,max(size(x))-1)
