clear all
n=500; %size
pfire=0.000005; %chance of 'fire' starting
pgrow=0.005; %chance of 'tree' growing
%initialize forest
empty=zeros(n);
forest=empty;
burn=empty;

figure(1);
imh = image(cat(3,empty,empty,empty));
drawnow
% The RG in RGB is useful
% Red=1--> Burning, Green=2--> Tree, and Black=0-->Empty
for i=1:1000
    r=rand(n);
    %creates adjcency matrix showing presence of fires with periodic
    %boundary conditions
    burn=(forest(:,[n 1:n-1])==1) + (forest(:,[2:n 1])==1) + (forest([n 1:n-1], 1:n)==1) + (forest([2:n 1],1:n)==1);
    %update the forest by subtracting burned trees and adding new trees.
    %Matlab converts logical arrays to numeric arrays on the fly
    forest = 2*(forest==2) - and((forest==2),(burn>0 | (r<pfire))) + 2*((forest==0) & r<pgrow);
    tempforest=int8(forest==2);
    pc=sum(tempforest(:))/n^2;
    set(imh, 'cdata', cat(3,(forest==1),(forest==2),empty))
    title(sprintf('pc=%.2f',pc))
    drawnow
    saveas(gcf, ['frame_',sprintf('%03d',i),'.png'])
end
