Let a user define two circles in a plane. Circle 1 will always occupy area. The user can
select an option for circle 2 to either occupy area, or it can represent a hole. The input is similar to the code below, but now, the only objective is to plot the circles and the border that is used to contain the two circles. In this part, graphically show how to optimize the procedure of area estimate by making the bordering rectangle as small as possible. This will be important for the case that circle 2 represents a hole.


function CalculateOverlapArea()
% Calculates the area of 2 circles in a plane that can overlap

% Inputs: R1, R2, d
% Procedure: Use Monte Carlo integration
% Output: Show a figure of overlapping circles according to input
    % parameters, and Show a graph of area VS distance of separation

% p-error = tolerance of the accuracy as percent error
% R1 = radius for circle 1
% R2 = radius for circle 2
% d = distance of separation b/t the two circles

R1 = input ('Enter radius of circle 1: ');
if (R1<0.1)
    error('1 = natural length scale: R1 should not be less than 0.01');
end
if (R1>10.0)
    error('1 = natural length scale: R1 should not be greater than 10.0');
end

R2 = input ('Enter radius of circle 2: ');
if (R2<0.1)
    error('1 = natural length scale: R2 should not be less than 0.01');
end
if (R2>10.0)
    error('1 = natural length scale: R2 should not be greater than 10.0');
end

d = input ('Enter separation distance: ');
if (d<0.0)
    error('The separation distance must be greater than 0');
end

% ————————— plot

GraphOverlappingCircles(R1,R2,d);

disp('The area of this object will be estimated');
disp('Method: Monte Carlo Integration');
disp('   ');
perror = input('Enter an acceptable percent error: ');
if (perror<0.1)
    error('Lowest reasonable value is 0.1, be less ambitious');
end
if (perror>20)
    error('20% is too crude of an estimate, be more ambitious');
end
% ————————— use method of Monte Carlo Integration
area = TotalArea(R1,R2,d,perror)
end


% ————————— calculate total area using Monte Carlo Integration
function A = TotalArea(R1,R2,d, perror)
[xmin, xmax, ymin, ymax] = GetBoxCoordinates(R1,R2,d); % rectangular box size
%           get area of box
Abox = (xmax-xmin)*(ymax-ymin);
%           generate uniformly random points
Npts = 10000;    % number of random points
points = xmin + (xmax-xmin)*rand(Npts,2);
%           check if each point is inside the overlapping area
countInside = sum(isInside(points(:,1), points(:,2), R1, R2, d));
%           calculate estimated area
A = Abox * countInside/Npts;
disp('Estimated area: ');
disp(A);
end

% ————————— check if a point is inside the overlapping area
function inside = isInside(x, y, R1, R2, d)
circle1 = ((x-R1).^2 + y.^2 <= R1^2);
circle2 = ((x-R2-d).^2 + y.^2 <= R2^2);
inside = circle1 & ~circle2;      % circle1 - circle2 (hole)
end

% ————————— get the coordinates of the rectangular box around the circles
function [xmin, xmax, ymin, ymax] = GetBoxCoordinates(R1, R2, d)
%           define rectangular box size
xmin = min(-R1,-R2-d);
xmax = max(R1,R2+d);
ymax = max(R1,R2);
ymin = -ymax;
end

% Plot the overlapping circles and the enclosing box
function GraphOverlappingCircles(R1, R2, d)
N = 200;
t = linspace(0, 2*pi, N);
figure;
hold on;
axis equal;
fill(R1*cos(t), R1*sin(t), 'k');
fill(d + R2*cos(t), R2*sin(t), 'w');
[xmin, xmax, ymin, ymax] = GetBoxCoordinates(R1, R2, d); % rectangular box size
% ————————— expand boundaries
xmin = xmin - 0.01;
xmax = xmax + 0.01;
ymin = ymin - 0.01;
ymax = ymax + 0.01;
% ————————— lot expanded box
plot([xmin, xmax, xmax, xmin, xmin], [ymin, ymin, ymax, ymax, ymin], 'r');
end