% ozone_10min_avg_ship.m
%{
    Averages and plots ozone data from the TEI 49C and chemiluminescent fast
    analyzer for specified range of days.
%}

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

%% adjust these parameters as desired, then run

% enter date range, year, leg#
% leg 1: jd 290:347, 2019
% leg 2: jd 348:365, 2019 & jd 001:054, 2020
% leg 3: jd 055:152, 2020
% leg 4: jd 160:225, 2020
% leg 5: jd 226:275, 2020

jd_start = 226;
jd_end = 275;
yyyy = 2020; YYYY = sprintf('%4i',yyyy); YY = YYYY(3:4);
leg = '5';

% root path to data files on this system
root_data_path = '/Users/byronb/Raw_Data_Files.noindex/MOSAIC_2019';
path_prog = fullfile(root_data_path,'_matlab');
slo3_path = fullfile(root_data_path,'ozone_slow_ship_bow','49C_Monitor');
fo3_path = fullfile(root_data_path,'ozone_fast_ship_bow');
fo3_LV_path = fullfile(root_data_path,'ozone_fast_ship_bow','FastO3#2_control_files');
plot_path = fullfile(root_data_path,'_scientific_analysis','raw_images','ozone_plots');
mkdir(plot_path)
save_path = fullfile(root_data_path,'_scientific_analysis','raw_1min');
mkdir(save_path)
save_path_10min = fullfile(root_data_path,'_scientific_analysis','flux_10min');
mkdir(save_path_10min)

plotit = 1; % set to 1 to save plot
N = jd_end - jd_start+ 1;

jd_bin = jd_start:1/144:jd_end+1; jd_bin = jd_bin';  % for 10-min average


%% set matlab script path
restoredefaultpath
cd(path_prog);
addpath(genpath(path_prog));
rehash toolboxcache;

