% b3_analyze_gases_metcity_8192.m
%{
        Loads trace gas data from the MetCity Picarro for entire cruise.
        Averages to flux window timebase.

        Run b1_hourly_met_metcity_8192.m first to define flux windows.
        Run concat_1min_metcity_co2.m to get 1-min gas concs for all legs.

        Table mc_flx will have the following columns:

        1     leg                       cruise leg number
        2     datenum                   matlab datenum, start of interval
        3     yyyy                      year
        4     mm                        month
        5     dd                        day of month
        6     hh                        hour of day, 0-23
        7     seg_begin                 segment begin point (36000 for 1 hr)
        8     seg_end                   segment end point (36000 for 1 hr)
        9     lat_tower                 columns 9-100 are met city tower data averaged over flux windows
        10    lon_tower
        11    heading_tower
        12    heading_mast
        13    lat_mast
        14    lon_mast
        15    zenith_true
        16    zenith_apparent
        17    azimuth
        18    ship_bearing
        19    ship_distance
        20    sr50_dist
        21    temp_2m
        22    temp_6m
        23    temp_10m
        24    temp_mast
        25    dew_point_2m
        26    dew_point_6m
        27    dew_point_10m
        28    dew_point_mast
        29    rh_2m
        30    rh_6m
        31    rh_10m
        32    rh_mast
        33    atmos_pressure_2m
        34    atmos_pressure_mast
        35    brightness_temp_surface
        36    subsurface_heat_flux_A
        37    subsurface_heat_flux_B
        38    snow_depth
        39    mixing_ratio_2m
        40    mixing_ratio_6m
        41    mixing_ratio_10m
        42    mixing_ratio_mast
        43    rhi_2m
        44    rhi_6m
        45    rhi_10m
        46    rhi_mast
        47    vapor_pressure_2m
        48    vapor_pressure_6m
        49    vapor_pressure_10m
        50    vapor_pressure_mast
        51    wspd_vec_mean_2m
        52    wspd_vec_mean_6m
        53    wspd_vec_mean_10m
        54    wspd_vec_mean_mast
        55    wdir_vec_mean_2m
        56    wdir_vec_mean_6m
        57    wdir_vec_mean_10m
        58    wdir_vec_mean_mast
        59    h2o_licor
        60    co2_licor
        61    temp_licor
        62    pressure_licor
        63    wspd_u_mean_2m
        64    wspd_u_mean_6m
        65    wspd_u_mean_10m
        66    wspd_u_mean_mast
        67    wspd_v_mean_2m
        68    wspd_v_mean_6m
        69    wspd_v_mean_10m
        70    wspd_v_mean_mast
        71    wspd_w_mean_2m
        72    wspd_w_mean_6m
        73    wspd_w_mean_10m
        74    wspd_w_mean_mast
        75    temp_acoustic_mean_2m
        76    temp_acoustic_mean_6m
        77    temp_acoustic_mean_10m
        78    temp_acoustic_mean_mast
        79    wspd_u_std_2m
        80    wspd_v_std_2m
        81    wspd_w_std_2m
        82    temp_acoustic_std_2m
        83    wspd_u_std_6m
        84    wspd_v_std_6m
        85    wspd_w_std_6m
        86    temp_acoustic_std_6m
        87    wspd_u_std_10m
        88    wspd_v_std_10m
        89    wspd_w_std_10m
        90    temp_acoustic_std_10m
        91    wspd_u_std_mast
        92    wspd_v_std_mast
        93    wspd_w_std_mast
        94    temp_acoustic_std_mast
        95    skin_temp_surface
        96    down_long_hemisp
        97    down_short_hemisp
        98    up_long_hemisp
        99    up_short_hemisp
        100   net_radiation
        101   wspd_10m                  10m true windspeed, m/s
        102   wdir_10m                  10m true wind direction, deg
        103   Pcav                      Picarro cavity pressure, torr
        104   Tcav                      Picarro cavity temperature, °C
        105   ch4w                      methane wet mole fraction, ppm
        106   ch4d                      methane dry mole fraction, ppm
        107   co2w                      carbon dioxide wet mole fraction, ppm
        108   co2d                      carbon dioxide wet mole fraction, ppm
        109   h2o                       water vapor mole fraction, %
        110   dco2_dt                   co2 rate of change, ppm/hr
        111   dch4_dt                   ch4 rate of change, ppb/hr
        112   sig_co2                   std dev co2 dry mole fraction, ppm
        113   sig_ch4                   std dev ch4 dry mole fraction, 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','metcity_flux_v2');
plot_path = fullfile(root_data_path,'_processed','metcity_flux_v2','flux_8192','plots');
mkdir(plot_path)
save_path = fullfile(root_data_path,'_processed','metcity_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_met_8192_all.mat'));  % load flux windows table 'mc_flx'
vars = mc_flx.Properties.VariableNames;  % cell array of variable names in flx
% add new variable names
%      1-102   103    104    105    106    107    108   109     110       111       112       113
vars = [vars,'Pcav','Tcav','ch4w','ch4d','co2w','co2d','h2o','dco2_dt','dch4_dt','sig_co2','sig_ch4'];

www = table2array(mc_flx);           % convert to array
www = [www,NaN(length(www),11)];     % add columns for new variables

% load 1 min Picarro
load(fullfile(data_path,'co2_ch4_metcity_1min_all_legs.mat')); % load co2_ch4_metcity table
co2_ch4_metcity = removevars(co2_ch4_metcity,'Valve1');

%% 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;  % flux windows per hour
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

for ii = 1:Nbins:length(www) % step through all flux windows, one hour at a time
    display(['window count ',sprintf('%05i',ii)]);
    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_metcity.datenum>=t1 & co2_ch4_metcity.datenum<t2);
        if isempty(kk)
            continue
        end
        % take mean of picarro variables over this window
        www(ii+jj-1,103:109) = nanmean1(table2array(co2_ch4_metcity(kk,3:9)));

        % compute dco2/dt and dch4/dt
        tt = co2_ch4_metcity.datenum(kk)*24 - hr; % time in fractional hours
        co2 = co2_ch4_metcity.co2dry(kk);
        ch4 = co2_ch4_metcity.ch4dry(kk);
        p = polyfit(tt,co2,1);
        www(ii+jj-1,110) = p(1);  % dco2/dt in ppm/hour
        p = polyfit(tt,ch4*1000,1);
        www(ii+jj-1,111) = p(1);  % dch4/dt in ppb/hour
        % add std dev of co2 / ch4
        www(ii+jj-1,112) = nanstd1(co2);         % std co2 ppm
        www(ii+jj-1,113) = nanstd1(ch4)*1000;    % std ch4 ppb
    end
end

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

%%  a couple plots
dt = datetime(mc_flx.datenum,'convertfrom','datenum') + seconds(819.2/2); % window midpoints

figure('position',[500,500,1000,600],'Name','co2_ch4');
yyaxis left; plot(dt,mc_flx.co2d,'.','markersize',6);grid minor;
ylabel('CO_2 ppm','fontsize',12); xlabel('date / time','fontsize',12); ylim([390,430]);
yyaxis right; plot(dt,mc_flx.ch4d,'.','markersize',6);
ylabel('CH_4 ppm'); ylim([1.92,2.1]);
title('MOSAiC MetCity 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
