clear all
close all
load('ozonebox_2020_02_15.mat')

% Just in case I forget to delete the old flagtimes
if exist('flagtimes_old','var')
    flagtimes = flagtimes_old;
    clear flagtimes_old
end

% Create timestamp for o3 data
out.timestamp = datenum(out.DateTimeUTC,'yyyy/mm/dd HH:MM:SS');

% Change timestamp to elapsed minutes
flagtimes_old = flagtimes;
first_flagtime = flagtimes{1}(1);
out.timestamp = (out.timestamp-first_flagtime).*1440;
for f = 1:length(flagtimes)
    for g = 1:length(flagtimes{f})
        flagtimes{f}(g) = (flagtimes{f}(g) - first_flagtime).*1440;
    end
end

% Compute the difference in o3 values
out.diff = out.o3_2adj - out.o3_1;

% Set plot variables
plot_vars = {'o3_1','o3_2adj';...
             'diffo3','diff'};
         
% Set plot labels
plot_labs = {'Ozone Mixing Ratio / ppbv' 'Differential Ozone Mixing Ratio / ppbv'};

% Set plot legends
plot_legs = {'Sonde1','Sonde2 (adjusted)';...
             'Differential TEI','Sonde Difference'};
fill_legs = {'Direct Comparison','Chamber Blank','Surface Measurement'};

% Set plot parameters
lw = 1.5;
l_cols = [0.000 0.447 0.741;...
          0.850 0.325 0.098];

% Set fill parameters
f_cols = [0.929 0.694 0.125;...
          0.494 0.184 0.556;...
          0.466 0.674 0.188];
alpha = 0.6;
edge_col = 'None';
f_marker = 'None';

% Plot data
figure
for s = 1:2
    % Create subplot
    subplot(2,1,s); hold on; box on;
    
    % Plot the data once and get y limits
    for p = 1:2
        plt(p) = plot(out.timestamp,out{:,strcmp(out.Properties.VariableNames,plot_vars{s,p})},...
                      'Color',l_cols(p,:),'LineWidth',lw);
    end
    yl = get(gca,'YLim');
    
    % Fill times
    for f = 1:3
        for g = 1:2:length(flagtimes{f})-1
            plt(f+2) = fill([flagtimes{f}(g) flagtimes{f}(g) flagtimes{f}(g+1) flagtimes{f}(g+1)],[yl,fliplr(yl)],f_cols(f,:),...
                             'EdgeColor',edge_col,'FaceAlpha',alpha,'Marker',f_marker);
        end
    end
    
    % Replot the data
    for p = 1:2
         plot(out.timestamp,out{:,strcmp(out.Properties.VariableNames,plot_vars{s,p})},...
              'Color',l_cols(p,:),'LineWidth',lw);
    end
    
    % Zero line for differential measurements
    if s == 2
        plot(get(gca,'XLim'),[0 0],'k--','LineWidth',lw)
    end
    
    % Set some plot parameters
    set(gca,'XLim',[flagtimes{1}(1) flagtimes{1}(end)])
    set(gca,'YLim',yl)
    xlabel('Elapsed Minutes')
    ylabel(plot_labs{s})
    legend(plt,[plot_legs(s,:),fill_legs])
end
    