function ws = wave_stats(vx,vy,zp,fs)
% WAVE_STATS -  Spectral processing of UVP data for wave statistics
%  ws = wave_stats(vx,vy,zp,fs)
%
% Input:
%   vx, vy - near-bottom velocity (m/s)
%   zp     - pressure converted to wave height (m) 
%   fs      - sampling frequency (Hz)
%
% Returned:
%   ws     - array of wave statistics containing:
%     peraray
%     twoA
%     Urms
%     alpham
%     alphap
%     alphas
%     thetas
%
% The spectrum resolution for the wave period found by spectral peak method 
% is 512 samples.

% Written by Doug Wilson, OHSU
% Modified by Chris Sherwood, USGS
% Last revised December 3, 2001

nfft = 512;
window = hanning(nfft);
noverlap = nfft/2;
lfc = 13;          % low-frequency cutoff when looking for peak
dfifty  = 0.00017; % Sand size used for Shields parameter 
                   % (nondimensionalized bottom stress)
varp = var( detrend(zp));
[ppres, freqs] = pwelch(detrend(zp),window,noverlap,nfft,fs);
varP = sum( ppres .* (freqs(3)-freqs(2)) );
if( (varp./varP) >1.1 | (varp./varP) < .9 ),
  fprintf('WARNING - Pressure variances: varp = %g, varP = %g.\n',varp,varP)
end
[Y,II] = max(ppres(lfc:length(ppres),1));
peraray = 1 ./ freqs(II+(lfc-1));

M = cov([zp(:) vx(:) vy(:)]);
M00 = M(1,1);
M20 = M(2,2);
vxstd = sqrt(M20);
vxmean = mean(vx);
M02 = M(3,3);
vystd = sqrt(M02);
vymean = mean(vy);
M01 = M(1,2);
M10 = M(1,3);
M11 = M(2,3);

% Wave directionality factors,  GAMLH is Longuet-Higgin's wave directionality
Gam2 = (M20 +M02); 
Gam3 = sqrt((M20-M02)^2 + 4*M11^2);
GamLH = sqrt((Gam2 - Gam3)/(Gam2 + Gam3));

alpham = atan(M10/M01).*(360/pi);            %  Mean wave direction
alphap = atan(2*M11/(M20-M02)).*(180/pi);    %  Principal wave direction
Urms = sqrt(2)*sqrt(vxstd.*vxstd + vystd.*vystd); % RMS
Usig = sqrt(2)*Urms;                 %  Significant near-bottom wave velocity
afit = polyfit(vx,vy,1);
alphas = atan(afit(1))*(180/pi);     %  Orientation of wave orbital ellipse.

% Now calculate the Jonsson friction factor from the orbital amplitude.
twoA = Usig*peraray/(2.*pi); 
fw = exp(5.213*(2.5*dfifty*2./twoA).^0.194 - 5.977);
thetas = fw*Usig.^2/(32.34*dfifty);

% copy local variables to structure
ws.perary = peraray;
ws.twoA = twoA;
ws.Urms = Urms;
ws.alpham = alpham;
ws.alphap = alphap;
ws.alphas = alphas;
ws.thetas = thetas;
ws.GamLH = GamLH;