function [C a] = prin_comp(A)
%               prin_com(A)
% Calculate the principle components of an array of two-dimensional
% data.
% Input: A - 2xn or nx2 dimensional array
% Output: C - Covariance matrix of the data
%         a - Mean of the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[m,n] = size(A);

if n==2
  A = A';
  n=m;
  m=2;
elseif (n~=2)&(m~=2)
    error('A is not a proper array');
end

clf
plot(A(1,:),A(2,:),'.');


X = mean(A(1,:));
Y = mean(A(2,:));

s = ones(1,n);
A(1,:) = A(1,:)-X*s;
A(2,:) = A(2,:)-Y*s;


A = A/sqrt(n);
hold on;
[U,S,V] = svd(A);
Sr = S(1:2,1:2);

C = U*Sr;
t = linspace(0,2*pi,200);
x = cos(t);y = sin(t);
V = [x;y];

% Compute ellipse, twice standard deviation
D = 2*C*V;


D(1,:) = D(1,:)+X;
D(2,:) = D(2,:)+Y;


plot(D(1,:),D(2,:));
hold on
% Plot principle axes
v1 = abs(U(:,1));
v2 = U(:,2);
plot([X X+6*v1(1)],[Y Y+6*v1(2)])
plot([X X+3*v2(1)],[Y Y+3*v2(2)])
axis equal; axis off;

% Covariance matrix
C = A*A';

% Mean
a = [X; Y];