function [licm,dT_hist_ddd,dT_bins_ddd,missing_ddd] = read_licor_day_FATIMA_2022(path_working_ddd,prtit,path_raw_images,read_raw)
%{
    Reads hourly licor files for one day and makes diagnostic plots
    Loops through 24 hours, calling read_licor_PISTON_MISOBOB_2018 for each hour
    Raw 10Hz data is averaged to 1Hz output
    
    Inputs: path_working_ddd: path to daily folder containing hourly files
            prtit: set true to make plots
            path_raw_images: path to root folder for saving plots
            read_raw: set true to force read of raw txt file
    
    Output: licm: 86400 x 8 array of 1-Hz avg gps data
            Returns NaN for data gaps
    
            1   dnum        matlab datenum 10 Hz timestamp
            2   h2o         H2O, mmol/m3
            3   co2         CO2, approx umol/mol or ppm
            4   Ta_lic      Licor box temperature, °C
            5   Pmb_lic     Licor box pressure, mb
            6   AGC         Licor diagnostic
            7   Rho_air     Licor air density, kg/m3
            8   q_lic       H2O, g/kg
%}

fclose all;

warning ('off','MATLAB:MKDIR:DirectoryExists');
mkdir(fullfile(path_raw_images,'qSpectra'));
mkdir(fullfile(path_raw_images,'QC_licor'));
mkdir(fullfile(path_raw_images,'co2_licor'));
mkdir(fullfile(path_raw_images,'q_licor'));

[~,folder_name,~] = fileparts(path_working_ddd);
YY = folder_name(1:2); yyyy = 2000 + str2double(YY);
DDD = folder_name(3:5); ddd = str2double(DDD);
disp(['Licor 7500 for day ',DDD]);

start = datenum(yyyy-1,12,31) + ddd;
delta = double(1.0/86400);
last = start + 86399*delta;
dnum = start:delta:last; dnum = dnum';	% ref 1 Hz timestamp

licm = NaN(86400,8);
licm(:,1) = dnum;

Rgas = 287.1;
Rgas_universal = 8.314472; % Pa m3 /K /mol
Mw = 18.01528;  % molar mass of H20 g/mol
Md = 28.964;    % molar mass dry air
epsilon = Mw/Md; % mass of water to mass of dry air = 0.622
dT_bins_ddd = [];
missing_ddd = [];

for hh = 0:23
    HH = sprintf('%02i',hh);
    dfl = fullfile(path_working_ddd,['lic0' YY,DDD,HH,'_raw.txt']);
    if exist(dfl,'file')==2
        [lic,dT_hist,dT_edges,missing] = read_licor_FATIMA_2022(dfl,read_raw);
        if isempty(dT_bins_ddd)
            dT_hist_ddd = NaN(24,length(dT_hist));
            dT_hist_ddd(hh+1,:) = dT_hist;
            dT_bins_ddd = dT_edges(1:end-1) + mean(diff(dT_edges))/2;
            missing_ddd = [missing_ddd; missing];
        else
            dT_hist_ddd(hh+1,:) = dT_hist;
            missing_ddd = [missing_ddd; missing];
        end
    else
        missing_ddd = [missing_ddd; ones(6,1)*6000];
        continue % no data this hour
    end

    % convert to mol/mol and specific humidity
    h2o_mr = (lic.h2o*1e-3).*Rgas_universal.*(lic.Ta_lic+273.15)./(lic.Pmb_lic*1e2);  % mmol/m3 to mol/mol
    rhoa_lic = (lic.Pmb_lic*1e2).*(1-h2o_mr*(1-epsilon))./(Rgas.*(lic.Ta_lic+273.15));% air density [kg/m3]
    q_lic = (lic.h2o*1e-3).*Mw./rhoa_lic;   % mmol/m3 to g/kg
    temp = [lic{:,2:6},rhoa_lic,q_lic];     % add columns to lic array

    % average this hour's data into 1Hz licm array
    begin = start + hh/24.0; % start time this hour
    dnum_bin = begin:1/86400:begin + 1/24; % 3600 averaging bins this hour
    dnum_bin = dnum_bin(:);
    temp = interval_avg(lic.dnum,[lic{:,2:6},rhoa_lic,q_lic],dnum_bin);
    ii = hh*3600 + 1; % index in hourly array for start of this hour
    licm(ii:ii+3599,2:8) = temp(:,2:8);

    if prtit && mod(hh,6)==0 % spectra plot
        Npts = length(q_lic);
        [Pqq,~] = psd2(detrend(despike2(q_lic)),Npts,10,hamming(Npts));
        Pqq = real(Pqq); [Pqq,Fxx] = specsmoo(Pqq,10);

        figure;loglog(Fxx,Fxx.*Pqq,'b-','LineWidth',1); hold on;
        loglog(Fxx(17:48),((Fxx(17:48)).^(-2/3))./(Fxx(35).^(-2/3)).*Pqq(35).*Fxx(35),...
               'r','LineWidth',1); hold off;
        title([DDD,' ',HH,' Spectrum of Specific Humidity']); grid;
        ppath = fullfile(path_raw_images,'qSpectra',['qSpectrum_',DDD,'_',HH,'.png']);
        print('-dpng',ppath);

        close all;
    end
