% a3_analyze_lag_times_ship_8192.m

%{
    Reads wind and trace gas data files analyzes lag times and inlet
    attenuation time constants using the from the N2 'puff'.

    Uses the 'flx' flux window table produced by a2_analyze_gases_ship_all_legs_8192.m
    Only computes lags for hours where QC for the first flux window is good.

    Appends the following variables for each window to the flx table:
        28       puff_lag_co2  lag between valve and start of Picarro gas puff, points
        29       co2_ch4_tau   inlet attn time const for co2 & ch4
        30       co2_lag       'optimal' lag for co2
        31       co2_std       std dev co2 (despiked), ppm
        32       ch4_lag       'optimal' lag for ch4
        33       ch4_std       std dev ch4 (despiked), ppm
        34       o3_puff_lag   number of data points to start of  o3 pulse
        35       o3_tau        inlet attn time const for o3
        36       o3_lag        'optimal' lag for o3
        37       o3_std        std dev o3 (despiked), ppb
        38       max_Corr_co2  maximum correlation coef, co2
        39       max_Corr_ch4  maximum correlation coef, ch4
        40       max_Corr_o3   maximum correlation coef, o3
        41       q_pic         specific humidity, Picarro, g/kg
        42       q_pic_tau     inlet attn time const for Picarro q
        43       q_pic_lag     'optimal' lag for Picarro q
        44       q_pic_std     std dev Picarro q (despiked), g/kg
        45       q_lic         specific humidity, Licor 7500, g/kg
        46       q_lic_std     std dev Licor q, g/kg
        47       q_lic_qc      Licor QC parameter (0-100) no units
        48       son_bad       number of NaN values in each flux window
        49       son_missing   number of missing points in each flux window
        50       son_azm       azimuth angle in each flux window (rel wdir)
        51       son_tilt      tilt angle in each flux window
        52       puff_stpt     array index for start of puff segment
        53       puff_endpt    array index for end of puff segment
        54       tair          air temperature, °C (D-ship)
        55       rh            relative humidity, % (D-ship)
        56       P_mb          air pressure, mb (D-ship)
        57       lat           latitude (D-ship)
        58       lon           longitude (D-ship)
        59       sog           speed-over-ground (D-ship)
        60       cog           course-over-ground (D-ship)
        61       hed           heading (D-ship)
        62       sog_std       std dev sog (D-ship)
        63       hed_std       std dev cog(D-ship)
%}

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

%% define fit function obtaining puff time constants
%     c(t) = fos(t,t0,tau,c0,c1) function should be in
%     the mosaic_subs folder as fos.m
fit_fos = fittype('fos(t,t0,tau,c0,c1)','independent','t','dependent','c');

%% adjust these parameters as desired, then run

plotit = 1;     % set to 1 to save plots
len = 8192;     % number of points in each flux window
wq_cutoff = datenum(2020,04,16);  % last datenum for picarro wq calcs - dryer installed
dms_puff_cutoff = datenum(2019,11,13); % puff control switched from DMS APIMS to Picarro

%% paths
root_data_path = '/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019';
path_prog = fullfile(root_data_path,'_matlab');
picarro_path = fullfile(root_data_path,'crds_ghg_ship_bow');
% picarro_path = fullfile(root_data_path,'crds_ghg_ice_station_tower');
ozone_path = fullfile(root_data_path,'ozone_fast_ship_bow');
lic_path = fullfile(root_data_path,'licor_7500ds_ship_bow');
wind_path = fullfile(root_data_path,'sonic_anemometer_3D_ship_bow');
dms_path = fullfile(root_data_path,'dms_apims_ship_bow','SIM');
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');

%% set matlab script path
path(pathdef);
cd(path_prog);
addpath(genpath(path_prog));
cd(fullfile(path_prog,'co2_ch4_o3_flux'));

%% load flx table and append new columns
load(fullfile(save_path,'wins_gases_8192_all.mat'));  % load flux windows table 'flx'
www = table2array(flx);
www = [www,NaN(length(www),36)];      % add columns for new variables
vars = flx.Properties.VariableNames;  % cell array of variable names in flx
% add new variable names for columns 28-63
vars = [vars,'co2_puff_lag','co2_ch4_tau','co2_lag','co2_std','ch4_lag','ch4_std','o3_puff_lag',...
             'o3_tau','o3_lag','o3_std','max_corr_co2','max_corr_ch4','max_Corr_o3','q_pic',...
             'q_pic_tau','q_pic_lag','q_pic_std','q_lic','q_lic_std','q_lic_qc','son_bad',...
             'son_missing','son_azm','son_tilt','puff_stpt','puff_endpt','tair','rh','P_mb',...
             'lat','lon','sog','cog','hed','sog_std','hed_std'];

% load fo3 calibration coefs for all legs
load(fullfile(save_path,'fo3_cal_coefs_all_legs.mat'));

% load TEI ozone for entire cruise - table slo3_ship: col1 is datenum, col 2 is O3_ppb
load(fullfile(root_data_path,'_processed','ship_flux_V1','o3slow_ship_1min_all_legs.mat'));

%% build list of hourly files to process
N = length(www);
flist_pic0  = cell(N,1);
flist_son0  = cell(N,1);
flist_fo3 = cell(N,1);
flist_lic0 = cell(N,1);

idx = 1;  % cell array index
dnum_start = floor(flx.datenum(1));
dnum_end = floor(flx.datenum(end));

for ii = dnum_start:dnum_end
    dv = datevec(ii);
    ddd = md2yd_vn(dv(1),dv(2),dv(3));
    DDD = sprintf('%03i',ddd);
    YYYY = sprintf('%04i',dv(1)); YY = YYYY(3:4);
    wind_fldr = fullfile(wind_path,[YY,DDD]);
    co2_fldr = fullfile(picarro_path,[YY,DDD]);
    o3_fldr = fullfile(ozone_path,[YY,DDD]);
    lic_fldr = fullfile(lic_path,[YY,DDD]);
    % check for existence of wind data for this day
    if ~exist(wind_fldr,'dir')
        continue    % go to next day
    end
    % if we have wind data, process hours for this day
    for hh = 0:23
        HH = sprintf('%02i',hh);
        wind_dfl = fullfile(wind_fldr,['son0',YY,DDD,HH,'_raw.txt']);
        co2_dfl = fullfile(co2_fldr,['pic0',YY,DDD,HH,'_raw.txt']);
        o3_dfl = fullfile(o3_fldr,['fo30',YY,DDD,HH,'_raw.txt']);
        lic_dfl = fullfile(lic_fldr,['lic0',YY,DDD,HH,'_raw.txt']);
        % check for existence of data for this hour
        if exist(wind_dfl,'file'); flist_son0{idx} = wind_dfl; end
        if exist(co2_dfl,'file'); flist_pic0{idx} = co2_dfl; end
        if exist(o3_dfl,'file'); flist_fo3{idx} = o3_dfl; end
        if exist(lic_dfl,'file'); flist_lic0{idx} = lic_dfl; end
        idx = idx+1;
    end
end

% trim off empty cells in arrays
flist_pic0 = flist_pic0(1:idx-1);
flist_son0 = flist_son0(1:idx-1);
flist_fo3 = flist_fo3(1:idx-1);
flist_lic0 = flist_lic0(1:idx-1);

% mark data gaps for the trace gas files
gap_pic = ones(length(flist_pic0),1);
gap_fo3 = ones(length(flist_fo3),1);
for ii=2:length(flist_pic0)
    if ~isempty(cell2mat(flist_pic0(ii))) && ~isempty(cell2mat(flist_pic0(ii-1)))
        gap_pic(ii) = 0;  % There is data in the prior hour, so set gap to 0
    end
    if ~isempty(cell2mat(flist_fo3(ii))) && ~isempty(cell2mat(flist_fo3(ii-1)))
        gap_fo3(ii) = 0;  % There is data in the prior hour, so set gap to 0
    end
end

%% process files and lag times
%{
    Loop through hourly files in consecutive order for each gas analyzer.
    Lags and tau are only computed for Picarro and fast ozone.
    Basic QC stats are computed for the Licor.

    Because timing for the N2 pulse is not exact, the beginning of the pulse may occur
    near the end of the previous hour rather than after the start of the hour, so
    it's necessary to process the data in two hour blocks and search for the complete
    pulse over a time range of +/- 30 sec or so at the top of the hour.

    Hourly files at the beginning of a data gap are processed as special cases, and we
    look for a puff at the start of that hour and compute lags if found, otherwise
    the first hour file is skipped and lag processing begins with the next hour in the series.

    We process the Picarro, fast ozone and Licor files in three separate loops because
    gaps in each time series occurr at different times, and it's too complex to try and
    account for all of that in one single loop.
%}

%% ------------------------------ begin loop to analyze Picarro puffs -----------------------------
ddd0 = 0; % set current day flag
clear h1 h2 h3 h4 h5

