%% skrypt realizuje automat komorkowy sandPile

close all;
clear all;
clc;

% bez tych trzech linijek animacja nie bedzie dzialac pod Windowsem:
opengl software;
figure(1)
set(gcf,'Renderer','OpenGL');
% dla animacji robionych pod Linuxem mozna te powyzsze linijki wykasowac


rozmiar = 50; % rozmiar planszy: 50x50

zycie = zeros(rozmiar,rozmiar);      % biezaca plansza zycia
liczbaCykli = 300;                   % ile iteracji zrobic

for n=1:liczbaCykli

    %--- KROK 1 ----------------------------------------------------
    % losowe ziarenka piasku 
    % z prawdopodobienstwem 0.3
    % gdzies na srodku u gory, +/- 3 komorki w prawo/lewo
    for k=(rozmiar/2-3):(rozmiar/2+3)
        if rand() < 0.3
            zycie(1,k) = 1;
        end
    end
    
    %--- KROK 2 ----------------------------------------------------
    % wyswietl cykl z tablicy zycie
    
    [I,J] = find(zycie);
    plot(J,rozmiar+2-I,'.','Marker', 's', 'MarkerFaceColor', 'r', 'MarkerSize', 4);
    axis([1 rozmiar+1 1 rozmiar+1]);
    daspect([1,1,1])    
    text(2,2.5,['cykl ' num2str(n)]);
        
    %--- KROK 3 ----------------------------------------------------
    % oblicz stany komorek do nastepnego cyklu
    % wykorzystamy zmienna n - w zaleznosci od tego, czy cykl ma wartosc
    % parzysta, czy tez nie, bedziemy odpowiednio dzielic siatke na kwadraty
    % wg sasiedztwa Margolusa
    
    % tworzymy plansze do odtwarzania zycia "za chwile"
    zyciePom= zeros(rozmiar,rozmiar); 

    for w=mod(n,2)+1:2:rozmiar-mod(n,2)
        
       if w==2
           % przepisac pierwszy rzadek i ostatni!!!!          
           for k=1:rozmiar
               zyciePom(1,k) = zycie(1,k);
               zyciePom(k,1) = zycie(k,1);
               zyciePom(rozmiar,k) = zycie(rozmiar,k);
               zyciePom(k,rozmiar) = zycie(k,rozmiar);
           end
       end
              
       for k=mod(n,2)+1:2:rozmiar-mod(n,2)
           % sytuacja nr 5 - zaklinowanie lub jednoczesne opadniecie
           if zycie(w,k)==1 && zycie(w,k+1)==1 && zycie(w+1,k)==0 && zycie(w+1,k+1)==0
              if rand()< 0.5    % jednoczesnie opadniecie
                  zyciePom(w,k)=0;
                  zyciePom(w,k+1)=0;
                  zyciePom(w+1,k)=1;
                  zyciePom(w+1,k+1)=1;
              else              % zaklinowanie
                  zyciePom(w,k)=1;
                  zyciePom(w,k+1)=1;
                  zyciePom(w+1,k)=0;
                  zyciePom(w+1,k+1)=0;                  
              end
           % sytuacja nr 3 - ziarenko na ziarenku, a obok pusto
           elseif zycie(w,k)==1 && zycie(w+1,k)==1 && zycie(w,k+1)==0 && zycie(w+1,k+1)==0
               zyciePom(w,k)=0;
               zyciePom(w,k+1)=0;
               zyciePom(w+1,k)=1;
               zyciePom(w+1,k+1)=1;           
           % sytuacja nr 4 - ziarenko na ziarenku, a obok pusto
           elseif zycie(w,k)==0 && zycie(w+1,k)==0 && zycie(w,k+1)==1 && zycie(w+1,k+1)==1
               zyciePom(w,k)=0;
               zyciePom(w,k+1)=0;
               zyciePom(w+1,k)=1;
               zyciePom(w+1,k+1)=1;                                        
           % sytuacja nr 1 - ziarenko ma pod spodem pusto, a obok i na ukos
           % cos jest 
           elseif zycie(w,k)==1 && zycie(w+1,k)==0
               zyciePom(w,k)=0;
               zyciePom(w+1,k)=1;
               zyciePom(w,k+1)=zycie(w,k+1);
               zyciePom(w+1,k+1)=zycie(w+1,k+1);
           % sytuacja nr 2 - ziarenko ma pod spodem pusto, a obok i na ukos
           % cos jest                
           elseif zycie(w,k+1)==1 && zycie(w+1,k+1)==0               
               zyciePom(w,k+1)=0;                   
               zyciePom(w+1,k+1)=1;
               zyciePom(w,k)=zycie(w,k);
               zyciePom(w+1,k)=zycie(w+1,k);               
           % w pozostalych przypadkach zostawiamy tak jak bylo
           else
               zyciePom(w,k)=zycie(w,k);
               zyciePom(w+1,k)=zycie(w+1,k);
               zyciePom(w,k+1)=zycie(w,k+1);
               zyciePom(w+1,k+1)=zycie(w+1,k+1);
           end
           
       end
    end
       
    %--- KROK 4 ----------------------------------------------------
    % przepisujemy wartosci z zyciePom do zycie
    
    % zerujemy znowu tablice zycie
    zycie= zeros(rozmiar,rozmiar);
    % gdzie w zyciePom sa jedynki?
    [I,J] = find(zyciePom);
    liczbaNieZer = length(I);           
    % przepisujemy jedynki z zycie do zyciePom
    for m=1:liczbaNieZer
        zycie(I(m),J(m)) = zyciePom(I(m),J(m));
    end    
    
    F(n) = getframe;
end

movie2avi(F, 'sandpile.avi', 'compression', 'None','fps',3,'quality',100)