%% Para ilustar o método:
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), '.r', Xs(D2s>1), Ys(D2s>1), '.b')

hold on
ang = 0:.01:pi/2;
plot( sin(ang), cos(ang), '-k' )
xlabel( 'X' )
ylabel( 'Y' )


%% Para fazer uma figura para ilustração do método de maneira didática:

clear
Ns = [100 1000 10000];
figure
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), '.r', Xs(D2s>1), Ys(D2s>1), '.b')
    hold on
    ang = 0:.01:pi/2;
    plot( sin(ang), cos(ang), '-k' )
    xlabel( 'X' )
    ylabel( 'Y' )
    
    title( sprintf( 'N = %d', N)  )    
end

%% AVANÇADO: Figura de ilustração didática com aparência melhorada

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), '.r', Xs(D2s>1), Ys(D2s>1), '.b')
    hold on
    ang = 0:.01:pi/2;
    plot( sin(ang), cos(ang), '-k', 'linewidth', 3 )
    
    xlabel( 'X', 'fontsize', 16, 'fontweight', 'bold'  )
    ylabel( 'Y', 'fontsize', 16, 'fontweight', 'bold'  )
    title( sprintf( '(%c) N = %d', char('a'+qN-1), N), 'fontsize', 16, 'fontweight', 'bold'  )
    
    xlim([0 1])
    ylim([0 1])
    set( gca, 'xtick', [0  1], 'ytick', [0  1], 'fontsize', 12, 'fontweight', 'bold', ...
        'tickdir', 'out', 'linewidth', 2 )
end
