function buoySWH()
fname = '../../waveData/xy029p101_201903020000-201903032359';
[r,stepError] = readBuoyFile(fname); % read data file
fs  = 1.28;                     % the sample rate
N30 = 60*30*fs;                 % number of samples in 30 minutes
Ns  = size(r,1);                % number of samples.
Ng  = floor(Ns/N30);            % number of 30 minute groups of samples.
eta   = r(1:Ng*N30,3);              % Truncate data to full group.
eta   = eta/100;                    % Convert data from cm to meters.
eta   = reshape(eta,N30,Ng);        % Group data.
error = stepError(1:Ng*N30);        % Truncate error data to full group.
error = sum(reshape(error,N30,Ng)); % Group error.
eta   = eta(:,error==0);            % Keep just the error free groups. 
eta   = eta - mean(eta,1);          % remove mean from each group.
m0    = mean(eta.^2,1);             % Compute the mean squared amplitude.
H_m0 = 4*sqrt(m0);                  % Compute H_m0.

Ng      = length(m0);               % Get number of groups.
H_third = NaN(1,Ng);                % Allocate an array to store H_1/3.
for ng = 1:Ng                               % Loop over groups.
   H_third(ng) = comp_H_third(eta(:,ng));   % Compute H_1/3.
end
slope = sum(H_m0.*H_third)/sum(H_m0.^2);    % Compute least squares slope.
figure(1)                                   % Graph H_1/3 versus H_m0.
    x = [0,3];                              % Set range of H_m0 to graph.
    plot(H_m0,H_third,'o',x,x,'r',x,slope*x,'g')
    texText('$H_{m_0}$ [m]','x')
    texText('$H_{1/3}$ [m]','y')
    str = sprintf('$%0.3f~H_{m_0}$',slope);
    legend({'$H_{1/3}$','$H_{m_0}$',str},'interpreter','latex')
    legend('location','northwest')
    graphToPDF('swh.pdf',4,4)
end

function H_third = comp_H_third(eta)
    ind_cross = find(diff(eta>0)==(-1));% Find downward crossing points.
    N = length(ind_cross)-1;            % The number of waves.    
    H = NaN(1,N);                       % Allocate memory for wave heights.
    for n = 1:N
        rng = (ind_cross(n)+1):ind_cross(n+1);  % One wave section of data.
        H(n) = max(eta(rng)) - min(eta(rng));   % Compute wave height.
    end  
    [~,ind] = sort(H,'descend');  
    H = H(ind);                                 % Put in descending order.
    H_third = mean(H(1:round(N/3)));            % Compute H_1/3.
end
