%  
%this code determines the distribution of decay positions of
%diffusing/decaying particles
clear all
 
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);
cmykflg=0;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;


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
for j = 1:N % for each particle
    
    x = 0;
   decay = 0;  % this particle has not yet decayed
       
    while (decay ==0) % step until particle decays
        x =  x+diff*randn(1,1); % move by diffusion
        if(alpha*dt>rand(1,1))
            decay = 1; %this particle has decayed
            % record its position
            X(j) = x;
           
        end
        
    end
end

 % make a histogram of the resting positions
[ND,Y]= hist(X,50);
dY = Y(2)-Y(1);
figure(1)
 
scl = sqrt(alpha/D);
bar(Y,ND/(N*dY))
% to add a comparison with the analytical answer do the following
hold on
plot(Y,scl*exp(-scl*abs(Y))/2,'--','linewidth',2)
xlabel('x','fontsize',20)
ylabel('Distribution of Particles','fontsize',20)
set(gca,'Xtick',-6:2:6)
 
 hold off
 
  if(cmykflg==1)
print('../../figs_c/chapt_6/prtcle_dk_hist','-deps','-cmyk')
  end

 