% Class Example - Tracking to a unit sphere, from a random start within 
%                 minimum bounding box. Includes lots of diagnostic coding.

clear all
close all

format compact

Nh = input('Number of histories: ');

Nd = input('Number of histories with individual output: ');

Inside  = 0; % Tally these
Outside = 0;
Exact   = 0; % Numerically on the surface
Exact2  = 0; % Numerically on the surface and tangent to it

xhit = [];
yhit = [];
zhit = [];
rhit = [];
Hits = 0;

teststart = 0;

for n = 1:Nh

   % pick a starting position
   if (teststart == 0)
      % The following code picks the initial position randomly with a 2x2 box
      % containing the sphere
      x = [1 - 2*rand,1 - 2*rand,1 - 2*rand];
   elseif (teststart == 1)
      % The following code puts it randomly on the surface
      costhe = 1 - 2*rand;
      sinthe = sqrt(1 - costhe^2);
      phi = 2*pi*rand;
      x = [sinthe*cos(phi),sinthe*sin(phi),costhe];
   elseif (teststart == 2)
      % The option puts it tangent to the surface
      x = [1,0,0];
   end

   if (n <= Nd)
      disp(sprintf('Starting [x,y,z] = [%g, %g, %g], x^2 = %g',x, x*x'))
   end
      
   % pick a starting direction vector
   costhe = 1 - 2*rand;
   sinthe = sqrt(1 - costhe^2);
   phi = 2*pi*rand;
   u = [sinthe*cos(phi),sinthe*sin(phi),costhe];
   if (teststart == 2)
      % The option puts it tangent to the surface
      u = [0,1,0];
   end
   if (n <= Nd)
      disp(sprintf('Starting [u,v,w] = [%g, %g, %g]',u))
   end

   % Determine inside or outside
   if (x*x' <  1)
      inside = true;
      Inside = Inside + 1;
      if (n <= Nd)
         disp('Starting inside')
      end
   elseif (x*x' >  1)
      inside = false;
      Outside = Outside + 1;
      if (n <= Nd)
         disp('Starting outside')
      end
   else
      Exact = Exact + 1;
      xdotu = x*u';
      disp(sprintf('x dot u = %g',xdotu))
      
      if (xdotu > 0) % On surface headed outside 
         inside = false; % Place it outside
         Outside = Outside + 1;
         if (n <= Nd)
            disp('Starting on the surface headed outside')
         end
      elseif (xdotu < 0) % On surface headed inside
         inside = true; % Place it outside
         Inside = Inside + 1;
         if (n <= Nd)
            disp('Starting on the surface headed inside')
         end
      else % On surface and tangent to it
         Exact2 = Exact2 + 1;
         inside = false; % Place it outside
         Outside = Outside + 1;
         if (n <= Nd)
            disp('Starting on the surface, tangent to it, headed outside')
         end
      end
   end
   if (teststart > 0)
      disp(' ')
      pause
   end
   
   [hit,s] = sphere(1,[0,0,0],x,u,inside);
   if (hit) 
      Hits = Hits + 1;
      xhit(Hits) = x(1) + s*u(1);
      yhit(Hits) = x(2) + s*u(2);
      zhit(Hits) = x(3) + s*u(3);
      rhit(Hits) = xhit(Hits)^2 + yhit(Hits)^2 + zhit(Hits)^2;
   end
   if (n <= Nd)
      if (hit)
         disp(sprintf('Hits the sphere after a step = %g',s))
      else
         disp('Particle misses the sphere')
      end
   end
    
end % End of loop over histories

plot3(xhit, yhit, zhit,'k.')