% -------------------------------------------------------------------------
% Zuse's Digitalteilchen, zweidimensional
% -------------------------------------------------------------------------
% Programmiersprache: MATLAB, getestet mit MATLAB 7.5.0 (R2007b)
%
% Autor: Dr. Uwe Renner
%        Arnold-Sommerfeld-Gesellschaft e.V., 2009
%
% für "Konrad-Zuse-Symposium 2009 - 40 Jahre Rechnender Raum"
% Beispiele stammen aus dem Buch von Konrad Zuse: Rechnender Raum.
% Vieweg, 1969
%
clear all;
 
bAvi = 0; %1: mit Avi-Video
 
aviobj = [];
if(1 == bAvi)
    aviobj = avifile('demo.avi', 'compression', 'rle', 'fps', 10);
end    
 
Zeitschritte = 103;  % Anzahl darzustellender Zeitschritte
N = 40;
 
NZustM = 4; % Anzahl der Monopod-Zustände
NZustB = 4; % Anzahl der Bipod-Zustände
 
%Definition der Nachbarschaft
%c = 1:N;        % Zentrum
s = [2:N 1];    % Süden
n = [N 1:N-1];  % Norden
e = [2:N 1];    % Osten
w = [N 1:N-1];  % Westen
 
% Monopods (nur die Pfeilspitzen dargestellt)
% 0: > (gehe nach Osten)
% 1: ^ (gehe nach Norden)
% 2: < (gehe nach Westen)
% 3: v (gehe nach Süden)
% NZustM*l + k + 1, k: Klasse (0...(NZustM-1)), l: Ordnung, 0: leer
M  = zeros(N, N, 'uint8'); 
M1 = false(N, N); M2 = false(N, N); M3 = false(N, N); M4 = false(N, N);
 
% Bipods (nur die Pfeilspitzen dargestellt)
% 0: ^>
% 1: <^
% 2: <v
% 3: v>
% NZustB*l + k + 1, k: Klasse (0...(NZustB-1)), l: Ordnung, 0: leer
B  = zeros(N, N, 'uint8'); 
B1 = false(N, N); B2 = false(N, N); B3 = false(N, N); B4 = false(N, N);
 
Nh = floor(N / 2);
 
% Anfangskonfigurationen
%
% Bild 43: Bipod
%B(Nh, Nh) = 1 + 3;
%
% Bild 44: Blinker, Nest 
%B(Nh, Nh) = 1 + 3; B(Nh + 1, Nh + 1) = 1 + 1;
%
% Bild 45: Nest
%B(Nh, Nh) = 1 + 3; B(Nh + 1, Nh + 1) = 1 + 1; 
%B(Nh + 1, Nh) = 1 + 0; B(Nh, Nh + 1) = 1 + 2;
%
% Bild 46: reaktive Kollision zweier Monopods zu einem Bipod
%M(Nh, Nh - 10) = 1 + 0; M(Nh - 10, Nh) = 1 + 3;
%
% Bild 52: reaktive Kollision zweier Bipods zu einem Monopod
%B(Nh - 10, Nh - 10) = 1 + 3; B(Nh - 10, Nh + 10) = 1 + 2;
%
% Bild 53: reaktionsfreie Kollision, Kreuzung zweier Bipods
%B(Nh - 10, Nh - 10) = 1 + 3; B(Nh - 10, Nh + 9) = 1 + 2;
%
%Bild 59: Bipod mit unterschiedlich langen Beinen
M(Nh,Nh) = 1 + 8 + 0; B(Nh,Nh) = 1 + 4 + 3;  % k = 3, l = 0
 
hfig = figure(1); clf(hfig);
set(hfig, 'DoubleBuffer','on', 'MenuBar', 'none', 'NumberTitle', 'off');
set(hfig, 'Name', 'Rechnender Raum und Digitalteilchen', 'Units', 'normalized');
ax = axes('Parent', hfig);
 
pos = get(ax, 'Position');
set(hfig, 'Position', [(1 - pos(3))/2,  (1 - pos(4))/2, pos(3), pos(4)]);
 
