Rewrite the following MatLab code for two spheres and a 3D plane:

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