% zplanes3a.m Class Example - 3 parallel planes, 1 cylinder, interacting beam
% a very simple implementation of multiple planes
% This is an adaptation of zplane3.m adding a cylinder and interactions

clear all, close all

% Source definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This example is a beam incident from the left, normal incidence
Rbeam = 1; % Radius of collimation of the beam

% Geometry definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nr = 2; % Number of elemental subvolumes

% Region 0 (escape zone) definition, bounded by:
z1 = 0; % z = z1; z-plane

% Region 1 definition, bounded by z1 defined above plus
z2 = 1; % z = z2; z-plane at z2
r1 = 2* Rbeam; % Cylinder with radius 2x that of the beam

% Region 2 definition, bounded by z2 and r1 defined above plus
z3 = 2; % z = z3; z-plane at z3

% Region 3 definition (escape zone), bounded by z3 defined above

% Region 4 definition, bounded by (escape zone)
% r1 defined above

% Physics set up %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

S = linspace(0.1,10,Nr); % Macroscopic Sigma in the exponential interaction law
% Note that it is region dependent and hence we are modelling an inhomogeneous
% object. Region 1 will have little scattering, Region 2 will have lots.
% For Nr = 2, S = [0.1,10];

% Tally set up %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xspot = []; % Tally vector used for plotting later

% Simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nh = 500; % Number of histories to simulate
for n = 1:Nh

   % Sample the source
   r = 1; % Source particle starts in region 1
   % Define the starting position vector 
   repeat = true; % This is a rejection loop for circular collimation
   while (repeat)
      x(1) = Rbeam*(2*rand - 1);
      x(2) = Rbeam*(2*rand - 1);
      if (x(1)^2 + x(2)^2 <= Rbeam^2) repeat = false; end
   end
   x(3) = z1;            % Starts on leftmost plane
   u = [0, 0, 1];        % Normal incidence
   xspot(end + 1,:) = x; % Record the starting point
   
   % Transport through the geometry
   
   while (r > 0 & r <= Nr) % True of particle is inside the geometry
   t0 = -log(rand)/S(r); % Distance to an interaction, note region dependence
   ri = r; % Current region
   
      % zplane() and ccylz() are NERS544 library routines
      
      if (r == 1) % Particle is in region 1
         [hit,t] = zplane(z1,x(3),u(3),false); % Outside z1
         if (hit & (t <= t0)) r = 0; t0 = t; end
         [hit,t] = zplane(z2,x(3),u(3),true); % Inside z2
         if (hit & (t <= t0)) r = 2; t0 = t; end
         [hit,t] = ccylz(r1,x,u,true); % Inside r1
         if (hit & (t <= t0)) r = 4; t0 = t; end

      elseif (r == 2) % Particle is in region 2
         [hit,t] = zplane(z2,x(3),u(3),false); % Outside z2
         if (hit & (t <= t0)) r = 1; t0 = t; end
         [hit,t] = zplane(z3,x(3),u(3),true); % Inside z3
         if (hit & (t <= t0)) r = 3; t0 = t; end
         [hit,t] = ccylz(r1,x,u,true); % Inside r1
         if (hit & (t <= t0)) r = 4; t0 = t; end
      end
      
      x = x + u*t0; % Transport the particle
      if (r == ri) % The particle stays in the same region thus must interact
         ws = 2*rand - 1; % Isotropic scattering
         u = rotate(u,ws); % NERS544 library routine
      end
         
      xspot(end + 1,:) = x; % Record the intercept or interaction point
     
   end 
      
   pause(0.001) % Matlab has a multithread error causing me to do this
              % Plotting starts before the data is written!
   plot3(xspot(:,1),xspot(:,2),xspot(:,3),'k.') % View results
   
end % End of simulation