cmap = zeros(5,3);
cmap(2,:) = [1, 1, 1]; % Bipod
cmap(3,:) = [1, 0, 0]; % Monopod
cmap(4,:) = [0, 1, 0]; % zusammengesetzt
cmap(5,:) = [0, 0, 1];
colormap(cmap);
 
if(~isempty(aviobj))    
    aviobj = set(aviobj, 'colormap', cmap);
end
 
CData = uint8(B > 0) + 2 * uint8(M > 0); % Farbkodierung
 
hImage = image(CData, 'parent', ax);
axis(ax, 'off');
axis(ax, 'image')
 
drawnow;
 
Q = ones(10, 10, 'uint8'); % benötigt zur Pixelentwicklung im Kroneckerprodukt
 
% aktualisiere Video
if(~isempty(aviobj))
    F = im2frame(kron(CData, Q), cmap);
    aviobj = addframe(aviobj, F);
end
 
for i=1:Zeitschritte       
    % Komposition aus Monopods und Bipods
    
    % Monopods
    Id = (M == 0);
    k = mod(M - 1, NZustM);
    k(Id) = NZustM;
    l = 1 + int8((M - 1 - k) / NZustM);
    l(Id) = 0;
    
    % Entwickelung in  X, Y
    dXI = int8(k(:,w) == 0) .* l(:,w) - int8(k(:,e) == 2) .* l(:,e);
    dYI = int8(k(s,:) == 1) .* l(s,:) - int8(k(n,:) == 3) .* l(n,:);
    
    % Bipods
    Id = (B == 0);
    k = mod(B - 1, NZustB);
    k(Id) = NZustB;
    l = 1 + int8((B - 1 - k) / NZustB);
    l(Id) = 0;
    
    % Entwicklung in X, Y mit Zuse-Regeln für Bipods
    dXB = (int8(k(s,:) == 0) - int8(k(s,:) == 1)) .* l(s,:) + ...
          (int8(k(n,:) == 3) - int8(k(n,:) == 2)) .* l(n,:);
    dYB = (int8(k(:,w) == 0) - int8(k(:,w) == 3)) .* l(:,w) + ...
          (int8(k(:,e) == 1) - int8(k(:,e) == 2)) .* l(:,e);     
     
    dX = dXI + dXB;
    dY = dYI + dYB;
        
    % Normalisierung der Digitalteilchen in M, B
         
    % Kodierung der Bipods
    LB = min(abs(dX), abs(dY));
        
    B1 = (dX > 0) & (dY > 0);    
    B2 = (dX < 0) & (dY > 0);
    B3 = (dX < 0) & (dY < 0);
    B4 = (dX > 0) & (dY < 0);
    
    B = 4*uint8(LB) - (4*uint8(B1) + 3*uint8(B2) + 2*uint8(B3) + uint8(B4));
    B = B + uint8(B1 | B2 | B3 | B4);
        
    % Subtraktion Bipod
    dX = dX - LB .* (int8(B1 | B4) - int8(B2 | B3));
    dY = dY - LB .* (int8(B1 | B2) - int8(B3 | B4));
    
    % Kodierung Monopods   
    M1 = (dX > 0)  & (dY == 0); 
    M2 = (dX == 0) & (dY > 0);
    M3 = (dX < 0)  & (dY == 0);
    M4 = (dX == 0) & (dY < 0);
    
    R  = (int8(M1) - int8(M3)) .* dX  + (int8(M2) - int8(M4)) .* dY;
        
    M  = 4*uint8(R) - (4*uint8(M1) + 3*uint8(M2) + 2*uint8(M3) + uint8(M4));
    M  = M + uint8(M1 | M2 | M3 | M4);
                
    % Aktualisierung für Bild
    CData = uint8(B > 0) + 2 * uint8(M > 0); % Farbkodierung
    set(hImage, 'CData', CData);
    drawnow;
    
    % optionale Aktualisierung für Video
    if(~isempty(aviobj))
        F = im2frame(kron(CData, Q), cmap);        
        aviobj = addframe(aviobj, F);
    end
    
end
 
if(~isempty(aviobj))
    aviobj = close(aviobj);
end
