%
%  w08rand_2d.m - 2d random walk (djm: 24 oct 2010)
%   

clear

%  domain
Nd = 40;  nn = -Nd:2:Nd;
prob = zeros(2*Nd+1,2*Nd+1);
prob(Nd+1,Nd+1) = 1.0;

%  iterations
Ns = (Nd^2)/2;

fr = 1:2*Nd;  ba = fr+1;

for jj = 1:Ns
	temp = zeros(size(prob));
	temp(fr,fr) =               3*prob(ba,ba);
	temp(ba,ba) = temp(ba,ba) + 3*prob(fr,fr);
	temp(fr,ba) = temp(fr,ba) +   prob(ba,fr);
	temp(ba,fr) = temp(ba,fr) +   prob(fr,ba);

	prob = temp/8;
end

figure(1);  clf
surf(nn,nn,prob(1:2:end,1:2:end));  
axis([-Nd Nd -Nd Nd 0 1.1*max(prob(:))])

title('\bf probability density, p(x,y)')

figure(2);  clf
contour(nn,nn,prob(1:2:end,1:2:end));  hold on
axis([-Nd Nd -Nd Nd]);  axis square

title('\bf probability density, p(x,y)')

%  a few random walks
Nw = 100;
rw = zeros(Nw,2);

for jj = 1:Ns
	step1 = sign(rand(Nw,1)-0.25);
	step2 = sign(rand(Nw,1)-0.50);

	step = [step2 , step2.*step1];
	
	rw = rw + step;
end

plot(rw(:,1),rw(:,2),'k.')