function [Guu, Gvv, hgt] = waveht_uv(udata, vdata, zdep, hgtsens, sampfreq)
% WAVEHT_UV - Determine wave heights from u, v velocity data
% 
% [Guu, Gvv, hgt] = waveht_uv(udata, vdata, zdep, hgtsens, sampfreq)
%
% Method of Thornton and Guza (1982)
%
% Input
%    udata		u velocity data
%    vdata		v velocity data
%    zdep	 	average water depth from pressure data
%    hgtsens	height of sensor off bottom
%    sampfreq	sampling freqency

% Meg Palmsten USF/USGS St Pete, 2/00
% Modified by Chris Sherwood, USGS\
% Last modified Sept 13, 2001

% Determine velocity spectra for u and v
% [Guu,junk1,freq,junk2] = powsub3(udata,2, zdep, hgtsens, sampfreq);
% [Gvv,junk1,junk,junk2] = powsub3(vdata,2, zdep, hgtsens, sampfreq);

[Guu, freq] = p_welch( detrend(udata), 512, sampfreq, 512, 256 );
[Gvv, freq] = p_welch( detrend(vdata), 512, sampfreq, 512, 256 );


% determine wave number
htot = zdep+hgtsens
ztosns = -1.*zdep
omega = 2*pi .* freq;
kreal = qkhf(omega, htot)./htot;

%compute depth correction
depcoruv = depthcor_uv(htot, ztosns, kreal, omega)

% create cut off freqency, so noise is not magnified
firstfreq = min(find(freq>=.05));
lastfreq = max(find(freq<=0.20));

% combine horizontal velocity spectra
Guv = Guu + Gvv;

% Determine wave height for combined spectra
Gnn = Guv./(depcoruv.^2);

% Determine significant wave height
area = inter(Gnn(firstfreq:lastfreq), freq(2)-freq(1));
hgt = (2.*(sqrt(2)))*(sqrt(area));

%fini

function [area] = inter(psd, dt)

% determine the area under the curve of a psd (integral) using the psd curve 
% and freqency resolution.
% Meg Palmsten, USF/USGS, St. Pete, 5/99.

area1 = psd * dt;
line = length(psd);
area = sum(area1(2:line));

% done