function [Evec, Eval, freqs, UD_power, EW_power, NS_power, First_eigenv_power, deg_of_pol, ThetaH, ThetaV, phiHH, phiVH] = FDPA_Park_Attributes(Z,E,N,Fs,Nwin);
%
% This code prepares eigenvalues and eigenvectors, evaluates power, and
% conducts polarization analysis. The timeseries are differentiated once
% in the code (to turn input velocity to output quantities for acceleration)
%
% To compute values in time windows, you can invoke this function for each
% time segment of data.
%
% This code also computes a polarization attribute not in Park et al (1987):
% the degree of polarization (Samson 1983).
%
% The INPUTS for the code are:
% Z     = Vertical Component velocity timeseries
% E     = East Component velocity timeseries
% N     = North Component velocity timeseries
% Fs    = Sampling Frequency (Hz)
% Nwin  = Number of subwindows
%
% The OUTPUTS of the code are:
% Evec  =
% Eval  =
% freq  =
% UD_power  = 10*log10(P), where P is the power (acceleration squared) in
%                         the up-down component
% EW_power  = 10*log10(P), where P is the power (acceleration squared) in
%                         the east-west component
% NS_power  = 10*log10(P), where P is the power (acceleration squared) in
%                         the nort-south component
% First_eigenv_power    = 10*log10(P), where P is the power (acceleration squared) in
%                         the dominant eigenvector of the spectral
%                         covariance matrix C
% deg_of_pol    = degree of polarization as defined by Samson (1983), which
%                         varies from zero, when all three eigenvalues are
%                         equal, to one, when the data are described by a
%                         single non-zero eigenvalue
% ThetaH    = The azimuth of the dominant eigenvector describing the particle
%                         motion ellipsoid, in degrees (clockwise from north)
% ThetaV    = The angle between the vertical and horizontal components of
%                         the eigenvector describing the particle motion
%                         ellipsoid, in degrees (from vertical)
% phiHH     = The horizontal component phase difference, in degrees
% phiHH     = The horizontal component phase difference, in degrees
%
%
% The polarization analysis closely follows a python implementation at this page:
% https://seiscode.iris.washington.edu/svn/noise-toolkit-microseism-energy-me-bundle/trunk/lib/polarLib.py
% https://ds.iris.edu/ds/products/noise-toolkit-polarization/
%
% Please cite the following when using this code:
%
% Goodling, P. Lekic, V. and K. Prestegaard (2017) Seismic signature of
% turbulence during the 2017 Oroville Dam spillway erosion crisis,
% Earth Surface Dynamics.
%
% Koper K.D, and V.L. Hawley (2010), Frequency dependent polarization
% analysis of ambient seismic noise recorded at a broadband seismometer
% in the Central United States, Earthquake Science, 23, 439-447.%
%
% and
%
% Park, J., Vernon, F. L., and Lindberg, C. R. (1987)
% Frequency-Dependent Polarization Analysis Of High-Frequency Seismograms,
% Journal of Geophysical Research-Solid Earth and Planets, 92, 12664-12674,
% doi:10.1029/JB092iB12p12664.


%% Prepare the Data for Polarization Analysis
%Step 1: Detrend the data within the sample window
sample_Z = detrend(Z);
sample_E = detrend(E);
sample_N = detrend(N);

clear Z E N

% Step 2: Divide each window into Nwin subwindows, each of which overlaps 
% by 50% Each window is a column in the hour_(Z,N,E)_w matrix

Npts = 2*(length(sample_Z)/(Nwin+1));

% Check if the total duration can be divided into Nwin equal size integer
% chunks
if(Npts ~= round(Npts)), error('Try a different Nwin'); end

% If Npts is an integer, then only the first window will be zero-padded and
% must be thrown out. Otherwise, both the first and last do.
YZ = buffer(sample_Z,Npts,Npts/2); YZ(:,1) = [];
YE = buffer(sample_E,Npts,Npts/2); YE(:,1) = [];
YN = buffer(sample_N,Npts,Npts/2); YN(:,1) = [];

% Step 4&5: Taper using a hanning function and calculate Fourier Transforms

win = hann(length(YZ(:,1))); %Create a hanning function the length of each window
%clear sample_Z sample_E sample_N;

%Create Tapers
for i = 1:Nwin
    % Divide by sqrt(Npts) to get the correct power, see fft documentation.
    sample_Z_t(:,i) = YZ(:,i) .* win; sample_Z_f(:,i) = fft(sample_Z_t(:,i))./sqrt(Npts);
    sample_E_t(:,i) = YE(:,i) .* win; sample_E_f(:,i) = fft(sample_E_t(:,i))./sqrt(Npts);
    sample_N_t(:,i) = YN(:,i) .* win; sample_N_f(:,i) = fft(sample_N_t(:,i))./sqrt(Npts);
end

% fft handles points differently if they are odd or even. We need the
% frequencies to reflect this to go to acceleration.

