function [Q,T,V,D,l,qf] = ellipse(A)
%
%  [Q,T,V,D,l,qf] = ellipse(T)
%
%   Draws the ellipse 
%      x' * Q * x = 1
%    corresponding to the positive definite symmetric 2 by 2 matrix 
%    where Q = A if A is symmetric; otherwise Q = A * A'
%
%  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 ELLIPSOID
%

if size(A) ~= [2 2], error('argument must be 2 by 2 matrix'); end

if A == A' 
   Q = A;
   T = chol(Q)';
   else
   Q = A * A';
   T = A;
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

q12 = abs(2*Q(1,2)); q12s = num2str(q12);
s12 = sign(Q(1,2));
if s12 > 0, sg = [' + ',q12s,' x y']; 
else if s12 < 0, sg = [' - ',q12s,' x y'];
else sg = '';
end 
end 
qf = [q1s,'x^2',sg,' + ',q2s,'y^2'];

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

disp(' ')

disp('Q = ')
disp(Q)

t = linspace(0,2*pi,100);

plot(T(1,1) .* cos(t) + T(1,2) .* sin(t), T(2,1) .* cos(t) + T(2,2) .* sin(t),'m')

[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,:);
N = diag([1 + 1/norm(W(:,1)),1 + 1/norm(W(:,2))]);
WN = W * N;
WN1 = WN(1,:); WN2 = WN(2,:);

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

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

axis equal;

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

