% Monte Carlo approximation of pi

L= 10;      % length of square (diameter of circle)
r= 2;       % Radius of inner circle
tol= 1e-2;  % tolerance

figure      % Open new plotting window
hold on     % don't erase old points
axis equal on

hits= 0;    % number of darts in the circle so far
myPi= 0;    % Current estimate of pi
k= 0;       % number of darts thrown so far
while abs(myPi - pi) > tol
    % Throw kth dart
    x= rand()*L - L/2;
    y= rand()*L - L/2;
    k= k + 1;

    if (x^2 + y^2 <= (L/2)^2) && (x^2 + y^2 >= r^2)
        hits= hits + 1;
        plot(x, y, 'b*')
    else
        plot(x, y, 'r*')
    end
    
    myPi= (hits/k) / (1/4 - (r/L)^2);
end
fprintf('Number of darts thrown: %f\n', k)
fprintf('Estimate of pi: %f\n', myPi)

hold off    % restore default plotting behavior
