%{
    hw_windsea_HiWinGS_2013.m

    computes estimate for significant wave height of the wind sea
    from smoothed wave spectra and U10 wind speed.  Spectra are
    from the fixed beam Riegl rangefinder.

    riegl data for julian dates 258 - 301 only

    output array [fw, u10n, hw] is appended to the stats array.
    fw is separation frequency, and hw is estimated sig wave ht
    for the wind sea component.

    spectrum data file format:
    columns     variable                description
    1           Decimal Day-of-Year     timestamp
    2:79        Freq                    Frequency bin array, Hz
    80:157      Df                      Delta F for bins, Hz
    158:235     Sw                      Binned spectral data, m^2/Hz

    BWB Feb 2016
%}

clear
close all
root_path = '/Users/byronb/Raw_Data_Files.noindex/HiWinGS_2013/Knorr/flux/Processed_v3/';

% load 10 minute flux result array, flux10, 7776 x 68 array
load(fullfile(root_path,'HiWinGS_2013_flux_10_v3_decorr.mat'));
xx = find(flux10(:,1)>=258 & flux10(:,1)<302);
flux10 = flux10(xx,:);  % trim to jd 258 to 301 only
% load wave spectra, jd 258 to 301, 5472 x 241 array
sp = load(fullfile(root_path,'riegl_wave_data','HiWinGS_2013_WaveSpectra.txt'));
% load wave stats, jd 258 to 301, 5472 x 17 array
st = load(fullfile(root_path,'riegl_wave_data','HiWinGS_2013_WaveStats.txt'));

% define variables of interest
jd = sp(:,1);
lat = flux10(:,34);     % latitude
u10n = flux10(:,68);    % u10n, COARE 3.5
fwv = sp(:,2:79);       % bin frequencies, Hz
dfw = sp(:,80:157);     % bin delta f
Swv = sp(:,158:235);    % spectra of high-pass filtered time series, Swv(f)

fw = 0.8*grv(lat)./(2*pi*u10n);             % wind sea separation frequency

hw = NaN(length(fw),1);                     % compute sig wave ht for wind sea
for jj = 1:length(fw)
    ii = find(fwv(jj,:)>=fw(jj),1,'first'); % start bin for integrations
    hw(jj) = 4*sqrt(sum(Swv(jj,ii:end).*dfw(jj,ii:end)));
end;

hw(hw>15) = NaN;        % kill bad values
out = [st,fw,u10n,hw];  % append results to st array

% save data file
save_path = fullfile(root_path,'riegl_wave_data','HiWinGS_2013_WaveStats_Hw.txt');
save(save_path,'out','-ascii');
save_path = fullfile(root_path,'riegl_wave_data','HiWinGS_2013_WaveStats_Hw.mat');
save(save_path,'out');

% average to 1 hr
jdhr = 258:1/24:302; jdhr = jdhr(:);
out_hr = interval_avg(out(:,1),out(:,2:end),jdhr);
out_hr(457:600,2:end) = NaN;  % in port, no data

% save data file
save_path = fullfile(root_path,'riegl_wave_data','HiWinGS_2013_WaveStats_Hw_hr.txt');
save(save_path,'out_hr','-ascii','-double');
save_path = fullfile(root_path,'riegl_wave_data','HiWinGS_2013_WaveStats_Hw_hr.mat');
save(save_path,'out_hr');

% save spectra as mat
save_path = fullfile(root_path,'riegl_wave_data','HiWinGS_2013_WaveSpectra.mat');
save(save_path,'sp');




