clc
clear all
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];
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
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.2;
y1range1 = -0.9;
y1range2 = 2.9;
z1range1 = -0.9;
z1range2 = 2.9;
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=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(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);
