%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% diffuse2d_exact_sol.m
%
% demonstration: Exact solution psi(x) of the density of random walkers
%                in 2D domain.  Solving the integral equation by power method.
%
% 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
Niter = 1000;

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);

%compute new transition probability in importance sampling
dr = [1 0; -1 0; 0 1; 0 -1];
K=ones(NGRID,NGRID,4)/4;   %old transition probability matrix
for i=1:NGRID,
   for j=1:NGRID,
      K(j,i,:)=K(j,i,:)*(1-q(j,i));
   end
end

%solve for psi(x,y)
psi=s;
for iter=1:Niter,
   psiK=zeros(size(psi));
   psiK1=psiK; psiK1(:,:)=psi.*K(:,:,1);
   psiK2=psiK; psiK2(:,:)=psi.*K(:,:,2);
   psiK3=psiK; psiK3(:,:)=psi.*K(:,:,3);
   psiK4=psiK; psiK4(:,:)=psi.*K(:,:,4);
   psiK=psiK1(:,[NGRID,1:NGRID-1])+psiK2(:,[2:NGRID,1]) ...
       +psiK3([NGRID,1:NGRID-1],:)+psiK4([2:NGRID,1],:) ;
   psi=s+psiK;
end

G=sum(sum(psi.*g));
disp(sprintf('exact solution G=psi*g= %e',G));

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

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

figure(2);
mesh(psi);

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

