%% Para ilustar o método (BEM SIMPLES):
clear
N = 1000;
Xs = rand(N,1);
Ys = rand(N,1);
D2s = ( Xs.^2 + Ys.^2 );
figure
plot( Xs(D2s<1), Ys(D2s<1), '.r', Xs(D2s>1), Ys(D2s>1), '.b' )


%% Para ilustar o método (UM POUCO MELHOR):
clear
N = 1000;
Xs = rand(N,1);
Ys = rand(N,1);
D2s = ( Xs.^2 + Ys.^2 );
figure
plot( Xs(D2s<1), Ys(D2s<1), 'xr', Xs(D2s>1), Ys(D2s>1), 'xb')
hold on
ang = 0:.01:pi/2;
plot( sin(ang), cos(ang), ':k' )
xlabel( 'X' )
ylabel( 'Y' )
set( gca, 'xtick', [0 1], 'ytick', [0 1] ) 


%% Para ilustar o método (DIDATICAMENTE):
clear
Ns = [100 1000 10000];
figure( 'position', [50 100 1200 380] )
for qN = 1 : length(Ns)
    N = Ns(qN);
    subplot(1, length(Ns), qN )
    Xs = rand(N,1);
    Ys = rand(N,1);
    D2s = ( Xs.^2 + Ys.^2 );
    plot( Xs(D2s<1), Ys(D2s<1), 'xr', Xs(D2s>1), Ys(D2s>1), 'xb')
    hold on
    ang = 0:.01:pi/2;
    plot( sin(ang), cos(ang), '--k', 'linewidth', 3 )
    title( sprintf( '(%c) N=%d, \\sigma=%.1e', char('a'+qN-1), N, sqrt(3/N) ), ...
        'fontsize', 16, 'fontweight', 'bold'  )
    xlabel( 'X', 'fontsize', 16, 'fontweight', 'bold'  )
    ylabel( 'Y', 'fontsize', 16, 'fontweight', 'bold'  )
    set( gca, 'xtick', [0 1], 'ytick', [0 1], 'fontsize', 12, 'fontweight', 'bold', ...
        'tickdir', 'out', 'linewidth', 2 )
end


%% Para ilustar o método em 3Dimensões
clear
N = 1e4;
Xs = rand(N,1);
Ys = rand(N,1);
Zs = rand(N,1);
D2s = ( Xs.^2 + Ys.^2 + Zs.^2 );
figure
plot3( Xs(D2s<1), Ys(D2s<1), Zs(D2s<1), 'xr', Xs(D2s>1), Ys(D2s>1), Zs(D2s>1), 'xb')
title( sprintf( 'N=%d, \\sigma_\\pi=%.1e', N, 3/sqrt(N) ) )
xlabel( 'X' )
ylabel( 'Y' )
zlabel( 'Z' )
set( gca, 'xtick', [0 1], 'ytick', [0 1], 'ztick', [0 1] )
theta = 0:.01:pi/2;
fi = 0:.1:pi/2;
[ Ts, Fs ] = meshgrid( theta, fi );
hold on
plot3( sin(Fs).*sin(Ts), sin(Fs).*cos(Ts), cos(Fs), '.-k' )


%% Comparação 2D com 3D
clear
home
N=1e4;
M=400;
disp( 'Incerteza no PI por Monte Carlo 2D usando 10.000 pontos' )
disp( std(4*mean( sum( rand(2,N,M).^2 )<1 )) )

disp( 'Incerteza no PI por Monte Carlo 3D usando 10.000 pontos' )
disp( std(6*mean( sum( rand(3,N,M).^2 )<1 )) )

