MATLAB Answers

Monte Carlo Two Overlapping Spheres

22 views (last 30 days)
Kelly McGuire
Kelly McGuire on 12 Jan 2019
Commented: M.S. Khan on 25 Feb 2020
I am trying to estimate the volume of two overlapping spheres using the Monte Carlo method. One sphere (white) is centered at 1,1,1. The other sphere is centered at 2.293,1,1. The length from one sphere center to the other sphere center is 1.293. The white sphere's radius is 1.2 and the green sphere's radius is 1.8. I am trying to fill the white sphere with white '+' marks and the green sphere with green '+' marks. Any volume outside of these overlapping spheres is filled with red '+' marks. So far, my code is very close to working the way I want. However, red '+' marks are still getting in the green sphere. It doesn't look like they are placed in the white sphere. I have four questions:
1) Why are red '+' marks still getting in based on the code I am pasting here?
2) If you comment out line 57 (no red marks plotted), and then choose around 400 points, you'll notice in the plot that the marks do not fill up the green sphere, but they do fill the volume of the white sphere. Why aren't marks being placed throughout all of the green sphere?
3) What should I do about the overlapping volume region?
4) How do I then calculate the volume of the overlapping spheres compared to the non-sphere volume? Including the overlapping part of the spheres (question 2).
Code:
clc
clear
N=input('Number of points: ');
[X, Y, Z]=sphere;
a = [1 1 1 1.2; 2.293,1,1,1.8];
%Draw sphere #1
s1=surf(X*a(1,4)+a(1,1), Y*a(1,4)+a(1,2), Z*a(1,4)+a(1,3),'FaceColor', [1 1 1],'edgecolor','none','FaceAlpha',0.6);
light
axis equal
set(gca,'Color','k')
hold on
s2=surf(X*a(2,4)+a(2,1), Y*a(2,4)+a(2,2), Z*a(2,4)+a(2,3),'FaceColor', [0 1 0],'edgecolor','none','FaceAlpha',0.5);
daspect([1 1 1])
view(30,10)
xlabel('x')
ylabel('y')
zlabel('z')
x1range1 = -0.2;
x1range2 = 4.093;
y1range1 = -0.8;
y1range2 = 2.8;
z1range1 = -0.8;
z1range2 = 2.8;
hit1 = 0;
hit2 = 0;
for i = 1:N
x = (x1range2-x1range1).*rand(N,1)+x1range1;
y = (y1range2-y1range1).*rand(N,1)+y1range1;
z = (z1range2-z1range1).*rand(N,1)+z1range1;
x1=x-1;
y1=y-1;
z1=z-1;
r1=x1.^2+y1.^2+z1.^2;
x2=x-2.293;
y2=y-1;
z2=z-1;
r2=x2.^2+y2.^2+z2.^2;
ii = r1<=1.2;
jj = r2<=1.8;
hit1 = sum(ii);
hit2 = sum(jj);
plot3(x(ii), y(ii), z(ii), 'w+');
plot3(x(~ii), y(~ii), z(~ii), 'r+');
plot3(x(jj), y(jj), z(jj), 'g+');
%plot3(x(~jj), y(~jj), z(~jj), 'r+');
drawnow
end

  0 Comments

Sign in to comment.

Accepted Answer

Kelly McGuire
Kelly McGuire on 13 Jan 2019
Here is the complete working code:
clc
clear
S=input('Number of Simulations: ');
N=input('Number of points: ');
fileID = fopen('Volume.txt','w');
for k = 1:S
[X, Y, Z]=sphere;
a = [1 1 1 1.2; 2.293,1,1,1.8];
%Draw sphere #1
figure('color','white');
s1=surf(X*a(1,4)+a(1,1), Y*a(1,4)+a(1,2), Z*a(1,4)+a(1,3),'FaceColor', [1 1 1],'edgecolor','none','FaceAlpha',0.6);
light
axis equal
set(gca,'Color','k')
hold on
%Draw sphere #2
s2=surf(X*a(2,4)+a(2,1), Y*a(2,4)+a(2,2), Z*a(2,4)+a(2,3),'FaceColor', [0 1 0],'edgecolor','none','FaceAlpha',0.5);
daspect([1 1 1])
view(30,10)
xlabel('x')
ylabel('y')
zlabel('z')
% x1range1 = -0.2;
% x1range2 = 4.093;
% y1range1 = -0.8;
% y1range2 = 2.8;
% z1range1 = -0.8;
% z1range2 = 2.8;
x1range1 = -1.0;
x1range2 = 5.0;
y1range1 = -1.0;
y1range2 = 3.0;
z1range1 = -1.0;
z1range2 = 3.0;
hit1 = 0; %Points inside spheres
hit2 = 0;
for i = 1:N
%Define Box Range and Randomize Points
x = (x1range2-x1range1).*rand(N,1)+x1range1;
y = (y1range2-y1range1).*rand(N,1)+y1range1;
z = (z1range2-z1range1).*rand(N,1)+z1range1;
x1=x-1;
y1=y-1;
z1=z-1;
r1=sqrt(x1.^2+y1.^2+z1.^2);
x2=x-2.293;
y2=y-1;
z2=z-1;
r2=sqrt(x2.^2+y2.^2+z2.^2);
ii = r1<=1.2;
jj = r2<=1.8;
hit1 = sum(ii);
hit2 = sum(jj);
plot3(x(ii), y(ii), z(ii), 'w+');
%plot3(x(~ii), y(~ii), z(~ii), 'r+');
plot3(x(jj), y(jj), z(jj), 'g+');
drawnow
end
SpherePercentVolume = (hit1+hit2)/N;
Volume = ((x1range2-x1range1)*(y1range2-y1range1)*(z1range2-z1range1))*SpherePercentVolume
fprintf(fileID,'%6.3f\r\n',Volume);
end
fclose(fileID);
load Volume.txt
average = mean(Volume);
stdev = std(Volume);
fileID = fopen('AvgVolume.txt','w');
fprintf(fileID,'%6s','Average Volume: ');
fprintf(fileID,'%6.3f\r\n',average);
fprintf(fileID,'%6s','Standard Deviation: ');
fprintf(fileID,'%6.3f',stdev);

  2 Comments

Image Analyst
Image Analyst on 13 Jan 2019
Please highlight the code and click the icon to format it as code to make it easy for people to help you by clicking the copy button.
M.S. Khan
M.S. Khan on 25 Feb 2020
could you please explain this code. Its very hard to understand. Could we simply measure the volume of these intersected spheres.

Sign in to comment.

More Answers (0)

Sign in to answer this question.