森林火灾和气体扩散的matlab元胞自动机模拟代码.zip

%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