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

dt = (0.1)^2;  dx = sqrt(dt);
Nw = 10000;  Ns = 100*(10)^2;

home = 5.0;

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

%  take Ns timesteps, but check for out-of-bounds
for jj = 1:Ns
	rw  = rw + dx*sign(rand(size(rw))-0.5);
	rwh = find(rw>home);
	rw(rwh) = -2*Ns*dx;
	th(jj) = size(rwh,1);
end
	
%  histogram
db= 0.05;  dbb = 4*sqrt(Ns*dt)*db;
bins = 4*sqrt(Ns*dt)*(-1:db:1);

figure(1);  clf
subplot(2,1,1)
[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'])

text(home,max(num)/2,['  \bf absorption wall at x = ' num2str(home)])
plot([1 1]*home,[0 max(num)],'g--')

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);
gaussR = exp(-((xb-2*home).^2)/(4*sig*tt)) / sqrt(4*pi*sig*tt) * (Nw*dbb);

gaussH = max([gauss-gaussR ; 0*gauss]);

plot(xb,gaussH,'r')

subplot(2,1,2)
plot((1:Ns)*dt,cumsum(th),'.');  hold on
xlabel(['\bf time'])
ylabel(['\bf # at home'])

ti = (1:Ns)*dt;
thy = Nw*(1-erf(home./sqrt(4*sig*ti)));
plot(ti,thy,'r')