%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% mcdiffuse2d.m
%
% demonstration: Monte Carlo simulation of reaction-diffusion in 2D
%                Compute the probability G of a walker (starting from source)
%                reaching the domain boundary before being absorbed.
%
% Wei Cai
% for LLNL CCMS Summer Institute Lectures Aug.1-3, 2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NGRID = 40; %number of grid points in both directions
r = 0.05;   %absorption rate inside domain
NWALKER = 1000;
NROUNDS = 10;

x=[1:NGRID]; y=[1:NGRID];
xx=ones(NGRID,1)*[1:NGRID];
yy=[1:NGRID]'*ones(1,NGRID);

xs=ceil(NGRID*0.25);
ys=ceil(NGRID*0.4);

%absorption function
q=max((abs(xx/NGRID-0.5)>0.42)|(abs(yy/NGRID-0.5)>0.35), r);

%source function
s=zeros(NGRID,NGRID);
s(ys,xs)=1;

%gain function (at domain boundary)
g=(q==1);

figure(1);
p1=contour(x,y,q,[0.5 0.5]*(r+1));
axis equal
hold on
plot(xs,ys,'ro');
hold off
xlim([0 NGRID]);
ylim([0 NGRID]);
xlabel('x');
ylabel('y');

%Start simulation

%number of visits
dr = [1 0; -1 0; 0 1; 0 -1];

G=zeros(NROUNDS,1);
for iter=1:NROUNDS,
 nvisit=zeros(NGRID,NGRID);
 for m=1:NWALKER,
    pos=[xs,ys];
    while 1,
      nvisit(pos(2),pos(1))=nvisit(pos(2),pos(1))+1;
      termprob = q(pos(2),pos(1));
      if(rand<termprob)
         %path terminated
         break;
      end
      pos=pos+dr(ceil(rand*4),:);
    end
 end
 G(iter)=sum(sum(nvisit.*g))/NWALKER;
 if(G(iter)>1)
    break;
 end
 
end
G
disp(sprintf('G = %e +/- %e',mean(G),std(G)/sqrt(NROUNDS)));

figure(1);
hold on
contour(nvisit/NWALKER,[0:0.01:0.2])
hold off

fs=18;
set(gca,'FontSize',fs);
xlabel('x');
ylabel('y');

figure(2);
mesh(nvisit/NWALKER);

set(gca,'FontSize',fs);
xlabel('x');
ylabel('y');
zlabel('\psi');

