close all

figure
hold
monthNum = 9;
index = 1;
spect4OneMonth = [];
lastwarn('');
for i=1:length(filename)
    a = datestr(WSStruct(i).timeStamp,5);
    b = str2num(a);
    mRows = find(b == monthNum);
    if length(mRows) > 0
        %waterfall(WSStruct(i).freqVector,WSStruct(i).timeStamp(mRows),WSStruct(i).spectrum(mRows,:))
    end
    for j=1:length(mRows)
        if (sum(isnan(WSStruct(i).spectrum(mRows(j),:)) == 0))
            spect4OneMonth(index,:) = interp1(WSStruct(i).freqVector, WSStruct(i).spectrum(mRows(j),:), x2,'linear', 'extrap');
            lw = lastwarn;
            if ( strcmp(lastwarn, '') ~= 1 )
                fprintf('i = %d   mRows(j) = %d  index = %d \n', i, mRows(j), index);
                lastwarn('');
            end
            MSV4OneMonth(index) = WSStruct(i).specMeanValue(mRows(j));
        end
        index = index + 1;
    end
end
% datetick('y',2)
% titleStr = sprintf('Month %d Spectrums',monthNum);
% title(titleStr)
% xlabel('Freq (Hz)')
% ylabel('Date')
% zlabel('Power (m^2/Hz)')

windSpeed_kts = 19
windSpeed_mps = windSpeed_kts * 1852.0 / 3600

g = 9.8;
U195 = windSpeed_mps;
T1 = 10;

f =  0.02:.001:0.5;
w = 2 * pi * f;
alpha = 8.1e-3;
beta = 0.74;
w0 = g / U195;
Spm = (2*pi)*((alpha * g^2) ./ (w.^5)) .* exp((-1*beta)*((w0./w).^4));

Hchar = 3.0;
Tavg = 9.0;
A = 172.75 * (Hchar^2 / Tavg^4);
B = 691 / Tavg^4
Sittc = (2*pi) * (A ./ w.^5) .* exp(-B./w.^4);

freqMeans = mean(spect4OneMonth);
freqMaxs = max(spect4OneMonth);

figure
plot(x2,Spm, x2,freqMeans, x2,freqMaxs)
titleStr = sprintf('Month %d Max of all Freqs, PM Spectrum (%4.1f kts), Mean of all Freqs',monthNum, windSpeed_kts);
title(titleStr)
xlabel('Freq (Hz)')
ylabel('Power (m^2/Hz)')

figure
plot(x2, 20*log10(Spm),x2,20*log10(freqMeans),x2,20*log10(freqMaxs))
titleStr = sprintf('Month %d Max of all Freqs, PM Spectrum (%4.1f kts), Mean of all Freqs',monthNum, windSpeed_kts);
title(titleStr)
xlabel('Freq (Hz)')
ylabel('Power dB(m^2/Hz)')

[sortedMSV, index] = sort(MSV4OneMonth);
numIncr = 100.0
incr = ceil(length(sortedMSV) / numIncr);
j = 1;
for i=1:incr:length(MSV4OneMonth)
    %fprintf(1,'%d %d\n',i, i+incr-1);
    temp = spect4OneMonth(index(i:min(length(MSV4OneMonth),i+incr-1)),:);
    sortedSpec(j,:) = mean(temp);
    j = j + 1;
end
figure
plot(x2,sortedSpec(end-5:end,:)',x2,Spm,'-.',x2,Sittc,':')
titleStr = sprintf('Mean of Month %d Spectrums sorted into %5.0f bins by Mean PSD Value\n(-. = PM Spectrum %5.0f m/s   .. = ITTC Spectrum %5.1f m Hchar %5.1f sec Tavg)' ...
    ,monthNum, numIncr, windSpeed_mps, Hchar, Tavg);
title(titleStr)
xlabel('Frequency (Hz)')
ylabel('PSD (m^2/Hz)')

figure
plot(x2,20.*log10(sortedSpec(end-3:end,:))',x2,20.*log10(Spm),'-.',x2,20.*log10(Sittc),':')
titleStr = sprintf('Mean of Month %d Spectrums sorted into %5.0f bins by Mean PSD Value\n(-. = PM Spectrum %5.0f m/s   .. = ITTC Spectrum %5.1f m Hchar %5.1f sec Tavg)' ...
    ,monthNum, numIncr, windSpeed_mps, Hchar, Tavg);
title(titleStr)
xlabel('Frequency (Hz)')
ylabel('PSD dB (m^2/Hz)')

medianIncr = ceil(incr/2):incr:length(MSV4OneMonth);
medianIncr = 
medianSpec = spect4OneMonth(index(medianIncr),:);

figure
plot(x2,medianSpec(end-5:end,:)',x2,Spm,'-.',x2,Sittc,':')
titleStr = sprintf('Median of Month %d Spectrums sorted into %5.0f bins by Mean PSD Value\n(-. = PM Spectrum %5.0f m/s   .. = ITTC Spectrum %5.1f m Hchar %5.1f sec Tavg)' ...
    ,monthNum, numIncr, windSpeed_mps, Hchar, Tavg);
title(titleStr)
xlabel('Frequency (Hz)')
ylabel('PSD (m^2/Hz)')

figure
plot(x2,20.*log10(medianSpec(end-3:end,:))',x2,20.*log10(Spm),'-.',x2,20.*log10(Sittc),':')
titleStr = sprintf('Median of Month %d Spectrums sorted into %5.0f bins by Mean PSD Value\n(-. = PM Spectrum %5.0f m/s   .. = ITTC Spectrum %5.1f m Hchar %5.1f sec Tavg)' ...
    ,monthNum, numIncr, windSpeed_mps, Hchar, Tavg);
title(titleStr)
xlabel('Frequency (Hz)')
ylabel('PSD dB (m^2/Hz)')