function [sonm,dT_hist_ddd,dT_bins_ddd,QC] = read_sonic_day_FATIMA_2022(path_working_ddd,...
         rotationsonic,read_raw,prtit,path_raw_images)
%{
    Reads hourly sonic files for one day and plots diagnostics
    Loops through 24 hours, calling read_sonic_FATIMA_2022 for each hour
    Raw 10Hz data is averaged to 1Hz output sonm table
    
    Inputs: path_working_ddd: string path to daily folder with hourly files
            rotationsonic: set 'true' for 30 deg rotation on some models
            read_raw = set 'true' to force read of raw txt file
            prtit: set 'true' to save spectra plots
            path_raw_images: string path to root daily folder for saving plots
    
    Output: sonm: 86400 x 5 table of 1Hz avg data
            Returns NaN for data gaps
    
            1   dnum        matlab datenum 10 Hz timestamp
            2   U           bow-stern axis wind velocity, m/s
            3   V           port-starboard axis wind velocity, m/s
            4   W           z axis wind velocity, m/s
            5   Tsonic      sonic temperature, C
    
            QC table: badPnts, missing per 10 min raw data
%}

fclose all;

warning ('off','MATLAB:MKDIR:DirectoryExists');
mkdir(fullfile(path_raw_images,'uSpectra'));
mkdir(fullfile(path_raw_images,'wSpectra'));
mkdir(fullfile(path_raw_images,'tSpectra'));
mkdir(fullfile(path_raw_images,'QC_sonic'));
mkdir(fullfile(path_raw_images,'Rwdir'));
mkdir(fullfile(path_raw_images,'Rwspd'));
mkdir(fullfile(path_raw_images,'UVWTs'));

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

fsons = 20;  % 20 Hz data for FATIMA

sonm = zeros(86400,5)*NaN;
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
sonm(:,1) = dnum;
QC = zeros(144,2)*NaN;      % for saving count of bad/missing points
dT_bins_ddd = [];
for hh = 0:23              % cycle thru 24 hourly sonm files
    HH = sprintf('%02i',hh);
    dfl = fullfile(path_working_ddd,['son0',YY,DDD,HH,'_raw.txt']);
    [son,dT_hist,dT_edges,badPnts,missing] = read_sonic_FATIMA_2022(dfl,rotationsonic,read_raw,fsons);
    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;
    else
        dT_hist_ddd(hh+1,:) = dT_hist;
    end

    % average this hour's data into 1Hz 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(son{:,1}, son{:,2:5}, dnum_bin);
    ii = hh*3600 + 1; % index in hourly array for start of this hour
    sonm(ii:ii+3599,2:5) = temp(:,2:5);
    QC(1+hh*6:1+hh*6+5,:) = [badPnts,missing];

    if prtit && mod(hh,6)==0 && all(isfinite(son.u)) % a few spectra plots
        Npts = height(son.u);
        [Puu,~] = psd2(detrend(son.u),Npts,10,hamming(Npts));
        [Puu,Fxx] = specsmoo(Puu,10);
        [Pww,~] = psd2(detrend(son.w),Npts,10,hamming(Npts));
        [Pww,~] = specsmoo(Pww,10);
        [Ptt,~] = psd2(detrend(son.Tsonic),Npts,10,hamming(Npts));
        [Ptt,~] = specsmoo(Ptt,10);

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

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

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

        close all;
    end
end

dnum10 = start:10/1440:start+1-10/1440; dnum10 = dnum10';
dt10 = datetime(dnum10,'ConvertFrom','datenum');
dnum10bin = start:10/1440:start+1; dnum10bin = dnum10bin';
dt10bin = datetime(dnum10bin,'ConvertFrom','datenum');

vars = {'dnum','u','v','w','Tsonic'};
units = {'none','m/s','m/s','m/s','degC'};
sonm = array2table(sonm,'VariableNames',vars);
sonm.Properties.VariableUnits = units;
dt1s = datetime(sonm.dnum,'ConvertFrom','datenum');

QC = [dnum10,QC];
QC = array2table(QC,'VariableNames',{'dnum10','badPnts','missing'});

if prtit
    figure; plot(dt10,QC.badPnts,'b.-',dt10,QC.missing,'r.-','markersize',10,'LineWidth',1.5);
    ylabel('Bad or Missing Raw Points per 10min');
    title(['FATIMA Day ',DDD,' Sonic Bad and Missing']);
    legend('Bad Pnts','Missing','location','best'); grid on; grid minor;
    ppath = fullfile(path_raw_images,'QC_sonic',['bad_missing_',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,3000]);
    title(['FATIMA Day ',DDD,' Daily Mean Sonic Delta-T Distribution']);
    ppath = fullfile(path_raw_images,'QC_sonic',['dT_hist_sonic_',DDD,'.png']);
    print('-dpng',ppath);

    [rwspd,rwdir] = uv_to_sd(sonm.u,sonm.v);
    temp = interval_avg(sonm.dnum, rwspd, dnum10bin);
    figure; plot(dt10,temp(:,2),'b.-','markersize',10,'LineWidth',1.5);
    title(['FATIMA Day ',DDD,' Sonic Relative Wind Speed']);
    ylabel('Rel Wind Spd, m/s'); grid on; grid minor;
    ppath = fullfile(path_raw_images,'Rwspd',['Rwspd_',DDD,'.png']);
    print('-dpng',ppath);
    
    hr = 12;
    ii = hour(dt1s)==hr;  % pick data for 1200 UTC
    figure('position',[300,300,800,600]);
    subtightplot(2,1,1,[0.08,0.08],[0.05,0.05],[0.08,0.05]);
        plot(dt1s(ii),sonm.Tsonic(ii),'.','MarkerSize',8); grid on; grid minor;
        ylabel('Tsonic, °C');
        title(['FATIMA Day ',DDD,', 1Hz Sonic Data for ',sprintf('%02i',hr),'00 UTC']);
    subtightplot(2,1,2,[0.08,0.08],[0.05,0.05],[0.08,0.05]);
        plot(dt1s(ii),zeros(size(dt1s(ii))),'k-','LineWidth',1); hold on;
        h1 =plot(dt1s(ii),sonm.u(ii),'.',dt1s(ii),sonm.v(ii),'.',...
             dt1s(ii),sonm.w(ii),'.','MarkerSize',8); hold off;
        ylabel('U,V,W m/s'); grid on; grid minor;
        legend(h1,'u','v','w','Location','best');
    ppath = fullfile(path_raw_images,'UVWTs',['UVWTs_',DDD,'.png']);
    print('-dpng',ppath);
    
    rwdir(rwdir>180) = rwdir(rwdir>180) - 360;
    temp = interval_avg(sonm.dnum, rwdir, dnum10bin);
    figure; plot(dt10bin,zeros(size(dt10bin)),'k-','LineWidth',1); hold on;
    plot(dt10,temp(:,2),'b.-','markersize',10,'LineWidth',1.5); hold off;
    title(['FATIMA Day ',DDD,' Sonic Relative Wind Direction']);
    ylabel('Rel Wind Dir, +/- 180°'); grid on; grid minor; ylim([-182,182]);
    ppath = fullfile(path_raw_images,'Rwdir',['Rwdir_',DDD,'.png']);
    print('-dpng',ppath);

    close all;
end

end

