% 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=10; % length of the slab
Lz=LzMax;
R=20000; %transverse size
ell=1; %scattering mean free path
Tdiff=0;


Nph=10000; % number of photons (use Nph=10000 - 6 secondes computing - or more for a nice curve)
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

if (z>Lz) 
    Tdiff=Tdiff+1; %just to count diffuse transmission 
    RTtime(iph)=Nsc; % register number of events, which we put to time
end


end;



Tdiff
ps=hist(RTtime, 100); %histrgram of transit times
ps(1)=0; %remove points which have not been scattered
plot(ps(1,1:100),'DisplayName','ps(1,1:100)','YDataSource','ps(1,1:100)');figure(gcf)

