求大佬解答，想建个球体，为什么出来是长方体？

Question

clear all
D1=input('请输入纳米球核的半径nm：');	%纳米球模型的直径（与自带相同）
D2=input('请输入纳米球壳的厚度nm：');
a=1;%	这个是用来记录当前行
d=5; %偶极子边长
l=(2*D1+2*D2)/(2*d);
D=2*(D1+D2);
if mod(D,2)==0
	xoff=0.5;
	yoff=0.5;
	zoff=0.5;
else
	xoff=0;
	yoff=0;
	zoff=0;
end
r=(D1+D2)/d;
xmin=-r-xoff+0.5;
xmax=r-xoff-0.5;
ymin=-r-l-yoff+0.5;
ymax=r-l-yoff-0.5;
zmin=-r-zoff+0.5;
zmax=r-zoff-0.5;
for x=xmin:1:xmax
    for y=ymin:1:ymax
	    for z=zmin:1:zmax
			if((x+xoff)^2+(y+l+yoff)^2+(z+zoff)^2-(D1)^2<=0) 	%条件
				I(a,1)=a;	%这个是点的编号
				I(a,2)=x;	%偶极子的x值
				I(a,3)=round(y);	%偶极子的y值
				I(a,4)=z;	%偶极子的z值
				I(a,5)=1;	%偶极子的组成
				I(a,6)=1;
				I(a,7)=1;
				a=a+1;		%换到下一行
            elseif(((x+xoff)^2+(y+l+yoff)^2+(z+zoff)^2<=(r)^2)&&((x+xoff)^2+(y+l+yoff)^2+(z+zoff)^2>(D1)^2))
                I(a,1)=a;	%这个是点的编号
				I(a,2)=x;	%偶极子的x值
				I(a,3)=round(y);	%偶极子的y值
				I(a,4)=z;	%偶极子的z值
				I(a,5)=2;	%偶极子的组成
				I(a,6)=2;
				I(a,7)=2;
				a=a+1;		%换到下一行
			end
		end
	end
end
xmin=-r-xoff+0.5;
xmax=r-xoff-0.5;
ymin=-r+l-yoff+0.5;
ymax=r+l-yoff-0.5;
zmin=-r-zoff+0.5;
zmax=r-zoff-0.5;
for x=xmin:1:xmax
    for y=ymin:1:ymax
	    for z=zmin:1:zmax
			if((x+xoff)^2+(y-l+yoff)^2+(z+zoff)^2-(D1)^2<=0) 	%条件
				I(a,1)=a;	%这个是点的编号
				I(a,2)=x;	%偶极子的x值
				I(a,3)=round(y);	%偶极子的y值
				I(a,4)=z;	%偶极子的z值
				I(a,5)=1;	%偶极子的组成
				I(a,6)=1;
				I(a,7)=1;
				a=a+1;		%换到下一行
            elseif(((x+xoff)^2+(y-l+yoff)^2+(z+zoff)^2<=(r)^2)&&((x+xoff)^2+(y-l+yoff)^2+(z+zoff)^2>(D1)^2))
                I(a,1)=a;	%这个是点的编号
				I(a,2)=x;	%偶极子的x值
				I(a,3)=round(y);	%偶极子的y值
				I(a,4)=z;	%偶极子的z值
				I(a,5)=2;	%偶极子的组成
				I(a,6)=2;
				I(a,7)=2;
				a=a+1;		%换到下一行
			end
		end
	end
end
aeff=d*((3*(a-1)/(4*pi))^(1/3))
aeff1=(2*(D/2)^3)^(1/3)

plot3(I(:,2),I(:,3),I(:,4),'b.','MarkerSize',0.5)	%绘制图形

fid=fopen('C:shapecoreshell503.dat','w');	%此处选择要保存文件的位置
[a b]=size(I);

for i=1:a
    for j=1:b
       fprintf(fid,'%10d',I(i,j)); %将I数组保存到文件中
    end

   fprintf(fid,'\n');
end
fclose(fid);

Answer 1

没有很仔细看，但是你建立出来一个长方体，再根据这三句话：

for x=xmin:1:xmax

for y=ymin:1:ymax

for z=zmin:1:zmax

我猜你应该是把非球体区域的也给赋值了，所以看起来就是个长方体了。假设球内核半径r，球壳的厚度为w，那么你球壳的部分应该是

r^2<=x^2+y^2+z^2<=(r+w)^2的这部分赋值，其他部分设置空白或者不赋值。