function [hit,s] = ellipsoid(xc,a,x0,u,inside)

% Calculates the distance to an ellipsoid centered at (xs,ys,zs)
% with the mathematical form
% (a(x-x0))^2 + (b(y-y0))^2 + (c(z-z0))^2 = 1

% Inputs
%%%%%%%%
%  double
%     xc:  coordinates of the center of the ellipsoid (row vector)
%     a:   axis distortion factors (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)
%     
% Internal variables
%  double
%     xr: relative coordinates of the particle wrt the ellipsoid center (row vector)
%     A:  Quadratic coefficient A (communicated to quadric) (scalar)
%     B:  Quadratic coefficient B (communicated to quadric) (scalar)
%     C:  Quadratic coefficient C (communicated to quadric) (scalar)



xr = x0 - xc; % Compute relative coordinates
xr = xr.*a;   % Scale the relative coordinates
u  =  u.*a;   % Scale the direction vector

% Compute the quadratic coefficients
A = u*u'; % i.e. [a(1)*u]^2 + [a(2)*v]^2 + [a(3)*w]^2 = 1
B = u*xr';
C = xr*xr' - 1;

% Call the quadratic solver 
[hit,s] = quadric(A,B,C,inside);

return