%
% Cellular automata, Greenberg-Hastings Automata
% MATH 499, Jae-Hun Jung
%


clear all, close all
N = 100;
W = zeros(N);
e =4;
W(N/2:N/2+2,N/2:N/2+2) = 1;

W(1,1) = 1;
W(1,end)=1;
W(end,1) = 1;
W(end,end)=1;

% for ix =1:N
%     for iy =1:N
%         if rand(1)>=0.5, W(ix,iy) = 1; end
%     end
% end
M = W; s = M;
imagesc(W),pause 
M = ones(N);

for t = 1:5000
    for iy=2:N-1
        s(1,iy) = sum(sum(W(1:2,iy-1:iy+1)))-W(1,iy);
        s(N,iy) = sum(sum(W(N-1:N,iy-1:iy+1)))-W(N,iy);
        s(iy,1) = sum(sum(W(iy-1:iy+1,1:2)))-W(iy,1);
        s(iy,N) = sum(sum(W(iy-1:iy+1,N-1:N)))-W(iy,N);
        for iz=2:N-1
            s(iy,iz) = sum(sum(W(iy-1:iy+1,iz-1:iz+1)))-W(iy,iz);
        end
    end
    s(1,1) = sum(sum(W(1:2,1:2)))-W(1,1);
    s(1,N) = sum(sum(W(1:2,N-1:N)))-W(1,N);
    s(N,1) = sum(sum(W(N-1:N,1:2)))-W(N,1);
    s(N,N) = sum(sum(W(N-1:N,N-1:N)))-W(N,N);

    p = 1;%0.4;
    
    fixindex=1; fixi(1) = 1;fixj(1) = 1;
    for ix = 1:N
        for iy = 1:N
            p = rand(1);
            prob = 1-(1-p)^s(ix,iy);
            pq=rand(1);
            if and(pq<=prob,W(ix,iy)==0) , W(ix,iy) = 1;
            elseif and(W(ix,iy)>0, W(ix,iy)<e), W(ix,iy) = W(ix,iy)+1;
            elseif W(ix,iy)==e ;
                W(ix,iy) = 0;
                fixindex = fixindex+1;
                fixi(fixindex) = ix; fixj(fixindex)=iy;
            else, W(ix,iy) = 0;
            end
        end
    end
    %M = W;
    for it = 1:fixindex
       W(fixi(it),fixj(it))=0;
       M(fixi(it),fixj(it))=0;
    end
%    W = W.*M; %------> Condition for "Once recovered, never infeced again"
    
    
    imagesc(W)
    pause(0.003)
end