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 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); 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'); %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%