%  function [deg,xs,ys]=prinaxis(u,v,iplot);
%
%  deg   is the angle rotated 
%   xs   is twice the std deviation of u
%   ys   is twice the std deviation of v
%  u&v   is the original vector time series
% iplot  is zero (0) for no plot and one (1) for plot
%
%
%  original code developed by Beardsley & Limeburner
%  with modifications/documentation by manning
%
function [deg,xs,ys]=prinaxis(u,v,iplot);
%  first, demean
%
xp=u-mean(u);
yp=v-mean(v);


%Find correlation matrix
xxms=mean(xp.*xp);
yyms=mean(yp.*yp);
xyms=mean(xp.*yp);


%Uncorrelated if uvms=0 so rotate axis to make uvms=0
t=.5*atan((2*xyms)/(xxms-yyms));
deg=t*180/pi;


%rotate data
for i=1:length(xp);
  [xr(i),yr(i)]=rot(-deg,xp(i),yp(i));
end;

%define axis length components
xs=std(xr);
ys=std(yr);

if iplot==1;
%draw ellipse
ellipse(0,0,xs,ys,-deg);
axis('equal');
end;





