% Self-avoiding Random walk 
% Illustrates polymer grow in 2D on square lattice
% January 2021 by AG.

clear  % Clears command window
clc    % Clears command history
clf    % Removes anything in the figure window before simulation. 
%tic    % set timer for the program

%rng('default')      % reset RNG
rng('shuffle')      % use seed as current time

Nmap   = 1000;      % area Nmap*Nmap 

L = 15;             % size of the plot ais -L to + L in x and y
amax = 0 + L;
amin = 0 - L;

% preparation
x = Nmap/2;         % start in the middle of the map
y = Nmap/2;

xi(1) = 0;          % origin for the plot
yi(1) = 0;
x0(1) = 0;          % mark the origin
y0(1) = 0;

Map = zeros(Nmap);  % initializing the map
Map(x,y) = 1;       % first point is already taken

out = 0;            % possible to walk 
monomer = 0;        % initial size

ip = 1;
k  = 1;
% end of preparation

while out == 0
    xold = x;
    yold = y;
    direction = randi(4);
    switch direction
        case 1
            x = x + 1;
        case 2
            x = x - 1;
        case 3
            y = y + 1;
        case 4
            y = y - 1;
    end
%   check if the place is vacant
    if Map(x,y) == 0
        Map(x,y) = 1;
        monomer = monomer +1;
        k = k+1;
        xi(k) = x-Nmap/2;
        yi(k) = y-Nmap/2;
        x0(k) = 0;
        y0(k) = 0;
    else
        x = xold;
        y = yold;
    end
%   check if we must stop (nowhere to go) or out of boundaries
    if x==1 || x==Nmap || y==1 || y==Nmap
        break
    end
%   condition 2: trapped - no more places around new x,y where to go
    if Map(x+1,y)*Map(x,y+1)*Map(x-1,y)*Map(x,y-1) == 1
        out = 1;
    end
    
    plot(xi,yi,'r-o',x0,y0,'b-x')
    axis([amin amax amin amax]);
    xlabel('x position')
    ylabel('y position')
    title ('Polymer simulation')
    grid;
    FR(ip)=getframe; 
    ip = ip +1;    
end

fprintf(' Lenght     = %4i   \n',monomer)
%toc