function [hit,s] = plane(xp,np,x0,u,inside)

% Calculates the distance of a particle to a planar surface normal to the
% an arbitrary axis

% Inputs
%%%%%%%%
%  double
%     xp:  coordinate that the plane passes through (row vector)
%     np:  normal direction of the plane (row vector)
%     x0:  coordinates of the particle (row vector)
%     u :  direction cosines of the particle (row vector)
%  logical
%     inside: == true   particle thinks it is inside
%           : == false  particle thinks it is outside
%      

% Outputs
%%%%%%%%%
%  logical
%     hit: == true   particle would hit  the surface
%        : == false  particle would miss the surface
%
%  double
%     s: distance to the surface (if hit)

npdotu = np*u'; % ' is the transpose operator, and this is the dot product of
                % the direction vector with the normal to the plane
if ((inside & npdotu > 0) | (~inside & npdotu < 0)) % headed towards the surface
   hit = true;
   s = max(0,np*(xp-x0)'/npdotu); % Note: equation was rearranged for efficiency
                                  % The max function corrects for 
                                  % delta ambiguities, signaling a region change
   return
else
   % If processing enters here, it is moving away from the plane in question
   hit = false;
   s = inf; % inf is a Matlab keyword for floating positive infinity
   return
end