% track2quadric.n Class Example
% Tracking to all the 9 quadrics, plus parallel and crossing
% planes, coincident planes, and even a single plan 
% from a random starting position within a 2x2x2x box.

clear all
close all
format compact

Ns = input('Input the surface option [1...13] (out of range value will give help): ');
if (Ns < 1 || Ns > 13)
   disp('1 is x^2 + y^2 + z^ 2 - 1 = 0, a sphere')
   disp('2 is x^2 + y^2 - z^ 2 = 0, a cone')
   disp('3 is x^2 + y^2 - 1 = 0, a cylinder')
   disp('4 is x^2 + y^2 - z^ 2 - 1 = 0, a hyperboloid of one sheet')
   disp('5 is x^2 + y^2 - z^ 2 + 1 = 0, a hyperboloid of two sheets')
   disp('6 is x^2 + y^2 + z = 0, an elliptic paraboloid')
   disp('7 is x^2 - y^2 + z = 0, a hyperbolic paraboloid')
   disp('8 is x^2 - y^2 + 1 = 0, a hyperbolic cylinder')
   disp('9 is x^2 + z = 0, a parabolic cylinder')
   disp('10 is z*(z-1) = 0, two parallel planes')
   disp('11 is xy = 0, two crossed planes')
   disp('12 is z^2 = 0, two coincident planes at z = 0')
   disp('13 is z = 0, a single plane at z = 0')
   return
end

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

Hits = 0; % Count the number of hits

xin  = []; % Tally for starting points that are inside
xout = []; % Tally for starting points that are outside

for n = 1:Nh

   % pick a starting position
   % The following code picks the initial position randomly with a 2x2x2 box
   % containing the sphere
   x = [1 - 2*rand, 1 - 2*rand, 1 - 2*rand];

   % 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 (Ns == 1)
      if (n == 1) disp('Surface: x^2 + y^2 + z^ 2 - 1 = 0, a sphere'); end
      A = 1;
      B = x*u';
      C = x*x' - 1;
   elseif (Ns == 2)
      if (n == 1) disp('Surface: x^2 + y^2 - z^ 2 = 0, a cone'); end
      A = u(1)*u(1) + u(2)*u(2) - u(3)*u(3);
      B = x(1)*u(1) + x(2)*u(2) - x(3)*u(3);
      C = x(1)*x(1) + x(2)*x(2) - x(3)*x(3);
   elseif (Ns == 3)
      if (n == 1) disp('Surface: x^2 + y^2 - 1 = 0, a cylinder'); end
      A = u(1)*u(1) + u(2)*u(2);
      B = x(1)*u(1) + x(2)*u(2);
      C = x(1)*x(1) + x(2)*x(2) - 1;
   elseif (Ns == 4)
      if (n == 1) disp('Surface: x^2 + y^2 - z^ 2 - 1 = 0, a hyperboloid of one sheet'); end
      A = u(1)*u(1) + u(2)*u(2) - u(3)*u(3);
      B = x(1)*u(1) + x(2)*u(2) - x(3)*u(3);
      C = x(1)*x(1) + x(2)*x(2) - x(3)*x(3) - 1;
   elseif (Ns == 5)
      if (n == 1) disp('Surface: x^2 + y^2 - z^ 2 + 1 = 0, a hyperboloid of two sheets'); end
      A = u(1)*u(1) + u(2)*u(2) - u(3)*u(3);
      B = x(1)*u(1) + x(2)*u(2) - x(3)*u(3);
      C = x(1)*x(1) + x(2)*x(2) - x(3)*x(3) + 1;
   elseif (Ns == 6)
      if (n == 1) disp('Surface: x^2 + y^2 + z = 0, an elliptic paraboloid'); end
      A = u(1)*u(1) + u(2)*u(2);
      B = x(1)*u(1) + x(2)*u(2) + 0.5*u(3);
      C = x(1)*x(1) + x(2)*x(2) + x(3);
   elseif (Ns == 7)
      if (n == 1) disp('Surface: x^2 - y^2 + z = 0, a hyperbolic paraboloid'); end
      A = u(1)*u(1) - u(2)*u(2);
      B = x(1)*u(1) - x(2)*u(2) + 0.5*u(3);
      C = x(1)*x(1) - x(2)*x(2) + x(3);
   elseif (Ns == 8)
      if (n == 1) disp('Surface: x^2 - y^2 + 1 = 0, a hyperbolic cylinder'); end
      A = u(1)*u(1) - u(2)*u(2);
      B = x(1)*u(1) - x(2)*u(2);
      C = x(1)*x(1) - x(2)*x(2) + 1;
   elseif (Ns == 9)
      if (n == 1) disp('Surface: x^2 + z = 0, a parabolic cylinder'); end
      A = u(1)*u(1);
      B = x(1)*u(1) + 0.5*u(3);
      C = x(1)*x(1) + x(3);
   elseif (Ns == 10)
      if (n == 1) disp('Surface: z*(z-1) = 0, two parallel planes'); end
      A = u(3)*u(3);
      B = x(3)*u(3) - 0.5*u(3);
      C = x(3)*(x(3) - 1);
   elseif (Ns == 11)
      if (n == 1) disp('Surface: xy = 0, two crossed planes'); end
      A = u(1)*u(2);
      B = 0.5*(x(1)*u(2) + x(2)*u(1));
      C = x(1)*x(2);
   elseif (Ns == 12)
      if (n == 1) disp('Surface: z^2 = 0, two coincident planes'); end
      A = u(3)*u(3);
      B = x(3)*u(3);
      C = x(3)*x(3);
   elseif (Ns == 13)
      if (n == 1) disp('Surface: z = 0, single plane'); end
      A = 0;
      B = 0.5*u(3);
      C = x(3);
   end
   % Determine inside or outside
      if  (C < 0)
      inside = true;
   elseif (C > 0)
      inside = false;
   else
      if (B > 0) % On surface headed outside 
         inside = false; % Place it outside
      elseif (B < 0) % On surface headed inside
         inside = true; % Place it inside
      else % On surface and tangent to it
         if (A > 0)
            inside = false; % Place it outside
         elseif (A < 0)
            inside = true; % Place it inside
         else
            % It is going along a ruling. Do it the MC way
            if (rand < 0.5)
               inside = false; % Place it outside
            else
               inside = true; % Place it inside
            end
         end
      end
   end
      
   [hit,s] = quadric(A,B,C,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);
      if      (inside) xin(end + 1,:)  = x;
      elseif (~inside) xout(end + 1,:) = x;
      end
   end
    
end % End of loop over histories



figure(1)
plot3(xhit, yhit, zhit,'k.') % Plot the surface intercepts
xlabel('x')
ylabel('y')
zlabel('z')
title('Where particles tracked to')
if (Ns == 1)
   axis equal
elseif(Ns == 5)
   axis([-3,3,-3,3,-3,3])
elseif(Ns == 10)
   axis([-2,2,-2,2,-0.5,1.5])
else
   axis([-2,2,-2,2,-2,2])
end

figure(2)
plot3(xout(:,1),xout(:,2),xout(:,3),'r.') % View results
axis equal
xlabel('x')
ylabel('y')
zlabel('z')
title('Where particles outside started')

figure(3)
plot3(xin(:,1),xin(:,2),xin(:,3),'b.') % View results
axis equal
xlabel('x')
ylabel('y')
zlabel('z')
title('Where particles inside started')