Monte Carlo Pi

clear
clc
tic
format long

k = round([exp(3:0.05:15)]);

%Use for better Pi estimation:
%k = [10:10:1000 1000:500:100000 100000:10000:1000000];

i = 1;
c = zeros(1,length(k));

hold on

for f = k
    figure(1);
    clf;

    x = rand(1,f);
    y = rand(1,f);


%     Using a non-random distribution does not return accurate
%     values for pi:
%      [x,y] = meshgrid(0:0.01+10/f:1);
%      x = x(:);
%      y = y(:);

    b = sqrt( x.^2+y.^2) <= 1;
    c(i) = 4*mean(b);


    clc
    subplot(211);
    axis square
    plot(x(b),y(b),'.r')
    hold on
    plot(x(~b),y(~b),'.g')
    pbaspect([1 1 1])

    P = polyfit(k(2:i-1),c(2:i-1),1);
    clc
    yfit = P(1)*k(2:i-1)+P(2);

    subplot(212);
    semilogx(k(2:i-1),c(2:i-1),'-b')
    title(['\pi = ' num2str(mean(c(1:i)))])
    hold on
    semilogx(k(2:i-1),yfit,'-r');
    pause(0.00000001)
    pbaspect([1 1 1])

    i = i+1;

end

Pi = mean(c)

lin_fit_gradient = P(1)

time_elapsed = toc