function forest();
%%% mcc -x forest for compile (on Linux)
clear all;
close all;

N=100;  % the grid size = N x N
th=0;

%a: 1
%q: 0
%r: 100

T = input('Enter # of interation (an integer:)');
if isempty(T)
T=100
end
p = input('Enter the growth rate (p):');
if isempty(p) 
p=0.01 
end
f = input('Enter the probability of the spontanious ignition(f):');
if isempty(f); 
f=0.01/500 
end

A=100*(rand(N,N)>0.5);

%A1=A;
%A=zeros(N,N);
%A(N/2,N/2)=1;

Ac=(A==1);
Qu=(A==0);
Rs=(A==100);

%plot([sum(sum(Ac)) sum(sum(Qu)) sum(sum(Rs)) p*sum(sum(Rs)) ]/N/N);
%axis([1 4 0 1])
%colormap(gray)
%image(Qu*50+Ac*100);

C=[0 1 0; 
    1 0 1; 
    0 1 0]; % 4 neighbor

for iii=1:T;

    A0=A;
Ac=(A0==1);
Ac=[Ac(N,N) Ac(N,:) Ac(N,1); Ac(:,N) Ac Ac(:,1) ;Ac(1, N) Ac(1,:) Ac(1,1)];
CoA=conv2(Ac,C,'full');
CoA=CoA(3:N+2,3:N+2);

for ii=1:N
    for jj=1:N
       
%A(N/2,N/2)=1; % fire in the center
%th=0+rand*2.2; % random fluctuation of th
        
        if A0(ii,jj)==1;   %active
            A(ii,jj)=100;
        elseif A0(ii,jj)==100; %rest
            A(ii,jj)=100*(rand>p); % 1 if rand >p
                               % 0 if rand <p 
        elseif A0(ii,jj)==0;      %quite
            if (CoA(ii,jj)>th)&(CoA(ii,jj)<5)
                A(ii,jj)=1;              

             else
                A(ii,jj)=(rand<f);
            end
            %else
            %pause;
         end
    end
 end
 
Ac=(A==1);
Qu=(A==0);
Rs=(A==100);

%plot(sum(sum(Qu))/N/N, sum(sum(Rs))/N/N);
%hold on

%axis([0 1 0 1])

%grid on;

    colormap(gray)
    image(Qu*50+Ac*100);
    pause(.1)
end
