如何用matlab画烧杯,DLA模型的Matlab程序.doc

大学里写了一些关于DLA模型的Matlab程序，由于不懂编程，而且也不懂matlab，在程序算法的编制上走了很多弯路，这里将陆续贴出关于DLA、DLCA、先DLA后DLCA模型的Matlab下实现和分维分析，请大家指导。

程序按gpl发布(:(不明白gpl什么意思，只知道共享)，如果你能有所改进也希望你能与我分享。

二维dla模型程序：

文件名：dla2d.m

功能：生成所需粒子数的二维dla絮体，绘制回转半径与粒子数的对数关系图，空隙率与絮体半径的关系图。

代码：

%(1) 设置初始凝结半径 ,总凝集粒子数nsum,绘图区域;

%(2)根据绘图区域确定凝集中心在屏幕上的坐标，并在中心位置画点，置n=0;

%(3)根据凝结半径 确定自由离子释放半径r，，最大可凝集半径 ;

%(4)在半径r0的圆周的任意位置释放一粒子，让其作布朗运动;

%(5)粒子每移动一次，都要检测其上下左右的紧邻位置是否为凝集体，若是则凝集下来，n=n+1.执行第(6)步,否则，执行第(7)步;

%(6)每凝结一个粒子，都要计算此粒子与凝集中心的距离是否大于凝结半径 r0,

%若是则将凝结半径 置为当前距离，然后执行第4步;

%(7)若粒子与凝集中心的距离小于rmax，继续作布朗运动。否则抛弃此粒子，执行第4步;

%循环执行第3、4、5、6、7步，直到n=nsum为止.

function dla2d(nsum)

r0=2;

r=8;

Rr=0.5;

L=2000;

Xhalf=floor(L/2);

Yhalf=Xhalf;

szPoints=sparse(L,L);

szPoints(Xhalf,Xhalf)=1;

nx=Xhalf;

ny=Yhalf;

n=0;

theta=2*pi*rand;

x=nx+r*cos(theta);

y=ny+r*sin(theta);

nx=round(x);

ny=round(y);

i=0;

j=1;

p=-1:1;

N2(Xhalf/2)=0;

N2(1)=1;

R=2; %%用于计算生长维的一个变量

while nr+20

%theta=2*pi*rand;

nx=Xhalf+r*cos(theta);

ny=Yhalf+r*sin(theta);

elseif szPoints(nx-1,ny)+szPoints(nx+1,ny)+szPoints(nx,ny+1)+szPoints(nx,ny-1)>0&szPoints(nx,ny)==0

szPoints(nx,ny)=1;

n=n+1;

s=sqrt((nx-Xhalf)^2+(ny-Yhalf)^2);

k=round(s)+1;

N2(k)=N2(k)+1;

if s>R

R=s;

R1(j)=log(R/Rr);

N1(j)=n;

I(j)=i;

j=j+1;

end

if s>r

r=s+5;

end

elseif szPoints(nx,ny)==1

%theta=2*pi*rand;

nx=Xhalf+r*cos(theta);

ny=Yhalf+r*sin(theta);

end

end

[u,v]=find(szPoints);

figure(1);

plot(u,v,*);

%%%%%%%%%%%%%%%%%%%%%%%%%

%%分维采样数据分析整理

%%%%%%%%%%%%%%%%%%%%%%%%%

num=find(N2);

b=size(num,2);

for i=2:b

N2(i)=N2(i)+N2(i-1);

num(i)=N2(i);

R2(i)=log(2*i);

end

R2(1)=[];

num(1)=[];

save num2 num;

save R2 R2;

figure(2);

plot(R2,log(num),*);

%%%%%%%%%%%%%%%%%%%%%%%

%%生长维曲线绘制

%%%%%%%%%%%%%%%%%%%%%%%

figure(3); %画回旋半径与其内粒子数的对数关系图

plot(R1,log(N1),*);

%%%%%%%%%%%%%%%%%%%%%%%

%%粒子数量与运行步数的关系图绘制

%%%%%%%%%%%%%%%%%%%%%%%

%figure(4);%画凝聚粒子数量与运行步数的关系图

%%plot(N1,I,o);

展开阅读全文