% Random walk (RW) in a finite slab 
%  All lengths (Lz, R, ell, x0, y0) have to be given in the same unit: e.g in 1/k. 
% The outputs are: the Trajectory, the number of scattering events Nsc (0 for transmitted photons), the output angles theta, phi, the output positions x, y

clear all;
% Parameters of the slab
LzMax=1000; % length of the slab
R=2000; %transverse size
ell=20; %scattering mean free path


for iLz=1:1:round(LzMax/10) % loop on slab thickness in large step sizes
Lz=iLz*10; 
    
Tcoh(iLz)=0; % put counter for number of photons without scattering to 0
Tdiff(iLz)=0; % put counter for number of photons AROUND forward direction to 0
Rdiff(iLz)=0; % put counter for number of photons in backfard direction to 0  

Nph=1000; % number of photons (use Nph=10000 for a nice curve; 1000 photons takes 23seconds)
for iph=1:1:Nph
clear X Y Z RR S; %clear the parameters of the loop
% Initial Conditions 
Nsc=0;  % number of scattering events already done
cos_phi = 1; theta=0; % direction of inicident propagation
x=0; %initial position along x
y=0; %initial position along y
z=0; % Initial impact with the vapor (NOT the 1st scattering event!)
r=sqrt(x^2+y^2);

% length of the 1st step using an exponential distribution (Beer-Lambert)
step = -ell*log(rand);  % use the inverse function f^(-1) to draw the random number following a distribution f: exp -> log

% The first step is ballistic in forward direction   
z = z+step;

while all(z<=Lz  && z>0  && r<R) % Next scattering event inside the slab
  
    Nsc = Nsc+1;  % add 1 to number of scattering events
      
    % draw a new random direction
    theta = 2*pi*rand;
    cos_phi = -1+2*rand;  % uniform between -1 and 1 (always  check that the distribution is uniform: mistakes often happen with cos, sin etc)
    sin_phi = sqrt(1-cos_phi^2);  % always >0 since phi is between 0 and pi
    
    % length of the next step in an exponential ditribution (Beer-Lambert)
    step = -ell*log(rand);  

    % new position
    x = x + step*sin_phi*cos(theta);   
    y = y + step*sin_phi*sin(theta);
    z = z + step*cos_phi;
    r = sqrt(x^2+y^2);
  
end

if (z>Lz) 
    if (Nsc>1)
    Tdiff(iLz)=Tdiff(iLz)+1; 
    else
    Tcoh(iLz)=Tcoh(iLz)+1;
    end
else
    Rdiff(iLz)=Rdiff(iLz)+1;
end




end;

b(iLz)=Lz/ell;
Tcoh(iLz);
Tdiff(iLz);
Rdiff(iLz);
end;

figure(1),
plot(b, Tcoh, 'Color','black'); 
hold on; 
plot(b, Tdiff, 'Color', 'red'); 
hold on; 
plot(b, Rdiff, 'Color', 'blue');
xlabel('b=L/l');
h = legend('Tcoh','Tdiff','Rdiff',3, 'Location', 'East');
ylabel('Tcoh (black), Tdiff (red), Rdiff (blue)'); 

