function jforest2();

clear all;

N=100;  % forest is NxN
it=100; % # iterations

B=zeros(N);  %burning matrix: 1 = firing right now; 0 = not firing right now
A=zeros(N);  %activation matrix: nonnegative real numbers
        %represents the pre-firing membrane potential, or, if it's a tree, temperature
T=0.5i*ones(N);  %thresholds matrix: positive real numbers
R=zeros(N);  %random activations matrix: this gets filled in (randomly!) and added to A each epoch

t=5; %number turns spent in refractory period

d=.5; %decay constant: the amount by which A decays each turn, equal to exp(-1/tau)


%SET UP THE CONNECTIVITY MATRIX

c=5; %radius of the matrix
ci=3; %radius of inner activation
co=5; %radius of outer inhibition
ni=0; %noise on inner radius
no=0; %noise on outer radius

%Defines C, the connectivity matrix.
for ii=1:(2*c+1); for jj=1:(2*c+1); C(ii,jj)=2*(((ii-c-1)^2+(jj-c-1)^2)^0.5<ci+ni*rand(1))-(((ii-c-1)^2+(jj-c-1)^2)^0.5<co+no*rand(1)); end; end;
C=0.5*C;

figure(4)
colormap(gray)
image(50+50*C)

tic

%MAIN LOOP

for ii=1:it;

	%DECIDE WHO'S GOING TO FIRE THIS TIME
	
	%take the existing activation level divided by threshold
	%A./abs(T)
	B=ceil(0.5.*(  ones(N)  +  sign( (A./abs(T))-ones(N) )  ));
	%This albegra is equivalent to a heaviside function centerd at 1.


	%CALCULATE NEW ACTIVATION LEVELS
	
	%make new random inputs
	R=0.3*rand(N);
	
	%map new activations onto the activity matrix...
	A=d.*A + conv2(B,C,'same') + R;
	% decay the old,
	% add the incoming connections of cells that fired,
	% and add in random activations
	
	
	%UPDATE THRESHOLDS
	
	%update thresholds to place those who just fired in the refractory state
	T=T+100*(2^(t)).*B;
	%the threshold is shot into the stratosphere
	
	%remove one layer of delay
	T=T-0.5.*real(T)-0.5.*(rem(T,100));
	%and it descends by one half each turn
	
	%reset the activation of anyone who is in the refractory period
	A=A.*(ones(N)-sign(abs(real(T))));
	%sodium won't go in if potassium's already out!
	
	%rectify
	A=0.5*(sign(A).*A+A);
	%negative activity cannot be tolerated!  This algebra makes any negative cells equal to 0.


	%DRAW IT
	
	figure(1)
	colormap(gray)
	image(100*B);
	figure(2)
	colormap(gray)
	image(real(T));
	figure(3)
	colormap(gray)
	image(50*A+100*B);
	pause(1)

end

toc