% This is to illustrate diffusion of particles (i.e., Brownian motion)
% in nd dimensions, and capture by a small inner region of radius rho
clear

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

nd = 2;  % set the dimension

N = 1000;  %number of particles
D = 1; % the diffusion coefficient
dt = .0005; % set the time step
L = 1; % radius of the outer  bounding circle
 
diff = sqrt(2*D*dt);
% find the capture time for a circle of radius rho, for several values of
% rho
K = 9;  % number of capture radii
dr = L/(K+1);

for k = 1:K
    rho = k*dr;
for j = 1:N %number of particles
   clear X 
  
   x = [1;zeros(nd-1,1)]; % initial position of the particles
X = x;

t = 0;
absx = x'*x;
    while(absx>rho^2 )
    dx = diff*randn(nd,1);
    xnew = x+dx;
    % check to make sure it is inside radius L
    xnR = xnew'*xnew;
    if (xnR > L^2) % project it back inside 
        scl=2*L/sqrt(xnR)-1;
         
        xnew = scl*xnew;
    end
    x=xnew;
    X = [X,x];
     t = t+dt;
     absx = x'*x;
    end
if (j==1)% make a plot showing one of the particles trajectory
    th = 2*[0:.01:1]*pi;
    figure(1)  % this is a 2d plot even though this might be a 3d simulation
    plot(X(1,:),X(2,:),L*sin(th),L*cos(th),rho*sin(th),rho*cos(th),'--','linewidth',2)
    axis([-L L -L L])
end

T(j) = t;
end
 

figure(2)
hist(T,40)
 xlabel('t','fontsize',20)
 ylabel('capture time distribution','fontsize',20)
 mT(k) = mean(T);
rhok(k) = rho;
end
% Question:  How does this compare to the analytical answer?
%
if(nd == 2)
mT_theor = L^2*(log(L./rhok)+(rhok.^2/L^2-1)/2)/(2*D);
end
if(nd==3)
    mT_theor = L^2*(L./rhok+rhok.^2/(2*L^2)-3/2)/(3*D);
end

figure(3)
plot(rhok,mT,'*',rhok,mT_theor,'linewidth',2)
ylabel('Mean Capture Time','fontsize',20)
xlabel('\rho','fontsize',20)