function [x xa xb] = get_ellipse(P, N)

% [x xa xb] = get_ellipse(P, N)
%
% DESC:
% computes the points that satisfy the quadratic equation x'*P*x=1, as long
% as P is symmetric positive definite
%
% AUTHOR
% Marco Zuliani - zuliani@ece.ucsb.edu
%
% VERSION
% 1.3
%
% INPUT:
% P             = symmetric positive definite matrix
% N             = number of points
%
% OUTPUT:
% x             = points that satisfy x'*P*x=1
% xa, xb        = points that identify the semiaxis

temp = P - P';
if ~all(temp < eps)
    error('P is not symmetric');
end;

% This can be numerically unstable !!!!!
% [V Lambda] = eig(P);
[V Lambda] = svd(P);
lambda = diag(Lambda);

if ~all(lambda > 0)
    error('P is not positive definte');
end;

theta = linspace(0, 2*pi, N);
y = [cos(theta); sin(theta)];
Q = V*diag(1./sqrt(lambda));

x = Q*y;
xa = Q(:, 1);
xb = Q(:, 2);

return