% ### EXrandomWalk2D.m ###    11.16.14    {C. Bergevin}
% Calculates mean-squared distance for a group of independent 2-D random
% walkers, each moving with equal probablity  a step left or right AND up or down as:
% o Method 1 - simply equal probability one step up/down (U/D) and left/right (L/R)
% o Method 2 - U/D and L/R steps are sampled from a uniform distribution over [-1,1]
% o Method 3 - U/D and L/R steps are sampled from a Gaussian distribution
% [Source (modified): Kevin Berwick (re 'Computational Physics', Giordano & Nakanishi)]
clear;
% -------------
N= 200;      % Total # of (independent) walkers (each starts at x=0)
M= 100;      % Total # of steps for each walker
K= 3;       % # of walkers to show individual traces for [3]
method= 3;  % see comments above
% -------------
% +++
step_number= zeros(1,M);        %
posAvg= zeros(1,M);              % allocate array to stored (suquentially averaged) MSD
step_number_array= [1:1:M];     %
% +++
% NOTE: the loop is set up in such a way to average posAvg across walkers
for r= 1:N
    x=0; y=0;   % initialize positions for r'th walker
    positionX(r,1)= 0;  positionY(r,1)= 0;
    % loop to go through M steps for r'th walker
    for nn=1:M;
        if method==1
            if (rand<0.5),  x=x+1;  % conditional for left or right
            else    x=x-1;  end;
            if (rand<0.5),  y=y+1;  % conditional for up or down
            else    y=y-1;  end;
        elseif method==2
            x= x+2*rand(1)-1;   y= y+2*rand(1)-1;
        elseif method==3
            x= x+randn(1);   y= y+randn(1);
        end
        posAvg(nn)= posAvg(nn)+ (x^2+y^2);    % store squared displacement (handles averging across r)
        positionX(r,nn+1)= x;          % store displacments for each walker and step
        positionY(r,nn+1)= y;
    end;
end;
posAvg= posAvg/N;     % Divide by number of walkers
% plot MSD
figure(1);
plot(step_number_array,posAvg, 'k'); hold on;
title('MSD for 2-D random walk');
xlabel('Step number'); ylabel('Mean-Squared Distance (x^2+y^2)');
% plot a subset of individual traces
figure(2); clf; hold on; grid on;
for nn=1:K
    shade= 1-(nn-1)/K;
    plot(positionX(nn,:),positionY(nn,:),'Color',[1 1 1]-shade);
end
xlabel('x'); ylabel('y'); title('Representative traces');
% also plot MSD (which is a circular arc in this case)
rBND= sqrt(posAvg(end));  % radius MSD
xBND= linspace(-rBND,rBND,100); yBND= sqrt(rBND^2-xBND.^2);
plot(xBND,yBND,'r--','LineWidth',2); plot(xBND,-yBND,'r--','LineWidth',2);