diffusion.m

森林火灾和气体扩散的matlab元胞自动机模拟(cellular Automata simulation)程序。

%diffusion + dla
clear all
clf

nx=200; %must be divisible by 4
ny=200;

z=zeros(nx,ny);
o=ones(nx,ny);
sand = z ;
sandNew = z;
sum = z;
gnd = z;
gnd(nx/2,ny/2) = 1 ;

sand = rand(nx,ny)<.1;

imh = image(cat(3,z,sand,gnd));
set(imh, 'erasemode', 'none')
axis equal
axis tight
 


for i=1:10000
    p=mod(i,2); %margolis neighborhood
   
    %upper left cell update
    xind = [1+p:2:nx-2+p];
    yind = [1+p:2:ny-2+p];
    %random velocity choice
    vary = rand(nx,ny)<.5 ;
    vary1 = 1-vary;
    
    %diffusion rule -- margolus neighborhood
    %rotate the 4 cells to randomize velocity
    sandNew(xind,yind) = ...
        vary(xind,yind).*sand(xind+1,yind) + ... %cw
        vary1(xind,yind).*sand(xind,yind+1) ;    %ccw
        
    sandNew(xind+1,yind) = ...
        vary(xind,yind).*sand(xind+1,yind+1) + ...
        vary1(xind,yind).*sand(xind,yind) ;
        
    sandNew(xind,yind+1) = ...    
        vary(xind,yind).*sand(xind,yind) + ...
        vary1(xind,yind).*sand(xind+1,yind+1) ;
        
     sandNew(xind+1,yind+1) = ...    
        vary(xind,yind).*sand(xind,yind+1) + ...
        vary1(xind,yind).*sand(xind+1,yind) ;
    
    sand = sandNew;
    
    %for every sand grain see if it near the fixed, sticky cluster
    sum(2:nx-1,2:ny-1) = gnd(2:nx-1,1:ny-2) + gnd(2:nx-1,3:ny) + ...
                        gnd(1:nx-2, 2:ny-1) + gnd(3:nx,2:ny-1) + ...
                        gnd(1:nx-2,1:ny-2) + gnd(1:nx-2,3:ny) + ...
                        gnd(3:nx,1:ny-2) + gnd(3:nx,3:ny);
    
    %add to the cluster
    gnd = ((sum>0) & (sand==1)) | gnd ;
    %and eliminate the moving particle
    sand(find(gnd==1)) = 0;
    
    set(imh, 'cdata', cat(3,gnd,gnd,(sand==1)) )
    drawnow
end