matlab下三维dla模型模拟,Matlab下三维DLA模型模拟.doc

Matlab下三维DLA模型模拟　2007-01-11 19:18

分类：science

字号： 大 中 小小

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(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))+szpoints(T(1),(T(2)+1),T(3))+szpoints(T(1),T(2),(T(3)-1))+szpoints(T(1),T(2),(T(3)+1))>0&szpoints(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);

for i=1:j3

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

hold on

end

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

%x=Z(:,1);

%y=Z(:,2);

%plot(x,y,+);

shading interp;

colormap(gray);

%%为了与dlacon连用而进行存储数据，不用时可以注销

%save datap p;

%save datarelease release;

%save dataszp szpoints;

%save dataM M;

%save dataNsum Nsum;

%save datan n;

%save datahalf Xhalf;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%分维分析计算

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

num=find(N);

R(1)=0.5;

for i=2:size(num,2)

N(i)=N(i)+N(i-1);

num(i)=N(i);

R(i)=(i-0.5)/0.5;

end

num(1)=[];

R(1)=[];

%save dataN num;

%save dataR R;

figure(5);

plot(log(R),log(num),*);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%figure(3);

%kxl=1-num.*(radius./R).^3;

%plot(R,kxl,*);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%图形绘制 %%

%figure; %绘制回转半径与其内粒子数的对数关系图%%

%plot(R, N,*); %%

%%figure(3); %绘制凝聚粒子数目与步长的关系图 %%

%%plot(Nn,I,o); %%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%