% p1p2c1es3.m short for Plane 1, Plane2, Cylinder 1, escape zone strategy 3
%             using a minimally encapsulating sphere
% Class Example: See the hand-written notes L8 pages 15 and beyond

clear all, close all, format compact

diagnostics = true;  
diagnostics = false; % Comment out this statement to enable diagnostics
if (diagnostics)  disp('Diagnostic mode is on: <return> to continue.'), pause, end

% Geometry definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nr = 4; % Number of elemental subvolumes
z1 = 0; % P1, z = z1
z2 = 1; % P2; z = z2
R  = 1; % C1 circular cylinder x^2 + y^2 -R^2 = 0
L  = z2 - z1; 
xs = [0,0,L/2]; % center of sphere in 3 space
Rs = sqrt((L/2)^2 + R^2); % center of sphere S1: x^2 + y^2 + (z - L/2)^2 - Rs^2
rplane = R/(R + L); % probability that a particle emanating from the surface
                    % originates on one of the planes 
% Region 0: escape zone    : Outside S1
% Region 1: physical volume: Outside P1, inside P2, inside C1
% Region 2: boundary zone  : Inside  P1, inside S1
% Region 3: boundary zone  : Outside P1, inside P2, Outside C1
% Region 4: boundary zone  : Outside P2
% All zones except 1 should be considered to be filled with vacuum and not allow
% any particle deflections

% Tally set up %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xspot = []; % Tally vector used for plotting later
esc = zeros(Nr,1);; % Tallies escapes from any region

% Diagnostic tallies
vsrhofails = 0; % Volume source rho failure counter
vszfails   = 0; % Volume source z   failure counter
ssrhofails = 0; % Surface source rho failure counter
  
% Source type %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
choice = input('0 for surface source, 1 for volume source: ' );
if     (choice == 0) surfaceSource = true;  volumeSource  = false;
elseif (choice == 1) surfaceSource = false; volumeSource  = true;
else   disp('No choice for source type, exiting'), return, end

if     (surfaceSource) disp('Surface source is modeled')
elseif (volumeSource)  disp('Volume source is modeled'), end

% Scattering type %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
choice = input('0 for forward, 1 for backward, 2 for isotropic scattering: ' );
if (choice == 0)
   scatteringForward   = true;
   scatteringBackward  = false;
   scatteringIsotropic = false;
elseif (choice == 1)
   scatteringForward   = false;
   scatteringBackward  = true;
   scatteringIsotropic = false;
elseif (choice == 2)
   scatteringForward   = false;
   scatteringBackward  = false;
   scatteringIsotropic = true;
else disp('No choice for scattering type, exiting'), return, end
   
if (scatteringBackward)  disp('Backward scattering is modeled' ), end
if (scatteringForward)   disp('Forward scattering is modeled'  ), end
if (scatteringIsotropic) disp('Isotropic scattering is modeled'), end

% Simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nh = 1000; % Number of histories to simulate
for n = 1:Nh
   if (diagnostics) n, end
   % Sample the source =========================================================
   r = 1; % Source particle starts in region 1
   if (volumeSource)
      % Define the starting position vector -----------------------------------
      % Samping loop for x and y
      repeat = true; % This is a rejection loop for circular collimation
      while (repeat)
         x(1) = 1 - 2*rand;
         x(2) = 1 - 2*rand;
         if (x(1)^2 + x(2)^2 < 1) repeat = false; % Note the "<" logical operator 
         else vsrhofails = vsrhofails + 1; end % Tally the fails
      end
      x = R*x;
      % Samping loop for z
      repeat = true;
      while (repeat)
         x(3) = z1 + (z2 - z1)*rand;
         if (x(3) > z1 && x(3) < z2) repeat = false; % Note the "<" logical operator
         else vszfails = vszfails + 1; end % Tally the fails
      end
   elseif (surfaceSource)
      % Define the starting position vector -----------------------------------
      if (rand < rplane) % Particle emerges from a plane
         % Samping loop for x and y
         repeat = true; % This is a rejection loop for circular collimation
         while (repeat)
            x(1) = 1 - 2*rand;
            x(2) = 1 - 2*rand;
            if (x(1)^2 + x(2)^2 < 1) repeat = false; % Note the "<" logical operator               
            else ssrhofails = ssrhofails + 1; end
         end
         x = R*x;
         % Pick a plane, left or right
         if (rand < 0.5) % Left plane is chosen
            x(3) = z1;
            if (diagnostics) disp(sprintf('Starting on plane at z = %g', z1)), end
         else
            x(3) = z2;
            if (diagnostics) disp(sprintf('Starting on plane at z = %g', z2)), end
         end
      else % particle emerges from the surface
         phi = 2*pi*rand; x(1) = R*cos(phi); x(2) = R*sin(phi); x(3) = z1 + (z2 - z1)*rand;
         if (diagnostics)
            disp(sprintf('Starting on cylinder at x = [%g,%g,%g]',x(1),x(2),x(3)))
         end
      end
   end % Starting position now given

   xspot(end + 1,:) = x; % Tally the starting point

   % Define a random starting direction vector
   u = rotate([0,0,1],1 - 2*rand); % Compact use of the rotate function

   % Transport through the geometry ============================================
   while (r > 0 & r <= Nr) % True if particle is inside the geometry
      t0 = inf;            % Distance to interaction (in vacuum)

      if (diagnostics)
         str = 'Before transport r:x:u,t0: = %g:[%g,%g,%g]:[%g,%g,%g],%g';
         disp(sprintf(str,r,x(1),x(2),x(3),u(1),u(2),u(3),t0))
      end
      
      ri = r; % Record the starting region number for later
      if (r == 1) % Particle is in region 1
         [hit,t] = zplane(z1,x(3),u(3),false); if (hit & (t <= t0)) r = 2; t0 = t; end
         [hit,t] = zplane(z2,x(3),u(3),true ); if (hit & (t <= t0)) r = 4; t0 = t; end
         [hit,t] = ccylz(R,x,u,true         ); if (hit & (t <= t0)) r = 3; t0 = t; end
      elseif (r == 2) % Particle is in region 2
         [hit,t] = zplane(z1,x(3),u(3),true);
         if (hit & (t <= t0)) r = 1; t0 = t; end
         [hit,t] = sphere(Rs,xs,x,u,true);
         if (hit & (t <= t0)) r = 0; t0 = t; end
         if (r == 0)
            esc(ri,1) = esc(ri,1) + 1;   
            if (diagnostics) disp(sprintf('Particle escape from region %g',ri)); end
         end
      elseif (r == 4) % Particle is in region 4
         [hit,t] = zplane(z2,x(3),u(3),false); 
         if (hit & (t <= t0)) r = 1; t0 = t; end
         [hit,t] = sphere(Rs,xs,x,u,true);
         if (hit & (t <= t0)) r = 0; t0 = t; end
         if (r == 0)
            esc(ri,1) = esc(ri,1) + 1;   
            if (diagnostics) disp(sprintf('Particle escape from region %g',ri)); end
         end
      elseif (r == 3) % Particle is in region 3
         [hit,t] = zplane(z1,x(3),u(3),false); if (hit & (t <= t0)) r = 2; t0 = t; end
         [hit,t] = zplane(z2,x(3),u(3),true ); if (hit & (t <= t0)) r = 4; t0 = t; end
         [hit,t] = ccylz(R,x,u,false        ); if (hit & (t <= t0)) r = 1; t0 = t; end
         [hit,t] = sphere(Rs,xs,x,u,true);
         if (hit & (t <= t0)) r = 0; t0 = t; end
         if (r == 0)
            esc(ri,1) = esc(ri,1) + 1;   
            if (diagnostics) disp(sprintf('Particle escape from region %g',ri)); end
         end
      end
      x = x + u*t0; % Transport the particle
      if (diagnostics)
         str = 'After  transport r:x:u,t0: = %g:[%g,%g,%g]:[%g,%g,%g],%g';
         disp(sprintf(str,r,x(1),x(2),x(3),u(1),u(2),u(3),t0))
      end

      % Deflect the particle at the end of the transport step (should have no effect)
      if (ri == 1) % Deflections only in physical regions
         if     (scatteringForward  ) w =  1;
         elseif (scatteringBackward ) w = -1;
         elseif (scatteringIsotropic) w =  1 - 2*rand; end
         u = rotate(u,w); % Scatter the particle
         if (diagnostics)
            str = 'After  rotation  r:x:u =     %g:[%g,%g,%g]:[%g,%g,%g]';
            disp(sprintf(str,r,x(1),x(2),x(3),u(1),u(2),u(3)))
         end
      end
      xspot(end + 1,:) = x; % Record the intercept point

   end % End of the transport loop ============================================

   if (diagnostics) % Graphical output each history in diagnostic mode only
      figure(1), plot3(xspot(:,1),xspot(:,2),xspot(:,3),'k.') % View results
      xlabel('x'), ylabel('y'), zlabel('z')
      axis equal
      pause
   end

end % End of simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Outputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1), plot3(xspot(:,1),xspot(:,2),xspot(:,3),'k.') % View results
xlabel('x'), ylabel('y'), zlabel('z'), title('Full geometry')
axis equal

% Diagnostic tallies (always) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (volumeSource)
   f = vsrhofails/Nh;
   s = sqrt(f*(1-f)/Nh);
   t = 1-pi/4;
   t = t/(1-t); % Account for multiple misses
   disp(sprintf('Fraction of rho rejections = %g +/- %g compared to theory %g',f, s, t))
   disp(sprintf('Fraction of z rejections = %g compared to theory %g',vszfails/Nh, 0))
elseif (surfaceSource)
   f = ssrhofails/Nh;
   s = sqrt(f*(1-f)/Nh);
   t = 1-pi/4;
   t = t/(1-t)*rplane; % Account for multiple misses, note the rplane fraction
   disp(sprintf('Fraction of rho rejections = %g +/- %g compared to theory %g',f, s, t))
end

for r = 1:Nr disp(sprintf('Region %g had %g escapees',r,esc(r,1))), end

pause, close all