%{
    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;

% gap_pic(157) = 1;  157:168 % test sequence for Nov 1 12-23, flist(157:168)
for kk = 1:length(flist_pic0)

    % get file names for this hour, if file missing set gap and go to next hour
    wind_dfl = cell2mat(flist_son0(kk));
    if isempty(wind_dfl)
        gap_pic(kk+1) = 1;
        continue
    end
    co2_dfl = cell2mat(flist_pic0(kk));
    if isempty(co2_dfl)
        gap_pic(kk+1) = 1;
        continue
    end

    % set up the various date variables and strings from wind file name
    temp = strsplit(wind_dfl,'/');  fname = cell2mat(temp(8));
    YY = fname(5:6); YYYY = ['20',YY]; DDD = fname(7:9); HH = fname(10:11);
    yyyy = str2double(YYYY); ddd = str2double(DDD); hh = str2double(HH);
    [mm,dd] = yd2md(yyyy,ddd); MM = sprintf('%02i',mm); DD = sprintf('%02i',dd);
    dnum_this_hour = datenum(yyyy,mm,dd,hh,0,0);

    % find leg number and www indices for this hour
    zz = find(www(:,3)==yyyy & www(:,4)==mm & www(:,5)==dd & www(:,6)==hh);
    if isempty(zz)         % no flux windows for this hour
        gap_pic(kk+1) = 1; % so set next hour gap flag
        continue           % and skip to next hour
    end
    LEG = sprintf('%01i',www(zz(1),1));

    % load met data
    %{
    metnav is 1440 x 32 array of 1-min met/nav data
    1   jd_min       1-min decimal day of year
    2   wspd0        ship 20m metek true wspd, m/s
    3   wdir0        ship 20m metek true wdir, deg
    4   tson0        ship 20m metek Tsonic, °C
    5   wspd2        met city 10m metek true wspd, m/s
    6   wdir2        met city 10m metek true wdir, deg
    7   tson2        met city 10m metek Tsonic, °C
    8   rwspd1       ship rm young rwspd, m/s
    9   rwdir1       ship rm young rwdir, deg
    10  wspd1        ship rm young true wspd, m/s
    11  wdir1        ship rm young true wdir, deg
    12  yr           year
    13  lat          latitude, decimal deg
    14  lon          longitude, decimal deg
    15  sog          gps speed over ground, m/s
    16  cog          gps course over ground, deg
    17  hed          heading, deg
    18  pitch        pitch angle, deg
    19  roll         roll angle, deg
    20  tair         air temperature, °C
    21  rh           relative humidity, %
    22  dpt          dewpoint, °C
    23  p_mb         air pressure, mb
    24  dir_rad      direct radiation, W/m2
    25  glb_rad      global radiation, W/m2
    26  rwspd        relative wind speed, m/s
    27  rwdir        relative wind direction, deg
    28  wspd         true wind speed, m/s
    29  wdir         true wind direction, deg
    30  vis          visibility, m
    31  tsea         sea surface temperature, °C
    32  ceil         ceiling, m
    %}

    if ddd~=ddd0  % if new day
        % load 1 min ship met for this day and compute hourly avg
        dfl_met = fullfile(root_data_path,'_scientific_analysis','raw_1Hz',...
            ['met_nav_Leg_',LEG],[YYYY,'_',MM,'_',DD,'_',DDD,'_Met_Nav.mat']);
        metnav = [];
        load(dfl_met,'metnav');

        % load ship heading for this day
        % this is used to read the raw ship sonic data
        path_dsp = fullfile(root_data_path,'ship_data','parsed',['dsp0',YY,MM,DD,'_raw.txt']);
        dsp = read_dsp_mosaic_2019(path_dsp,ddd);
        hed0 = [dsp(:,1),dsp(:,7)]; % 1-min ship timestamp and heading for this day

        ddd0 = ddd;  % remember current day

        metnav_dnum = metnav(:,1) + datenum(yyyy-1,12,31);
    end

    % compute mean met/nav data for flux windows this hour
    for zzz = 1:length(zz)
        tm0 = find(metnav_dnum>=www(zz(zzz),2) & metnav_dnum<www(zz(zzz)+1,2));
        www(zz(zzz),54) = nanmean1(metnav(tm0,20)); % mean air temperature
        www(zz(zzz),55) = nanmean1(metnav(tm0,21)); % mean rh
        www(zz(zzz),56) = nanmean1(metnav(tm0,23)); % mean pressure, mb
        www(zz(zzz),57) = nanmean1(metnav(tm0,13)); % mean latitude
        www(zz(zzz),58) = nanmean1(metnav(tm0,14)); % mean longitude
        www(zz(zzz),59) = nanmean1(metnav(tm0,15)); % mean speed-over-ground
        www(zz(zzz),60) = nanmean1(metnav(tm0,16)); % mean course-over-ground
        www(zz(zzz),61) = nanmean1(metnav(tm0,17)); % mean heading
        www(zz(zzz),62) = nanstd1(metnav(tm0,15));  % std dev speed-over-ground
        www(zz(zzz),63) = nanstd1(metnav(tm0,17));  % std dev heading
    end

    % we will process picarro & wind data in 2-hour sequential segments, so we decide here
    % if we will load a new hour in the current series or start a new series after a gap
    if gap_pic(kk)==1 % this is the first hour in a new series, so check for a puff in first 60 sec
        % and process if found; if no puff then load second hour and continue processing in 2 hour
        % chunks until another gap occurs in the file series.
        [son,~,~,~,~,~,~,~] = read_sonic3D_ship_mosaic_2019_bins(wind_dfl,hed0,ddd,hh,bins,len);
        picA = read_picarro_mosaic_2019(co2_dfl,yyyy,ddd,hh);
        % if this file doesn't go to end of hour it's an orphan, so set gap on next file & contine
        if all(isnan(picA(35800:36000,2))) % the last 200 data points this hour are all NaN
            gap_pic(kk+1) = 1;
            continue
        end

        w = son(:,8);  % streamline vertical wind
        badw = length(find(isnan(w)));
        badp = length(find(isnan(picA(:,6))));
        % if too many bad points, skip this hour and set gap flag
        if badw>100 || badp>100
            gap_pic(kk+1) = 1;  % set begin gap for next file
            continue  % skip this hour if more than 10 sec of NaN in wind data
        end
        w = replace_NaN_nearest_neighbor(w);  % fix a few NaNs if present

        % load Picarro variables
        q = picA(:,10)*0.01*0.622*1000;      % q, currnet hour segment
        q = replace_NaN_nearest_neighbor(q);
        co2 = picA(:,9);                     % dry co2, ppm, current hour segment
        co2 = replace_NaN_nearest_neighbor(co2);
        ch4 = picA(:,7);
        ch4 = replace_NaN_nearest_neighbor(ch4);
        puff = picA(:,5)-32;                % valve signal, 1 = valve open
        puff = replace_NaN_nearest_neighbor(puff);

        % if before Nov 13, 2019 puff is actually controlled by DMS APIMS...
        if dnum_this_hour<dms_puff_cutoff
            dms_dfl = fullfile(dms_path,[YY,DDD],['DMS_2019_',DDD,'_',HH,'00.txt']);
            dms_rawA = read_dms_mosaic_2019(dms_dfl,ddd,hh);
            puff = dms_rawA(:,4) - 0.1;  % define puff variable for this hour
        end

        % Look for N2 puff valve signal in the first minute of this hour
        ii = find(puff(1:600)>0); % check valve signal for puff
        if isempty(ii)            % if no puff is found at start of hour
            continue              % load next hour
        else
            % define puff window
            puff_stpt = ii(1)-30;      % seg start is 3 sec before start of puff
            if puff_stpt<1             % if puff is too close to start of hour...
                puff_stpt=1;           % set to 1
            end
            % and seg end is 9 sec past end of puff
            puff_endpt = ii(end)+(ii(end)-ii(1))+90;

            % save puff begin and end points in this hour
            www(zz,52) = puff_stpt;
            www(zz,53) = puff_endpt;

            plot_path_hh = fullfile(plot_path,'lags',[YY,DDD],HH);
            mkdir(plot_path_hh)

            if plotit
                if exist('h1','var')
                    set(0,'CurrentFigure', h1);
                else
                    h1 = figure('position',[500,100,700,800],'Name','pic_puff_plot');
                end
                xdat = (picA(puff_stpt:puff_endpt,1) - picA(ii(1),1))*86400; % time stamp in seconds
                subplot(4,1,1);
                plot(xdat,puff(puff_stpt:puff_endpt),'linewidth',2);
                grid on; grid minor;
                xlim([xdat(1),xdat(end)]); ylabel('N_2 pulse');
                title(['Nitrogen Sync Pulse, DOY ',DDD,', HR, ',HH]);
                subplot(4,1,2);
                plot(xdat,co2(puff_stpt:puff_endpt),'linewidth',2);
                grid on; grid minor;
                xlim([xdat(1),xdat(end)]); ylabel('CO_2 ppm');
                subplot(4,1,3);
                plot(xdat,ch4(puff_stpt:puff_endpt),'linewidth',2);
                grid on; grid minor;
                xlim([xdat(1),xdat(end)]); ylabel('CH_4 ppm');
                subplot(4,1,4);
                plot(xdat,q(puff_stpt:puff_endpt),'linewidth',2);
                grid on; grid minor;
                xlim([xdat(1),xdat(end)]); ylabel('H_2O g/kg');
                xlabel('seconds');
                ppath = fullfile(plot_path_hh,[DDD,'_',HH,'a_PicarroN2puff.png']);
                print('-dpng',ppath);
            end

            % Select indices for a data segment in current hour, excluding first minute and
            % last 2.5 min (to exclude puff)
            % we will compute the correlation over this segment
            seg = 600:34500;
            w_seg = w(seg);

            % get rough lag time from the CO2 puff response
            jj = find(diff(ch4(1:600))<-0.02,1,'first');
            if isempty(jj)  % if no puff is found in the gas data, so go to next hour
                continue
            end
            shift = jj-ii(1);         % lag between valve signal and gas response
            www(zz,28) = shift;       % save puff_lag_pic for all window segs this hour
            iStart = jj+30;           % start tau fit 3 sec into puff event
            iEnd = iStart + 120;      % and use 12 seconds of data
            xdat = (picA(iStart:iEnd,1) - picA(iStart,1))*86400;

            % synchronize w and q if < 90 NaNs and dryer NOT installed
            if dnum_this_hour<wq_cutoff && length(find(isnan(q)))<90
                q = replace_NaN_nearest_neighbor(q);  % fix NaNs

                % compute inlet attenuation time constant, tau, for q
                % choose rising edge at end of 6 sec puff event
                q_puff = q(iStart:iEnd);
                xx = isfinite(q_puff);
                q_puff = q_puff(xx);
                xdat_q = xdat(xx);
                % first guess at fit coefficients [c0,c1,t0,tau]
                guess = [q_puff(1),q_puff(end),2,0.5];
                lower = [0.1,0.1,0.2,0.08];  % lower/upper limits on coeffs
                upper = [10,10,4,2.0];
                ft = fit(xdat_q,q_puff,fit_fos,'StartPoint',guess,'Lower',lower,'Upper',upper);
                www(zz,42) = ft.tau;  % save time constant from fit

                % plot fit
                if plotit
                    if exist('h2','var')
                        set(0,'CurrentFigure', h2);
                    else
                        h2 = figure('Name','tau_fit_q');
                    end
                    plot(ft); hold on; plot(xdat_q,q_puff,'o'); hold off;
                    grid on; grid minor;
                    xlabel('seconds'); ylabel('q g/kg');
                    title(['Water Vapor Puff, DOY ',DDD,' HR ',HH]);
                    legend(['tau = ',sprintf('%5.4f',ft.tau)],'puff data','location','northwest');
                    ppath = fullfile(plot_path_hh,[DDD,'_',HH,'b_q_puff_fit.png']);
                    print('-dpng',ppath);
                end

                % compute lag correlation using best lag estimate
                q_seg = detrend(despike2(q(seg+shift)));
                if dnum_this_hour<dms_puff_cutoff
                    ixx = 175:290; % timing was more uncertain when APIMS ran puff
                else
                    % just search for lags up to 2 to 8 sec from puff response
                    ixx = 220:280; % index range in XCorr where peak should appear
                end
                [XLag,~,~,bestCorrq] = lagcorrection2(w_seg,q_seg,200,0,ixx);
                if abs(bestCorrq)>0.02  % if there's a peak above threshold
                    q_seg = detrend(despike2(q(seg+shift+XLag)));  % use adjusted lag
                else
                    %%%%% ADJUST THIS DEFAULT LAG AS NEEDED %%%%%
                    XLag = 35; % default additional shift
                    q_seg = detrend(despike2(q(seg+shift+XLag)));  % use approximate inlet lag
                end
                www(zz,43) = shift + XLag;  % save q lag
                % compute lag correlation on shifted data and plot
                [~,XYcorr,xlags,~] = lagcorrection2(w_seg,q_seg,200,0,ixx);
                www(zz,41) = nanmean1(q_seg);   % save mean & std dev q
                www(zz,44) = nanstd1(q_seg);
                if plotit
                    if exist('h3','var')
                        set(0,'CurrentFigure', h3);
                    else
                        h3 = figure('Name','lagcorr_q');
                    end
                    if max(abs(XYcorr))>0.05
                        lims = [-max(abs(XYcorr))-0.01, max(abs(XYcorr)+0.01)];
                    else
                        lims = [-0.06,0.06];
                    end
                    % shade region where peak should appear (shifted by XLag)
                    fill([ixx(1)-200-XLag,ixx(end)-200-XLag,ixx(end)-200-XLag,ixx(1)-200-XLag],...
                        [lims(1),lims(1),lims(2),lims(2)],[0.6,0.8,1]); hold on;
                    plot([xlags(1),xlags(end)],[0,0],'k-',xlags,XYcorr,'.','markersize',10);
                    hold off; grid on; grid minor; ylim(lims); ax = gca; ax.Layer = 'top';
                    ylabel('Detrended Corr.Coef.'); xlabel('Lag points');
                    title(['Shifted q:w Lag Correlation, DOY ',DDD,' HR ',HH,';  Shift = ',...
                            sprintf('%3g',www(zz(1),43)/10),' sec']);
                    ppath = fullfile(plot_path_hh,[DDD,'_',HH,'c_qw_lagCorr.png']);
                    print('-dpng',ppath);
                end
            end

            % synchronize CO2/CH4
            % lag is approximately 1 to 4 sec after puff
            co2_seg = detrend(despike2(co2(seg+shift)));
%             co2_seg = filt_lowpass(co2_seg,10,0.8);   % low-pass filter, remove Picarro oscillation
            ch4_seg = detrend(despike2(ch4(seg+shift)));
%             ch4_seg = filt_lowpass(ch4_seg,10,0.8);   % low-pass filter, remove Picarro oscillation
            % XLag is adjustment to shift for best correlation
            % just search for lags up to 1.5 to 3 sec longer than puff response
            ixx = 215:237; % index range in XCorr where peak should appear
            [Xlagco2,XYcorrco2,xlagsco2,bestCorrco2] = lagcorrection2(w_seg,co2_seg,200,0,ixx);
            [Xlagch4,XYcorrch4,xlagsch4,bestCorrch4] = lagcorrection2(w_seg,ch4_seg,200,0,ixx);

            % if both 'peak heights' are > 0.025, we use the mean lag adjustment for XLag,
            % or whichever is greater than 0.025, or none if correlation is below 0.025
            if abs(bestCorrco2)>0.025 && abs(bestCorrch4)>0.025
                XLag = round((Xlagco2 + Xlagch4)/2);
            elseif abs(bestCorrco2)>0.025 && abs(bestCorrch4)<0.025
                XLag = Xlagco2;
            elseif abs(bestCorrco2)<0.025 && abs(bestCorrch4)>0.025
                XLag = Xlagch4;
            else
                %%%%% ADJUST THIS DEFAULT LAG AS NEEDED %%%%%
                XLag = 25;  % default lag is puff lag + 2.5 sec
            end

            co2_seg = detrend(despike2(co2(seg+shift+XLag)));  % use adjusted lag shift
%             co2_seg = filt_lowpass(co2_seg,10,0.8);   % low-pass filter, remove Picarro oscillation
            ch4_seg = detrend(despike2(ch4(seg+shift+XLag)));
%             ch4_seg = filt_lowpass(ch4_seg,10,0.8);   % low-pass filter, remove Picarro oscillation
            
            www(zz,30) = shift+XLag; www(zz,32) = shift+XLag; % save adjusted 'optimal' lag

            www(zz,31) = std(co2_seg); % save co2 stdev and corr coeff
            www(zz,38) = bestCorrco2;
            www(zz,33) = std(ch4_seg);
            www(zz,39) = bestCorrch4;

            % compute inlet attenuation time constant, tau, for co2/ch4
            % assume tau is the same for co2 & ch4
            % start/end points for puff segment same as for q
            % xdat is also the same
            ch4_puff = ch4(iStart:iEnd);
            xx = isfinite(ch4_puff);
            ch4_puff = ch4_puff(xx);
            xdat_ch4 = xdat(xx);
            % first guess at fit coefficients [c0,c1,t0,tau]
            guess = [ch4_puff(1),ch4_puff(end),2,0.5];
            lower = [0.1,0.1,0.2,0.08];  % lower/upper limits on coeffs
            upper = [3,3,4,1.2];
            ft = fit(xdat_ch4,ch4_puff,fit_fos,'StartPoint',guess,'Lower',lower,'Upper',upper);
            www(zz,29) = ft.tau;  % save time constant from fit

            % plot fit
            if plotit
                if exist('h4','var')
                    set(0,'CurrentFigure', h4);
                else
                    h4 = figure('Name','tau_fit_ch4');
                end
                plot(ft); hold on;
                plot(xdat_ch4,ch4_puff,'o'); hold off;
                hold off; grid on; grid minor;
                xlabel('seconds'); ylabel('CH_4 ppm');
                title(['CH_4 Puff, DOY ',DDD,' HR ',HH]);
                legend(['tau = ',sprintf('%5.4f',ft.tau)],'puff data','location','northwest');
                ppath = fullfile(plot_path_hh,[DDD,'_',HH,'d_ch4_puff_fit.png']);
                print('-dpng',ppath);
            end

            % plot lag correlation plots for CH4 & CO2
            if plotit
                [~,XYcorrco2,xlagsco2,~] = lagcorrection2(w_seg,co2_seg,200,0,ixx);
                [~,XYcorrch4,xlagsch4,~] = lagcorrection2(w_seg,ch4_seg,200,0,ixx);
                if exist('h5','var')
                    set(0,'CurrentFigure', h5);
                else
                    h5 = figure('position',[500,100,700,800],'Name','lagcorr_co2_ch4_final');
                end
                subplot(2,1,1)
                if max(abs(XYcorrco2))>0.05
                    lims = [-max(abs(XYcorrco2)),max(abs(XYcorrco2))];
                else
                    lims = [-0.05,0.05];
                end
                fill([15-XLag,37-XLag,37-XLag,15-XLag],[lims(1),lims(1),lims(2),lims(2)],...
                         [0.6,0.8,1]); hold on; % shade region where peak should appear (shifted by XLag)
                plot([xlagsco2(1),xlagsco2(end)],[0,0],'k-',xlagsco2,XYcorrco2,'.','markersize',10);
                hold off; grid on; grid minor; ylabel('Detrended Corr.Coef.'); xlabel('Lag points');
                ylim(lims);  ax = gca; ax.Layer = 'top';
                title(['Final CO_2 Lag Correlation, DOY ',DDD,' HR ',HH,';  Shift = ',...
                        sprintf('%3g',www(zz(1),30)/10),' sec']);
                subplot(2,1,2)
                if max(abs(XYcorrch4))>0.05
                    lims = [-max(abs(XYcorrch4)),max(abs(XYcorrch4))];
                else
                    lims = [-0.05,0.05];
                end
                fill([15-XLag,37-XLag,37-XLag,15-XLag],[lims(1),lims(1),lims(2),lims(2)],...
                         [0.6,0.8,1]); hold on; % shade region where peak should appear (shifted by XLag)
                plot([xlagsch4(1),xlagsch4(end)],[0,0],'k-',xlagsch4,XYcorrch4,'.','markersize',10);
                hold off; grid on; grid minor; ylabel('Detrended Corr.Coef.'); xlabel('Lag points');
                ylim(lims);  ax = gca; ax.Layer = 'top';
                title(['Final CH_4 Lag Correlation, DOY ',DDD,' HR ',HH,';  Shift = ',...
                        sprintf('%3g',www(zz(1),32)/10),' sec']);
                ppath = fullfile(plot_path_hh,[DDD,'_',HH,'e_co2_ch4_lagCorr.png']);
                print('-dpng',ppath);
            end
        end

        continue % we have completed processing the initial hour in a series, so go back and load
        % the next hour of data, concatenate, and continue processing in 2-hour chunks until a
        % gap is found in the file series.

    else  % gap==0 for this hour; there is existing data for prior hour; load this hour
        % & concatenate with the prior hour and look for the puff in a segment +/-30 sec from
        % the top of the hour.

        [son,~,~,~,~,~,~,~] = read_sonic3D_ship_mosaic_2019_bins(wind_dfl,hed0,ddd,hh,bins,len);
        pic = read_picarro_mosaic_2019(co2_dfl,yyyy,ddd,hh);
        % if this file doesn't go to end of hour set gap on next file but process normally
        if all(isnan(picA(35800:36000,2))) % the last 200 data points this hour are all NaN
            gap_pic(kk+1) = 1;
        end

        pic2 = [picA;pic];  % 2 hour segment for analyzing N2 pulse
        picA = pic;         % replace prior hour data with new hour

        w = son(:,8);  % streamline vertical wind
        badw = length(find(isnan(w)));
        badp = length(find(isnan(pic(:,6))));
        if badw>100 || badp>100
            gap_pic(kk+1) = 1;  % set begin gap for next file
            continue  % skip this hour if more than 10 sec of NaN in wind data
        end
        w = replace_NaN_nearest_neighbor(w);  % fix a few NaNs if present

        % load Picarro variables
        q_2 = pic2(:,10)*0.01*0.622*1000;   % h2o from pct mole fraction to g/kg
        q = pic(:,10)*0.01*0.622*1000;      % q, currnet hour segment
        q_2 = replace_NaN_nearest_neighbor(q_2);
        q = replace_NaN_nearest_neighbor(q);
        co2_2 = pic2(:,9);                  % dry co2, ppm, 2 hour segment
        co2 = pic(:,9);                     % dry co2, ppm, current hour segment
        co2_2 = replace_NaN_nearest_neighbor(co2_2);
        co2 = replace_NaN_nearest_neighbor(co2);
        ch4_2 = pic2(:,7);                  % dry ch4, ppm
        ch4 = pic(:,7);
        ch4_2 = replace_NaN_nearest_neighbor(ch4_2);
        ch4 = replace_NaN_nearest_neighbor(ch4);
        puff = pic2(:,5)-32;                % valve signal, 1 = valve open
        puff = replace_NaN_nearest_neighbor(puff);

        % if before Nov 13, 2019 puff is actually controlled by DMS APIMS...
        if dnum_this_hour<dms_puff_cutoff
            dms_dfl = fullfile(dms_path,[YY,DDD],['DMS_2019_',DDD,'_',HH,'00.txt']);
            dms_raw = read_dms_mosaic_2019(dms_dfl,ddd,hh);  % load this hour
            dms_raw2 = [dms_rawA;dms_raw];  % concat with last hour for 2-hour segment
            dms_rawA = dms_raw;             % replace prior hour data with new file
            puff = dms_raw2(:,4) - 0.1;     % define puff variable
        end

        % Look for N2 puff around the top of the hour
        ii = find(puff(34000:37000)>0);
        if isempty(ii)
            continue  % skip this hour if no N2 puff
        else          % else, define the puff segment
            ii = ii + 33999;  % ii index marks start of puff in 2 hr pic2 array
            % define puff window
            puff_stpt = ii(1)-30;                    % back up 3 sec from start of puff
            puff_endpt = ii(end)+(ii(end)-ii(1))+90; % and go to 9 sec after end of puff

            % save puff begin and end points
            % if index is prior to beginning of current hour, it's a negative index
            www(zz,52) = puff_stpt - 36000;
            www(zz,53) = puff_endpt - 36000;

            plot_path_hh = fullfile(plot_path,'lags',[YY,DDD],HH);
            mkdir(plot_path_hh)

            if plotit
                xdat = (pic2(puff_stpt:puff_endpt,1) - pic2(ii(1),1))*86400; % time stamp in seconds
                if exist('h1','var')
                    set(0,'CurrentFigure', h1);
                else
                    h1 = figure('position',[500,100,700,800],'Name','pic_puff_plot');
                end
                subplot(4,1,1);
                plot(xdat,puff(puff_stpt:puff_endpt),'linewidth',2);
                grid on; grid minor;
                xlim([xdat(1),xdat(end)]); ylabel('N_2 pulse');
                title(['Nitrogen Sync Pulse, DOY ',DDD,', HR, ',HH]);
                subplot(4,1,2);
                plot(xdat,co2_2(puff_stpt:puff_endpt),'linewidth',2);
                grid on; grid minor;
                xlim([xdat(1),xdat(end)]); ylabel('CO_2 ppm');
                subplot(4,1,3);
                plot(xdat,ch4_2(puff_stpt:puff_endpt),'linewidth',2);
                grid on; grid minor;
                xlim([xdat(1),xdat(end)]); ylabel('CH_4 ppm');
                subplot(4,1,4);
                plot(xdat,q_2(puff_stpt:puff_endpt),'linewidth',2);
                grid on; grid minor; xlabel('seconds');
                xlim([xdat(1),xdat(end)]); ylabel('H_2O g/kg');
                ppath = fullfile(plot_path_hh,[DDD,'_',HH,'a_PicarroN2puff.png']);
                print('-dpng',ppath);
            end

            % Select indices for a data segment in current hour, excluding first & last minutes
            % we will compute the correlation over this segment
            if puff_endpt>36000                     % valve puff came after top of hour
                seg = (puff_endpt-36000):34500;     % so define indices based on current hour
            else                               % valve puff came before top of hour
                seg = 600:34500;               % so just take default segment from current hour
            end
            w_seg = w(seg);

            % get rough lag time from the CH4 puff response
            jj = find(diff(ch4_2(34000:37000))<-0.02,1,'first');  % search for gas puff near top of hour
            if isempty(jj)  % if no puff is found in the gas data, go to next hour
                continue
            end
            jj = jj+33999;            % index in co2_2 for start of gas puff
            shift = jj-ii(1);         % lag between valve action and gas response, in points

            www(zz,28) = shift;       % save puff_lag_pic for all window segs this hour
            iStart = jj+30;           % start tau fit 3 sec into puff event
            iEnd = iStart + 120;      % and use 12 seconds of data
            xdat = (pic2(iStart:iEnd,1) - pic2(iStart,1))*86400;

            % synchronize w and q if < 90 NaNs and date is prior to dryer installation
            if dnum_this_hour<wq_cutoff
                % compute inlet attenuation time constant, tau, for q
                % choose rising edge at end of 6 sec puff event
                q_puff = q_2(iStart:iEnd);
                xx = isfinite(q_puff);
                q_puff = q_puff(xx);
                xdat_q = xdat(xx);
                % first guess at fit coefficients [c0,c1,t0,tau]
                guess = [q_puff(1),q_puff(end),2,0.5];
                lower = [0.1,0.1,0.2,0.08];  % lower/upper limits on coeffs
                upper = [10,10,4,2.0];
                ft = fit(xdat_q,q_puff,fit_fos,'StartPoint',guess,'Lower',lower,'Upper',upper);
                www(zz,42) = ft.tau;  % save time constant from fit

                % plot fit
                if plotit
                    if exist('h2','var')
                        set(0,'CurrentFigure', h2);
                    else
                        h2 = figure('Name','tau_fit_q');
                    end
                    plot(ft); hold on;
                    plot(xdat_q,q_puff,'o'); hold off;
                    grid on; grid minor;
                    xlabel('seconds'); ylabel('q g/kg');
                    title(['Water Vapor Puff, DOY ',DDD,' HR ',HH]);
                    legend(['tau = ',sprintf('%5.4f',ft.tau)],'puff data','location','northwest');
                    ppath = fullfile(plot_path_hh,[DDD,'_',HH,'b_q_puff_fit.png']);
                    print('-dpng',ppath);
                end

                % compute lag correlation using best lag estimate
                q_seg = detrend(despike2(q(seg+shift)));
                if dnum_this_hour<dms_puff_cutoff
                    ixx = 175:290; % timing was more uncertain when APIMS ran puff
                else
                    % just search for lags up to 2 to 8 sec from puff response
                    ixx = 220:280; % index range in XYCorr where peak should appear
                end
                [XLag,~,~,bestCorrq] = lagcorrection2(w_seg,q_seg,200,0,ixx);
                if abs(bestCorrq)>0.02  % if there's a peak above threshold
                    q_seg = detrend(despike2(q(seg+shift+XLag)));  % use adjusted lag
                else
                    %%%%% ADJUST THIS DEFAULT LAG AS NEEDED %%%%%
                    XLag = 35; % default additional shift
                    q_seg = detrend(despike2(q(seg+shift+XLag)));  % approximate inlet lag
                end
                www(zz,43) = shift+XLag;  % save q lag
                % compute lag correlation on shifted data and plot
                [~,XYcorr,xlags,~] = lagcorrection2(w_seg,q_seg,200,0,ixx);
                www(zz,41) = nanmean1(q_seg);   % save mean & std dev q
                www(zz,44) = nanstd1(q_seg);
                if plotit
                    if exist('h3','var')
                        set(0,'CurrentFigure', h3);
                    else
                        h3 = figure('Name','lagcorr_q');
                    end
                    if max(abs(XYcorr))>0.05
                        lims = [-max(abs(XYcorr))-0.01, max(abs(XYcorr)+0.01)];
                    else
                        lims = [-0.06,0.06];
                    end
                    % shade region where peak should appear (shifted by XLag)
                    fill([ixx(1)-200-XLag,ixx(end)-200-XLag,ixx(end)-200-XLag,ixx(1)-200-XLag],...
                        [lims(1),lims(1),lims(2),lims(2)],[0.6,0.8,1]); hold on;
                    plot([xlags(1),xlags(end)],[0,0],'k-',xlags,XYcorr,'.','markersize',10);
                    hold off; grid on; grid minor; ylim(lims); ax = gca; ax.Layer = 'top';
                    ylabel('Detrended Corr.Coef.'); xlabel('Lag points');
                    title(['Shifted q:w Lag Correlation, DOY ',DDD,' HR ',HH,';  Shift = ',...
                            sprintf('%3g',www(zz(1),43)/10),' sec']);
                    ppath = fullfile(plot_path_hh,[DDD,'_',HH,'c_qw_lagCorr.png']);
                    print('-dpng',ppath);
                end
            end
            
            % synchronize CO2/CH4
            % lag is approximately 1 to 4 sec after puff
            co2_seg = detrend(despike2(co2(seg+shift)));
%             co2_seg = filt_lowpass(co2_seg,10,0.8);   % low-pass filter, remove Picarro oscillation
            ch4_seg = detrend(despike2(ch4(seg+shift)));
%             ch4_seg = filt_lowpass(ch4_seg,10,0.8);   % low-pass filter, remove Picarro oscillation
            
            % XLag is adjustment to shift for best correlation
            % just search for lags up to 1.5 to 3 sec longer than puff response
            ixx = 215:237; % index range in XYCorr where peak should appear
            [Xlagco2,XYcorrco2,xlagsco2,bestCorrco2] = lagcorrection2(w_seg,co2_seg,200,0,ixx);
            [Xlagch4,XYcorrch4,xlagsch4,bestCorrch4] = lagcorrection2(w_seg,ch4_seg,200,0,ixx);

            % if both 'peak heights' are > 0.025, we use the mean lag adjustment for XLag,
            % or whichever is greater than 0.025, or none if correlation is below 0.025
            if abs(bestCorrco2)>0.025 && abs(bestCorrch4)>0.025
                XLag = round((Xlagco2 + Xlagch4)/2);
            elseif abs(bestCorrco2)>0.025 && abs(bestCorrch4)<0.025
                XLag = Xlagco2;
            elseif abs(bestCorrco2)<0.025 && abs(bestCorrch4)>0.025
                XLag = Xlagch4;
            else
                %%%%% ADJUST THIS DEFAULT LAG AS NEEDED %%%%%
                XLag = 25;  % default lag is puff lag + 2.5 sec
            end

            co2_seg = detrend(despike2(co2(seg+shift+XLag)));  % use adjusted lag shift
%             co2_seg = filt_lowpass(co2_seg,10,0.8);   % low-pass filter, remove Picarro oscillation
            ch4_seg = detrend(despike2(ch4(seg+shift+XLag)));
%             ch4_seg = filt_lowpass(ch4_seg,10,0.8);   % low-pass filter, remove Picarro oscillation
            
            www(zz,30) = shift+XLag; www(zz,32) = shift+XLag; % save adjusted 'optimal' lags

            www(zz,31) = std(co2_seg); % save co2 stdev and corr coeff
            www(zz,38) = bestCorrco2;
            www(zz,33) = std(ch4_seg);
            www(zz,39) = bestCorrch4;

            % compute inlet attenuation time constant, tau, for co2/ch4
            % assume tau is the same for co2 & ch4
            % start/end points for puff segment same as for q
            % xdat is also the same
            ch4_puff = ch4_2(iStart:iEnd);
            xx = isfinite(ch4_puff);
            ch4_puff = ch4_puff(xx);
            xdat_ch4 = xdat(xx);
            % first guess at fit coefficients [c0,c1,t0,tau]
            guess = [ch4_puff(1),ch4_puff(end),2,0.5];
            lower = [0.1,0.1,0.2,0.08];  % lower/upper limits on coeffs
            upper = [3,3,4,1.2];
            ft = fit(xdat_ch4,ch4_puff,fit_fos,'StartPoint',guess,'Lower',lower,'Upper',upper);
            www(zz,29) = ft.tau;  % save time constant from fit

            % plot fit
            if plotit
                if exist('h4','var')
                    set(0,'CurrentFigure', h4);
                else
                    h4 = figure('Name','tau_fit_ch4');
                end
                plot(ft); hold on;
                plot(xdat_ch4,ch4_puff,'o'); hold off;
                grid on; grid minor;
                xlabel('seconds'); ylabel('CH_4 ppm');
                title(['CH_4 Puff, DOY ',DDD,' HR ',HH]);
                legend(['tau = ',sprintf('%5.4f',ft.tau)],'puff data','location','northwest');
                ppath = fullfile(plot_path_hh,[DDD,'_',HH,'d_ch4_puff_fit.png']);
                print('-dpng',ppath);
            end

            % plot lag correlation plots for CH4 & CO2
            if plotit
                [~,XYcorrco2,xlagsco2,~] = lagcorrection2(w_seg,co2_seg,200,0,ixx);
                [~,XYcorrch4,xlagsch4,~] = lagcorrection2(w_seg,ch4_seg,200,0,ixx);
                if exist('h5','var')
                    set(0,'CurrentFigure', h5);
                else
                    h5 = figure('position',[500,100,700,800],'Name','lagcorr_co2_ch4_final');
                end
                subplot(2,1,1)
                if max(abs(XYcorrco2))>0.05
                    lims = [-max(abs(XYcorrco2)),max(abs(XYcorrco2))];
                else
                    lims = [-0.05,0.05];
                end
                fill([15-XLag,37-XLag,37-XLag,15-XLag],[lims(1),lims(1),lims(2),lims(2)],...
                         [0.6,0.8,1]); hold on; % shade region where peak should appear (shifted by XLag)
                plot([xlagsco2(1),xlagsco2(end)],[0,0],'k-',xlagsco2,XYcorrco2,'.','markersize',10);
                hold off; grid on; grid minor; ylabel('Detrended Corr.Coef.'); xlabel('Lag points');
                ylim(lims);  ax = gca; ax.Layer = 'top';
                title(['Final CO_2 Lag Correlation, DOY ',DDD,' HR ',HH,';  Shift = ',...
                        sprintf('%3g',www(zz(1),30)/10),' sec']);
                subplot(2,1,2)
                if max(abs(XYcorrch4))>0.05
                    lims = [-max(abs(XYcorrch4)),max(abs(XYcorrch4))];
                else
                    lims = [-0.05,0.05];
                end
                fill([15-XLag,37-XLag,37-XLag,15-XLag],[lims(1),lims(1),lims(2),lims(2)],...
                         [0.6,0.8,1]); hold on; % shade region where peak should appear (shifted by XLag)
                plot([xlagsch4(1),xlagsch4(end)],[0,0],'k-',xlagsch4,XYcorrch4,'.','markersize',10);
                hold off; grid on; grid minor; ylabel('Detrended Corr.Coef.'); xlabel('Lag points');
                ylim(lims);  ax = gca; ax.Layer = 'top';
                title(['Final CH_4 Lag Correlation, DOY ',DDD,' HR ',HH,';  Shift = ',...
                        sprintf('%3g',www(zz(1),32)/10),' sec']);
                ppath = fullfile(plot_path_hh,[DDD,'_',HH,'e_co2_ch4_lagCorr.png']);
                print('-dpng',ppath);
            end
        end  % finished processing this 2-hour chunk, so load next hour and repeat until
        % a gap is detected in the file series
    end
end     % end of picarro loop
close all
clear ch4 ch4_2 ch4_seg co2 co2_seg co2_2 q q_2 q_seg pic pic2 picA h1 h2 h3 h4 h5

%% ------------------------------ begin loop to analyze fast ozone puffs -----------------------------
% this will be just like the Picarro loop above
ddd0 = 0; % set current day flag
clear h6 h7 h8
www(:,34:36) = NaN;

% gap_fo3(157) = 1; kk = 157:168
for kk = 1:length(flist_fo3)

    % get file names for this hour - if any file missing, set gap and go to next hour
    wind_dfl = cell2mat(flist_son0(kk));
    if isempty(wind_dfl)
        gap_fo3(kk+1) = 1;
        continue
    end
    % picarro is controlling puff, so we need these files to plot puff signal...
    co2_dfl = cell2mat(flist_pic0(kk));
    if isempty(co2_dfl)
        gap_fo3(kk+1) = 1;
        continue
    end
    fo3_dfl = cell2mat(flist_fo3(kk));
    if isempty(fo3_dfl)
        gap_fo3(kk+1) = 1;
        continue
    end

    % set up the various date variables and strings from wind file name
    temp = strsplit(wind_dfl,'/');  fname = cell2mat(temp(8));
    YY = fname(5:6); YYYY = ['20',YY]; DDD = fname(7:9); HH = fname(10:11);
    yyyy = str2double(YYYY); ddd = str2double(DDD); hh = str2double(HH);
    [mm,dd] = yd2md(yyyy,ddd); MM = sprintf('%02i',mm); DD = sprintf('%02i',dd);
    dnum_this_hour = datenum(yyyy,mm,dd,hh,0,0);

    % find leg number and www indices for this hour
    zz = find(www(:,3)==yyyy & www(:,4)==mm & www(:,5)==dd & www(:,6)==hh);
    if isempty(zz) % no flux windows for this hour
        gap_fo3(kk+1) = 1; % set next hour gap flag
        continue           % skip to next hour
    end

    if any(www(zz,12)<25)  % ozone inlet flow shut down, so skip this hour
        gap_fo3(kk+1) = 1; % set next hour gap flag
        continue           % skip to next hour
    end

    LEG = sprintf('%01i',www(zz(1),1));

    if isnan(www(zz(1),28)) % check for a co2/ch4 puff lag in this hour
        gap_fo3(kk+1) = 1;  % no puff, so set next hour gap flag
        continue            % and skip to next hour
    else                    % use puff indices from CO2/CH4 Picarro analysis
        puff_stpt = www(zz(1),52);
        puff_endpt = www(zz(1),53);
    end

    % load ship met/nav data
    if ddd>ddd0  % if new day
        % load ship heading for this day
        % this is required to read the raw sonic data
        path_dsp = fullfile(root_data_path,'ship_data','parsed',['dsp0',YY,MM,DD,'_raw.txt']);
        dsp = read_dsp_mosaic_2019(path_dsp,ddd);
        hed0 = [dsp(:,1),dsp(:,7)]; % 1-min ship timestamp and heading for this day

        ddd0 = ddd;  % remember current day
    end

    % we will process data in 2-hour sequential segments, so we decide here
    % if we will load a new hour in the current series or start a new series following a gap
    if gap_fo3(kk)==1 % this is the first hour in a new series, so check for a puff in first 60 sec
        % and process if found; if no puff then load second hour and continue processing in 2 hour
        % chunks until another gap occurs in the file series.
        [son,~,~,~,~,~,~,~] = read_sonic3D_ship_mosaic_2019_bins(wind_dfl,hed0,ddd,hh,bins,len);
        picA = read_picarro_mosaic_2019(co2_dfl,yyyy,ddd,hh);  % need this file for the puff signal
        fo3A = read_fo3_cnts_mosaic_2019(fo3_dfl,ddd,hh);
        % recompute 10Hz O3 from raw count rate and cal coefs for each leg
        % this is necessary because concentration recorded in the file is LP filtered
        temp = fo3A(:,2)*fo3_coef(www(kk,1),1) + fo3_coef(www(kk,1),2);
        fo3A(:,3) = temp;

        % if this file doesn't go to end of hour it's an orphan, so set gap on next file & contine
        if all(isnan(fo3A(35800:36000,3))) % the last 200 data points this hour are all NaN
            gap_fo3(kk+1) = 1;
            continue
        end

        % if puff occurred prior to start of hour, don't analyze this hour for the puff,
        % continue to load next hour and process a 2 hour segment
        if puff_stpt<1
            continue
        end

        w = son(:,8);  % streamline vertical wind
        badw = length(find(isnan(w)));
        bado3 = length(find(isnan(fo3A(:,3))));
        % if too many bad points, skip this hour and set gap flag
        if badw>100 || bado3>100
            gap_fo3(kk+1) = 1;  % set begin gap for next file
            continue  % skip this hour if more than 10 sec of NaN in wind data
        end
        w = replace_NaN_nearest_neighbor(w);  % fix a few NaNs if present

        % load ozone and puff variables
        o3 = fo3A(:,3);
        o3 = replace_NaN_nearest_neighbor(o3);
        puff = picA(:,5)-32;                  % valve signal, 1 = valve open
        puff = replace_NaN_nearest_neighbor(puff);

        % if before Nov 13, 2019 puff is actually controlled by DMS APIMS...
        if dnum_this_hour<dms_puff_cutoff
            dms_dfl = fullfile(dms_path,[YY,DDD],['DMS_2019_',DDD,'_',HH,'00.txt']);
            dms_rawA = read_dms_mosaic_2019(dms_dfl,ddd,hh);
            puff = dms_rawA(:,4) - 0.1;  % define puff variable for this hour
            puff = replace_NaN_nearest_neighbor(puff);
        end

        jjj = find(puff(puff_stpt:puff_endpt)>0.2,1,'first');  % start of valve pulse

        % Select indices for a data segment in current hour, excluding first & last minutes
        % we will compute the correlation for this hour over this segment
        seg = 600:34500;
        w_seg = w(seg);

        % analyze puff, using puff_stpt/puff_endpt indices from CO2/CH4 run
        % plot ozone puff
        if plotit
            plot_path_hh = fullfile(plot_path,'lags',[YY,DDD],HH);
            if exist('h6','var')
                set(0,'CurrentFigure', h6);
            else
                h6 = figure('Name','O3_puff_plot');
            end
            xdat = (fo3A(puff_stpt:puff_endpt,1) - fo3A(puff_stpt+jjj,1))*86400; % time stamp in seconds
            subplot(2,1,1);
            plot(xdat,puff(puff_stpt:puff_endpt),'linewidth',2); grid on; grid minor;
            xlim([xdat(1),xdat(end)]); ylabel('N_2 pulse');
            title(['Ozone Sync Pulse, DOY ',DDD,', HR, ',HH]);
            subplot(2,1,2);
            plot(xdat,o3(puff_stpt:puff_endpt),'linewidth',2); grid on; grid minor;
            xlim([xdat(1),xdat(end)]); ylabel('O_3 ppb');
            xlabel('seconds');
            ppath = fullfile(plot_path_hh,[DDD,'_',HH,'f_OzoneN2puff.png']);
            print('-dpng',ppath);
        end

        % find puff lag time and synchronize o3_seg
        o3_filt = smooth(o3,10/36000,'moving'); % smoothed ozone

        % need to use different sensitivities to detect the puff in the ozone signal,
        % depending on the mean ozone concentration, and wait at least 4 seconds after
        % valve injection before searching for o3 puff
        if dnum_this_hour<dms_puff_cutoff
            delay = 0;
        else
            delay = 40;
        end
        if mean(o3_filt)>15
            jj = find(abs(diff(o3_filt(puff_stpt+jjj+delay:puff_endpt)))>0.50,1,'first');
            jj = jj+jjj+delay;
        elseif mean(o3_filt)>10
            jj = find(abs(diff(o3_filt(puff_stpt+jjj+delay:puff_endpt)))>0.10,1,'first');
            jj = jj+jjj+delay;
        elseif mean(o3_filt)>5
            jj = find(abs(diff(o3_filt(puff_stpt+jjj+delay:puff_endpt)))>0.06,1,'first');
            jj = jj+jjj+delay;
        elseif mean(o3_filt)>2
            jj = find(abs(diff(o3_filt(puff_stpt+jjj+jjj+delay:puff_endpt)))>0.03,1,'first');
            jj = jj+jjj+delay;
        else
            jj = [];   % ozone concentration is too low to detect puff
        end
        if isempty(jj)                            % can't find an O3 puff in this hour
            continue
        else                                      % start of ozone puff is found
            shift = jj-jjj;                       % lag points between valve opening and ozone response
            www(zz,34) = shift;                   % save lag time to start of pulse
            o3_seg = o3(seg+shift);               % shift O3 by observed puff lag
            [o3_seg,~] = despike2(o3_seg);        % despike
            www(zz,37) = std(o3_seg);             % save O3 std dev
            % just search for lags up to 1.8 to 4.6 sec longer than puff response
            ixx = 218:246; % index range in XCorr where peak should appear
            [XLag,~,~,bestCorro3] = lagcorrection2(w_seg,o3_seg,200,0,ixx);
            www(zz,40) = bestCorro3;              % save default 'optimal' o3 corr coeff
            if abs(bestCorro3)>0.015              % if bestCorro3 is above threshold
                o3_seg = detrend(despike2(o3(seg+shift+XLag)));  % use adjusted lag
                www(zz,36) = shift + XLag;        % save 'optimal' o3 lag
            else                                  % use default 'optimal' o3 lag = puff lag + 32
                %%%%%%%%%  readjust this default value  %%%%%%%%%%
                XLag = 32;
                o3_seg = detrend(despike2(o3(seg+shift+XLag)));
                www(zz,36) = shift + XLag;
            end

            % compute inlet attenuation time constant, tau, for ozone
            % start/end points for puff segment same as for q
            % xdat is also the same
            iStart = jj+40;      % start tau fit 4 sec into puff event
            iEnd = iStart + 100; % and use 10 seconds of data
            o3_puff = o3(iStart:iEnd);
            xdat = (fo3A(iStart:iEnd,1) - fo3A(iStart,1))*86400;
            xx = isfinite(o3_puff);
            o3_puff = o3_puff(xx);
            xdat_o3 = xdat(xx);
            % first guess at fit coefficients [c0,c1,t0,tau]
            guess = [o3_puff(1),o3_puff(end),2,0.5];
            lower = [0.1,0.1,0.2,0.08];  % lower/upper limits on coeffs
            upper = [50,50,4,2.0];
            ft = fit(xdat_o3,o3_puff,fit_fos,'StartPoint',guess,'Lower',lower,'Upper',upper);
            www(zz,35) = ft.tau;  % save time constant from fit

            % plot fit
            if plotit
                if exist('h7','var')
                    set(0,'CurrentFigure', h7);
                else
                    h7 = figure('Name','O3_fit_tau');
                end
                plot(ft);  grid on; grid minor; hold on;
                plot(xdat_o3,o3_puff,'o'); hold off;
                xlabel('seconds'); ylabel('ozone, ppb');
                title(['Ozone Puff, DOY ',DDD,' HR ',HH]);
                legend(['tau = ',sprintf('%5.4f',ft.tau)],'puff data','location','northwest');
                ppath = fullfile(plot_path_hh,[DDD,'_',HH,'g_o3_puff_fit.png']);
                print('-dpng',ppath);
            end
        end

        if plotit
            % plot o3:w lag correlation
            [~,XYcorro3,xlagso3,~] = lagcorrection2(w_seg,o3_seg,200,0,ixx);
            if exist('h8','var')
                set(0,'CurrentFigure', h8);
            else
                h8 = figure('Name','O3_lagCorr_final');
            end
            if max(abs(XYcorro3))>0.05
                lims = [-max(abs(XYcorro3)),max(abs(XYcorro3))];
            else
                lims = [-0.05,0.05];
            end
            fill([18-XLag,46-XLag,46-XLag,18-XLag],[lims(1),lims(1),lims(2),lims(2)],...
                         [0.6,0.8,1]); hold on; % shade region where peak should appear (shifted by XLag)
            plot(xlagso3,XYcorro3,'r.','markersize',10); hold off; grid on; grid minor;
            ylabel('Detrended Corr.Coef.'); xlabel('Lag points');
            ylim(lims);  ax = gca; ax.Layer = 'top';
            title(['Final O_3 Lag Correlation, DOY ',DDD,' HR ',HH,';  Shift = ',...
                sprintf('%3g',www(zz(1),36)/10),' sec']);
            ppath = fullfile(plot_path_hh,[DDD,'_',HH,'h_o3w_lagCorr.png']);
            print('-dpng',ppath);
        end


    else  % gap==0 for this hour; there is existing data for prior hour; load this hour
        % & concatenate with the prior hour and look for the puff in a segment +/-30 sec from
        % the top of the hour.
        [son,~,~,~,~,~,~,~] = read_sonic3D_ship_mosaic_2019_bins(wind_dfl,hed0,ddd,hh,bins,len);
        pic = read_picarro_mosaic_2019(co2_dfl,yyyy,ddd,hh);  % need this file for the puff signal
        pic2 = [picA;pic];  % concat 2 hours of Picarro data
        picA = pic;         % replace prior hour data with new hour

        fo3 = read_fo3_cnts_mosaic_2019(fo3_dfl,ddd,hh);
        % recompute 10Hz O3 from raw count rate and cal coefs for each leg
        % this is necessary because concentration recorded in the file is LP filtered
        temp = fo3(:,2)*fo3_coef(www(kk,1),1) + fo3_coef(www(kk,1),2);
        fo3(:,3) = temp;

        % if this file doesn't go to end of hour it's an orphan, so set gap on next file & contine
        if all(isnan(fo3(35800:36000,3))) % the last 200 data points this hour are all NaN
            gap_fo3(kk+1) = 1;
            continue
        end

        fo3_2 = [fo3A;fo3];  % concat 2 hours of ozone data
        fo3A = fo3;          % save this hour's data for next cycle

        w = son(:,8);  % streamline vertical wind
        badw = length(find(isnan(w)));
        bado3 = length(find(isnan(fo3(:,3))));
        if badw>100 || bado3>100 % skip this hour if more than 10 sec of NaN data
            gap_fo3(kk+1) = 1;  % set begin gap for next file
            continue
        end
        w = replace_NaN_nearest_neighbor(w);  % fix a few NaNs if present

        % load ozone variables
        o3_2 = fo3_2(:,3);      % 2 hour segment
        o3 = fo3(:,3);          % current hour only
        o3_2 = replace_NaN_nearest_neighbor(o3_2);
        o3 = replace_NaN_nearest_neighbor(o3);
        puff = pic2(:,5)-32;    % 2 hour segment of valve signal, 1 = valve open
        puff = replace_NaN_nearest_neighbor(puff);

        % if before Nov 13, 2019 puff is actually controlled by DMS APIMS...
        if dnum_this_hour<dms_puff_cutoff
            dms_dfl = fullfile(dms_path,[YY,DDD],['DMS_2019_',DDD,'_',HH,'00.txt']);
            dms_raw = read_dms_mosaic_2019(dms_dfl,ddd,hh);  % load this hour
            dms_raw2 = [dms_rawA;dms_raw];  % concat with last hour for 2-hour segment
            dms_rawA = dms_raw;             % replace prior hour data with this hour
            puff = dms_raw2(:,4) - 0.1;     % define puff variable
            puff = replace_NaN_nearest_neighbor(puff);
        end

        % analyze puff, using puff_stpt/puff_endpt indices from CO2/CH4 run
        puff_stpt = puff_stpt + 36000;    % adjust indices for 2hr time series
        puff_endpt = puff_endpt + 36000;
        jjj = find(puff(puff_stpt:puff_endpt)>0.2,1,'first');  % start of valve pulse

        % Select indices for a data segment in current hour, excluding first & last minutes
        % we will compute the correlation for this hour over this segment
        seg = 600:34500;
        w_seg = w(seg);

        % plot ozone puff
        if plotit
            plot_path_hh = fullfile(plot_path,'lags',[YY,DDD],HH);
            if exist('h6','var')
                set(0,'CurrentFigure', h6);
            else
                h6 = figure('Name','O3_puff_plot');
            end
            xdat = (fo3_2(puff_stpt:puff_endpt,1) - fo3_2(puff_stpt+jjj,1))*86400; % time stamp in seconds
            subplot(2,1,1);
            plot(xdat,puff(puff_stpt:puff_endpt),'linewidth',2); grid on; grid minor;
            xlim([xdat(1),xdat(end)]); ylabel('N_2 pulse');
            title(['Ozone Sync Pulse, DOY ',DDD,', HR, ',HH]);
            subplot(2,1,2);
            plot(xdat,o3_2(puff_stpt:puff_endpt),'linewidth',2); grid on; grid minor;
            xlim([xdat(1),xdat(end)]); ylabel('O_3 ppb');
            xlabel('seconds');
            ppath = fullfile(plot_path_hh,[DDD,'_',HH,'f_OzoneN2puff.png']);
            print('-dpng',ppath);
        end

        % find puff lag time and synchronize o3_seg
        o3_filt = smooth(o3_2,10/36000,'moving'); % smoothed ozone
        jjj = find(puff(puff_stpt:puff_endpt)>0.2,1,'first');  % start of valve pulse
        % need to use different sensitivities to detect the puff in the ozone signal,
        % depending on the mean ozone concentration, and wait at least 4 seconds after
        % valve injection before searching for o3 puff
        if dnum_this_hour<dms_puff_cutoff
            delay = 0;
        else
            delay = 40;
        end
        if mean(o3_filt)>15
            jj = find(abs(diff(o3_filt(puff_stpt+jjj+delay:puff_endpt)))>0.50,1,'first');
            jj = jj+jjj+delay;
        elseif mean(o3_filt)>10
            jj = find(abs(diff(o3_filt(puff_stpt+jjj+delay:puff_endpt)))>0.10,1,'first');
            jj = jj+jjj+delay;
        elseif mean(o3_filt)>5
            jj = find(abs(diff(o3_filt(puff_stpt+jjj+delay:puff_endpt)))>0.06,1,'first');
            jj = jj+jjj+delay;
        elseif mean(o3_filt)>2
            jj = find(abs(diff(o3_filt(puff_stpt+jjj+delay:puff_endpt)))>0.03,1,'first');
            jj = jj+jjj+delay;
        else
            jj = [];   % ozone concentration is too low to detect puff
        end
        if isempty(jj)                            % can't find an O3 puff in this hour
            continue
        else                                      % start of ozone puff is found
            shift = jj-jjj;                       % lag points between valve opening and ozone response
            www(zz,34) = shift;                   % save lag time to start of pulse
            o3_seg = o3(seg+shift);               % shift O3 by observed puff lag
            [o3_seg,~] = despike2(o3_seg);        % despike
            www(zz,37) = std(o3_seg);             % save O3 std dev
            % just search for lags up to 1.8 to 4.6 sec longer than puff response
            ixx = 218:246; % index range in XCorr where peak should appear
            [XLag,~,~,bestCorro3] = lagcorrection2(w_seg,o3_seg,200,0,ixx);
            www(zz,40) = bestCorro3;              % save default 'optimal' o3 corr coeff
            if abs(bestCorro3)>0.015               % if bestCorro3 is above threshold
                o3_seg = detrend(despike2(o3(seg+shift+XLag)));  % use adjusted lag
                www(zz,36) = shift+XLag;          % save 'optimal' o3 lag
            else                                  % use default 'optimal' o3 lag = puff lag + 32
                %%%%%%%%%  readjust this default value  %%%%%%%%%%
                XLag = 32;
                o3_seg = detrend(despike2(o3(seg+shift+XLag)));
                www(zz,36) = shift + XLag;
            end

            % compute inlet attenuation time constant, tau, for ozone
            % start/end points for puff segment same as for q
            % xdat is also the same
            iStart = puff_stpt+jj+40;     % start tau fit 4 sec into puff event
            iEnd = iStart + 100;          % and use 10 seconds of data
            o3_puff = o3_2(iStart:iEnd);
            xdat = (fo3_2(iStart:iEnd,1) - fo3_2(iStart,1))*86400;
            xx = isfinite(o3_puff);
            o3_puff = o3_puff(xx);
            xdat_o3 = xdat(xx);
            % first guess at fit coefficients [c0,c1,t0,tau]
            guess = [o3_puff(1),o3_puff(end),2,0.5];
            lower = [0.1,0.1,0.2,0.08];  % lower/upper limits on coeffs
            upper = [50,50,4,2.0];
            ft = fit(xdat_o3,o3_puff,fit_fos,'StartPoint',guess,'Lower',lower,'Upper',upper);
            www(zz,35) = ft.tau;  % save time constant from fit

            % plot fit
            if plotit
                if exist('h7','var')
                    set(0,'CurrentFigure', h7);
                else
                    h7 = figure('Name','O3_fit_tau');
                end
                plot(ft);  grid on; grid minor; hold on;
                plot(xdat_o3,o3_puff,'o'); hold off;
                xlabel('seconds'); ylabel('ozone, ppb');
                title(['Ozone Puff, DOY ',DDD,' HR ',HH]);
                legend(['tau = ',sprintf('%5.4f',ft.tau)],'puff data','location','northwest');
                ppath = fullfile(plot_path_hh,[DDD,'_',HH,'g_o3_puff_fit.png']);
                print('-dpng',ppath);
            end
        end

        if plotit
            % plot o3:w lag correlation
            [~,XYcorro3,xlagso3,~] = lagcorrection2(w_seg,o3_seg,200,0,[ixx]);
            if exist('h8','var')
                set(0,'CurrentFigure', h8);
            else
                h8 = figure('Name','O3_lagCorr_final');
            end
            if max(abs(XYcorro3))>0.05
                lims = [-max(abs(XYcorro3)),max(abs(XYcorro3))];
            else
                lims = [-0.05,0.05];
            end
            fill([18-XLag,46-XLag,46-XLag,18-XLag],[lims(1),lims(1),lims(2),lims(2)],...
                         [0.6,0.8,1]); hold on; % shade region where peak should appear (shifted by XLag)
            plot(xlagso3,XYcorro3,'r.','markersize',10); hold off; grid on; grid minor;
            ylabel('Detrended Corr.Coef.'); xlabel('Lag points');
            ylim(lims);  ax = gca; ax.Layer = 'top';
            title(['Final O_3 Lag Correlation, DOY ',DDD,' HR ',HH,';  Shift = ',...
                sprintf('%3g',www(zz(1),36)/10),' sec']);
            ppath = fullfile(plot_path_hh,[DDD,'_',HH,'h_o3w_lagCorr.png']);
            print('-dpng',ppath);
        end
    end
end

close all
clear o3 o3_2 o3_seg o3_filt h6 h7 h8 fo3 fo32 fo3A pic pic2 picA son temp w w_seg
clear dms_raw dms_raw2 dms_rawA

%% save the flx data table with additional columns
%{
    Flux table variables for overlapping 8192 point flux windows (7 per hour)
    1        leg           cruise Leg number (1-5)
    2        datenum       matlab datenum @ start of flux window
    3        yyyy          year
    4        mm            month
    5        dd            day
    6        hh            hour
    7        seg_begin     flux window start data point in 1 hr 10Hz array (36000 points/hr)
    8        seg_end       flux window end data point in 1 hr 10Hz array (36000 points/hr)
    9        rwdir         mean rel wind direction over flux window
    10       rwspd         mean rel wind speed over flux window
    11       co2_flo       Picarro main inlet line flow rate (LPM)
    12       o3_flo        fast ozone main inlet line flow rate (LPM)
    13       Pcav          Picarro cavity pressure, torr
    14       Tcav          Picarro 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
    28       puff_lag_co2  lag between valve and start of Picarro gas puff, points
    29       co2_ch4_tau   inlet attn time const for co2 & ch4
    30       co2_lag       'optimal' lag for co2
    31       co2_std       std dev co2 (despiked), ppm
    32       ch4_lag       'optimal' lag for ch4
    33       ch4_std       std dev ch4 (despiked), ppm
    34       o3_puff_lag   number of data points to start of  o3 pulse
    35       o3_tau        inlet attn time const for o3
    36       o3_lag        'optimal' lag for o3
    37       o3_std        std dev o3 (despiked), ppb
    38       max_Corr_co2  maximum correlation coef, co2
    39       max_Corr_ch4  maximum correlation coef, ch4
    40       max_Corr_o3   maximum correlation coef, o3
    41       q_pic         specific humidity, Picarro, g/kg
    42       q_pic_tau     inlet attn time const for Picarro q
    43       q_pic_lag     'optimal' lag for Picarro q
    44       q_pic_std     std dev Picarro q (despiked), g/kg
    45       q_lic         specific humidity, Licor 7500, g/kg
    46       q_lic_std     std dev Licor q, g/kg
    47       q_lic_qc      Licor QC parameter (0-100) no units
    48       son_bad       number of NaN values in each flux window
    49       son_missing   number of missing points in each flux window
    50       son_azm       azimuth angle in each flux window (rel wdir)
    51       son_tilt      tilt angle in each flux window
    52       puff_stpt     array index for start of puff segment
    53       puff_endpt    array index for end of puff segment
%}

% first kill known bad ozone segments - no puff due to icing during Leg 5
ii = www(:,2)>datenum(2019,12,31)+262 & www(:,2)<datenum(2019,12,31)+263.3;
www(ii,[35,36,40]) = NaN;
ii = www(:,2)>datenum(2019,12,31)+268 & www(:,2)<datenum(2019,12,31)+271;
www(ii,[35,36,40]) = NaN;

clear flx
flx = array2table(www,'VariableNames',vars);
sfl = fullfile(save_path,'wins_gases_lags_8192_all.mat');
save(sfl,'flx');

%% plot lags and std deviations

dt = datetime(flx.datenum,'convertfrom','datenum');

figure('position',[500,100,1500,1100]);
subplot(5,1,1);
    yyaxis left; plot(dt,flx.q_pic_lag/10,'.','markersize',15); grid on; grid minor;
    ylabel('q lag seconds'); xlabel('day-of-year'); xlim([dt(1),dt(end)]);
    ax=gca; set(ax,'fontsize',14);
    yyaxis right; plot(dt,flx.q_pic_std,'.','markersize',15);
    ylabel('q std dev, g/kg'); ylim([0,3e-2]);
    ax=gca; set(ax,'fontsize',14);
    if any(isfinite(flx.q_pic_std))
        med_sigma = nanmedian1(flx.q_pic_std);
        ylim([med_sigma*0.25,med_sigma*1.75]);
    end
    title('Inlet Lag Times & Time Constants, Conc Std Dev, Ship','fontweight','bold','fontsize',18);
subplot(5,1,2);
    yyaxis left; plot(dt,flx.co2_lag/10,'.','markersize',15); grid on; grid minor;
    ylabel('co2 lag seconds'); xlabel('day-of-year'); xlim([dt(1),dt(end)]);
    ax=gca; set(ax,'fontsize',14);
    yyaxis right; plot(dt,flx.co2_std,'.','markersize',15);
    ylabel('co2 std dev, ppm'); ax=gca; set(ax,'fontsize',14);
    if any(isfinite(flx.co2_std))
        med_sigma = nanmedian1(flx.co2_std);
        ylim([med_sigma*0.25,med_sigma*1.75]);
    end
subplot(5,1,3);
    yyaxis left; plot(dt,flx.ch4_lag/10,'.','markersize',15); grid on; grid minor;
    ylabel('ch4 lag seconds'); xlabel('day-of-year'); xlim([dt(1),dt(end)]);
    ax=gca; set(ax,'fontsize',14);
    yyaxis right; plot(dt,flx.ch4_std,'.','markersize',15); ax=gca; set(ax,'fontsize',14);
    ylabel('ch4 std dev, ppm');
    if any(isfinite(flx.ch4_std))
        med_sigma = nanmedian1(flx.ch4_std);
        ylim([med_sigma*0.25,med_sigma*1.75]);
    end
subplot(5,1,4);
    yyaxis left; plot(dt,flx.o3_lag/10,'.','markersize',15); grid on; grid minor;
    ylabel('o3 lag seconds'); xlabel('day-of-year'); xlim([dt(1),dt(end)]);
    ax=gca; set(ax,'fontsize',14);
    yyaxis right; plot(dt,flx.o3_std,'.','markersize',15);
    ylabel('o3 std dev, ppb'); ylim([0,10]); ax=gca; set(ax,'fontsize',14);
    if any(isfinite(flx.o3_std))
        med_sigma = nanmedian1(flx.o3_std);
        ylim([med_sigma*0.25,med_sigma*1.75]);
    end
subplot(5,1,5);
    plot(dt,flx.co2_ch4_tau,'.',dt,flx.o3_tau,'.',dt,flx.q_pic_tau,'.','markersize',15);
    ax=gca; set(ax,'fontsize',14);
    ylabel('Tau, sec'); xlabel('day-of-year'); xlim([dt(1),dt(end)]); grid on; grid minor;
    legend('co2/ch4','o3','q','location','northeast');
if plotit
    ppath = fullfile(plot_path,'ship_inlet_lags_all_gases_8192.png');
    print('-dpng',ppath);
    ppath = fullfile(plot_path,'ship_inlet_lags_all_gases_8192.fig');
    savefig(ppath);
end

med_lag = nanmedian1(flx.puff_lag_co2/10);
figure; plot(dt,flx.puff_lag_co2/10,'.','markersize',10); grid;
ylabel('Puff Lag, sec'); xlabel('day-of-year'); ax=gca; set(ax,'fontsize',14);
title(['Raw N_2 Lag Time, Ship CO_2/CH_4, median = ',sprintf('%2.1f',med_lag),' sec']);
if plotit
    ppath = fullfile(plot_path,'ship_picarro_raw_N2_lag_8192.png');
    print('-dpng',ppath);
    ppath = fullfile(plot_path,'ship_picarro_raw_N2_lag_8192.fig');
    savefig(ppath);
end

figure; plot(dt,flx.max_corr_co2,'.',dt,flx.max_corr_ch4,'.',[dt(1),dt(end)],[0,0],'k-',...
             'markersize',10); grid;
ylabel('Correlation Coeff'); xlabel('day-of-year'); ax=gca; set(ax,'fontsize',14);
title('CO_2/CH_4 : W Correlation, Ship'); ylim([-0.2,0.2]);
legend('co2','ch4','location','north');
if plotit
    ppath = fullfile(plot_path,'ship_picarro_corr_coeffs_8192.png');
    print('-dpng',ppath);
    ppath = fullfile(plot_path,'ship_picarro_corr_coeffs_8192.fig');
    savefig(ppath);
end
