% Cellular Automaton for 2D Precipitation Patterns with Aging
% Remove Comment Sign in Line 93 to Obtain Basic Model (No Aging)
% Oliver Steinbock and Bruno C. Batista (FSU);  February 9, 2023

clear all
close all

N=150;         % 2D lattice contains N x N cells
KAGE=10^-0.5;  % aging rate constant for C cells
DELTA=0;       % spatial slack when inserting deleted C cell
BMIN=2;        % number of B's needed to become target (1, 2 or 3)
ATTCH=0;       % new C required to be attached to old C (0=NO/1=YES)

r0=10;         % radius of initial B disk
rmaxC=0;       % distance of farthest C cell from domain center

rng(123);      % controlled random numbers by use of a seed

figure(1), set(gcf,'color','w')

A=zeros(N,N)+1; newA=A;         % A initially in entire domain
B=zeros(N,N);   newB=B;         % B is injected at constant rate
C=zeros(N,N);   newC=C;         % C is the precipitation product

t=0;                            % time step counter

[X,Y]=meshgrid(1:N,1:N);        % coordinates for next line use only
rC=(X-N/2).^2+(Y-N/2).^2;       % center distances to find farthest C cell

% --- create injection site at domain center ---
for y=1:N
    for x=1:N
        [theta,rho]=cart2pol(x-N/2,y-N/2);
        if rho<=r0
           A(y,x)=0; B(y,x)=1;  % small disk in middle contains B only
        end
    end
end

% --------------------------- iterate ---------------------------

while rmaxC<(N/2-2)^2
    
    % --- create composite color image and display pattern ---
    if mod(t,50)==0
        minC=min(C(C>0)); maxC=max(C(C>0)); dum1=zeros(size(C));
        dum1(C>0)=1-(C(C>0)-minC)/(maxC-minC);                 % scaled C>0  
        img1=cat(3,0.8*A+dum1,0.8*A+dum1*0.6+0.6*B,0.8*A+B);   % RGB
        imagesc((0.5:N-0.5),(0.5:N-0.5),img1)
        axis('square'), xlabel('x'), ylabel('y')      
        title(sprintf('t=%d, k=%2.3e, \\delta=%d',t,KAGE,DELTA))
        drawnow
    end

    t=t+1;

    % --- execute basic rule of C production when A and B close-by ---
    for y=2:N-1
        for x=2:N-1
            newA(y,x)=A(y,x);
            newB(y,x)=B(y,x);
            newC(y,x)=C(y,x);            
            sumA=sum(sum(A(y-1:y+1,x-1:x+1)));  % A count in 3x3 neighborh.
            sumB=sum(sum(B(y-1:y+1,x-1:x+1)));  % B count in 3x3 neighborh.
            if sumA>0 && sumB>0 && C(y,x)==0    % reaction to C
                newA(y,x)=0;
                newB(y,x)=0;
                newC(y,x)=1;
            end
            if C(y,x)>0
                newA(y,x)=0;
                newB(y,x)=0;
                newC(y,x)=C(y,x)+1;             % product ages
            end
        end
    end
    
    % --- randomly replace one fresh C cell with new B (injection) ---
    clear xtarget ytarget Ctarget
    m=0;
    for y=2:N-1
        for x=2:N-1
            if C(y,x)>0 && sum(sum(B(y-1:y+1,x-1:x+1)))>=BMIN  
                m=m+1;
                xtarget(m)=x;
                ytarget(m)=y;
                Ctarget(m)=C(y,x);
            end
        end
    end
    if m>0
        Cprobs=exp(-KAGE*(Ctarget-1.0));           % exponential aging
%       Cprobs=1.0+Ctarget*0.0;                    % no ageing (basic model)
%       Cprobs=(-atan(Ctarget-TAGE)+pi/2)/pi;      % sigmoidal aging (TAGE not defined)

        probs=Cprobs/sum(Cprobs);       % compute probabilities proportional to the values in Ctarget
        xtmp=mnrnd(1, probs);           % generate a random sample from the multinomial distribution with probabilities probs
        [~, ti]=max(xtmp);              % find the index of the nonzero element in xtmp
        newB(ytarget(ti),xtarget(ti))=1;
        newC(ytarget(ti),xtarget(ti))=0;
    
        
        % --- place the extracted C element into nearest available A site ---
        
        % the following is three-step process:
        
        % i. find the shortest distance dsqfinal
        dsqfinal=inf;
        for y=2:N-1
            for x=2:N-1
                if A(y,x)==1                                
                    dsq=sqrt((x-xtarget(ti))^2+(y-ytarget(ti))^2);
                    if dsq<dsqfinal
                       dsqfinal=dsq;
                    end
                end
            end
        end

        
        % ii. find all A cells with distances of dsqfinal+DELTA or less
        p=0;
        clear xp yp
        for y=2:N-1
            for x=2:N-1
                if A(y,x)==1 && sum(sum(C(y-1:y+1,x-1:x+1)))>=ATTCH  % attachment requirement
                   dsq=sqrt((x-xtarget(ti))^2+(y-ytarget(ti))^2);
                   if dsq<=dsqfinal+DELTA
                      p=p+1;
                      xp(p)=x;
                      yp(p)=y;
                    end
                end
            end
        end
        
        % iii. pick one of those cells and set it to C=1
        ppick=randi([1 p]);
        newA(yp(ppick),xp(ppick))=0;
        newB(yp(ppick),xp(ppick))=0; 
        newC(yp(ppick),xp(ppick))=1;
    end
    
    % --- refresh variables ---
    
    A=newA; B=newB; C=newC;
 
    % --- find most distant C to determine possible termination ---

    rmaxC = max(rC(C>0));
    
end


