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=2000;
p = input('Enter the growth rate (p): ');
if isempty(p) 
	p=0.06 
end
f = input('Enter the probability of the spontanious ignition(f): ');
if isempty(f)
	f=.00006 
end
I = input('Enter the # of iterations between display updates: ');
if isempty(I) 
	I=50 
end

%D records forest density over time
D=zeros(1,T);

%E records the number of fires
E=zeros(1,T);

%A=forest;  1=empty, 0=green tree, -1=fire
%S=starting up fires;  0=no lightning strike, 1=lightning strike
%B=burning trees;  0=no fire, 1=fire;  (extrapolates fire spread)
%G=new growth

A=(rand(N,N)>.35);
B=zeros(N);  %Initially nothing is burning

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

figure;
figure;

g=get(2,'Position');
set(1,'Position',get(2,'Position')+[0 -g(4)-50 0 0]);

%Begin Loop
for iii=1:T

	G=rand(N)>(1-p);  %Grow new trees at random locations
	S=rand(N)>(1-f);  %Add lightning strikes at random locations

	%eliminate burned trees from forest
	A=abs(A);  %changes fires to empty squares

	%propagate burning and add in new fires (independently of where there are trees!)
	%b is a temporary variable.  0 = no fire; n>0 = fire
	b=(conv2(B,C,'full'));b=b(2:N+1,2:N+1)+S;
	B=fix(  b ./ (b-.000001*ones(N))  );     %This makes all entries either 1 or 0.

	%Add fires and growth.  Right now forest has only 0's (trees) and 1's (empties)
	A= A - (ones(N)-A).*B  - A.*G;
		%subtract 1 from entries with 0 in A (a tree) and 1 in B (fire)
		%and subtract 1 from entries with 1 in A (empty) and 1 in G (new growth)

	%Update the burning tree matrix
	B=0.5*(abs(A)-A);  %1 -> 0; 0 -> 0; -1 -> 1

	%output density, number of fires
	D(1,iii)=1-nnz(A)/(N*N);
	E(1,iii)=30*(nnz(B)/(N*N));

	if rem(iii,I)==0

		%Draw it
		%figure(1)
	  colormap(gray)
		image(10*(A.*ones(N))-20*(A-abs(A))+60*S);
 	 pause(.01)

		%figure(2)
		%plot(1:T,D,1:T,E)

	end

%End Loop
end

