function [AMAJ,BMIN,thetam,wrr]=var_el_params(xt,yt)   
%COMP_VAR_EL_PARAMS compute variance ellipse parameters
% COMP_VAR_EL_PARAMS computes variance ellipse parameters 
% for observation pairs.  All NaN's are removed from the
% input data before computation.
% 
% Input: x,y - 2 vectors of observations, considered pairs
% Output: a  - variance ellipse major axis length
%         b  - variance ellipse minor axis length
%         c  - variance ellipse rotation in radians, 
%              COUNTERclockwise 
%         wr - rotated data, in complex form (xr+i*yr)
% 
% Call as: 
%         >>  [a,b,c,wr]=comp_var_el_params(x,y); 
%
% This routine computes the variance ellipse parameters
% directly; i.e., by computing the variance of the
% data and rotating the original data into the 
% new reference frame, etc.  It would be shorter to
% use MATLAB's cov routine, but this shows explicitly
% what is happening, for educational purposes.  The 
% notation follows exactly that of Preisendorfer, Chapter 2,
% PCA in Met. and Ocean. 

% Find not NaN's in data
temp=[xt(:) yt(:)];
idx=find(~isnan(sum(temp.')));

xp=xt(idx);
yp=yt(idx);

% Demean
n=length(xp);
xbar=mean(xp);
ybar=mean(yp);
x=xp-xbar;y=yp-ybar;

% Find Variances
Sxx=sum(x.*x)/(n-1);	 % Eqn 2.7
Syy=sum(y.*y)/(n-1);  
 
% Find Covariance
Sxy=sum(x.*y)/(n-1);   

% Find Principal angle;  this is the angle along 
% which the variance is maximum
thetam=.5*atan2(2*Sxy,(Sxx-Syy));

AMAJ=sqrt(.5*(Sxx+Syy + sqrt((Sxx-Syy)^2 + 4*Sxy^2)));
BMIN=sqrt(.5*(Sxx+Syy - sqrt((Sxx-Syy)^2 + 4*Sxy^2)));

% The rotated data (zeta,eta), centered at (xbar,ybar), Eqn 2.3
wr=((x+xbar)+i*(y+ybar))*exp(i*thetam);

wrr=(NaN+i*NaN)*ones(size(xt));
wrr(idx)=wr;

%
%        Brian O. Blanton
%        Department of Marine Sciences
%        Ocean Processes Numerical Modeling Laboratory
%        12-7 Venable Hall
%        CB# 3300
%        University of North Carolina
%        Chapel Hill, NC
%                 27599-3300
%
%        919-962-4466
%        blanton@marine.unc.edu
%
%        Spring  1999
%