% a2_analyze_gases_ship_all_legs_8192.m
%{
    Loads trace gas data from the Picarro, TEI 49C and fast ozone
    analyzers for entire cruise.  Averages to flux window timebase.

    Run a1_flux_windows_8192.m first to define flux windows.
    Run concat_1min_co2_o3.m to 1-min gas concs for all legs.

    Adds the following columns to the flx table:
        
        13   Pcav     cavity pressure, torr
        14   Tcav     cavity temperature, °C
        15   ch4w     methane mole ratio, ppm, wet
        16   ch4d     methane mole ratio, ppm, dry
        17   co2w     carbon dioxide mole ratio, ppm, wet
        18   co2d     carbon dioxide mole ratio, ppm, dry
        19   h20      water vapor mole ratio, %
        20   o3slo    TEI 49C mean ozone mole ratio, ppb
        21   o3fst    mean fast ozone mole ratio, ppb
        22   dco2_dt  rate of change in co2, ppm/hr
        23   dch4_dt  rate of change in ch4, ppb/hr
        24   do3_dt   rate of change in o3, ppb/hr
        25   sig_co2  std dev co2, ppm
        26   sig_ch4  std dev ch4, ppb
        27   sig_o3   std dev o3, ppb
%}

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

%% adjust these parameters as desired, then run
root_data_path = '/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019';
path_prog = fullfile(root_data_path,'_matlab');
data_path = fullfile(root_data_path,'_processed','ship_flux_V2');
plot_path = fullfile(root_data_path,'_processed','ship_flux_V2','flux_8192','plots');
mkdir(plot_path)
save_path = fullfile(root_data_path,'_processed','ship_flux_V2','flux_8192');
mkdir(save_path)

plotit = 1; % set to 1 to save plot

%% set matlab script path
path(pathdef);
cd(path_prog);
addpath(genpath(path_prog));

%% load flux window table
load(fullfile(save_path,'wins_8192_all.mat'));  % load flux windows table 'flx'
vars = flx.Properties.VariableNames;  % cell array of variable names in flx
% add new variable names
%       1-12   13     14     15     16     17     18     19     20     21      22        23
vars = [vars,'Pcav','Tcav','ch4w','ch4d','co2w','co2d','h2o','o3slo','o3fst','dco2_dt','dch4_dt',...
             'do3_dt','sig_co2','sig_ch4','sig_o3'];
%               24       25        26        27
www = table2array(flx);               % convert to array
www = [www,NaN(length(www),15)];      % add columns for new variables

%% load 1 min Picarro and ozone data
load(fullfile(data_path,'co2_ch4_ship_1min_all_legs.mat')); % load co2_ch4_ship table
load(fullfile(data_path,'o3fast_ship_1min_all_legs.mat'));
load(fullfile(data_path,'o3slow_ship_1min_all_legs.mat'));

