%{
    hw_windsea_CAPRICORN_2016.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.

    REQUIRES 10MIN FINAL FLUX RESULTS FILE (da_red scripts must be
    run first...)

    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

% matlab script path
path_prog = fullfile('c:\Users\Sergio\Documents\MATLAB\2016_CAPRICORN\Riegl_Hs_fixed');
cd(path_prog);                  % cd to matlab directory
restoredefaultpath              % reset matlab path
addpath(genpath(path_prog));	% add scripts to path
savepath; rehash;

root_path = 'D:\Data\CAPRICORN_2016\Investigator\flux\Processed\';

% load 10 minute flux result array, flux10, N x 61 array
load(fullfile(root_path,'CAPRICORN_2016_flux_10_v4_decorr.mat'));
% load wave spectra, jd 258 to 301, 5472 x 241 array
sp = load(fullfile(root_path,'riegl_wave_data','CAPRICORN_2016_WaveSpectra.txt'));
% load wave stats, jd 258 to 301, 5472 x 17 array
st = load(fullfile(root_path,'riegl_wave_data','CAPRICORN_2016_WaveStats.txt'));

% define variables of interest
jd = sp(:,1);
lat = flux10(:,34);     % latitude
u10n = flux10(:,59);    % 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
hw(hw==0) = NaN;
Riegl_1D_10min = [st,fw,u10n,hw];  % append results to st array

figure; plot(jd,st(:,3),'r.-',jd,hw,'b.-'); grid;
xlim([min(jd)-0.5, max(jd)+0.5]);
xlabel('day-of-year'); ylabel('Sig Wave Ht, m');
legend('Sig Ht', 'Sig Ht Wind Sea','location','best');
title('Wave Height Time Series');
ppath = fullfile(root_path,'riegl_wave_data','riegl_wave_height.png');
print('-dpng',ppath);

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

% average to 1 hr
jdhr = floor(jd(1)):1/24:ceil(jd(end)); jdhr = jdhr(:);
Riegl_1D_hr = interval_avg(Riegl_1D_10min(:,1),Riegl_1D_10min(:,2:end),jdhr);

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

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




