%  
%this code determines the distribution of decay positions of
%diffusing/decaying particles in 2dimensions
clear all
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7); 
D = 1; % diffusion coefficient

alpha = 1;  % decay rate 
 
dt = .0001;  %time step
diff = sqrt(2*D*dt); % diffusion step size
 
N = 1000;  %number of particles
nd = 0;  % number of decayed particles so far
X = zeros(2,N);
for j = 1:N % for each particle
    
    x = [0;0];
 
   decay = 0;  % this particle has not yet decayed
       
    while (decay ==0) % step until particle decays
        x =  x+diff*randn(2,1); % move by 2 d diffusion
        if(alpha*dt>rand(1,1))
            decay = 1; %this particle has decayed
            % record its current position
            X(:,j) = x;
           
        end
        
    end
end

figure(1)
plot(X(1,:),X(2,:),'*')
 % make a histogram of the resting positions
R = sqrt(X(1,:).^2+X(2,:).^2);
rmx = max(R);
dr = rmx/50;
Rn = fix(R/dr)+1;
r= [1:51]'*dr;
Rhist = zeros(51,1);

for j = 1:N
    Rhist(Rn(j))= Rhist(Rn(j))+1;
end
mu = sqrt(alpha/D);
figure(2)
plot(r,Rhist./(r*N),'*',r,mu*besselk(0,mu*r)/(2*pi),'--','linewidth',2)
xlabel('r','fontsize',20)
ylabel('v(r)','fontsize',20)

 
 