% p1p2c1es1.m short for Plane 1, Plane2, Cylinder 1, escape zone strategy 1
% Class Example: See the hand-written notes L8 pages 1 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Region 1: physical volume: Outside P1, inside P2, inside C1
% Region 0: escape region  : everything not in Region 1
Nr = 1; % 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;
rplane = R/(R + L); % probability that a particle emanating from the surface
                    % originates on one of the planes.

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

% 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
      
      if (r == 1) % If particle is in region 1
         [hit,t] = zplane(z1,x(3),u(3),false); if (hit & (t <= t0)) r = 0; t0 = t; end
         [hit,t] = zplane(z2,x(3),u(3),true ); if (hit & (t <= t0)) r = 0; t0 = t; end
         [hit,t] = ccylz(R,x,u,true         ); if (hit & (t <= t0)) r = 0; t0 = t; 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     (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
      xspot(end + 1,:) = x; % Record the intercept point

      if (diagnostics) pause, end
      
   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([-2,2,-2,2,0,1])
      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([-2,2,-2,2,0,1])

figure(2)
plot3(xspot(:,1),xspot(:,2),xspot(:,3),'k.') % View results
xlabel('x'), ylabel('y'), zlabel('z'), title('Chopping the ends off')
title('Chopping the ends off')
axis([-2,2,-2,2,0.001,0.999])

% 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

pause, close all