function [sonm,QC] = read_sonic_day_SeaState_2015(path_working_ddd,ddd,yyyy,...
         fsonic,prtit,path_raw_images,graphformat,graphdevice)
%{
Reads hourly sonic files for one day
Loops through 24 hours, calling read_sonic_SeaState_2015 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
        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 9 array of 1Hz avg data from high and low sonics
        Returns NaN for data gaps

        1   jd_ref      1 Hz timestamp
        2   U_hi        m/s
        3   V_hi        m/s
        4   W_hi        m/s
        5   Ts_hi       rad/s
        6   U_lo        m/s
        7   V_lo        m/s
        8   W_lo        m/s
        9   Ts_lo       rad/s
%}

fclose all;

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

for hhh = 0:23              % cycle thru 24 hourly sonm files
    hr = sprintf('%02i',hhh);
    dfl0 = fullfile(path_working_ddd,['mtx0' yyyy(3:4),jd,hr,'_raw.txt']);
    dfl1 = fullfile(path_working_ddd,['mtx1' yyyy(3:4),jd,hr,'_raw.txt']);
    
    if exist(dfl0,'file')==0 && exist(dfl1,'file')==0
        continue % no data this hour
    end;
    
    if exist(dfl0,'file')==2
        [son0,bad0,missing0] = read_sonic_SeaState_2015(dfl0,ddd,hhh);    
        jdson = son0(:,1);
        % average this hour's data into 1Hz 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, son0(:,2:5), jd_ref(ii:ii+3600));
        sonm(ii:ii+3599,2:5) = temp(:,2:5);
        QC0(1+hhh*6:1+hhh*6+5,:) = [bad0',missing0'];
    end;
    
    if exist(dfl1,'file')==2
        [son1,bad1,missing1] = read_sonic_SeaState_2015(dfl1,ddd,hhh);    
        jdson = son1(:,1);
        % average this hour's data into 1Hz 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, son1(:,2:5), jd_ref(ii:ii+3600));
        sonm(ii:ii+3599,6:9) = temp(:,2:5);
        QC1(1+hhh*6:1+hhh*6+5,:) = [bad1',missing1'];
    end;
    
    if prtit && mod(hhh,6)==0 % a few spectra plots
        if exist(dfl0,'file')==2
            Npts = length(son0(:,2));
            [Puu0,Fxx0] = psd2(detrend(son0(:,2)),Npts,fsonic,hamming(Npts));
            [Puu0,Fxx0] = specsmoo(Puu0,fsonic);
            [Pww0,~] = psd2(detrend(son0(:,4)),Npts,10,hamming(Npts));
            [Pww0,~] = specsmoo(Pww0,fsonic);
            [Ptt0,~] = psd2(detrend(son0(:,5)),Npts,10,hamming(Npts));
            [Ptt0,~] = specsmoo(Ptt0,fsonic);
        else
            Puu0 = NaN(48,1); Fxx0 = Puu0; Pww0 = Puu0; Ptt0 = Puu0;
        end;
        
        if exist(dfl1,'file')==2
            Npts = length(son1(:,2));
            [Puu1,Fxx1] = psd2(detrend(son1(:,2)),Npts,fsonic,hamming(Npts));
            [Puu1,Fxx1] = specsmoo(Puu1,fsonic);
            [Pww1,~] = psd2(detrend(son1(:,4)),Npts,10,hamming(Npts));
            [Pww1,~] = specsmoo(Pww1,fsonic);
            [Ptt1,~] = psd2(detrend(son1(:,5)),Npts,10,hamming(Npts));
            [Ptt1,~] = specsmoo(Ptt1,fsonic);
        else
            Puu1 = NaN(48,1); Fxx1 = Puu1; Pww1 = Puu1; Ptt1 = Puu1;
        end;
        
        figure;loglog(Fxx0,Fxx0.*Puu0,'b-',Fxx1,Fxx1.*Puu1,'r-');
        hold on;
        loglog(Fxx0(17:48),((Fxx0(17:48)).^(-2/3))./(Fxx0(35).^(-2/3)).*Puu0(35).*Fxx0(35),'g-');
        title([jd,' ',hr,' Spectrum of Velocity U']); grid;
        legend('Tower','Bow','location','northwest');
        ppath = fullfile(path_raw_images,['uSpectrum_',jd,'_',hr,'_',graphformat]);
        print(graphdevice,ppath);

        figure;loglog(Fxx0,Fxx0.*Pww0,'b-',Fxx1,Fxx1.*Pww1,'r-');
        hold on;
        loglog(Fxx0(17:48),((Fxx0(17:48)).^(-2/3))./(Fxx0(35).^(-2/3)).*Pww0(35).*Fxx0(35),'g-');
        title([jd,' ',hr,' Spectrum of Velocity W']); grid;
        legend('Tower','Bow','location','northwest');
        ppath = fullfile(path_raw_images,['wSpectrum_',jd,'_',hr,'_',graphformat]);
        print(graphdevice,ppath);

        figure;loglog(Fxx0,Fxx0.*Ptt0,'b-',Fxx1,Fxx1.*Ptt1,'r-');
        hold on;
        loglog(Fxx0(17:48),((Fxx0(17:48)).^(-2/3))./(Fxx0(35).^(-2/3)).*Ptt0(35).*Fxx0(35),'g-');
        title([jd,' ',hr,' Spectrum of Sonic Temperature']); grid;
        legend('Tower','Bow','location','northwest');
        ppath = fullfile(path_raw_images,['tSpectrum_',jd,'_',hr,'_',graphformat]);
        print(graphdevice,ppath);
    end;
end;

QC = [QC0,QC1];
