princax

 PRINCAX Principal axis, rotation angle, principal ellipse

   [theta,maj,min,wr]=princax(w) 

   Input:  w   = complex vector time series (u+i*v)

   Output: theta = angle of maximum variance, math notation (east == 0, north=90)
           maj   = major axis of principal ellipse
           min   = minor axis of principal ellipse
           wr    = rotated time series, where real(wr) is aligned with 
                   the major axis.

 For derivation, see Emery and Thompson, "Data Analysis Methods
   in Oceanography", 1998, Pergamon, pages 325-327.  ISBN 0 08 0314341
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Version 1.0 (12/4/1996) Rich Signell (rsignell@usgs.gov)
 Version 1.1 (4/21/1999) Rich Signell (rsignell@usgs.gov) 
     fixed bug that sometimes caused the imaginary part 
     of the rotated time series to be aligned with major axis. 
     Also simplified the code.
 Version 1.2 (3/1/2000) Rich Signell (rsignell@usgs.gov) 
     Simplified maj and min axis computations and added reference
     to Emery and Thompson book
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [theta,maj,min,wr]=princax(w)
0002 % PRINCAX Principal axis, rotation angle, principal ellipse
0003 %
0004 %   [theta,maj,min,wr]=princax(w)
0005 %
0006 %   Input:  w   = complex vector time series (u+i*v)
0007 %
0008 %   Output: theta = angle of maximum variance, math notation (east == 0, north=90)
0009 %           maj   = major axis of principal ellipse
0010 %           min   = minor axis of principal ellipse
0011 %           wr    = rotated time series, where real(wr) is aligned with
0012 %                   the major axis.
0013 %
0014 % For derivation, see Emery and Thompson, "Data Analysis Methods
0015 %   in Oceanography", 1998, Pergamon, pages 325-327.  ISBN 0 08 0314341
0016 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0017 % Version 1.0 (12/4/1996) Rich Signell (rsignell@usgs.gov)
0018 % Version 1.1 (4/21/1999) Rich Signell (rsignell@usgs.gov)
0019 %     fixed bug that sometimes caused the imaginary part
0020 %     of the rotated time series to be aligned with major axis.
0021 %     Also simplified the code.
0022 % Version 1.2 (3/1/2000) Rich Signell (rsignell@usgs.gov)
0023 %     Simplified maj and min axis computations and added reference
0024 %     to Emery and Thompson book
0025 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0026 
0027 % use only the good (finite) points
0028 ind=find(finite(w));
0029 wr=w;
0030 w=w(ind);
0031 
0032 % find covariance matrix
0033 cv=cov([real(w(:)) imag(w(:))]); 
0034 
0035 % find direction of maximum variance
0036 theta=0.5*atan2(2.*cv(2,1),(cv(1,1)-cv(2,2)) );
0037 
0038 % find major and minor axis amplitudes
0039 
0040 term1=(cv(1,1)+cv(2,2));
0041 term2=sqrt((cv(1,1)-cv(2,2)).^2 + 4.*cv(2,1).^2);
0042 maj=sqrt(.5*(term1+term2));
0043 min=sqrt(.5*(term1-term2));
0044 
0045 % rotate into principal ellipse orientation
0046 
0047 wr(ind)=w.*exp(-i*theta);
0048 theta=theta*180./pi;