npoints = 10000;
t = (2*pi)*rand(1,npoints); % angulo
rho = rand(1,npoints); % modulo
xc=3;yc=3;
a=3;b=2;phi=atan(2/3);

% círculo de raio unitário
% x^2 + y^2 = 1^2

% círculo centrado em (xc,yc) de raio unitário
% (x-xc)^2 + (y-yc)^2 = 1

% elipse centrada em (xc,yc) e raios a e b.
% (x-xc)^2/a + (y-yc)^2/b = 1

% representação polar de uma elipse de raios a e b centrada em (xc,yc)
% parametrizada por t
% x(t) = xc + a*cos(t)
% y(t) = yc + b*cos(t)

% matriz de rotação por phi
% rot(phi) =  [cos(phi) -sin(phi); sin(phi) cos(phi)]

% nuvem de pontos determinados por uma elipse de raio a e b centrada em
% (xc,yc)
x = xc + (rho*a) .* cos(t) * cos(phi) - (rho*b) .* sin(t) * sin(phi);
y = yc + (rho*a) .* cos(t) * sin(phi) + (rho*b) .* sin(t) * cos(phi);

%data = [x'-repmat(mean(x),npoints,1) y'-repmat(mean(y),npoints,1)];
data = [x' y'];
covdata = cov(data);

[autovetor,score,autovalor] = princomp(data);
% [autovetor,autovalor,explained] = pcacov(cov(data));

% score_nonorm(:,1) = data * autovetor(:,1);
% score_nonorm(:,2) = data * autovetor(:,2);

for(i=1:2)
    %cov(data) * autovetor(:,i) % A * v(i)
    % lambda * v(i)
    fprintf(1,'v%d * lambda%d = (%7.4f,%7.4f) => |v%d * lambda%d| = %7.4f, sqrt(|v%d * lambda%d|) = %7.4f\n',i,i,autovalor(i) * autovetor(1,i),autovalor(i) * autovetor(2,i) ...
        ,i,i,sqrt( (autovalor(i) * autovetor(1,i))^2 + (autovalor(i) * autovetor(2,i))^2 ) ... 
        ,i,i,sqrt(sqrt( (autovalor(i) * autovetor(1,i))^2 + (autovalor(i) * autovetor(2,i))^2 )) );
end

%

figure;
subplot(2,2,1),hold
subplot(2,2,1),axis([0 6 0 6]);
subplot(2,2,1),plot(x,y,'o');
subplot(2,2,1),plot([xc xc+autovalor(1) * autovetor(1,1)],[yc yc+autovalor(1) * autovetor(2,1)],'r-d');
subplot(2,2,1),plot([xc xc+autovalor(2) * autovetor(1,2)],[yc yc+autovalor(2) * autovetor(2,2)],'r-d');

subplot(2,2,2),hold
subplot(2,2,2),axis([-3 3 -3 3]);
subplot(2,2,2),plot(score(:,1),score(:,2),'rx')
subplot(2,2,2),plot([0 autovalor(1) * sqrt(autovetor(1,1)^2+autovetor(2,1)^2) ],[0 0],'b-d');
subplot(2,2,2),plot([0 0],[0 autovalor(2) * sqrt(autovetor(1,2)^2+autovetor(2,2)^2)],'b-d');

subplot(2,2,3),hold
subplot(2,2,3),axis([-3 3 -3 3]);
subplot(2,2,3),plot(score(:,1),zeros(npoints,1),'rx')

subplot(2,2,4),hold
subplot(2,2,4),axis([-3 3 -3 3]);
subplot(2,2,4),plot(zeros(npoints,1),score(:,2),'rx')