function [hed,dT_hist,dT_edges,missing] = read_hed_roll_FATIMA_2022(dfl)
%{
    reads heading/roll file from PSD Hemisphere system
    
    outputs:
    1   dnum        matlab datenum timestamp
    2   hed         heading, deg
    3   roll        roll angle, deg
    
    input parameters: dfl = string, full path to hed0 file
%}

[~,fname,~] = fileparts(dfl);
YY = fname(5:6); yyyy = 2000 + str2double(YY);
DDD = fname(7:9); ddd = str2double(DDD);
HH = fname(10:11); hh = str2double(HH);

%% 10 Hz reference timestamp and averaging bins
start = datenum(yyyy-1,12,31) + ddd + hh/24;
delta = double(0.5/86400);
dnum = start:delta:start + 1/24 - delta;
dnum = dnum(:); % exact 2 Hz timestamp
dnum_bin_10min = start:10/1440:start+1/24;  dnum_bin_10min = dnum_bin_10min(:); % 10 min bins

hed = NaN(7200,3);
hed(:,1) = dnum;

missing = ones(6,1)*7200/6;  % initialize to all missing
dT_hist = NaN(1,32);
dT_edges = NaN(1,33);

%% read data if file exists, if not return NaNs
if exist(dfl,'file')==2
    disp(['Reading hed roll file for ',DDD,' ',HH]);
    %% read file
    flist = fopen(dfl);
    temp = {};    % empty arrays for textscan
    stx = [];

    while feof(flist)==0	% read entire file into cell array
        try
            % temp will be 1x6 cell array, cells will be ~36000x1 arrays
            % 0000106 $PSAT,HPR,000002.00,184.19,,-1.73,N*3D
            temp = textscan(flist,'%2f%2f%3f %*5c %*3c %f %f %*f %f %*[^\n]','delimiter',', ','headerlines',1,'emptyvalue', NaN);
            stx = [stx cell2mat(temp)'];  % 6x~36000 array
        catch
            for ii=1:6      % length of last cell reflects missing values from any/all fields
                 if length(temp{1,ii})~=length(temp{1,6})
                     temp{1,ii}(length(temp{1,6})+1)=[];  %truncate length of all cells to length(temp{1,11})
                 end
            end
            if ~isempty(temp{1,1}) % save good values from last read, if any
                stx = [stx cell2mat(temp)'];
            end
        end
    end
    fclose(flist);
    stx = stx';
    dnum0 = start + stx(:,1)/1440 + (stx(:,2) + stx(:,3)/1000)/86400;
    % histogram of delta T for the raw data, should all be 0.5 sec
    [dT_hist,dT_edges] = histcounts(diff(dnum0)*86400,-0.2:0.1:3);
    % count missing points in each 10-min interval
    % should be 1200 points per 10min
    for ii=1:6
        missing(ii) = 1200 - numel(dnum0(dnum0>=dnum_bin_10min(ii) & dnum0<dnum_bin_10min(ii+1)));
    end

    heading = despike(stx(:,5));
    roll = despike(stx(:,6));

    %% format to 3x7200 2Hz output array
    % copy over original points to the 10 Hz grid - this is slow
    temp = dnum;
    for ii = 1:length(dnum0)                         % loop over all data points
        if any(isfinite(temp))                       % if there are unused grid points
            [~,jj] = nanmin1(abs(temp - dnum0(ii))); % find closest grid point for this datum
            dT = (dnum0(ii)-temp(jj))*86400*1000;    % compute deltaT in mSec
            if abs(dT)<250                           % only if dT is within range of +/- 250 ms...
                hed(jj,2) = heading(ii);             % copy data over to gridded array point
                hed(jj,3) = roll(ii);
                temp(jj) = NaN;                      % then remove this grid point from further consideration
            end
        end
    end
else
    % just return empty hed  and missing arrays
end

% reformat to table for output
vars = {'dnum','hed','roll'};
units = {'none','deg','deg'};
hed = array2table(hed,'VariableNames',vars);
hed.Properties.VariableUnits = units;

end

