Find area of two overlapping circles using monte carlo method

Question

Actually i have two intersecting circles as specified in the figure

i want to find the area of each part separately using Monte carlo method in Matlab .

The code doesn't draw the rectangle or the circles correctly so i guess what is wrong is my calculation for the x and y and i am not much aware about the geometry equations for solving it so i need help about the equations.

this is my code so far :

n=1000;
%supposing that a rectangle will contain both circles so :
% the mid point of the distance between 2 circles will be (0,6) 
% then by adding the radius of the left and right circles the total distance 
% will be 27 , 11 from the left and 16 from the right 
% width of rectangle = 24

x=27.*rand(n-1)-11;
y=24.*rand(n-1)+2;
count=0;

for i=1:n

if((x(i))^2+(y(i))^2<=25 && (x(i))^2+(y(i)-12)^2<=100)
count=count+1;        
    figure(2);
    plot(x(i),y(i),'b+')
    hold on

elseif(~(x(i))^2+(y(i))^2<=25 &&(x(i))^2+(y(i)-12)^2<=100)  
    figure(2);
    plot(x(i),y(i),'y+')
    hold on

else 
     figure(2);
    plot(x(i),y(i),'r+')

end

end

Answer 1

Without using For loop.

    n = 100000;
    data = rand(2,n);
    data = data*2*30 - 30;
    x = data(1,:);
    y = data(2,:);
    plot(x,y,'ro');
    inside5 = find(x.^2 + y.^2 <=25);
    hold on 
    plot (x(inside5),y(inside5),'bo');
    hold on 
    inside12 = find(x.^2 + (y-12).^2<=144);
    plot (x(inside12),y(inside12),'g');
    hold on 
    insidefinal1 = find(x.^2 + y.^2 <=25 & x.^2 + (y-12).^2>=144);
    insidefinal2 = find(x.^2 + y.^2 >=25 & x.^2 + (y-12).^2<=144);
    % plot(x(insidefinal1),y(insidefinal1),'bo');
    hold on
    % plot(x(insidefinal2),y(insidefinal2),'ro');
    insidefinal3 = find(x.^2 + y.^2 <=25 & x.^2 + (y-12).^2<=144);
    % plot(x(insidefinal3),y(insidefinal3),'ro');
    area1=(60^2)*(length(insidefinal1)/n);
    area3=(60^2)*(length(insidefinal2)/n);
    area2= (60^2)*(length(insidefinal3)/n);

Answer 2

Here is my generic solution for any two circles (without any hardcoded value):

function [ P ] = circles_intersection_area( k1, k2, N )
%CIRCLES_INTERSECTION_AREA Summary...
%   Adnan A.
x1 = k1(1);
y1 = k1(2);
r1 = k1(3);

x2 = k2(1);
y2 = k2(2);
r2 = k2(3);

if sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)) >= (r1 + r2)
    % no intersection
    P = 0;
    return
end

% Wrapper rectangle config
a_min = x1 - r1 - 2*r2;
a_max = x1 + r1 + 2*r2;
b_min = y1 - r1 - 2*r2;
b_max = y1 + r1 + 2*r2;

% Monte Carlo algorithm
n = 0;
for i = 1:N
    rand_x = unifrnd(a_min, a_max);
    rand_y = unifrnd(b_min, b_max);

    if sqrt((rand_x - x1)^2 + (rand_y - y1)^2) < r1 && sqrt((rand_x - x2)^2 + (rand_y - y2)^2) < r2
        % is a point in the both of circles
        n = n + 1;
        plot(rand_x,rand_y, 'go-');
        hold on;
    else
        plot(rand_x,rand_y, 'ko-');
        hold on;
    end
end

P = (a_max - a_min) * (b_max - b_min) * n / N;

end

Call it like: circles_intersection_area([-0.4,0,1], [0.4,0,1], 10000) where the first param is the first circle (x,y,r) and the second param is the second circle.

Answer 3

Here are the errors I found:

x = 27*rand(n,1)-5
y = 24*rand(n,1)-12

The rectangle extents were incorrect, and if you use rand(n-1) will give you a (n-1) by (n-1) matrix.

and

first If:

(x(i))^2+(y(i))^2<=25 && (x(i)-12)^2+(y(i))^2<=100

the center of the large circle is at x=12 not y=12

Second If:

~(x(i))^2+(y(i))^2<=25 &&(x(i)-12)^2+(y(i))^2<=100

This code can be improved by using logical indexing.

For example, using R, you could do (Matlab code is left as an excercise):

n = 10000
x = 27*runif(n)-5
y = 24*runif(n)-12
plot(x,y)

r = (x^2 + y^2)<=25 & ((x-12)^2 + y^2)<=100
g = (x^2 + y^2)<=25
b = ((x-12)^2 + y^2)<=100
points(x[g],y[g],col="green")
points(x[b],y[b],col="blue")
points(x[r],y[r],col="red")

which gives: