function di2 = demo_int2D()
%
% Integrace 2D funkce jednoduchymi formulemi
% Integrujeme na intervalu <0,1> x <0,1>
%
% - 2D obdelnikove pravidlo (midpoint)
% - 1D obdelnikove pravidlo postupne po smerech
%      (na teto pravouhle siti je to zcela ekvivalentni pouziti 2D obd)
% - 1D lichobeznikove pravidlo postupne po smerech
%      (srovnani naivni implementace po bunkach a slozeneho vzorce)
% 
% POZN.: Da se dale (mnohonasobne) urychlit, napr. operacemi s poli
%        misto s prvky, t.j. sum(pole(:)) misto cyklu pres pole(i), atd.
%


test = 2;

%% integrovana funkce a presne reseni (jen pro vypocet chyb)
if(test==1)
    f = inline('abs(sin(pi^2.*((x-0.5).^2+(y-0.5).^2)))','x','y');
    presny = 0.655172456; % by Matlab. By Maple: 0.655172424
elseif(test==2)
    f = inline('x+y','x','y');
    presny = 1.0;   
end;

%% integracni oblast
xmin = 0.0;
xmax = 1.0;
ymin = 0.0;
ymax = 1.0;

%% prostorovy krok (stejny v obou smerech)
h = 5.0e-2;
%h = 1.0e-2;
%h = 5.0e-3;

%% 2D Obdelnikove pravidlo (midpoint)
%% Do xArr,yArr ulozime rovnou polohy stredu bunek
xArr = (xmin+h/2:h:xmax-h/2);
yArr = (ymin+h/2:h:ymax-h/2);

integral=0;

tic % zacni merit cas
for xi=xArr       % cyklus pres x
    for yi=yArr   % cyklus pres y
        % obdelnikove pravidlo v pulce podintervalu (midpoint rule)
        integral = integral + f(xi,yi);
    end;
end;
integral = integral*h^2;
cas = toc; % skonci mereni casu

err_rel = abs(integral-presny)/max(abs(presny),1e-30);
disp(['2D obdelnikove pravidlo (midpoint):']);
%sprintf('  Vypocteny integral je %.12f',integral)
disp(['  Vypocteny integral je ' num2str(integral) ...
      ',  rel.ch. je ' num2str(err_rel) ]);
disp(['  Vypocet trval ' num2str(cas) ' sekund']);




%% 1D integrace ve smeru y, pak ve smeru x, obdelnikove pravidlo
%% V xArr,yArr mame stale polohy stredu bunek
inty     = zeros(size(xArr,2));  % pole spocitanych integralu dy
i        = 0;                    % budeme potrebovat index do pole
integral = 0;

tic % zacni merit cas
%% spocti integraly po stredech bunek ve smeru y
for xi = xArr
    i       = i+1;
    % obdelnikove pravidlo ve smeru y
    % ulozime pro kazde x integral v dy
    inty(i) = 0;
    for yi=yArr
       inty(i) = inty(i) + f(xi,yi);
    end;
end;
%% zintegruj integraly pres smer x
% obdelnikove pravidlo ve smeru x
for i=1:size(inty,2)
    integral = integral + inty(i);
end
integral = integral*h^2;
cas = toc; % skonci mereni casu

err_rel = abs(integral-presny)/max(abs(presny),1e-30);
disp(['1D Obdelnikove pravidlo (midpoint) po smerech:']);
%sprintf('  Vypocteny integral je %.12f',integral)
disp(['  Vypocteny integral je ' num2str(integral) ...
      ',  rel.ch. je ' num2str(err_rel) ]);
disp(['  Vypocet trval ' num2str(cas) ' sekund']);





%% 1D integrace ve smeru y, pak ve smeru x, lichobeznikove pravidlo
%% Tentokrat do xArr,yArr ulozime polohy uzlu
xArr = (xmin:h:xmax);   
yArr = (ymin:h:ymax); 

inty     = zeros(size(xArr,2)); % pole spocitanych integralu dy
i        = 0;                   % budeme potrebovat index do pole
integral = 0;

tic % zacni merit cas
%% spocti integraly po hranach ve smeru y
for xi = xArr
    i       = i+1;
    % lichobeznikove pravidlo ve smeru y
    % (slo by urychlit pouzitim slozeneho vzorce)
    % ulozime pro kazde x integral v dy
    inty(i) = 0;
    for yi=yArr(1:size(yArr,2)-1)
       inty(i)=inty(i) + 0.5*( f(xi,yi) + f(xi,yi+h) );
    end;
end;
%% zintegruj integraly pres smer x
% lichobeznikove pravidlo ve smeru x
% (slo by urychlit pouzitim slozeneho vzorce)
for i=1:size(inty,2)-1
    integral = integral + 0.5*( inty(i) + inty(i+1) );
end
integral = integral*h^2;
cas = toc; % skonci mereni casu

err_rel = abs(integral-presny)/max(abs(presny),1e-30);
disp(['1D Lichobeznikove pravidlo po smerech - scitame po bunkach:']);
%sprintf('  Vypocteny integral je %.12f',integral)
disp(['  Vypocteny integral je ' num2str(integral) ...
      ',  rel.ch. je ' num2str(err_rel) ]);
disp(['  Vypocet trval ' num2str(cas) ' sekund']);





%% 1D integrace ve smeru y, pak ve smeru x, lichobeznikove pravidlo
%% V xArr,yArr mame stale polohy uzlu

inty     = zeros(size(xArr,2)); % pole spocitanych integralu dy
i        = 0;                   % budeme potrebovat index do pole

tic % zacni merit cas
%% spocti integraly po hranach ve smeru y
for xi = xArr
    i       = i+1;
    % lichobeznikove pravidlo ve smeru y, slozeny vzorec
    % ulozime pro kazde x integral v dy
    yimax = size(yArr,2);
    inty(i) = 0.5*f(xi,yArr(1));
    for yi=yArr(2:yimax-1)
       inty(i) = inty(i) + f(xi,yi);
    end;
    inty(i) = inty(i) + 0.5*f(xi,yArr(yimax));
end;
%% zintegruj integraly pres smer x
% lichobeznikove pravidlo ve smeru x, slozeny vzorec
imax = size(inty,2);
integral = 0.5*inty(1);
for i=2:imax-1
    integral = integral + inty(i);
end
integral = integral + 0.5*inty(imax);
integral = integral*h^2;
cas = toc; % skonci mereni casu

err_rel = abs(integral-presny)/max(abs(presny),1e-30);
disp(['1D Lichobeznikove pravidlo po smerech - slozeny vzorec:']);
%sprintf('  Vypocteny integral je %.12f',integral)
disp(['  Vypocteny integral je ' num2str(integral) ...
      ',  rel.ch. je ' num2str(err_rel) ]);
disp(['  Vypocet trval ' num2str(cas) ' sekund']);

end