%% load data
sfl = fullfile(save_path,[YYYY,'_Leg',leg,'_ozone_1min.mat']);
% if exist(sfl,'file')
%     load(sfl,'fo3','slo3');
% else

    %% load and average fast ozone data
    fo3 = NaN(N*1440,3); % empty array for 1-min data
    jj = 1;

    for ii = jd_start:jd_end
        JD = sprintf('%03i',ii);
        fldr = fullfile(fo3_path,[YY,JD]);
        if ~exist(fldr,'dir')
            disp([fldr,' not found'])
            continue
        end
        cd(fldr);
        files = dir('fo3*.txt');
        M = length(files);
        for kk = 1:M
            F = files(kk).name;
            HH = F(10:11); hh = str2num(HH);
            dfl = fullfile(fldr,F);
            fo3_raw = read_fo3_cnts_mosaic_2019(dfl,ii,hh);
            if all(isnan(fo3_raw(:,2)))
                continue
            end
            
            % find puff and remove (what a pain)
            o3_filt = smooth(fo3_raw(:,3),10/36000,'moving'); % smoothed ozone
            % look for puff at end of hour and delete data if found
            % need to use various criteria depending on mean value of ozone...
            if mean(o3_filt(32000:36000))<5
                xx = find(abs(diff(o3_filt(32000:36000)))>0.03,1,'first');
            elseif mean(o3_filt(32000:36000))<10
                xx = find(abs(diff(o3_filt(32000:36000)))>0.06,1,'first');
            elseif mean(o3_filt(32000:36000))<15
                xx = find(abs(diff(o3_filt(32000:36000)))>0.10,1,'first');
            else
                xx = find(abs(diff(o3_filt(32000:36000)))>0.20,1,'first');
            end
            if ~isempty(xx)
                if 32000+xx+200>36000
                    bad = 32000+xx-10:36000;
                else
                    bad = 32000+xx-10:32000+xx+200;
                end
                fo3_raw(bad,2:3) = NaN;
            end
            % delete first 20 sec no matter what
            fo3_raw(1:200,2:3) = NaN;
            % search remaining first 300 sec for pulse and remove if found
            if mean(o3_filt(1:3000))<5
                xx = find(abs(diff(o3_filt(1:3000)))>0.03,1,'first');
            elseif mean(o3_filt(1:3000))<10
                xx = find(abs(diff(o3_filt(1:3000)))>0.06,1,'first');
            elseif mean(o3_filt(1:3000))<15
                xx = find(abs(diff(o3_filt(1:3000)))>0.10,1,'first');
            else
                xx = find(abs(diff(o3_filt(1:3000)))>0.20,1,'first');
            end
            if ~isempty(xx)
                if xx-10<1
                    bad = 1:xx+200;
                else
                    bad = xx-10:xx+200;
                end
                fo3_raw(bad,2:3) = NaN;
            end
            
            % average 10 Hz to 1 min
            tref_1min_bin = ii+hh/24:1/1440:ii+(hh+1)/24;
            tref_1min_bin = tref_1min_bin';
            temp2 = interval_avg(fo3_raw(:,1),fo3_raw(:,2:end),tref_1min_bin);
            fo3(jj:jj+length(temp2)-1,:) = temp2;
            jj = find(isnan(fo3(:,1)),1,'first');
        end
    end

    jj = find(isnan(fo3(:,1)),1,'first');
    if ~isempty(jj)
        fo3 = fo3(1:jj-1,:);
    end

    %% load slow ozone
    
    % fractional day time indices for first and last minute of each hour where puff may occur
    puff_idx0 = (0.0:23.0)/24;         
    puff_idx1 = ((1.0:24.0)-1/60)/24;

    slo3 = NaN(N*1440,2); % empty array for 1-min data
    jj = 1;
    
    for ii = jd_start:jd_end
        [mm,dd] = yd2md(yyyy,ii); MM = sprintf('%02i',mm); DD = sprintf('%02i',dd);
        dfl = fullfile(slo3_path,['49C_Monitor_',YYYY,MM,DD,'.txt']);
        temp = read_TEI_O3_mosaic(dfl,ii,yyyy,mm,dd);
        if all(isnan(temp(:,2)))
            continue
        end
        bad = find(diff(temp(:,1))<0);  % remove data points out of sequence
        if bad==1
            temp(bad,:) = [];
        elseif bad == length(temp)-1
            temp(bad+1,:) = [];
        end
        for zz = 1:length(puff_idx0)
            xx = find(temp(:,1)==ii+puff_idx0(zz));
            temp(xx,2) = NaN; % remove data for first min of each hour
        end
        for zz = 1:length(puff_idx1)
            xx = find(temp(:,1)==ii+puff_idx1(zz));
            temp(xx,2) = NaN; % remove data for last min of each hour
        end
        slo3(jj:jj+length(temp)-1,:) = temp;
        jj = find(isnan(slo3(:,1)),1,'first');
    end
    
    jj = find(isnan(slo3(:,1)),1,'first');
    if ~isempty(jj)
        slo3 = slo3(1:jj-1,:);
    end
    
    cd(path_prog);

    % save 1-min O3 data
    sfl = fullfile(save_path,[YYYY,'_Leg',leg,'_ozone_1min.mat']);
    save(sfl,'fo3','slo3'); % save raw 1-min data in .mat
    
% end
    
%% load fast O3 1-min control files

sfl = fullfile(save_path,[YYYY,'_Leg',leg,'_fO3_control_1min.mat']);
% if exist(sfl,'file')
%     load(sfl,'fo3LV_1min');
% else
    fo3LV_1min = NaN(N*1440,22); % empty array for 1-min data
    jj = 1;
    
    for ii = jd_start:jd_end
        [mm,dd] = yd2md(yyyy,ii); MM = sprintf('%02i',mm); DD = sprintf('%02i',dd);
        cd(fo3_LV_path);
        files = dir(['control_',YY,'-',MM,'-',DD,'*.lvm']);
        for kk = 1:length(files)
            dfl = files(kk).name;
            temp = read_fo3_control_mosaic_2019(dfl,ii);
            [r,c] = size(temp);
            fo3LV_1min(jj:jj+r-1,:) = temp;
            jj = find(isnan(fo3LV_1min(:,1)),1,'first');
        end
    end
    
    jj = find(isnan(fo3LV_1min(:,1)),1,'first');
    if ~isempty(jj)
        fo3LV_1min = fo3LV_1min(1:jj-1,:);
    end
    
    % save 1 min control data
    save(sfl,'fo3LV_1min'); % save raw 1-min data in .mat
