function buoyCheckJS(day,nKeep)
% Check to see how well the Pierson-Moskowitz spectrum fits the data. 
march2 = '../../waveData/xy029p101_201903020000-201903032359';
march6 = '../../waveData/xy029p101_201903060000-201903062359';
march7 = '../../waveData/xy029p101_201903070000-201903072359';
march8 = '../../waveData/xy029p101_201903080000-201903082359';
switch day
    case 'march2'
        fname = march2;
    case 'march6'
        fname = march6;
    case 'march7'
        fname = march7;
    case 'march8'
        fname = march8;
end
[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);                % Set up array for spectral density.
df = fs/Nfft;                       % Compute the frequency spacing.
for ng = 1:Ng                       % Compute spectral density.
    if ~error(ng)
        H = fft(squeeze(eta(:,:,ng)));  % Compute FFT of sections in group.
        H = H(1:Nfft/2,:);              % Cut off the redundant half.
        epsilon = (2/(Nfft*fs))*abs(H).^2; 
        En(:,ng) = mean(epsilon,2);     % Density is the mean of epsilon.
    end
end 
En = En(2:end,:);
fn =  (1:(Nfft/2-1))'*df;
m0 = sum(       En,1)*df;
%m2 = sum(fn.^2.*En,1)*df;
logE = log(En./m0);
figure(1)
surf(logE)
zlim([-6,6])

if nargin == 1
	nKeep = 1:size(En,2);
end
E   = mean(En(:,nKeep),2,'omitnan');
mm0 = mean(m0(  nKeep)  ,'omitnan');
[~,ind] = max(E);

% initial guesses
fp_0    = fn(ind);
alpha   = 5*mm0*(2*pi*fp_0)^4/9.8^2 ;
sigma_a = 0.07;
sigma_b = 0.09;
gamma   = 3.3;
par_0 = [alpha,fp_0,sigma_a,sigma_b,gamma];
[par_fit,E_fit,f_fit] = fit_JONSWAP(fn,E,par_0);

figure(2)
plot(fn,E,f_fit,E_fit);
legend({'$E(f)$','$E_{\rm JS}(f)$'},'interpreter','latex')
texText('$f$','x')
texText(sprintf('%s: group %d to %d',day,min(nKeep),max(nKeep)),'t')
parName = {'\alpha','f_p','\sigma_a','\sigma_b','\gamma'};
str = '';
for n = 1:5
    str = sprintf('%s$%s=%0.3f$\n',str,parName{n},par_fit(n));
end
texText(str,'g',0.7*max(fn),0.5*max(E))
graphToPDF(sprintf('JS-fit-%s.pdf',day),4,3);

end

