分形维数的计算

function DLAmodel（NsumWstep）             %定义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（LLL）;            %%网格点阵
szpoints（XhalfYhalfZhalf）=1;    %%种子位置标志
%%释放初始粒子    
（NsumWstep）             %定义dla函数，Nsum为所生成絮体包含的颗粒数，Wstep为计算过程中所采取的步长                            
theta=2*pi*rand;
gama=pi*rand;
M=p（1:）+Wstep*[cos（theta） sin（theta） cos（gama）];
while n    theta=2*pi*rand（1）;         %%粒子随机移动
    gama=pi*rand（1）;
    %Wstep=Wstep*[2*rand（1）-1 2*rand（1）-1 2*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）];
    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 s>release     %%调整释放半径
              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）];