% 1D Cellular Automata for recrystallization
% Site saturation Random nucleation mechanism
% Uniform distribution of driving pressure and mobility
% 
% Developed by M. S. Salehi

clc
clear

%% Solution parameters
N = 5000;                        % Cell number
b = 1;                           % Boundary width
NN = 25;                         % Number of nuclei 
Time = 500;                      % Annealing Time

M=1;                             % Mobility   
p=1;                             % Driving pressure
lc=1;                            % Cell size
Neighbor = [-1 1];               % Neighbor indices 
%% Initial condition

X = zeros(1,N+2*b);              % State matrix (0: deformed and 1: rexed)
step = 0;                        % CA step number
t=0;                             % Solution time

%% Nucleation

I = [];                          % Indices of nuclei
n=0;
while n<NN
    i=int32(rand(1)*N+b);        % select a random cell
    
    check=true;
    for m=1:numel(I)
        if I(m)==i
            check = false;
            break;
        end;
    end;
    if check
        n = n+1;
        I=[I,i];
        X(i)=n;
    end;
end;

%%
scrsz=get(0,'screensize');
figure1 = figure('position',scrsz);

X_new = X;                       % state in current step
XX = X;                          % history of states
N_Rex = NN;                      % Number of rexed cells
x=[0;0];                         % rex fraction and time
while t<Time
    
    % Boundary condition (periodic)
    for n=1:b
        X_new(n) = X(N+n);
        X_new(N+b+n) = X(b+n);
    end;
    
    % Solution
    for i=b+1:N+b
        for n=1:numel(Neighbor) % check neighbors
            if X(i)==0&&X(i+Neighbor(n))~=0
                X_new(i)=X(i+Neighbor(n));
                N_Rex = N_Rex+1;
                break;
            end;
        end;
    end;
    
    X = X_new;                   % Update state
    XX=[XX;X];
    
    v = M*p;
    dt = lc/v;
    t=t+dt                       % Update time
    
    x=[x,[N_Rex/N;t]];           % Update rex fraction
    
    %% visualize the results
    % Create subplot
    subplot1 = subplot(1,2,1,'Parent',figure1);
    box(subplot1,'on');
    imagesc(XX,'Parent',subplot1)
    daspect([N Time 1])
    xlabel('position');
    ylabel('CA Steps');
    colormap('hot');
    title({'Black cell: Deformed, Colored: Rexed grains'},'Visible','on');

    % Create subplot
    subplot2 = subplot(1,2,2,'Parent',figure1);
    box(subplot2,'on');
    plot(log10(x(2,:)),x(1,:),'Marker','o','LineWidth',2,'Parent',subplot2)  
    daspect([log10(Time) 1 1])
    xlim(subplot2,[0 log10(Time)]);
    ylim(subplot2,[0 1.05]);
    xlabel('Log(time)');
    ylabel('recrystallization fraction');

    pause(0.05)
end