function [hit,s] = quadric(A,B,C,inside)
% Calculates the first positive distance to an arbitrary quadric surface
% Inputs:
%  double
%     A:  Quadratic coefficient A (communicated to quadric)
%     B:  Quadratic coefficient B (communicated to quadric)
%     C:  Quadratic coefficient C (communicated to quadric)
%  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)

% Internal variables
%  double
%     Q: quadratic coefficient

Q = B.^2 - A.*C;
if (Q < 0)
    hit = false;
    s = nan;
    return
end

if (inside)           % inside the surface
   if (B <= 0)        % headed away from surface
      if (A > 0)      % but, surface curving up
         hit = true;  % always a hit in this case
         s = (sqrt(Q) - B)./A;
         return
      else            % surface curving away and headed in
         hit = false; % never a hit in this case
         s = inf;
         return
      end
   else               % headed toward the surface
      hit = true;     % always a hit in this case
      s = max(0, -C./(sqrt(Q) + B));
      return
   end
end

% Must be outside the surface
if (B >= 0)          % headed away
    if (A >= 0)      % surface curves away
        hit = false; % never a hit in this case
        s = inf;
        return
    else             % surface curves up
       hit = true;   % always a hit in this case
       s = -(sqrt(Q) + B)./A;
       return
    end
else                 % headed toward the surface
    hit = true;      % always a hit in this case
    s = max(0, C./(sqrt(Q) - B));
    return
end
return