function integral_richardson(vzorec, f, a, b, N, pocet_korekci)
%
% Vypocte hodnotu urciteho integralu zadane funkce f na <a,b> 
% s kroky h, h/2, h/4, ... h/2^pocet_korekci (kde h=(b-a)/N)
% pomoci Newton-Cotesovych kvadraturnich vzorcu

if mod(N,2)
  N=N+1;
end
h=(b-a)/N;
clc;
disp('----------------------------------------------------------------');
fprintf('  Vypocte hodnotu urciteho integralu zadane funkce \n');
fprintf('  f=%s na <%f,%f> s kroky\n',f,a,b);
kroky=sprintf('  h=[');
for j=1:pocet_korekci+1
  kroky=[kroky,sprintf('%g',h/(2^(j-1)))];
  if j<=pocet_korekci
    kroky=[kroky,', '];
  else 
    kroky=[kroky,']'];
  end
end
disp(kroky);
disp(' ');

if strncmp(vzorec,'I_O',3)
  fprintf('  Pro vypocet se pouzije slozene obdelnikove pravidlo I_O.\n'); 
  fprintf('  Ke zpresneni se pouzije Richardsonova extrapolace.\n');
else 
  if strncmp(vzorec,'I_L',3)
    fprintf('  Pro vypocet se pouzije slozene lichobeznikove pravidlo I_L.\n'); 
    fprintf('  Ke zpresneni se pouzije Richardsonova extrapolace.\n');
  else
    if strncmp(vzorec,'I_S',3)
      fprintf('  Pro vypocet se pouzije slozene Simpsonovo pravidlo I_S.\n'); 
      fprintf('  Ke zpresneni se pouzije Richardsonova extrapolace.\n');
    else
      disp(' ');
      fprintf('  Neni zadan zadny vzorec ("I_O", "I_L", "I_S") !!!\n'); 
      disp('----------------------------------------------------------------');
      return;       
      % break;
    end
  end
end
disp('----------------------------------------------------------------');


for i=1:pocet_korekci+1
  [int_fce,I_O,I_L,I_S]=integral_vypocet_newton(f, a, b, N*(2^(i-1)), 0);
  if strncmp(vzorec,'I_O',3) 
    I(i)=I_O;
    m=2;%mm=0;
  else 
    if strncmp(vzorec,'I_L',3)
      I(i)=I_L;
      m=2;%mm=0;
    else 
      if strncmp(vzorec,'I_S',3)
        I(i)=I_S;
        m=4;%mm=2;
      end
    end
  end
end

TAB=zeros(pocet_korekci+1);
TAB(:,1)=I';
for i=2:pocet_korekci+1
  for j=i:pocet_korekci+1 
    TAB(j,i)=TAB(j,i-1)-(TAB(j,i-1)-TAB(j-1,i-1))/(1-2^((i-1)*m));
%    TAB(j,i)=TAB(j,i-1)-(TAB(j,i-1)-TAB(j-1,i-1))/(1-2^(mm+(i-1)*2));
  end
end

disp(' ');
tiskni=sprintf(' |     h     |  %s(f,a,b,h) | ',vzorec(1:3));
for i=1:pocet_korekci
  tiskni=[tiskni,sprintf('  %d.korekce   | ',i)];
end
disp(tiskni);
pom=' ';
for i=1:length(tiskni)-2
  pom=[pom,'='];
end
disp(pom)
for i=1:pocet_korekci+1
  tiskni=sprintf(' |  %5.5f  ',h/(2^(i-1)));
  for j=1:i
    tiskni=[tiskni,sprintf('| %12.9f  ',TAB(i,j))];
  end
  tiskni=[tiskni,sprintf('|')];
  disp(tiskni);
end
disp(' ');

if ~isa(int_fce,'sym')
  fprintf('  Presna hodnota integralu funkce f na <a,b> je  %15.12f\n',int_fce); 
end

disp(' ');





