clear
N = 1000;
x = rand(N,1);
y = rand(N,1);
D = (x.^2 + y.^2);
pi_mc = 4 * sum( D<=1 ) / N
ang = linspace(0,pi/2,1000);
figure
plot( x(D<=1), y(D<=1), 'xb' )
hold on
plot( x(D>1), y(D>1), 'xr' )
plot( cos(ang), sin(ang), '--k', 'linewidth', 3 )
