matlab下三维dla模型模拟

《matlab下三维dla模型模拟》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、Matlab下三维DLA模型模拟　2007-01-1119:18分类：science字号：大中小小functiondla3dv5(Nsum,Wstep)%定义dla函数，Nsum为所生成絮体包含的颗粒数，Wstep为计算过程中所采取的步长%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%本程序内变量的定义%%%%radius为颗粒半径，release为起始释放半径%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%radius=0.5;%%颗粒半径rel

2、ease=2;%%起始释放半径L=200;Xhalf=floor(L/2);Yhalf=Xhalf;Zhalf=Xhalf;n=1;%粒子计数N(1)=1;N(Xhalf)=0;p(1,:)=[XhalfYhalfZhalf];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)];whilen

3、i*rand(1);%%粒子随机移动gama=pi*rand(1);%Wstep=Wstep*[2*rand(1)-12*rand(1)-12*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)co

4、s(gama)];elseifszpoints((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=s

5、qrt((M(1)-Xhalf)^2+(M(2)-Yhalf)^2+(M(3)-Zhalf)^2);k=round(s)+1;N(k)=N(k)+1;ifs>release%%调整释放半径release=s+3;endelseifszpoints(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

6、=0;fori=1:nmaxifp(i,1)==Xhalfj1=j1+1;X(j1,:)=p(i,:);endifp(i,2)==Yhalfj2=j2+1;Y(j2,:)=p(i,:);endifp(i,3)==Zhalfj3=j3+1;Z(j3,:)=p(i,:);endend%绘制立体图figure(1);fori=1:nmaxssphere(p(i,:),radius);holdonendshadinginterp;colormap(gray);title('絮凝分形仿真模拟结果');%沿轴线切割图形绘制，分别为垂直于x,y,z轴的切割面图figure(2

7、);fori=1:j1ssphere(X(i,:),radius);holdonendshadinginterp;colormap(gray);title('过中心垂直于X轴的切割面');%y=X(:,2);%z=X(:,3);%plot(y,z,'+');figure(3);fori=1:j2ssphere(Y(i,:),radius);holdonendshadinginterp;colormap(gray);title('过中心垂直于Y轴的切割面');%x=Y(:,1);%z=Y(:,3);%plot(x,z,'+');figure(4);fori=1:j

8、3ssphere(Z(i