%
%  2-D random walk code
%

nw = 20000;  				% number of walkers
ns =  200;				% number of steps taken
dr   = 0.04;                            % histogram binning parameter

close all

%
%  random walk (unit step in random direction)
%
pw = zeros(nw,1);			% initial position of walkers
for j1 = 1:ns
   angles = 2*pi*rand(nw,1);		% vector of nw random angles
   pw = pw + exp(i*angles);          % unit step in angle direction
end 

%
%  optional scatterplot of walkers
%	- if less than 501 walkers
%
if(nw<501)
   figure(1);  clf;  plot(pw(:,1),'o')
   xlabel('x-axis');     ylabel('y-axis');
   title('scatterplot of walkers');
   maxdist = max(max(abs(pw)))
   axis(maxdist*[-1 1 -1 1]); axis('square')
end

rpw = abs(pw(:,1));

disp(' ');  
disp(['average coordinates of walkers:   ' num2str(mean(pw))])
disp(['longest walk:   ' num2str(max(max(rpw)))])
disp(['mean displacement from origin:   ' num2str(mean(rpw))])
disp(' ');  

figure(2);   clf

%
%  normalized histogram of radial displacement
%
rbin = [dr/2:dr:3.0];
hpw  = hist(rpw/sqrt(ns),rbin);

bar(rbin,hpw/nw/dr,1);
h = get(gca,'children');  set(h,'facecolor','w','edgecolor','k')
hold on

xlabel('r-axis    [ r^2 = (x^2 + y^2) / (# of steps) ]');     
ylabel('probability density');
title('normalized histogram of radial displacement')
axis([0 3.0 0 1.0]);

rbin1 = [0 rbin];  plot(rbin1,2*rbin1.*exp(-rbin1.^2)); 

text(1.7,0.8,['# of walkers = ' num2str(nw)])
text(1.7,0.75,['# of steps    = ' num2str(ns)])