function [int_fce,I_O,I_L,I_S]=integral_vypocet_newton(f, a, b, N, zobraz)
%
% Vypocte numericky urcity integral symbolicky zadane funkce f 
% na intervalu <a,b> s poctem deleni N pomoci Newton-Cotesovych vzorcu
% I_O ... obdelnikove pravidlo
% I_L ... lichobeznikove pravidlo
% I_S ... Sipmsonovo pravidlo


if mod(N,2)
  N=N+1;
end
krok=(b-a)/N;
deleni=krok/100;

if zobraz==1
 clc;
 disp('------------------------------------------------------------------');
 fprintf('  Numericky vypocet urciteho integralu funkce f(x)=%s\n',f);
 fprintf('  na intervalu <%f,%f> s poctem deleni N=%d\n',a,b,N);
 disp('------------------------------------------------------------------');
end

syms 'x';
%int_fce=eval(int(eval(f),a,b));
int_fce=double(int(eval(f),a,b));

I_O=0;
I_L=0;
I_S=0;

% obdelnikove pravidlo
for i=1:N
  x=a+(i-0.5)*krok;
  ff(i)=eval(f); 
  I_O=I_O+ff(i)*krok;  
end

% lichobeznikove pravidlo
for i=1:N
  x=a+(i-1)*krok;
  levy=eval(f);
  x=a+i*krok;
  pravy=eval(f);
  ff(i)=(pravy+levy)/2;
  I_L=I_L+ff(i)*krok;  
end

% Simpsonovo pravidlo
for i=1:N/2
  x=a+2*(i-1)*krok;
  levy=eval(f);
  x=a+(2*i-1)*krok;
  stred=eval(f);
  x=a+2*i*krok;
  pravy=eval(f);
  ff(i)=(pravy+4*stred+levy)/3;
  I_S=I_S+ff(i)*krok;  
end


if zobraz==1
 figure(1);
 clf;
 hold on;
 zoom on;
 xx=a:deleni:b;
 i=1;
 for x=a:deleni:b
   fcef(i)=eval(f);
   i=i+1;
 end

 max_y=max([max(fcef),0]);
 min_y=min([min(fcef),0]);
 sy=max_y-min_y;
 sx=b-a;

 axis([a-sx/10 b+sx/10 min_y-.1*sy max_y+.1*sy]);
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 texty(1,:)=sprintf('zadana funkce                   ');
 %legend(texty,'Location','NorthWest');
 legend(texty,'Location','Best');
 pause;
 clf;
 hold on;
 zoom on;
 axis([a-sx/10 b+sx/10 min_y-.1*sy max_y+.1*sy]);

 disp(' ');
 fprintf(' Presna hodnota integralu je  %f \n',int_fce);
 disp(' ');
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 barva=[0.75 0.75 0.85];
 fill([xx,b,a],[fcef,0,0],barva)
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 legend(texty,'Location','Best');
 title('Presna hodnota integralu');
 pause;

 disp(sprintf([' Priblizna hodnota integralu',...
                '                                  chyba ']));
 fprintf(...
 '     - pomoci obdelnikoveho pravidla       %f          %f\n',I_O,I_O-int_fce);
 clf;
 hold on;
 zoom on;
 axis([a-sx/10 b+sx/10 min_y-.1*sy max_y+.1*sy]);
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 barva=[0.75 0.85 0.75];
 for i=1:N
   levy=a+(i-1)*krok;
   pravy=a+i*krok;
   x=(pravy+levy)/2;
   pom=eval(f); 
   fill([levy,levy,pravy,pravy],[0 pom pom 0], barva);
 end
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 legend(texty,'Location','Best');
 title('Obdelnikove pravidlo');
 pause;              

 fprintf(...
 '     - pomoci lichobeznikoveho pravidla    %f          %f\n',I_L,I_L-int_fce);
 clf;
 hold on;
 zoom on;
 axis([a-sx/10 b+sx/10 min_y-.1*sy max_y+.1*sy]);
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 barva=[0.85 0.75 0.75];
 for i=1:N
   levy=a+(i-1)*krok;
   pravy=a+i*krok;
   x=levy;
   pom1=eval(f); 
   x=pravy;
   pom2=eval(f); 
   fill([levy,levy,pravy,pravy],[0 pom1 pom2 0], barva);
 end
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 legend(texty,'Location','Best');
 title('Lichobeznikove pravidlo');
 pause;              

 fprintf(...
 '     - pomoci Simpsonova pravidla          %f          %f\n',I_S,I_S-int_fce);
 clf;
 hold on;
 zoom on;
 axis([a-sx/10 b+sx/10 min_y-.1*sy max_y+.1*sy]);
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 barva=[0.85 0.65 0.75];
 for i=1:N/2
   x=a+2*(i-1)*krok;
   xlevy=x;
   levy=eval(f);
   x=a+(2*i-1)*krok;
   xstred=x;
   stred=eval(f);
   x=a+2*i*krok;
   xpravy=x;
   pravy=eval(f);
   [xparabola,parabola,~]=polynom([xlevy,xstred,xpravy],[levy,stred,pravy],deleni,0);
   fill([xparabola,xpravy,xlevy],[parabola,0,0], barva);
 end
 o=plot(xx,fcef,'r-');
 set(o,'LineWidth',2);
 legend(texty,'Location','Best');
 title('Simpsonovo pravidlo');
 disp(' ');
 pause; 
 close(1);
end









