function hed = read_hed_day_WHOTS_2022(path_working_ddd,ddd,yyyy)
%{
    Reads raw 2 Hz heading files for one day, averages to 10-min
    Loops through 24 hours, calling read_hed_WHOTS_2018 for each hour.

    Inputs: path_working_ddd: path to daily folder containing hourly files
            ddd: day-of-year variable
            yyyy: year variable

    Output: hedm: 144x5 table of heading and pitch/roll angle
            Returns NaN for data gaps

        1   jd_ref      10 min timestamp
        2   hed         heading
        3   pitch       pitch angle
        4   roll        roll angle
        5   hed_std     std deviation in heading
%}

fclose all;

hedmm = NaN(144,7);
delta = double(10/1440);
last = ddd + 144*delta;
jd_bin = ddd:delta:last; 
jd_bin = jd_bin'; % 10 min bin edges
hedmm(:,1) = jd_bin(1:end-1);
JD = sprintf('%03i',ddd);
YYYY = sprintf('%04i',yyyy);

for hhh = 0:23 % cycle thru 24 hourly hed files
    hr = sprintf('%02i',hhh);
    dfl = fullfile(path_working_ddd,['hed0' YYYY(3:4),JD,hr,'_raw.txt']);
    if exist(dfl,'file')==2
        hed = read_hed_WHOTS_2022(dfl,ddd,hhh);
        hedm = table2array(hed);
    else
        continue % no data this hour
    end
    
    [hed_n, hed_e] = pol2cart(hedm(:,2)*pi/180,1);  % N/E heading components
    hedm = [hedm,hed_n,hed_e]; % add N/E components as cols 5,6
    
    % average this hour's data into 10min hedm array
    start = ddd + hhh/24.0; % jd start time this hour
    diff = jd_bin - start;  % look for closest time stamp to start
    [~,ii] = min(abs(diff));
    temp = interval_avg(hedm(:,1), hedm(:,2:6), jd_bin(ii:ii+6));
    hedmm(ii:ii+5,2:6) = temp(:,2:6);
    
    % compute 10-min heading std dev and add to array in col 7
    hed_unwp = unwrap(hed{:,2}*pi/180)*180/pi;
    temp = interval_std(hedm(:,1), hed_unwp, jd_bin(ii:ii+6));
    hedmm(ii:ii+5,7) = temp(:,2);
end
vars = hed.Properties.VariableNames; % save var names from last hed table
vars = [vars,'hed_std'];

% recompute hed from averaged E/N components
hed_10 = atan2(hedmm(:,6), hedmm(:,5))*180/pi; % from avg components
hed_10 = mod(hed_10 + 360*10, 360);
hedmm(:,2) = hed_10;
hedmm = hedmm(:,[1:4,7]); % discard E/N components

% output table
hed = array2table(hedmm,'VariableNames',vars);

end



