% track2sphere.n Class Example
% Tracking to a unit sphere, radius = 1, centered at the origin 
% from a random starting position within a minimum bounding box.

clear all, close all

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

Hits = 0; % Count the number of hits

for n = 1:Nh

   % Pick a random starting position within the box that contains the sphere
   x = [1 - 2*rand,1 - 2*rand,1 - 2*rand];

   % Pick a random starting direction vector
   costhe = 1 - 2*rand;
   sinthe = sqrt(1 - costhe^2);
   phi = 2*pi*rand;
   u = [sinthe*cos(phi),sinthe*sin(phi),costhe];

   % Determine wherther inside or outside
   if     (x*x' <  1) inside = true;  % Starts inside
   elseif (x*x' >  1) inside = false; % Starts outside
   else % x*x' = 1, on the sphere, need another way to determine inside/outside
      xdotu = x*u';      % Dot product of x[] aand u[]
      if (xdotu > 0)     % On surface headed outside 
         inside = false; % Place it outside
      elseif (xdotu < 0) % On surface headed inside
         inside = true;  % Place it inside
      else % On surface and tangent to it
         inside = false; % Place it outside (this works for spheres)
      end
   end
   
   [hit,s] = sphere(1,[0,0,0],x,u,inside); % See sphere.m for explanation of I/O
   
   if (hit) 
      Hits = Hits + 1; % Count the hits
      xhit(Hits) = x(1) + s*u(1); % Track to the surface
      yhit(Hits) = x(2) + s*u(2);
      zhit(Hits) = x(3) + s*u(3);
      r2(Hits) = xhit(Hits)^2 + yhit(Hits)^2 + zhit(Hits)^2; % Diagnostic
   end
    
end % End of loop over histories

disp(sprintf('Minimum radius^2 - 1 = %g, Maxiumum radius^2 - 1 = %g', min(r2) - 1, max(r2) - 1))

plot3(xhit, yhit, zhit,'k.') % Plot the surface intercepts
axis equal