%% process picarro data
%{
    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(7,2);
bins(:,1) = 600:4096:34800-8192;
bins(:,2) = bins(:,1) + 8191;
bins = (bins-1)/36000;          % bin edges as fraction-of-hour

disp('processing picarro')
for ii = 1:Nbins:length(www) % step through all flux windows, one hour at a time
    hr = datenum([www(ii,3:6),0,0])*24;  % base hour-datenum for this hour
    for jj=1:Nbins                       % for 7 windows in 1 hour of data
        t1 = (hr+bins(jj,1))/24;         % start datenum for averaging
        t2 = (hr+bins(jj,2))/24;         % end datenum for averaging
        kk = find(co2_ch4_ship.datenum>=t1 & co2_ch4_ship.datenum<t2);
        if isempty(kk)
            continue
        end
        % take mean of picarro variables over this window
        www(ii+jj-1,13:19) = nanmean1(table2array(co2_ch4_ship(kk,[3,4,6:10])));
        
        % compute dco2/dt and dch4/dt
        tt = co2_ch4_ship.datenum(kk)*24 - hr; % time in fractional hours
        co2 = co2_ch4_ship.co2dry(kk);
        ch4 = co2_ch4_ship.ch4dry(kk);
        p = polyfit(tt,co2,1);
        www(ii+jj-1,22) = p(1);  % dco2/dt in ppm/hour
        p = polyfit(tt,ch4*1000,1);
        www(ii+jj-1,23) = p(1);  % dch4/dt in ppb/hour
        % add std dev of co2 / ch4
        www(ii+jj-1,25) = nanstd1(co2);         % std co2 ppm
        www(ii+jj-1,26) = nanstd1(ch4)*1000;    % std ch4 ppb
    end
end

%% process slow ozone
disp('processing TEI ozone')
for ii = 1:Nbins:length(www) % step through all flux windows, one hour at a time
    hr = datenum([www(ii,3:6),0,0])*24;  % base hour-datenum for this hour
    for jj=1:Nbins                       % for 7 windows in 1 hour of data
        t1 = (hr+bins(jj,1))/24;         % start datenum for averaging
        t2 = (hr+bins(jj,2))/24;         % end datenum for averaging
        kk = find(slo3_ship.datenum>=t1 & slo3_ship.datenum<t2);
        if isempty(kk)
            continue
        end
        % take mean of picarro variables over this window
        www(ii+jj-1,20) = nanmean1(table2array(slo3_ship(kk,2)));
    end
end

%% process fast ozone
disp('processing fast ozone')
for ii = 1:Nbins:length(www) % step through all flux windows, one hour at a time
    hr = datenum([www(ii,3:6),0,0])*24;  % base hour-datenum for this hour
    for jj=1:Nbins                       % for 7 windows in 1 hour of data
        t1 = (hr+bins(jj,1))/24;         % start datenum for averaging
        t2 = (hr+bins(jj,2))/24;         % end datenum for averaging
        kk = find(fo3_ship.datenum>=t1 & fo3_ship.datenum<t2);
        if isempty(kk)
            continue
        end
        % take mean of picarro variables over this window
        www(ii+jj-1,21) = nanmean1(table2array(fo3_ship(kk,3)));
        % compute do3/dt and dch4/dt
        tt = fo3_ship.datenum(kk)*24 - hr; % time in fractional hours
        o3 = fo3_ship.O3_ppb(kk);
        p = polyfit(tt,o3,1);
        www(ii+jj-1,24) = p(1);  % do3/dt in ppb/hour
        % add std dev of o3
        www(ii+jj-1,27) = nanstd1(o3);  % std o3 ppb
    end
end

%% save flx table with additional columns
clear flx
flx = array2table(www,'VariableNames',vars);
sfl = fullfile(save_path,'wins_gases_8192_all.mat');
save(sfl,'flx');

%%  a couple plots
ii = flx.co2_flo>25;
jj = flx.o3_flo>25;
dt = datetime(flx.datenum,'convertfrom','datenum') + seconds(819.2/2); % window midpoints

figure('position',[500,500,1000,600],'Name','co2_ch4');
yyaxis left; plot(dt(ii),flx.co2d(ii),'.','markersize',6);grid minor;
ylabel('CO_2 ppm','fontsize',12); xlabel('date / time','fontsize',12); ylim([390,430]);
yyaxis right; plot(dt(ii),flx.ch4d(ii),'.','markersize',6);
ylabel('CH_4 ppm'); ylim([1.92,2.1]);
title(['MOSAiC Ship CO_2 / CH_4'],'fontsize',18);
if plotit
    ppath = fullfile(plot_path,'co2_ch4_8192.png');
    print('-dpng',ppath);
    ppath = fullfile(plot_path,'co2_ch4_8192.fig');
    save(ppath);
end

figure('position',[500,500,1000,600],'Name','ozone');
plot(dt(jj),flx.o3fst(jj),'.','markersize',6);grid minor;
ylabel('O_3 ppb','fontsize',12); xlabel('date / time','fontsize',12); ylim([0,60]);
title(['MOSAiC Ship Ozone'],'fontsize',18);
if plotit
    ppath = fullfile(plot_path,'o3_8192.png');
    print('-dpng',ppath);
    ppath = fullfile(plot_path,'o3_8192.fig');
    save(ppath);
end
