MATLAB实现长方体容器内球形颗粒堆积
右手坐标系；容器的中心点为原点；
容器尺寸；球的体积分数；球的半径可以调整；
有过程的动态演示；
最后可以输出球的位置信息并保存。
PS：希望能写注释比较详细！
动态绘图，并存储。请采纳
xl=200;%长
yl=200;%宽
zl=100;%高
xmin=-xl/2;%计算x轴最负值，因为0是原点。下同
ymin=-yl/2;
zmin=-zl/2;
r=10;%球半径
n=40;%球数量
for i = 1:n%循环生成初始位置
if i==1%如果是第一个，直接随机生成
x(i)=rand()*(xl-2*r)+xmin+r;
y(i)=rand()*(yl-2*r)+ymin+r;
z(i)=rand()*(zl-2*r)+zmin+r;
else%如果不是第一个
while 1%进入死循环，并生成位置
x(i)=rand()*(xl-2*r)+xmin+r;
y(i)=rand()*(yl-2*r)+ymin+r;
z(i)=rand()*(zl-2*r)+zmin+r;
st=0;%标记
for j=1:i-1%遍历所有已生成的球心
if ((x(i)-x(j))^2+(y(i)-y(j))^2+(z(i)-z(j))^2)<(2*r)^2 %如果干涉
st=1;%更改标记
break %跳出for循环
end
end
if st==0%如果for循环完毕后，st依然为0，表明没有干涉
break%则推出死循环
end
end
end
end
while 1%循环，用于绘制动态球下落
st2=0;%标记，用于识别是否所有球都下落完
for i= 1:n%遍历球
st=0;%标记
for j=1:n%遍历球
if i~=j%如果ij不是一个球，则进行判断，如果两球干涉，且i球高于j球，则停止下落
if ((x(i)-x(j))^2+(y(i)-y(j))^2+(z(i)-z(j))^2)<(2*r)^2 && z(i)>z(j)
st=1;%修改下落标记
break
end
end
if z(i)-zmin<r %如果i球已经到底面
st=1;%修改下落标记
break
end
end
if st==0%如果标记位0，则下落
z(i)=z(i)-1;%下落计算，步长为1，可更改，越小越慢
st2=1;%更改st2标记，表明有球下落，循环继续
end
end
if st2==0%如果无球下落，则跳出循环
break
end
for i =1:n%遍历所有球绘图
[xx,yy,zz]=ellipsoid(x(i),y(i),z(i),r,r,r);%设定绘制的球心和半径
surf(xx,yy,zz)%绘制球
if i==1%当绘制第一个球的时候
hold on%打开hold on，否则只显示一个球
axis([xmin,xmin+xl,ymin,ymin+yl,zmin,zmin+zl])%设置坐标轴
end
end
hold off%每次绘制完一轮球，关闭hold，以清屏
pause(0.01)%延时
end
aa=[x;y;z];%把球的坐标输入到aa
t=table(aa);%把aa变成表t
tmp={'x','y','z'};%设置表头
t.Properties.RowNames=tmp;%设置表头
writetable(t,'1.xlsx','WriteRowNames',true);%输出到excel保存
采纳一下，谢谢
function spheres
clc,clear;
n=50; %球的个数
x1=0;x2=100; %x轴范围
y1=0;y2=100; %y轴范围
z1=0;z2=20; %z轴范围
r=2; %球的半径
A=[x2-x1;y2-y1;z2-z1];
A=A-2*r;
centre=zeros(3,n-1);%保存圆心坐标
for i=1:n
c=rand(3,1);
c=c.*A; %圆心坐标
c=c+[x1;y1;z1]+r;
if i==1
[x,y,z] = ellipsoid(c(1),c(2),c(3),r,r,r);
surf(x,y,z,ones(size(x))) %画出来球
hold on
else
s=0;
while(s==0)
c=rand(3,1);
c=c.*A; %圆心坐标
c=c+[x1;y1;z1]+2;
s=judge(i,c,centre,r); %判断圆心是否可以 即新画出的球不会与之前的重叠
end
[x,y,z] = ellipsoid(c(1),c(2),c(3),r,r,r);
surf(x,y,z,ones(size(x))) %画出来球
hold on
end
end
axis([0 100 0 100 0 20]);
function s=judge(i,c,centre,r)
for j=1:i-1
temp=centre(:,j);
f=c-temp;
d=norm(f,2);
if d<2*r
s=0; %表示新圆心不可以在c表示的点
return
end
end
s=1;%表示新圆心可以