%
%  w08rand_wk.m - 1d random walk (djm: 24 oct 2010)
%   

dt = (0.1/2)^2;  dx = sqrt(dt);
Nw = 8*10000;  Ns = 5*(2*10)^2;

%  initialize Nw random walkers at x=0
rw = zeros(Nw,1);

%  take Ns timesteps
for jj = 1:Ns
	rw = rw + dx*sign(rand(size(rw))-0.5);	
end
	
%  histogram
db= 0.05;  dbb = 4*sqrt(Ns*dt)*db;
bins = 4*sqrt(Ns*dt)*(-1:db:1);

figure(1);  clf
[num,cent] = histc(rw,bins-(dbb/2));
plot(bins,num,'b*');  hold on
title(['\bf position histogram of ' num2str(Nw) ' walkers'])
ylabel('\bf # of walkers')
xlabel(['\bf locations after ' num2str(Ns) ' steps'])

sig = dx^2/dt/2;
tt = Ns*dt;
xb = 4*sqrt(Ns*dt)*(-1:db/8:1);
gauss = exp(-(xb.^2)/(4*sig*tt)) / sqrt(4*pi*sig*tt) * (Nw*dbb);

plot(xb,gauss,'r')