function fpsm = read_paroP_day_SeaState_2015(path_working_ddd,ddd,yyyy,...
         freq,prtit,path_raw_images,graphformat,graphdevice)
%{
Reads hourly Paro fast pressure files for one day
Loops through 24 hours, calling read_paroP_SeaState_2015 for each hour
Raw 10Hz data is averaged to 1Hz output.
Creates spectrum plot for hours 0,6,12,18 and saves if desired.

Inputs: path_working_ddd: path to daily folder containing hourly files
        ddd: day-of-year variable
        yyyy: year string
        freq: 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: 36000 x 2 array of 1Hz avg data from high and low sonics
        Returns NaN for data gaps

        1   jd_ref      1 Hz timestamp
        2   fps0        pressure, mb or hPa, tower fast pressure
        3   fps1        pressure, mb or hPa, bow fast pressure
%}

fclose all;

fpsm = NaN(86400,3);
delta = double(1.0/86400);
last = ddd + 86400*delta;
jd_ref = [ddd:delta:last];	% ref 1 Hz timestamp
fpsm(:,1) = jd_ref(1:end-1)';
jd = sprintf('%03i',ddd);

for hhh = 0:23              % cycle thru 24 hourly fps files
    hr = sprintf('%02i',hhh);
    dfl0 = fullfile(path_working_ddd,['fps0' yyyy(3:4),jd,hr,'_raw.txt']);
    dfl1 = fullfile(path_working_ddd,['fps1' 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
        fps0 = read_paroP_SeaState_2015(dfl0,ddd,hhh);    
        jdfps = fps0(:,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(jdfps, fps0(:,2), jd_ref(ii:ii+3600));
        fpsm(ii:ii+3599,2) = temp(:,2);
    end;
    
    if exist(dfl1,'file')==2
        fps1 = read_paroP_SeaState_2015(dfl1,ddd,hhh);    
        jdfps = fps1(:,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(jdfps, fps1(:,2), jd_ref(ii:ii+3600));
        fpsm(ii:ii+3599,3) = temp(:,2);
    end;
    
    if prtit && mod(hhh,6)==0 % a few spectra plots
        if exist(dfl0,'file')==2
            Npts = length(fps0(:,2));
            [Pfp0,~] = psd2(detrend(fps0(:,2)),Npts,freq,hamming(Npts));
            [Pfp0,Fxx0] = specsmoo(Pfp0,freq);
        else
            Pfp0 = NaN(48,1); Fxx0 = Pfp0;
        end;
        
        if exist(dfl1,'file')==2
            Npts = length(fps1(:,2));
            [Pfp1,~] = psd2(detrend(fps1(:,2)),Npts,freq,hamming(Npts));
            [Pfp1,Fxx1] = specsmoo(Pfp1,freq);
        else
            Pfp1 = NaN(48,1); Fxx1 = Pfp1;
        end;

        figure;loglog(Fxx0,Fxx0.*Pfp0,'b-',Fxx1,Fxx1.*Pfp1,'r-');
%         hold on;
%         loglog(Fxx0(17:48),((Fxx0(17:48)).^(-2/3))./(Fxx0(35).^(-2/3)).*Pfp0(35).*Fxx0(35),'g-');
        title([jd,' ',hr,' Spectrum of Pressure']); grid;
        legend('Tower','Bow','location','northwest');
        ppath = fullfile(path_raw_images,['pSpectrum_',jd,'_',hr,'_',graphformat]);
        print(graphdevice,ppath);
        
        close all;
    end;
end;
