Path: news.mathworks.com!not-for-mail
From: "Bruno Luong" <b.luong@fogale.findmycountry>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Vectors on ellipsoid surface
Date: Fri, 21 May 2010 17:35:21 +0000 (UTC)
Organization: FOGALE nanotech
Lines: 47
Message-ID: <ht6g8p$5g7$1@fred.mathworks.com>
References: <ht5ldf$7tf$1@fred.mathworks.com> <ht6fr5$73l$1@fred.mathworks.com>
Reply-To: "Bruno Luong" <b.luong@fogale.findmycountry>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1274463321 5639 172.30.248.35 (21 May 2010 17:35:21 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Fri, 21 May 2010 17:35:21 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 390839
Xref: news.mathworks.com comp.soft-sys.matlab:638148

Sorry make, it:

% Data, generate ellipse centered at origin, x'*Q*x = c
L=randn(2);
Q = L'*L;
c = rand;

% Given Point (to be projected on ellipse)
P = randn(2,1);

% Engine
[U S V] = svd(Q);
s = diag(S);
A = U*diag(sqrt(c./s));
R = [0 1;
     -1 0];
B = R*Q*A;
alpha = -P.'*R*Q*A;
rho = sqrt(s(1)/s(2));
beta = c*(rho-1/rho);

Poly = beta*[0 -2 0 2 0] + ...
       alpha(1)*[-1 0 0 0 1] + ...
       alpha(2)*[0 2 0 2 0];
tprj = roots(Poly);
tprj = tprj(imag(tprj)==0);
tprj = reshape(tprj, 1, []);

% Candidate of the projections
xprj = A*[(1-tprj.^2); 2*tprj];
xprj = bsxfun(@rdivide, xprj, (1+tprj.^2));

% Graphical check
% Generate point on ellipse
theta = linspace(0,2*pi);
x = A*[cos(theta); sin(theta)];
clf;
plot(x(1,:),x(2,:));
axis equal;
hold on;
% Plot P projected to the ellipse
plot(P(1),P(2),'or');
for k=1:size(xprj,2)
    plot([xprj(1,k) P(1)],[xprj(2,k) P(2)],'-r');
end

% Bruno