end

vars = {'dnum','h20','co2','Ta_lic','Pmb_lic','agc','rhoa','q_lic'};
units = {'none','mmol/m3','ppm','degC','mb','none','kg/m3','g/kg'};
licm = array2table(licm,'VariableNames',vars);
licm.Properties.VariableUnits = units;

dt1 = datetime(dnum,'ConvertFrom','datenum');
dnum_10min = start:10/1440:start + 1 - 10/1440;
dnum_10min = dnum_10min(:); % 10 min timestamp

if prtit
    figure; plot(dt1,licm.agc,'LineWidth',1.5);
    ylabel('AGC'); grid on; grid minor;
    title(['FATIMA Day ',DDD,' Licor 7500 QC']);
    ppath = fullfile(path_raw_images,'QC_licor',['agc_licor_',DDD,'.png']);
    print('-dpng',ppath);

    figure; plot(dt1,despike(licm.co2),'LineWidth',1.5);
    ylabel('co2, ppm'); grid on; grid minor;
    title(['FATIMA Day ',DDD,' Licor CO2']);
    ppath = fullfile(path_raw_images,'co2_licor',['co2_licor_',DDD,'.png']);
    print('-dpng',ppath);

    figure; plot(dt1,despike(licm.q_lic),'LineWidth',1.5);
    ylabel('q, g/kg'); grid on; grid minor;
    title(['FATIMA Day ',DDD,' Licor Specific Humidity']);
    ppath = fullfile(path_raw_images,'q_licor',['q_licor_',DDD,'.png']);
    print('-dpng',ppath);
    
    figure; plot(datetime(dnum_10min,'ConvertFrom','datenum'),...
            missing_ddd,'o-','markersize',6,'LineWidth',1.5);
    ylabel('missing raw points per 10-min'); grid on; grid minor;
    title(['FATIMA Day ',DDD,' Licor Missing Data']); xlabel('hour of day');
    ppath = fullfile(path_raw_images,'QC_licor',['missing_licor_',DDD,'.png']);
    print('-dpng',ppath);

    figure; bar(dT_bins_ddd,nanmean1(dT_hist_ddd,1)); grid on;
    xlabel('Raw Delta Time, sec'); ylabel('counts'); ylim([0,3600]);
    title(['FATIMA Day ',DDD,' Daily Mean Licor Delta-T Distribution']);
    ppath = fullfile(path_raw_images,'QC_licor',['dT_hist_licor_',DDD,'.png']);
    print('-dpng',ppath);

    close all;
end

end
