Matlab下三维DLA模型模拟知识分享.doc

Matlab下三维DLA模型模拟

精品文档
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<Nsum
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))+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'); %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

收集于网络，如有侵权请联系管理员删除