% Uri Wilensky
% http://web.eecs.utk.edu/~mclennan/Classes/527/NetLogo/B-Z%20Reaction.html
%   http://web.eecs.utk.edu/courses/fall2010/cs527/NetLogo/B-Z%20Reaction.nlogo
%   http://web.eecs.utk.edu/~mclennan/Classes/420-594-F07/NetLogo/B-Z-Reaction.nlogo
% http://ccl.northwestern.edu/netlogo/models/B-ZReaction
%(1) If the cell is healthy (i.e., in state 0) then its new state is [a/k1] + [b/k2], where a is the number of infected cells
%among its eight neighbors, b is the number of ill cells among its neighbors, and k1 and k2 are constants. Here "[]" means
%the integer part of the number enclosed, so that, for example, [7/3] = [2+1/3] = 2.
%(2) If the cell is ill (i.e., in state n) then it miraculously becomes healthy (i.e., its state becomes 0).
%(3) If the cell is infected (i.e., in a state other than 0 and n) then its new state is [s/(a+b+1)] + g, where a and b are
%as above, s is the sum of the states of the cell and of its neighbors and g is a constant. 
% If resulted cell state > n then its state become n.

% cell state: [z;m]
% a - neighbors count (among 8) in state (z;m)
% b - neighbors count (among 8) in state m
% s - states sum of the current cell and of its 8 neighbors
%|current cell state | transform to state|
%|          m        |           z       |
%|          z        |     [a/k1]+[b/k2] |
%|        (z;m)      |     [s/(a+b+1)]+g |

clear all;
clf;
z=0; m=256;
k1=2; k2=2; g=100;
rN=256; cN=rN;
iter=500;
c1 = round((z-0.5)+(m-z+1)*rand(rN,cN));
c2 = zeros(rN,cN);
ih=imagesc(c1);
colormap(hot);
th=title('titileText');
axis image off;
set(gcf,'Position',[1300 500 rN+40 rN+40]);
set(gca,'Units','pixels','ActivePositionProperty','Position','Position',[20 20 rN rN]);
figure(1);
for cycle=1:iter
tic;
  for x=1:rN
    for y=1:cN
      if c1(x,y)==m  % (2)
          c2(x,y)=z; % (2)
      else
        s=0; a=0; b=0;
        for i=x-1:x+1
          for j=y-1:y+1
              cc=c1(mod(i+rN-1,rN)+1,mod(j+cN-1,cN)+1); % current cell
              s=s+cc;
              if (i~=x || j~=y) % exclude center cell
                  if  cc == m 
                      b=b+1;
                  elseif cc > z 
                      a=a+1;
                  end
              end
          end
        end
        if c1(x,y)==z 
            c2(x,y)=floor(a/k1)+floor(b/k2); % (1)
        else
            c2(x,y)=floor(s/(a+b+1)) + g; % (3)
            if c2(x,y)>m c2(x,y)=m; end   % (3)
        end
      end
     end
  end
  c1=c2; 
  set(ih,'Cdata',c1);
  set(th,'string',sprintf('i=%d dt=%f',cycle,toc));
  drawnow
end

% http://www.hermetic.ch/pca/algorithms.htm#bz1 ()
% slightly modified version: http://www.hermetic.ch/pca/algorithms.htm#bz1
% (i) Select an integer q in the range 2 through 255. Cells may be in any of the states 1 through q.
% (ii) Select two integers k1 and k2 in the range 1 through 8 and an integer g in the range 0 through 100.
% (iii) In the transition from one "step" to the next the state of each cell is changed once according to rules (iv) - (vii) below:
% (iv) A cell in state q changes to state 1.
% (v) A cell in state 1 changes to state a/k1 + b/k2 + 1 where a is the number of neighbors of the cell which are in states 2 through q-1 and b is the number of neighbors in state q.
% (vi) A cell in any of states 2 through q-1 changes to S/(9 - c) + g, where S is the sum of the states of the cell and its neighbors and c is the number of neighbors in state 1.
% (vii) If the application of rule (v) or rule (vi) would result in a cell having a state > q then the state of that cell becomes q. 
