function buoyStats()
fname = '../../waveData/xy029p101_201903020000-201903032359';
[r,stepError] = readBuoyFile(fname);
fs = 1.28;
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.
Nfft = 128;         % the rest assumes N30/Nfft is an integer.
Ne = N30/Nfft;      % number of FFT's in a 30 minute ensemble. 
eta = reshape(r(1:Ng*N30,3)/100,Nfft,Ne,Ng);
error = sum(reshape(stepError(1:Ng*N30),N30,Ng));
eta = eta(:,:,error==0); % keep just the error free groups. 
Ng = size(eta,3);
eta = eta - mean(eta,1); % remove mean from each column.
En = NaN(Nfft/2,Ng);
df = fs/Nfft;
for ng = 1:Ng  % compute spectral density and m0 and m2.
    H = fft(eta(:,:,ng)); % compute FFT of all sections in group
    H = H(1:Nfft/2,:);  % cut off the redundant half of the spectrum
    epsilon = (2/(Nfft*fs))*abs(H).^2;
    En(:,ng) = mean(epsilon,2);
end 
fn =  (0:(Nfft/2-1))'*df;
m0 = sum(       En,1)*df;
m2 = sum(fn.^2.*En,1)*df;
eta = reshape(eta,Nfft*Ne,Ng); % regroup by half hour rather than Nfft

% Compute the mean period for each value of eta threshold
    eta_lim = 0:0.01:1.6;
    meanT = NaN(length(eta_lim),Ng);
    for ng = 1:Ng
        if error(ng)==0
            meanT(:,ng) = meanPeriod(eta(:,ng),eta_lim,fs);
        end
    end
    upsilon = m0.*log(meanT.^2.*(m2./m0));

figure(1) % graph the pdf of eta
    histogram(eta./sqrt(m0),'normalization','pdf','displaystyle','stairs')
    hold on
    u = -4:0.01:4;
    plot(u,sqrt(1/(2*pi))*exp(-0.5*u.^2))
    hold off
    texText('$u = \eta/\sqrt{m_0}$','x')
    texText('probability density','y')
    legend({'P.D.','$\sqrt{1\over 2\pi}e^{-u^2/2}$'},'interpreter','latex')
    graphToPDF('etaPDF.pdf',4,3)

figure(2) % plot Upsilon verus eta
    plot(eta_lim,upsilon)
    hold on
    plot(eta_lim,eta_lim.^2,'g:','linewidth',2)
    hold off
    texText('$\eta$','x')
    texText('$\Upsilon$','y')
    xlim([0,1.6])
    ylim([0,2.6])
    graphToPDF('LHSvsEta.pdf',4,3)

figure(3) % a 2-D histogram of (eta,upsilon)
    eta_m = eta_lim'*ones(1,size(upsilon,2));
    eta_m = reshape(eta_m,[],1);
    upsilon   = reshape(upsilon    ,[],1);
    xEdges = eta_lim-0.001;
    yEdges = 0:0.01:2.5;
    histogram2(eta_m,upsilon,xEdges,yEdges,'displaystyle','tile')
    hold on
    plot(eta_lim,eta_lim.^2,'r:','linewidth',1)
    hold off
    texText('$\eta$','x')
    texText('$\Upsilon$','y')
    texText('Histogram','t')
    graphToPDF('LHS-PDF.pdf',4,4);

figure(4) % a 2-D histogram of (eta,sqrt(upsilon))
    yEdges = 0:0.005:1.6;
    histogram2(eta_m,sqrt(upsilon),xEdges,yEdges,'displaystyle','tile')
    hold on
    plot(eta_lim,eta_lim,':','linewidth',2)
    hold off
    texText('$\eta$','x')
    texText('$\sqrt{\Upsilon}$','y')
    texText('Histogram','t')

figure(5) % a surface plot of upsilon verus eta over time
    surf(upsilon)

end

function T = meanPeriod(eta,eta_lim,fs)
% computes the mean period of eta for each threshold value in the 
% vector eta_lim
    N = length(eta_lim);
    T = NaN(1,N);
    for n = 1:N
        ind = find(diff(eta>eta_lim(n))==(-1)); % find level crossings
        NT = length(ind)-1;  % the number of periods found
        if NT>6 % if there is enough periods for trusty mean
            T(n) = (1/fs)*(ind(end)-ind(1))/NT;
        end
    end
end