Modify the MatLab code below so that it is in three dimensions and the circles are now spheres. The code needs to plot the total occupied volume covered by the two spheres as a function of d, where d is the distance between the centers of the spheres. Must also include error bars in the plot.


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 (graph of area as function of separation distance),
    % including error bars. Also give the resulting area estimated by the
    % Monte Carlo Integration.

% 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
plotOption = input ('Enter plot option (1 or 0): ');
GraphOverlappingCircles(R1,R2,d,plotOption);

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,plotOption);
end

% ————————— calculate total area using Monte Carlo Integration
function A = TotalArea(R1,R2,d,perror,plotOption)
[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, plotOption));
%           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, plotOption)
circle1 = ((x-R1).^2 + y.^2 <= R1^2);
circle2 = ((x-R2-d).^2 + y.^2 <= R2^2);

if plotOption == 1
    inside = circle1 & ~circle2;    % circle1 - circle2 (hole)
else
    inside = circle1 & circle2;
end
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 overlapping circles and enclosing box
function GraphOverlappingCircles(R1,R2,d,plotOption)
N = 200;
t = linspace(0, 2*pi, N);
figure;
hold on;
axis equal;
fill(R1*cos(t), R1*sin(t), 'k');

if plotOption == 1
    fill(d+R2*cos(t), R2*sin(t), 'w');
else
    fill(d+R2*cos(t), R2*sin(t), 'k');
end

[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;
% ————————— plot expanded box
plot([xmin,xmax,xmax,xmin,xmin], [ymin,ymin,ymax,ymax,ymin], 'r');
end