[MATLAB数学相关] 求助：椭球随机分布程序

[fusijie1]

谁能看看这个程序啊，是生成随机球体的，可总是不对，希望得到帮助，谢谢

xyzr=[];
n=10000;
while n
x=100*rand;
y=100*rand;
z=100*rand;
r=3+2*rand;
if ~isempty(xyzr)
d1=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,3)-z).^2);
d2=xyzr(:,4)+r;
if all(d1>=d2)
if x>=r && x<=100-r && y>=r && y<=100-r && z>=r && z<=100-r
xyzr=[xyzr;x y z r];
end
if x<r && y>=r && y<=100-r%x<rµÄÇé¿ö
d1=sqrt((xyzr(:,1)-x-100).^2+(xyzr(:,2)-y).^2+(xyzr(:,3)-z).^2);
d2=xyzr(:,4)+r;
%% if y>=r && y<=100-r %x<rʱ,y>=r
if z>=r && z<=100-r %
if all(d1>=d2)
xyzr=[xyzr;x y z r;x+100 y z r];
end
else if z<r
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z-100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x+100 y z r;x y z r];
end
else
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z+100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x+100 y z r;x y z-100 r];
end
end
elseif x<r && y<r %x<r,y<r
if z>=r && z<=100-r %
if all(d1>=d2)
xyzr=[xyzr;x y z r;x+100 y+100 z r];
end
else if z<r
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y-100).^2+(xyzr(:,2)-z-100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x+100 y+100 z r;x y z r];
end
else
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y+100).^2+(xyzr(:,2)-z+100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x+100 y z r;x y-100 z-100 r];
end
end
end

if x>100-r %x>100-r
d1=sqrt((xyzr(:,1)-x+100).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z).^2);
d2=xyzr(:,4)+r;
if y>=r && y<=100-r
if z>=r && z<=100-r
if all(d1>=d2)
xzyr=[xyzr;x y z r;x-100 y z r];
end
elseif z<r
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z-100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x-100 y z r;x y z+100 r];
end
else
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z+100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x-100 y z r;x y z-100 r];
end
end
elseif y<r
if z>=r && z<=100-r
if all(d1>=d2)
xzyr=[xyzr;x y z r;x-100 y-100 z r];
end
elseif z<r
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y-100).^2+(xyzr(:,2)-z-100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x-100 y z r;x y+100 z+100 r];
end
else
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y+100).^2+(xyzr(:,2)-z+100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x-100 y z r;x y-100 z-100 r];
end
end
end

else
n=n-1;
end
else
xyzr=[x y z r];
if x<r
xyzr=[xyzr;x+100 y z r];
end
if x>100-r
xyzr=[xyzr;x-100 y z r];
end
if y<r
xyzr=[xyzr;x y+100 z r];
end
if y>100-r
xyzr=[xyzr;x y-100 z r];
end
if z<r
xyzr=[xyzr;x y z+100 r];
end
if z>100-r
xyzr=[xyzr;x y z-100 r];
end
end

%%t=[0:0.1:2*pi 0]';
delta=pi/12;
theta=0:deltai;
phi=0:delta:2*pi;

x=sin(theta)*cos(phi)*xyzr(:,[1 4])';
y=sin(theta)*sin(phi)*xyzr(:,[2 4])';
z=cos(theta)*xyzr(:,[3 4])';
h=plot(x,y,z);
hold on
plot(xyzr(:,1),xyzr(:,2),xyzr(:,3),'.')
set(h,'color','r')
axis equal
axis([0 100 0 100 0 100])
hold off