% mcPi
% Monte Carlo approximation of pi

close all % close all figure windows
figure % new figure window
axis equal % all axes: same scale
axis off % all axes: don't display
hold on % hold all calls to plot on current axes

L= 10; % length of square (diameter of circle)
R= L/2; % half length (radius)
axis([-R R -R R]) % axis limits ([xmin xmax ymin ymax])

N= 2000; % total number of darts
hits= 0; % number of darts in the circle so far
for k = 1:N
% Throw kth dart
x= rand()*L - R;
y= rand()*L - R;
if sqrt(x^2 + y^2) <= R
plot(x,y,'.r')
hits= hits + 1;
else
plot(x,y,'.b')
end
end

myPi= 4*hits/N;
err= abs(myPi-pi);
titleText = sprintf('Pi Estimate = %6.3f Error = %6.2e', ...
myPi, err);
title(titleText, 'Fontsize',14)

hold off % next call to plot will be on new axes
shg % Show Handle Graphics (bring figure window to the front)