% zplanes3.m Class Example - 3 parallel planes, non-interacting beam,
% a very simple implementation of multiple planes

clear all, close all

% 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

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

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

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

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

% Simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nh = 200; % 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 if particle is inside the geometry
      t0 = inf;   % Start with an impossibly big transport distance
      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
      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
      end
      x = x + u*t0;         % Transport the particle
      xspot(end + 1,:) = x; % Record the intercept point
   end % End of the transport loop
   plot3(xspot(:,1),xspot(:,2),xspot(:,3),'k.') % View results
   pause(0.001) % The pause enables a dynamic view as points are added

end % End of the history loop