% a1_flux_windows_8192.m

%{
    Sets up 8192 point windows for the ship flux systems over the entire cruise.
    Computes  mean rel wind dir and inlet flows for each window from the 1-min bow QC data.

    Run MOSAiC_bow_QC_1min.m before this script.
%}

%% initialize
clear
close all
warning('off','MATLAB:MKDIR:DirectoryExists');

root_data_path = '/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019';

path(pathdef);
path_prog = fullfile(root_data_path,'_matlab');
cd(fullfile(path_prog,'co2_ch4_o3_flux'));
addpath(genpath(path_prog));

save_path = fullfile(root_data_path,'_processed','ship_flux_V2','flux_8192');
mkdir(save_path);
QC_path = fullfile(root_data_path,'_processed','Bow_Flux_Inlet_QC_V2');
fig_path = fullfile(save_path,'plots'); % for saving plots
mkdir(fig_path);

%% load QC arrays and concatenate for legs 1-5

load(fullfile(QC_path,'Leg1-3_flux_inlet_QC_1min_V2.mat'));
temp = flux_inlet_QC; clear flux_inlet_QC
load(fullfile(QC_path,'Leg4-5_flux_inlet_QC_1min_V2.mat'));
QC = cat_tables(temp,flux_inlet_QC);
clear temp flux_inlet_QC
dt1 = datetime(QC.datenum,'convertfrom','datenum');  % 1min timestamp
time_vec = datevec(dt1);  % array of yyyy,mm,dd,hh,mn,ss for 1-min timestamp
% figure; plot(dt1,flux_inlet_QC.rwdir); grid minor;

%% define flux windows for entire cruise

%{
    Define start/stop indices for 7 overlapping 8192 point windows in 1 hr of 10Hz data (36000 total).
    Windows are 13.653 minutes (2^13 data points each).
             600: 8791
            4696:12887
            8792:16983
           12888:21079
           16984:25175
           21080:29271
           25176:33367
    Flux will be computed for these windows each hour.
    The windows specifically exclude the last two minutes and first minute of each hour
    due to the 'puff' which happens some time in this range.
%}
Nbins = 7;
bins = zeros(Nbins,2);
bins(:,1) = 600:4096:34800-8192;
bins(:,2) = bins(:,1) + 8191;

dnum_start = floor(datenum(time_vec(1,:)));
dnum_end = floor(datenum(time_vec(end,:)));
N = (dnum_end - dnum_start)*24*Nbins;
% empty output array for time stamps and bin edges
% for each window:[leg,dnum,yyyy,mm,dd,hh,bin_st,bin_end,rwdir,rwspd,co2_flo,o3_flo]
www = NaN(N,12);
ii = 1;

% build array of flux window variables for entire cruise
for jj=dnum_start:dnum_end
    for kk = 0:23
        ttt = NaN(Nbins,12);
        % compute timestamp columns for 8192 point windows
        % ttt = dnum timestamps for begin/end of each window this hour
        ttt(:,2) = jj + kk/24 + bins(:,1)/36000/24;
        % set leg number for this date/time
        if ttt(1,2)>=datenum(2019,10,17) && ttt(1,2)<datenum(2019,12,14)
            ttt(:,1) = 1; % leg 1
        elseif ttt(1,2)>=datenum(2019,12,14) && ttt(1,2)<datenum(2020,02,25)
            ttt(:,1) = 2; % leg 2
        elseif ttt(1,2)>=datenum(2020,02,25) && ttt(1,2)<datenum(2020,06,01)
            ttt(:,1) = 3; % leg 3
        elseif ttt(1,2)>=datenum(2020,06,01) && ttt(1,2)<datenum(2020,08,14)
            ttt(:,1) = 4; % leg 4
        elseif ttt(1,2)>=datenum(2020,08,14) && ttt(1,2)<datenum(2020,10,02)
            ttt(:,1) = 5; % leg 5
        end
        % insert yyyy,mm,dd,hh for these windows
        ddd = datevec(ttt(:,2));
        ttt(:,3:6) = ddd(:,1:4);
        % insert begin/end 10Hz indicies for these windows
        ttt(:,7:8) = bins;
        % evaluate QC data for each window
        for mm = 1:Nbins
            t0 = floor(ttt(mm,2)*1440)/1440;  % first minute of window
            xx = find(QC.datenum>=t0 & QC.datenum<=t0+13/1440);  % QC data for this window
            if isempty(xx)
                continue
            end
            rd = nanmean1(QC.rwdir_filt(xx));
            rs = nanmean1(QC.rwspd_filt(xx));
            co2f = nanmean1(QC.co2_flow(xx));
            o3f = nanmean1(QC.fo3_flow(xx));
            ttt(mm,9:12) = [rd,rs,co2f,o3f];
        end
        % save to flx array if QC is good
        if all(isnan(ttt(:,9)))
            continue
        else
            www(ii:ii+Nbins-1,:) = ttt;
            ii = ii+Nbins;
        end
    end
end

% remove NaN rows at the end of flx array
xx = find(isnan(www(:,1)),1,'first');
www = www(1:xx-1,:);

% save flx as mat file
vars = {'leg','datenum','yyyy','mm','dd','hh','seg_begin','seg_end','rwdir','rwspd','co2_flo','o3_flo'};
flx = array2table(www,'VariableNames',vars);
sfl = fullfile(save_path,'wins_8192_all.mat');
save(sfl,'flx');

% check relative wind direction plot
dt = datetime(www(:,2),'convertfrom','datenum');
figure; plot(dt,www(:,9)); grid minor;