% end

cd(path_prog);

%% 10-min average
fo3_10min = interval_avg(fo3(:,1),fo3(:,2:end),jd_bin);
slo3_10min = interval_avg(slo3(:,1),slo3(:,2),jd_bin);
fo3LV_10min = interval_avg(fo3LV_1min(:,1),fo3LV_1min(:,2:end),jd_bin);
sfl = fullfile(save_path_10min,[YYYY,'_Leg',leg,'_ozone_10min.mat']);
save(sfl,'fo3_10min','slo3_10min','fo3LV_10min');

%% plot timeseries

if yyyy==2019
    dtslo3 = datetime(slo3(:,1)+datenum(2018,12,31),'convertfrom','datenum');
    dtfo3 = datetime(fo3(:,1)+datenum(2018,12,31),'convertfrom','datenum');
    dtfo3LV = datetime(fo3LV_1min(:,1)+datenum(2018,12,31),'convertfrom','datenum');
elseif yyyy==2020
    dtslo3 = datetime(slo3(:,1)+datenum(2019,12,31),'convertfrom','datenum');
    dtfo3 = datetime(fo3(:,1)+datenum(2019,12,31),'convertfrom','datenum');
    dtfo3LV = datetime(fo3LV_1min(:,1)+datenum(2019,12,31),'convertfrom','datenum');
end

ii = fo3LV_1min(:,11)>20;
jj = fo3LV_1min(:,11)<20;
figure('position',[100,600,900,700]);
subplot(2,1,1);
    plot(dtslo3,slo3(:,2),'-',dtfo3,fo3(:,3),'-','linewidth',2); grid;
    ylabel('ozone, ppb'); xlabel('date / time');  ylim([-0.1,50]);
    legend('TEI','Fast','location','southwest','fontsize',12);
    title(['Surface Ozone Timeseries, 1-min, Leg',leg],'fontweight','bold','fontsize',14);
subplot(2,1,2)
    plot(dtfo3LV(ii),fo3LV_1min(ii,11),'g.',dtfo3LV(jj),fo3LV_1min(jj,11),'r.','markersize',6); grid;
    ylabel('Inlet Flow, LPM'); xlabel(['day of year, ',YYYY]); ylim([-0.1,32]);
if plotit
    ppath = fullfile(plot_path,[YYYY,'_Leg',leg,'_ozone.png']);
    print('-dpng',ppath);
    ppath = fullfile(plot_path,[YYYY,'_Leg',leg,'_ozone.fig']);
    savefig(ppath)
end

% plot counts vs. o3 ppb
ii = isfinite(fo3(:,3)) & isfinite(fo3(:,2));
p_coef = polyfit(fo3(ii,2),fo3(ii,3),1);
save(fullfile(save_path,[YYYY,'_Leg',leg,'_fo3_calib.mat']),'p_coef');
ft = polyval(p_coef,[min(fo3(ii,2)),max(fo3(ii,2))]);
figure; plot(fo3(:,2),fo3(:,3),'.',[min(fo3(:,2)),max(fo3(:,2))],ft,'-','markersize',5);
grid minor; ylabel('O_3 ppb'); xlabel('counts');
txt = ['O_3 ppb = ',sprintf('%1.5g',p_coef(2)),' + ',sprintf('%1.5e',p_coef(1)),'*cnts'];
legend(txt,'location','north','fontsize',12);
title(['Fast Ozone raw Calibration, Leg ',leg]);
if plotit
    ppath = fullfile(plot_path,[YYYY,'_Leg',leg,'_ozone_raw_calibration.png']);
    print('-dpng',ppath);
end


