Matlab下三维DLA模型模拟

function dla3dv5(Nsum,Wstep) %定义dla函数，Nsum为所生成絮体包含的颗粒数，Wstep为计算过程中所采取的步长%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%

%%本程序内变量的定义%%

%%radius为颗粒半径，release为起始释放半径%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%

radius=0.5; %%颗粒半径

release=2; %%起始释放半径

L=200;

Xhalf=floor(L/2);

Yhalf=Xhalf;

Zhalf=Xhalf;

n=1; %粒子计数

N(1)=1;

N(Xhalf)=0;

p(1,:)=[Xhalf Yhalf Zhalf];

szpoints=zeros(L,L,L); %%网格点阵

szpoints(Xhalf,Yhalf,Zhalf)=1; %%种子位置标志

%%释放初始粒子

theta=2*pi*rand;

gama=pi*rand;

M=p(1,:)+Wstep*[cos(theta) sin(theta) cos(gama)];

while n

theta=2*pi*rand(1); %%粒子随机移动

gama=pi*rand(1);

%Wstep=Wstep*[2*rand(1)-1 2*rand(1)-1 2*rand(1)-1];

step=Wstep*[sin(theta) cos(theta) cos(gama)];

M=M+step;

T=round(M);

if (M(1)-Xhalf)^2+(M(2)-Yhalf)^2+(M(3)-Zhalf)^2>(release+15)^2 %%判断是否逃逸

theta=2*pi*rand;

gama=pi*rand;

M=p(1,:)+release*[cos(theta) sin(theta) cos(gama)];

elseif

szpoints((T(1)-1),T(2),T(3))+szpoints((T(1)+1),T(2),T(3))+szpoints(T(1),(T(2)-1),T(3))+szp oints(T(1),(T(2)+1),T(3))+szpoints(T(1),T(2),(T(3)-1))+szpoints(T(1),T(2),(T(3)+1))>0&szp oints(T(1),T(2),T(3))~=1 %%判断是否凝结

n=n+1;

szpoints(T(1),T(2),T(3))=1;

p(n,:)=T; %存储凝聚颗粒的球心坐标。

s=sqrt((M(1)-Xhalf)^2+(M(2)-Yhalf)^2+(M(3)-Zhalf)^2);

k=round(s)+1;

N(k)=N(k)+1;

if s>release %%调整释放半径

release=s+3;

end

elseif szpoints(T(1),T(2),T(3))==1 %%检查是否出现漏检,即运动一步后进入粒子内部的情况.

theta=2*pi*rand(1);

gama=pi*rand(1);

M=p(1,:)+Wstep*[cos(theta) sin(theta) cos(gama)];

end

end

nmax=size(p,1);

j1=0;

j2=0;

j3=0;

for i=1:nmax

if p(i,1)==Xhalf

j1=j1+1;

X(j1,:)=p(i,:);

end

if p(i,2)==Yhalf

j2=j2+1;

Y(j2,:)=p(i,:);

end

if p(i,3)==Zhalf

j3=j3+1;

Z(j3,:)=p(i,:);

end

end

%绘制立体图

figure(1);

for i=1:nmax

ssphere(p(i,:),radius);

hold on

end

shading interp;

colormap(gray);

title('絮凝分形仿真模拟结果');

%沿轴线切割图形绘制，分别为垂直于x,y,z轴的切割面图figure(2);

for i=1:j1

ssphere(X(i,:),radius);

hold on

end

shading interp;

colormap(gray);

title('过中心垂直于X轴的切割面');

%y=X(:,2);

%z=X(:,3);

%plot(y,z,'+');

figure(3);

for i=1:j2

ssphere(Y(i,:),radius);

hold on

end

shading interp;

colormap(gray);

title('过中心垂直于Y轴的切割面');

%x=Y(:,1);

%z=Y(:,3);

%plot(x,z,'+');

figure(4);