function s = pca2d(u,v)
%s = pca2d(u,v) Principal Component Analysis of 2D datasets
%   
%Input:
%   u,v = 2D vector data
%
%Output:
%   s = Structure with the following fields:
%       center = Ellipse center (centroid)
%       area = Ellipse area
%       axes = Ellipse axes
%       eigen_values = Eigen values
%       eigen_vectors = Eigen vectors
%       total_variance = Total variance
%       percentage_variance = Percentage variance for each axes
%       angles = Angles for each axes in degrees
%       eccentricity = Ellipse eccentricity
%       slope = Slope of the linear fit which goes through the origin
%
% Disclaimer:
%  This program is provided as is without any guaranty or warranty.
%  The user assumes all risk involved in the application of this software.
%
% Version: 1.0
% Date: Novermber 1, 2012
% Author: Alex Sanchez, 2012
if nargin==0
    load wind
    u = u(:,20,:);
    v = -v(:,20,:);
end
u = u(:); v = v(:);
ind = (isfinite(u) & isfinite(v));
u = u(ind); v = v(ind);
s.center = [mean(u),mean(v)];
[eigvecmat,eigvalmat] = eig(cov(u,v));   
eigval = svd(eigvalmat)';
s.area = pi*prod(eigval);
s.axes = sqrt(eigval);
s.eigen_values = eigval;
s.eigen_vectors = eigvecmat;
s.total_variance = sum(eigval);
s.percentage_variance = 100*eigval/s.total_variance;
s.angles = atan2(eigvecmat(1,:),eigvecmat(2,:)); %note swich of axis for conversion
s.angles = mod(s.angles*180/pi-90,360); %note rotation for conversion
s.eccentricity = 1-(eigval(2)/eigval(1));
s.slope = (u-s.center(1))\(v-s.center(2));
if nargin>0
    return
end
s.handles = ones(4,1); %Init handles
clf
hold on
%scatter data
s.handles(1) = plot(u,v,'r.','markersize',5);
%ellipse
N=100;
t=linspace(0,2*pi,N);
ellip=(eigvecmat*sqrt(eigvalmat)*[cos(t);sin(t)]+repmat(s.center',1,N))';
s.handles(2)=line(ellip(:,1),ellip(:,2),'linewidth',2','color','b');
%linear fit
x = ([min(u),max(u)]-s.center(1))*0.8;
s.handles(3) = line(x+s.center(1),s.center(2)+x*s.slope,...
    'linewidth',2,'color','k','linestyle','--');
%axes
line1=eigvecmat*sqrt(eigvalmat)*[-1,1;0,0]+repmat(s.center',1,2);
line2=eigvecmat*sqrt(eigvalmat)*[0,0;-1,1]+repmat(s.center',1,2);
dat=[line1'; NaN,NaN; line2'];
s.handles(4) = line(dat(:,1),dat(:,2),'linewidth',2,'color','b');
%text
table= sprintf('Center = %4.2f, %4.2f',s.center(1),s.center(2));
table=char(table,sprintf('Area = %4.2f',s.area));
table=char(table,sprintf('Axes = %4.2f, %4.2f',s.axes(1),s.axes(2)));
table=char(table,sprintf('Angles = %4.2f, %4.2f',s.angles(1),s.angles(2)));
table=char(table,sprintf('Slope = %4.2f',s.slope));
axis equal
box on
limits=axis;
xtxt=limits(1)+(limits(2)-limits(1))*0.2;
ytxt=limits(3)+(limits(4)-limits(3))*0.2;
text(xtxt,ytxt,table,'FontSize',8);
xlabel('x-Component of Wind speed, m/s')
ylabel('y-Component of Wind speed, m/s')
legend(s.handles(1:3),'Scatterset','Ellipse','Linear Fit',0)
disp(s)
return
