%  
%this code determines the probability of decay vs. exit at r=R, starting
%from r=0, in a 3d spherical
clear all

R = 1; % absorb at r=R
D = 1; % diffusion coefficient

Ng = 50;% number of gamma values to simulate
for g = 1:Ng
gamma = 0.2*g; %decay rate
gm(g) = gamma;
dt = .001;  %time step
diff = sqrt(2*D*dt); % diffusion step size
 
N = 1000;  %number of particles
nd = 0;  % number of decayed particles
for j = 1:N % for each particle
    
    x = [0;0;0]; % a 3-vector
    r = sum(x.^2);
   decay = 0;  %this particle has not yet decayed
        
    while ((r<R^2)& (decay ==0)) % step until x reaches L or particle decays
        x =   x+diff*randn(3,1); % move by 3 d diffusion
        r = sum(x.^2);% the new radius
        if(gamma*dt>rand(1,1))
            decay = 1; %this particle has decayed
            nd = nd +1; % count the number that decay
        end
        
    end
end

 pd(g) =  1-nd/N ;  % the percentage that reached r=R
end
mu = gm*R^2/D;
plot(mu,pd,'*',mu,sqrt(mu)./sinh(sqrt(mu)),'--','linewidth',2)
xlabel('\alpha R^2/D','fontsize',20)
ylabel('Probability of reaching r = R','fontsize',20)