% CHIRAL WALKER WITH REFRACTORY TAIL ON SQUARE LATTICE
% CW.m (04/20/09)
% This script computes the trajectory for straight initial conditions and
% for a single (L,A) pair
%
% For further information see Dutta&Steinbock, Int. J. Bif. Chaos (2010) 
% or go to http://www.chem.fsu.edu/steinbock/hp1.htm

clear all
close all

figure(1)

L=8;                                                    % length of refractory tail (in bonds)
A=4;                                                    % length of absolute refractory tail (in bonds), i.e. no crossing
N=L+1;                                                  % n bonds involves N=n+1 lattice points 
M=A+1;                                                  % m bonds involves M=M+1 lattice points
TEND=5000;                                              % total number of iterations to be performed

TRANSIENTCUT=50;                                        % transients must have died off after this time

NDISP=10;                                               % displayed area is (-NDISP...+NDISP) x (-NDISP...+NDISP)
DELAY=0.05;                                             % time delay between iterations in seconds

x=(1:N)*0;  xmemory(1,1:N)=x;                           % initial condition (straight walker)
y=(1:N)-1;  ymemory(1,1:N)=y;
t=0;

while (t <= TEND)                                       % iterate
        
    information=sprintf('N=%d, M=%d, t=%d', N, M, t);
    subplot(1,2,1)
    plot(x,y,'ob',x,y,'b-',x(1),y(1),'g*',x(1:M),y(1:M),'r-',x(1:M),y(1:M),'ro')
    title(information), axis('square'), axis ([-NDISP NDISP -NDISP NDISP]), drawnow
    if (t>TRANSIENTCUT)
        subplot(1,2,2), plot(xmemory(TRANSIENTCUT-1:t),ymemory(TRANSIENTCUT-1:t),'k-')
        axis('square'), axis ([-NDISP NDISP -NDISP NDISP]), drawnow
    end

    %waitforbuttonpress;                                % slow things down                     
    pause(DELAY);                                       % slow things down       

    t=t+1;
    
    xnew(2:N)=x(1:N-1);        ynew(2:N)=y(1:N-1);      % copying of main tail info dropping the last bond
    xnew(1)=x(1)+y(1)-y(2);    ynew(1)=y(1)-x(1)+x(2);  % left turn for new leading bond
    direct(t)=0;                                        % direct=0 means 'turned'
    
    % --- test whether the new leading bond is already occupied
    
    test=0;                                                         
    for n=1:N-1                                         % CHECKING WHETHER NEW LEADING BOND IS ALREADY OCCUPIED
       if (xnew(1)==x(n) & ynew(1)==y(n) & xnew(2)==x(n+1) & ynew(2)==y(n+1))
           test=1;
       end
       if (xnew(2)==x(n) & ynew(2)==y(n) & xnew(1)==x(n+1) & ynew(1)==y(n+1))
           test=1;
       end
    end
    for n=1:M                                           % CHECKING WHETHER NEW LEADING BOND CROSSES ABSOLUTE REFRACTORY TAIL
       if (xnew(1)==x(n) & ynew(1)==y(n))
           test=1;
       end
    end

    if (test==1)                                        % STRAIGHT MOVE NECESSARY B/C BOND IS OCCUPIED
        xnew(1)=x(1)+x(1)-x(2);
        ynew(1)=y(1)+y(1)-y(2);
        direct(t)=1;                                    % direct=1 means 'straight'
    end
        
    test=0;                                             % test whether the straight move lead the walker onto its own tail
    for n=1:N-1
       if (xnew(1)==x(n) & ynew(1)==y(n) & xnew(2)==x(n+1) & ynew(2)==y(n+1))
           test=1;
       end
       if (xnew(2)==x(n) & ynew(2)==y(n) & xnew(1)==x(n+1) & ynew(1)==y(n+1))
           test=1;
       end
    end
    for n=1:M
       if (xnew(1)==x(n) & ynew(1)==y(n))
           test=1;
       end
    end

    count=0;
    while(test==1)                                        % change the earlier step(s) from turn to straight if applicable
          count=count+1;
          t=t-1;
          x=xmemory(t,1:N); y=ymemory(t,1:N);             % go back t-th iteration
          xnew(2:N)=x(1:N-1);        ynew(2:N)=y(1:N-1);  % copying of main tail info dropping the last bond
          if (direct(t+1)==0)                             % ... and change to 'straight' if this was not the original move 
              xnew(1)=x(1)+x(1)-x(2);
              ynew(1)=y(1)+y(1)-y(2);
              direct(t+1)=1;                              % direct=1 means 'straight'
              test=0;
              for n=1:N-1
                if (xnew(1)==x(n) & ynew(1)==y(n) & xnew(2)==x(n+1) & ynew(2)==y(n+1))
                   test=1;
                end
                if (xnew(2)==x(n) & ynew(2)==y(n) & xnew(1)==x(n+1) & ynew(1)==y(n+1))
                   test=1;
                end
              end
              for n=1:M
               if (xnew(1)==x(n) & ynew(1)==y(n))
                test=1;
               end
              end
          end 
    end

    if (count~=0)
        t=t+1;
    end
        
    x=xnew;  xmemory(t,1:N)=x;    
    y=ynew;  ymemory(t,1:N)=y;
    
    % DETERMINE PERIOD AND TERMINATE SCRIPT
    if (t>TRANSIENTCUT)             
        diff=abs(x-xmemory(TRANSIENTCUT,1:N))+abs(y-ymemory(TRANSIENTCUT,1:N));
        if (diff==0)
           period=t-TRANSIENTCUT;
           fprintf('***\nL=%d, A=%d: PERIOD=%d\n***\n', L, A, period);
           t=TEND+1;                                        % terminates script through while loop
        end
    end  

end

