runtime = datestr(now);
clear all
close all

plotLinear = 1; % 1 for linear amplitude plot anything else for amplitude in dB 

filename = {'WaveSpectra1996.txt'; 'WaveSpectra1997.txt'; 'WaveSpectra1998.txt'; 'WaveSpectra1999.txt'; 'WaveSpectra2000.txt'; ...
            'WaveSpectra2001.txt'; 'WaveSpectra2001b.txt'; 'WaveSpectra2002.txt'; 'WaveSpectra2003.txt'; ...
            'WaveSpectra2004.txt'; 'WaveSpectra2005.txt'; 'WaveSpectra2005b.txt'; 'WaveSpectra2006.txt'; ...
            'WaveSpectra2007.txt'; 'WaveSpectra2008.txt'; 'WaveSpectra2009.txt'}; 
        
offset1 = 0;
offset2 = 0;
offset3 = 0;

figure
hold
for i = 1:length(filename)    
    fprintf(1, 'Processing file %s\n', char(filename(i)));
    fileData = importdata(char(filename(i)), ' ', 1);
    
    headerString = char(fileData.colheaders);
    firstDataCol = 1;
    while length(str2num(headerString(firstDataCol,:))) == 0
        firstDataCol = firstDataCol + 1;
    end
    freqVector = str2num(char(fileData.colheaders(firstDataCol:end)));

    if fileData.data(1) < 100
        timestamp = datenum( (1900 + fileData.data(:,1)), fileData.data(:,2), (fileData.data(:,3)+ fileData.data(:,4)/24.0));
    else
        timestamp = datenum( fileData.data(:,1), fileData.data(:,2), (fileData.data(:,3)+ fileData.data(:,4)/24.0));    
    end
    
    spec = fileData.data(:,firstDataCol:end);
    [r, c] = find(spec == 999.00);
    spec(r,c) = NaN;
    
    %Matlab mean function doesn't deal with NaN so this section does
    for j=1:12
        rows = find(fileData.data(:,2) == j);
        monthlySpec = spec(rows,:);
        for k=1:length(freqVector)
           temp = monthlySpec(:,k);
           nanRows = find( isnan(temp) == 1 );
           temp(nanRows) = [];
           yearlyMeanSpec(j,k) = mean(temp);
        end
    end
    
    for j=1:length(spec)
        temp = spec(j,:);
        cols = isnan(temp);
        temp(cols) = [];
        specMeanValue(j) = mean(temp);
    end
        
    WSStruct(i) = struct('timeStamp',timestamp,'freqVector',freqVector,'spectrum',spec, ...
                         'specMeanValue',specMeanValue,'meanSpectrum',yearlyMeanSpec);
          
    clear fileData spec temp monthlySpec yearlyMeanSpec
    
    if plotLinear == 1
        %waterfall(WSStruct(i).freqVector, WSStruct(i).timeStamp, WSStruct(i).spectrum);
    else
        %waterfall(WSStruct(i).freqVector, WSStruct(i).timeStamp, 20*log10(WSStruct(i).spectrum));
    end
end
%datetick('y',2);
%title('Wave Spectral Values')
%xlabel('Freq (Hz)')
%ylabel('Date')
%zlabel('Power (m^2/Hz)')

clear meanSpec temp
x2 = 0.02:0.001:0.5;
for i=1:length(filename)
    for j = 1:12
        if (sum(isnan(WSStruct(i).meanSpectrum(j,:)) == 0))
            interpSpec(i,j,:) = interp1(WSStruct(i).freqVector, WSStruct(i).meanSpectrum(j,:),x2);
        end
    end
end

%Matlab mean function doesn't deal with NaN so this section does
for i = 1:12
    for j=1:length(interpSpec(1,1,:))
        temp = interpSpec(1:length(filename),i,j);
        nanRows = find( isnan(temp) == 1 );
        temp(nanRows) = [];
        meanSpec(i,j) = mean(temp);
        clear temp 
    end
end
figure
waterfall(x2,1:12,meanSpec)
title('Mean of all Years Spectrum by Month')
xlabel('Freq (Hz)')
ylabel('Month')
zlabel('Power (m^2/Hz)')

figure
hold
for i=1:length(filename)
    waterfall(WSStruct(i).freqVector, 1:12, WSStruct(i).meanSpectrum)
end
title('Every Years Mean Spectrums by Month')
xlabel('Freq (Hz)')
ylabel('Month')
zlabel('Power (m^2/Hz)')

mVectors = {[] [] [] [] [] [] [] [] [] [] [] [] };
for i = 1:12
    for j = 1:length(filename)
        a = datestr(WSStruct(j).timeStamp,5);
        b = str2num(a);
        mRows = find(b == i);
        mVectors{i} = [mVectors{i}; WSStruct(j).specMeanValue(mRows)'];
    end
    nans = isnan(mVectors{i});
    mVectors{i}(nans) = [];
    mMean(i) = mean(mVectors{i});
    mStdDev(i) = std(mVectors{i});
    mMin(i) = min(mVectors{i});
    mMax(i) = max(mVectors{i});
end
x = 1:12;
figure
errorbar(x, mMean, mMin, mMax)
title('Mean, Max, Min of the Mean Power Across All Frequencies')
xlabel('Month')
ylabel('Mean Power (m^2/Hz) Across All Frequencies')

figure
hold
monthNum = 9;
index = 1;
spect4OneMonth = [];
lw = 'none';
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(lw, 'none') ~= 1 )
                fprintf('i = %d   j = %d\n', i, j)
                lw = 'none';
            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 = 25
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 = 2.5;
Tavg = 15;
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',x2,Spm)
titleStr = sprintf('Mean of Month %d Spectrums sorted into %5.0f bins by Mean PSD Value (-. = PM Spectrum)',monthNum, numIncr);
title(titleStr)
xlabel('Frequency (Hz)')
ylabel('PSD (m^2/Hz)')