function [sonm,QC] = read_sonic_day_STRATUS_2017(path_working_ddd,ddd,yyyy,sonicmodel,...
         rotationsonic,fsonic,prtit,path_raw_images,graphformat,graphdevice)
%{
Reads hourly sonic files for one day
Loops through 24 hours, calling read_sonic_STRATUS_2017 for each hour
Raw 10Hz data is averaged to 1Hz output

Inputs: path_working_ddd: path to daily folder containing hourly files
        ddd: day-of-year variable
        yyyy: year string
        sonicmodel: string (WindMasterPro, R3, R3A, R2, R2A)
        rotationsonic: true for 30 deg rotation on some models
        fsonic: sampling frequency, Hz
        prtit: set 'true' to save spectra plots
        path_raw_images: path for saving plots
        graphformat: for plots, e.g. .png
        graphdevice: print driver, e.g. -dpng

Output: gprm: 86400 x 5 array of 1Hz avg data
        Returns NaN for data gaps

        1   jd_ref      10 Hz timestamp
        2   U           x axis acceleration, m/s
        3   V           y axis acceleration, m/s
        4   W           z axis acceleration, m/s
        5   Tsonic      x axis rotation rate, rad/s
        
        QC array, 144 x 2, # bad, # missing per 10 min

%}

fclose all;

sonm = zeros(86400,5)*NaN;
delta = double(1.0/86400);
last = ddd + 86400*delta;
jd_ref = ddd:delta:last; jd_ref=jd_ref';	% ref 1 Hz timestamp
sonm(:,1) = jd_ref(1:end-1)';
jd = sprintf('%03i',ddd);
QC = zeros(144,2)*NaN;      % for saving count of bad/missing points

for hhh = 0:23              % cycle thru 24 hourly sonm files
    hr = sprintf('%02i',hhh);
    dfl = fullfile(path_working_ddd,['son0' yyyy(3:4),jd,hr,'_raw.txt']);
    if exist(dfl,'file')==2
        [son,badPnts,missing] = read_sonic_STRATUS_2017(dfl,ddd,hhh,...
                                        sonicmodel,rotationsonic);
    else
        continue % no data this hour
    end
    jdson = son(:,1);
    % average this hour's data into 1Hz gprm array
    start = ddd + hhh/24.0; % jd start time this hour
    diff = jd_ref - start;  % look for closest time stamp to start
    [~,ii] = min(abs(diff));
    temp = interval_avg(jdson, son(:,2:5), jd_ref(ii:ii+3600));
    sonm(ii:ii+3599,2:5) = temp(:,2:5);
    QC(1+hhh*6:1+hhh*6+5,:) = [badPnts',missing'];

    if prtit && mod(hhh,6)==0 % a few spectra plots
        Npts = length(son(:,2));
        [Puu,~] = psd2(detrend(son(:,2)),Npts,fsonic,hamming(Npts));
        [Puu,Fxx] = specsmoo(Puu,fsonic);
        [Pww,~] = psd2(detrend(son(:,4)),Npts,10,hamming(Npts));
        [Pww,~] = specsmoo(Pww,fsonic);
        [Ptt,~] = psd2(detrend(son(:,5)),Npts,10,hamming(Npts));
        [Ptt,~] = specsmoo(Ptt,fsonic);

        figure;loglog(Fxx,Fxx.*Puu,'b-');
        hold on;
        loglog(Fxx(17:48),((Fxx(17:48)).^(-2/3))./(Fxx(35).^(-2/3)).*Puu(35).*Fxx(35),'r');
        title([jd,' ',hr,' Spectrum of Velocity U']); grid;
        ppath = fullfile(path_raw_images,['uSpectrum_',jd,'_',hr,'_',graphformat]);
        print(graphdevice,ppath);

        figure;loglog(Fxx,Fxx.*Pww,'g-');
        hold on;
        loglog(Fxx(17:48),((Fxx(17:48)).^(-2/3))./(Fxx(35).^(-2/3)).*Pww(35).*Fxx(35),'r');
        title([jd,' ',hr,' Spectrum of Velocity W']); grid;
        ppath = fullfile(path_raw_images,['wSpectrum_',jd,'_',hr,'_',graphformat]);
        print(graphdevice,ppath);

        figure;loglog(Fxx,Fxx.*Ptt,'r-');
        hold on;
        loglog(Fxx(17:48),((Fxx(17:48)).^(-2/3))./(Fxx(35).^(-2/3)).*Ptt(35).*Fxx(35),'r');
        title([jd,' ',hr,' Spectrum of Sonic Temperature']); grid;
        ppath = fullfile(path_raw_images,['tSpectrum_',jd,'_',hr,'_',graphformat]);
        print(graphdevice,ppath);
        
        close all;
    end
end

if prtit
    jd10 = ddd:10/1440:ddd+1-10/1440;
    figure; plot(jd10,QC(:,1),'bo-',jd10,QC(:,2),'r.-');
    title([jd,' ',hr,' Sonic Bad and Missing Counts / 10min']); grid;
    legend('Bad Pnts','Missing','location','best');
    ppath = fullfile(path_raw_images,['SonicQC_',jd,'_',hr,'_',graphformat]);
    print(graphdevice,ppath);
    close all;
end
end