if(mod(Npts,2)==1) % Odd number of points
    freqs = linspace(0,Fs/2,(Npts-1)/2+1);
    freqs = [freqs fliplr(freqs(1:end-1))];
else                % Even number of points
    freqs = linspace(0,Fs/2,Npts/2);
    freqs = [freqs fliplr(freqs)];
end

% Differentiate to acceleration by multiplying by i*angular frequency in 
% the Fourier domain
% freqs = linspace(0,Fs/2,length(sample_Z_f(:,1)));

for n = 1:Nwin;
    sample_Z_a(:,n) = sample_Z_f(:,n).*transpose((1i.*freqs.*2.*pi));
    sample_E_a(:,n) = sample_E_f(:,n).*transpose((1i.*freqs.*2.*pi));
    sample_N_a(:,n) = sample_N_f(:,n).*transpose((1i.*freqs.*2.*pi));
end

% Step 6: Compute a spectral covariance matrix at EACH frequency.
% Compute the elementwise mean.
% Diagonal values represent power spectra of each component, off-diagonal
% values represent cross spectra of each component.

% When we compute the covariance of each frequency at each time step,
% we create a 3 x 3 x [# of windows] x [# of freqs] matrix
clear C Eval Evec;
freqs = linspace(0,Fs/2,length(sample_Z_f(:,1))/2+1);


% First, we take the fourier coefficients and construct one 3x3 for each
% time window

for j = 1:length(sample_Z_f(:,1))/2+1% Go one frequency at a time for all windows
    
    M(:,1) = sample_E_a(j,:); % Component 1 is E
    M(:,2) = sample_N_a(j,:); % Component 2 is N
    M(:,3) = sample_Z_a(j,:); % Component 3 is Z
    
    C(:,:,j) = (M'*M)./length(M(:,1)); %Compute the cross-spetral matrix
end


%Compute eigenvectors/eigenvalues for each frequency
for j = 1:length(sample_Z_f(:,1))/2+1;
    [V,S] = eig(C(:,:,j)); Evec(:,:,j) = V; Eval(:,j) = diag(S);
end

clear V S M

%Compute polarization metric (not used in Park, 1987)
for j = 1:length(sample_Z_f(:,1))/2+1
    deg_of_pol(j) = (3*trace(C(:,:,j)^2) - trace(C(:,:,j))^2)/(2*trace(C(:,:,j))^2); %From Samson (1983), also used by Koper and Hawley (2010)
end

clear C

% Compute the power (from the accelerograms) of each component at each
% frequency, put into decibels, and compute average in each of Nwin windows.

UD_power_windows = 10*log10(sample_Z_a.*conj(sample_Z_a));
UD_power_mean = mean(UD_power_windows,2);
UD_power = transpose(UD_power_mean(1:(length(UD_power_mean(:,1))/2)+1));
clear UD_power_windows UD_power_mean

EW_power_windows = 10*log10(sample_E_a.*conj(sample_E_a));
EW_power_mean = mean(EW_power_windows,2);
EW_power = transpose(EW_power_mean(1:(length(EW_power_mean(:,1))/2)+1));
clear EW_power_windows EW_power_mean

NS_power_windows = 10*log10(sample_N_a.*conj(sample_N_a));
NS_power_mean = mean(NS_power_windows,2);
NS_power = transpose(NS_power_mean(1:(length(NS_power_mean(:,1))/2)+1));
clear NS_power_windows NS_power_mean

First_eigenv_power = 10*log10(Eval(3,:));

%% Perform Polarization analysis

% This code takes the dominant eigenvector XYZ components and computes the
% parameters described in Park 1987. Input the complex x, y, and z
% components of the dominant eigenvector of a 3x3 spectral covariance matrix.
% In the Park paper, Z1 is vertical, E is Z2, and N is Z3.
% Inconsistencies in the Park 1987 paper are handled in a similar way to
% the code by Keith Koper:

% phiHH is the phase difference between the horizontal components of motion.
% If phiHH is ~0 or 180, it is linear motion. If it is ~90, then the
% horizontal motion is elliptical, with motion aligned with the
% instruments. if zx = +/- i*zy, then the motion is circular and ThetaV is
% meaningless.
% Following Park, will be restricted to -180 and +180.

for j = 1: length(freqs);
    z1 = Evec(3,3,j);
    z2 = Evec(1,3,j);
    z3 = Evec(2,3,j);
    
    
    phiHH_temp(j) = (180/pi)*(angle(z3)-angle(z2));
    if phiHH_temp(j) >= 180;
        phiHH(j) = phiHH_temp(j) - 360;
    elseif phiHH_temp(j) <= -180;
        phiHH(j) = phiHH_temp(j) + 360;
    else
        phiHH(j) = phiHH_temp(j);
    end
    clear phiHH_temp
    
    % To find the degree of tilt of the polarization ellipsoid, we must find
    % little theta H. To do so, we must find the value for l that maximizes
    % equation 5 in Park 1987, since little_theta_l will be substituted for
    % 2pift in equation 5
    
    % We want to find the l closest to zero. We will look at values 0 through 4.
    % If there isn't any difference in values, we will select l to be zero.
    
    val = -0.5 *angle(z2*z2 + z3*z3);
    
    l = 0:1:4;
    tmp_angle = val+(l*(pi/2));
    parkeqn_5 = abs(z2).^2.*(cos(tmp_angle+angle(z2)).^2)  ...
        + abs(z3).^2.*(cos(tmp_angle+angle(z3)).^2);
    l_text = l(parkeqn_5==max(parkeqn_5));
    
    % Once we have l, we can identify the horizontal phase angle (little
    % thetaH) in equation 7 in Park 1987
    
    little_theta_h = -0.5*angle(z2^2 + z3^2) + (l_text(1)*pi/2);
    
    clear tmp_angle parkeqn_5 l_text l val
    
    % ThetaH is the horizontal azimuth of the dominant eigenvector. It is
    % big Theta H in equation 8 in Park 1987. Depending on the sign of the
    % real(z*conj(x)) value, it will be a positive or negative angle.
    % Note the Park et al. gives ThetaH as counterclockwise from East, wheras
    % for backazimuth we normally want clockwise from North.
    
    %Eqn 8
    ztmp = exp(1i*little_theta_h);
    
    zval3 = z3 * ztmp;
    zval2 = z2 * ztmp;
    ThetaH(j) = atan2(real(zval3),real(zval2));
    
    if (real(z1*conj(z3)) <  0) && (ThetaH(j) < 0)
        ThetaH(j) = ThetaH(j) + pi;
    end
    if (real(z1*conj(z3))) >= 0 && (ThetaH(j) > 0)
        ThetaH(j) = ThetaH(j) - pi;
    end
    ThetaH(j) = pi/2 - ThetaH(j);   % measure azimuth clockwise from north
    % rather than counterclockwise from east
    if(ThetaH(j) < 0)               % Restrict to 0 to 360 interval.
        ThetaH(j) = ThetaH(j) + 2*pi;
    end
    ThetaH(j) = 180/pi*ThetaH(j);
    
    clear ztmp zval2 zval2
    
    % Using equations 9 and 10 in Park 1987, we need to find the vertical
    % tilt of the ellipsoid, in this case maximizing equation 9 with by
    % selecting the appropriate value of m
    
    % We want to find the integer value of m that is closest to zero and
    % maximizes equation 9. If there isn't any difference, then set m to
    % zero.
    
    val2 = -0.5 * angle(z1^2 + z2^2 + z3^2);
    
    m = 0:1:4;
    tmp_angle2 = val2 + (m*(pi/2));
    parkeqn_9 = abs(z1).^2*cos(tmp_angle2+angle(z1)).^2 ...
        +abs(z2).^2*cos(tmp_angle2+angle(z2)).^2 ...
        +abs(z3).^2*cos(tmp_angle2+angle(z3)).^2;
    
    
    m_text = m(parkeqn_9==max(parkeqn_9));
    
    % Once we have m, we can identify the vertical phase angle (little
    % theta v) in equation 9
    little_theta_v = -0.5*angle(z1^2 + z2^2 + z3^2) + (m_text(1)*pi/2);
    
    clear val2 m tmp_angle2 parkeqn_9 m_text
    
    % With this, we can determine the phase difference between the
    % horizontal and vertical components (phiVH in text after equation 10;
    phiVH(j) = (180/pi)*(little_theta_h - angle(z1));
    
    % Restrict phiVH between -90 and 90
    while phiVH(j) > 90
        phiVH(j) = phiVH(j) - 180;
    end
    
    while phiVH(j) < -90
        phiVH(j) = phiVH(j) + 180;
    end
    
    % The vertcial angle of incidence (tilt angle) is restricted between
    % 0 and 90 degrees. The angle is measured with respect to the vertical.
    zh = sqrt(z2^2 + z3^2);
    
    if imag(zh)<0
        zh = -zh;
    end
    
    ThetaV(j) = atand(abs(real(z1*exp(-1i*little_theta_v))/real(zh*exp(-1i*little_theta_v))));
    ThetaV(j) = 90- ThetaV(j); %Measure from the vertical
    
    clear zh
    
    % plotting
    % subplot(2,2,1);semilogx(freqs,ThetaH);xlabel('f/Hz');ylabel('Degrees from E (counterclockwise)');
    % ylim([0 360]);
    % subplot(2,2,2);semilogx(freqs,ThetaV);xlabel('f/Hz');ylabel('Degrees from Vertical');
    % ylim([0 90]);
    % subplot(2,2,3);semilogx(freqs,phiVH);xlabel('f/Hz');ylabel('Vertical to Horizontal Phase Difference');
    % ylim([-90 90]);
    % subplot(2,2,4);semilogx(freqs,phiHH);xlabel('f/Hz');ylabel('Horizontal Component Phase Difference');
    % ylim([-180 180]);
    
end
end
