function buoyCheckPM(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);
fp_0 = fn(ind);
fp_fit = fit_Epm(fn,E,mm0,fp_0);
Epm_fit = PM_spec(fn,fp_fit,mm0);

figure(2)
plot(fn,E,fn,Epm_fit);
legend({'$E(f)$','$E_{\rm PM}(f)$'},'interpreter','latex')
texText('$f$','x')
texText(sprintf('%s: group %d to %d',day,min(nKeep),max(nKeep)),'t')
graphToPDF(sprintf('PM-fit-%s.pdf',day),4,3);

end

function Epm = PM_spec(f,fp,m0)
    u = fp./f;
    Epm = (5*m0/fp)*u.^5.*exp(-1.25*u.^4);
end

function fp_fit = fit_Epm(f,E,m0,fp_0)
% This function finds the Pierson-Moskowitz spectrum that minimizes the
% mean squared error between the observed data (f,E) and the model. This is
% accompished by finding that parameter fp that minimizes the error.
    fp_fit = fminsearch(@compMSE,fp_0);
    function mse = compMSE(fp)
        E_fit = PM_spec(f,fp,m0);
        mse = sum((E-E_fit).^2);
    end
end