function [Q,T,V,D,l,qf] = ellipsoid(varargin)
%
%  Draws the ellipsoid 
%      x' * Q * x = 1
%    corresponding to the positive definite symmetric 3 by 3 matrix Q 
%
%  ellipsoid  with no arguments sets Q = eye(3) - a sphere
%  ellipsoid(A)    sets Q = A if A is symmetric; otherwise Q = A * A'
%  ellipsoid(A,n)  same as ellipsoid(A) with n the number of mesh points
%                    default is n = 25
%  ellipsoid(a,b,c)  sets Q = diag([a b c])
%  ellipsoid(a,b,c,n)  sets Q = diag([a b c]) and n = # mesh points
%
%  Output arguments:  [Q,T,V,D,l,qf] = ellipsoid(...)
%
%  Q - symmetric positive definite matrix used
%  T - Cholesky decomposition of Q = T * T'
%  V - columns are eigenvectors or principal axes of ellipse
%  D - diagonal matrix with eigenvalues, so Q = V' * D * V
%  l = sqrt(diag(D)) - lengths of axes of ellipse
%  qf - string giving quadratic form x'*Q*x
%
%  See also ELLIPSE, TR3
%

switch length(varargin)
  case 0
     T = eye(3); n = 25;
  case 1 
     T = varargin{1}; n = 25;
  case 2 
     T = varargin{1}; n = varargin{2};
  case 3 
     T = diag([varargin{1},varargin{2},varargin{3}]); n = 25;
  case 4 
     T = diag([varargin{1},varargin{2},varargin{3}]); n = varargin{4};
end

if T == T' 
   Q = T;
   T = chol(Q)';
   else
   Q = T * T';
end

if Q(1,1) == 1, q1s = ''; else q1s = [num2str(Q(1,1)),' ']; end
if Q(2,2) == 1, q2s = ''; else q2s = [num2str(Q(2,2)),' ']; end
if Q(3,3) == 1, q3s = ''; else q3s = [num2str(Q(3,3)),' ']; end

q12 = abs(2*Q(1,2)); q12x = num2str(q12);
s12 = sign(Q(1,2));
if s12 > 0, q12s = [' + ',q12x,' x y']; 
else if s12 < 0, q12s = [' - ',q12x,' x y'];
else q12s = '';
end 
end 
q13 = abs(2*Q(1,3)); q13x = num2str(q13);
s13 = sign(Q(1,3));
if s13 > 0, q13s = [' + ',q13x,' x z']; 
else if s13 < 0, q13s = [' - ',q13x,' x z'];
else q13s = '';
end 
end 
q23 = abs(2*Q(2,3)); q23x = num2str(q23);
s23 = sign(Q(2,3));
if s23 > 0, q23s = [' + ',q23x,' y z']; 
else if s23 < 0, q23s = [' - ',q23x,' y z'];
else q23s = '';
end 
end 


qf = [q1s,'x^2',q12s,' + ',q2s,'y^2',q13s,q23s,' + ',q3s,'z^2'];

disp(' ')
disp('Quadratic form')
disp(' ')
disp(qf)

disp(' ')

disp('Q = ')
disp(Q)

df = linspace(0,pi,n); 
dt = [0:2*pi/n:2*pi]; 

x0 = sin(df)' * cos(dt); 
y0 = sin(df)' * sin(dt);
z0 = cos(df)' * ones(size(dt));

[x,y,z] = tr3(T,x0,y0,z0);

surf(x,y,z);

[V, D] = eig(Q);

disp('Eigenvectors - principal axes')
disp(' ')
disp(V)
disp('Eigenvalues')
ev = diag(D);
disp(' ')
disp(ev)
disp('Lengths of principal axes')
l = sqrt(ev);
disp(' ')
disp(l)

L = sqrt(D);
W = V * L;

W1 = W(1,:); W2 = W(2,:); W3 = W(3,:); 
N = diag([1 + 1/norm(W(:,1)),1 + 1/norm(W(:,2)),1 + 1/norm(W(:,3))]);
WN = W * N;
WN1 = WN(1,:); WN2 = WN(2,:);  WN3 = WN(3,:);

xa = [W1;WN1];
ya = [W2;WN2];
za = [W3;WN3];
line(xa,ya,za);

xa = - xa;
ya = - ya;
za = - za;
line(xa,ya,za);

axis equal;
shading interp;
colormap(bone);

title(['Ellipsoid:   ',qf,' = 1']);
