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 srelease %调整释放半径 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); endendnmax=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,:); endend%绘制立体图figure(1);for i=1:nmax ssphere(p(i,:),radius); hold onendshading interp;colormap(gray);title(絮凝分形仿真模拟结果);%沿轴线切割图形绘制，分别为垂直于x,y,z轴的切割面图figure(2);for i=1:j1 ssphere(X(i,:),radius); hold onendshading 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 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;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;endnum(1)=;R(1)=;%save dataN num;%save dataR R;figure(5);plot(log(R),log(num),*);%figure(3);%kxl=1-num.*(radius./R