%Estimate PI Monte Carlo Simulation
%By Xiaoping Du, 02/19/2014
%Center of the circle = (0,0)
%Radius of the circle = 1
clc;close all; clear all;
% n=input('Enter number of samples:');   %# of samples
n=1000000;
rand('seed',108);                                   %Initialize the random generator
x=2*rand(n,1)-1;                                     %Samples of x 
y=2*rand(n,1)-1;                                     %Samples of y
m=sum(x.^2+y.^2<1);                             %# of samples in the circle
pi_est=4*m/n;                                         %Estimate of PI 


%Display
rectangle('Position',[-1,-1,2,2]);             %Draw the square
hold on;
theta=linspace(0,2*pi,1e3);                    %Draw the circle
rho=ones(1,1e3);
[xc,yc] = pol2cart(theta,rho);
plot(xc,yc,'k-','linewidth',2);                     
axis square;

inside=find(x.^2+y.^2<1);                       %Get samples in the circle
outside=find(x.^2+y.^2>=1);                   %Get samples outside the circle
plot(x(inside),y(inside),'r.');                     %Plot samples in the circle
plot(x(outside),y(outside),'b.');               %Plot samples outside the circle

text(-0.8,-1.2,['m/n = ',num2str(m),'/',num2str(n),...
    ', PI = ',num2str(pi_est,'%8.6f')],'FontSize',14); %Display the result
