function [aprox_pi, sigma] = ukol2_reseni(itermax, delta)
% [aprox_pi, sigma] = ukol2_reseni(itermax)
%
% funkce vypocte metodou Monte Carlo s body generovanymi pomoci nahodne prochazky obsah plochy ohranicene grafem funkce sqrt(1 - x^2), 
% kladnym smerem osy x a osou y a smerodatnou odchylku sigma 
% pomoci zmineneho obsahu aproximuje cislo pi
%
% itermax ... pocet pokusu
% delta ... maximalni delka hodu

x = rand(); % generuje x-ovou souradnici bodu ve ctverci <0,1> x <0,1>
y = rand(); % generuje y-ovou souradnici bodu ve ctverci <0,1> x <0,1>

q = zeros(1,itermax);

for i = 1 : itermax
    dx = -delta + (2*delta)*rand(); % generuje x-ovou slozku hodu o velikosti maximalne delta 
    dy = -delta + (2*delta)*rand(); % generuje y-ovou slozku hodu o velikosti maximalne delta

    if ((x+dx) < 1.0 && (y+dy) < 1.0) && ((x+dx) > 0.0 && (y+dy) > 0.0) % hod miri do ctverce <0,1> x <0,1> 
        x = x+dx;
        y = y+dy;
    end

    if x^2 + y^2 < 1.0
        q(i) = 1; % bod lezi v plose ohranicene grafem funkce sqrt(1 - x^2), kladnym smerem osy x a osou y
    end

end

aprox_pi = 4.0 * sum(q) / double(itermax);

% vypis vysledku
disp(['Aproximace pi: ', num2str(aprox_pi, '%.15f')])
disp(['Absolutni chyba: ', num2str(abs(aprox_pi - pi), '%.15f')])
disp(['Relativni chyba: ', num2str(100*abs(aprox_pi - pi)/pi, '%.3f'),'%'])

mu = sum(q) / double(itermax);

sigma2 = 1/itermax * (q - mu) * (q - mu)';
sigma = 2 * sqrt(sigma2)/sqrt(itermax);

% konec
