function theta=princax(u)
%
%     u is vector velocity, E then N, theta is angle clockwise from N
%	  of principal axis (+-180 )
%
cv=cov(u)
theta=0.5*atan(2.*cv(2,1)/(cv(2,2)-cv(1,1)) );
theta=theta*180./pi;
if cv(2,1)>0. and theta<0.  
   theta=theta+90.;
else if cv(2,1)<0. and theta >0. 
   theta=theta-90.;
end
theta
