《matlab下三维dla模型模拟》由会员分享，可在线阅读，更多相关《matlab下三维dla模型模拟（6页珍藏版）》请在金锄头文库上搜索。
1、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 Y
2、half 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
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
4、)+1;if srelease %调整释放半径release=s+3;endelseif 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);endendnmax=size(p,1);j1=0;j2=0;j3=0;for i=1:nmaxif p(i,1)=Xhalfj1=j1+1;X(j1,:)=p(i,:);endif p(i,2)=Yhalfj2=j2+1;Y(j2,:)=
5、p(i,:);endif p(i,3)=Zhalfj3=j3+1;Z(j3,:)=p(i,:);endend%绘制立体图figure(1);for i=1:nmaxssphere(p(i,:),radius);hold onendshading interp;colormap(gray);title(絮凝分形仿真模拟结果);%沿轴线切割图形绘制，分别为垂直于 x,y,z 轴的切割面图figure(2);for i=1:j1ssphere(X(i,:),radius);hold onendshading interp;colormap(gray);title(过中心垂直于 X 轴的切割面);%y
6、=X(:,2);%z=X(:,3);%plot(y,z,+);figure(3);for i=1:j2ssphere(Y(i,:),radius);hold onendshading 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 onendtitle(过中心垂直于 Z 轴的切割面);%x=Z(:,1);%y=Z(:,2);%plot(x,y,+);shading interp;colorma
7、p(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;endnum(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); %
