Matlab下三维DLA模型模拟

Matlab下三维DLA模型模拟2007-01-1119:18分类：science字号：大中小functiondla3dv5(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,:)=[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)];whilenNsumtheta=2*pi*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)cos(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=sqrt((M(1)-Xhalf)^2+(M(2)-Yhalf)^2+(M(3)-Zhalf)^2);k=round(s)+1;N(k)=N(k)+1;ifsrelease%%调整释放半径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=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);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:j3ssphere(Z(i,:),radius);holdonendtitle('过中心垂直于Z轴的切割面');%x=Z(:,1);%y=Z(:,2);%plot(x,y,'+');shadinginterp;colormap(gray);%%为了与dlacon连用而进行存储数据，不用时可以注销%savedatapp;%savedatareleaserelease;%savedataszpszpoints;%savedataMM;%savedataNsumNsum;%savedatann;%savedatahalfXhalf;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%分维分析计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%num=find(N);R(1)=0.5;fori=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)=[];%savedataNnum;%savedataRR;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');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%