function [nav,dT_hist_ddd,dT_bins_ddd,missing_ddd] = read_nav_day_FATIMA_2022(path_gps0_ddd,path_gps1_ddd,...
                path_hed0_ddd,path_raw_images,prtit,PosLims)
%{
    Reads hourly gprm and heading files for one day and plots diagnostics.
    Loops through 24 hours, calling read_gprm_FATIMA_2022 and
        read_hed_roll_FATIMA_2022 for each hour.
    Raw 10Hz data is averaged to 1Hz output.
    
    Inputs: path_gps0_ddd: path to daily folder of Hemisphere gps hourly files
            path_gps1_ddd: path to daily folder of MRU gps hourly files
            path_hed0_ddd: path to daily folder of Hemisphere hed hourly files
            path_raw_images: string path to root daily folder for saving plots
            prtit: set true for plots
            PosLims: [lon_lo,lon_hi,lat_lo,lat_hi] limits
    
    Output: gprm: 86400 x 15 table of 1-Hz avg nav data
            Returns NaN for data gaps
    
            1   dnum        matlab datenum 1 Hz timestamp
            2   cog0        Hemisphere cog, deg
            3   sog0        Hemisphere sog, m/s
            4   Ngps0       Hemisphere N component of ship speed, m/s
            5   Egps0       Hemisphere E component of ship speed, m/s
            6   lat0        Hemisphere N latitude (+/- 90°)
            7   lon0        Hemisphere E longitude (0-360°)
            8   cog1        MRU cog, deg
            9   sog1        MRU sog, m/s
            10  Ngps1       MRU N component of ship speed, m/s
            11  Egps1       MRU E component of ship speed, m/s
            12  lat1        MRU latitude (+/- 90°)
            13  lon1        MRU E longitude (0-360°)
            14  hed0        Hemisphere heading, deg
            15  roll0       Hemisphere roll angle, deg
%}

fclose all;

warning ('off','MATLAB:MKDIR:DirectoryExists');
mkdir(fullfile(path_raw_images,'heading'));
mkdir(fullfile(path_raw_images,'sog'));
mkdir(fullfile(path_raw_images,'cog'));
mkdir(fullfile(path_raw_images,'track_map'));
mkdir(fullfile(path_raw_images,'QC_heading'));

[~,folder_name,~] = fileparts(path_gps0_ddd);
YY = folder_name(1:2); yyyy = 2000 + str2double(YY);
DDD = folder_name(3:5); ddd = str2double(DDD);
disp(['GPS/HED for day ',DDD]);

nav = NaN(86400,15);
start = datenum(yyyy-1,12,31) + ddd;
delta = double(1.0/86400);
last = start + 86399*delta;
dnum = start:delta:last; dnum = dnum';	% ref 1 Hz timestamp
nav(:,1) = dnum;
dT_bins_ddd = [];
missing_ddd = [];

%% load and average all nav data
for hh = 0:23
    HH = sprintf('%02i',hh);
    % average this hour's gps0 data into 1Hz nav array
    dfl0 = fullfile(path_gps0_ddd,['gprm' YY,DDD,HH,'_raw.txt']);
    gps0 = read_gprm_FATIMA_2022(dfl0,PosLims);
    begin = start + hh/24.0; % dnum start time this hour
    dnum_bin = begin:1/86400:begin + 1/24; % 3600 averaging bins this hour
    dnum_bin = dnum_bin(:);
    temp = interval_avg(gps0.dnum, gps0{:,2:7}, dnum_bin);
    ii = hh*3600 + 1; % index in hourly array for start of this hour
    nav(ii:ii+3599,2:7) = temp(:,2:7);

    % average this hour's gps1 data into 1Hz nav array
    dfl1 = fullfile(path_gps1_ddd,['gprm' YY,DDD,HH,'_raw.txt']);
    gps1 = read_gprm_FATIMA_2022(dfl1,PosLims);
    temp = interval_avg(gps1.dnum, gps1{:,2:7}, dnum_bin);
    nav(ii:ii+3599,8:13) = temp(:,2:7);
    
    % average this hour's hed0 data into 1 Hz nav array
    dfl3 = fullfile(path_hed0_ddd,['hed0' YY,DDD,HH,'_raw.txt']);
    [hed,dT_hist,dT_edges,missing] = read_hed_roll_FATIMA_2022(dfl3);
    if isempty(dT_bins_ddd)
        dT_hist_ddd = NaN(24,length(dT_hist));
        dT_hist_ddd(hh+1,:) = dT_hist;
        dT_bins_ddd = dT_edges(1:end-1) + mean(diff(dT_edges))/2;
        missing_ddd = [missing_ddd; missing];
    else
        dT_hist_ddd(hh+1,:) = dT_hist;
        missing_ddd = [missing_ddd; missing];
    end
    temp = interval_avg(hed.dnum, hed{:,2:3}, dnum_bin);
    nav(ii:ii+3599,14:15) = temp(:,2:3);
end

%% reformat to table for output
vars = {'dnum','cog0','sog0','Ngps0','Egps0','lat0','lon0','cog1','sog1','Ngps1','Egps1','lat1','lon1','hed0','roll0'};
units = {'none','deg','m/s','m/s','m/s','degN','degE','deg','m/s','m/s','m/s','degN','degE','deg','deg'};
nav = array2table(nav,'VariableNames',vars);
nav.Properties.VariableUnits = units;

%% plots
dt1 = datetime(nav.dnum,'ConvertFrom','datenum');
dnum_10min = start:10/1440:start + 1 - 10/1440;
dnum_10min = dnum_10min(:); % 10 min timestamp
if prtit
    figure; plot(dt1,nav.hed0,'LineWidth',1);
    ylabel('Heading, deg'); grid on; grid minor;
    title(['FATIMA Day ',DDD,' Ship Heading']);
    ppath = fullfile(path_raw_images,'heading',['hed_',DDD,'.png']);
    print('-dpng',ppath);

    figure; plot(datetime(dnum_10min,'ConvertFrom','datenum'),...
            missing_ddd,'o-','markersize',6,'LineWidth',1.5);
    ylabel('missing raw points per 10 min'); grid on; grid minor;
    title(['FATIMA Day ',DDD,' Heading Missing Data']); xlabel('hour of day');
    ppath = fullfile(path_raw_images,'QC_heading',['missing_hed_',DDD,'.png']);
    print('-dpng',ppath);

    figure; bar(dT_bins_ddd,nanmean1(dT_hist_ddd,1)); grid on;
    xlabel('Raw Delta Time, sec'); ylabel('counts'); ylim([0,1200]);
    title(['FATIMA Day ',DDD,' Daily Mean Heading Delta-T Distribution']);
    ppath = fullfile(path_raw_images,'QC_heading',['dT_hist_hed_',DDD,'.png']);
    print('-dpng',ppath);

    figure; plot(dt1,nav.cog1,dt1,nav.cog0,'LineWidth',1);
    ylabel('COG, deg'); grid on; grid minor;
    title(['FATIMA Day ',DDD,' GPS Course Over Ground']);
    legend('MRU','Hemisphere','location','best');
    ppath = fullfile(path_raw_images,'cog',['cog',DDD,'.png']);
    print('-dpng',ppath);

    figure; plot(dt1,nav.sog1,dt1,nav.sog0,'LineWidth',1);
    ylabel('SOG, m/s'); grid on; grid minor;
    title(['FATIMA Day ',DDD,' GPS Speed Over Ground']);
    legend('MRU','Hemisphere','location','best');
    ppath = fullfile(path_raw_images,'sog',['sog',DDD,'.png']);
    print('-dpng',ppath);

    map_FATIMA_2022(nav.lon0,nav.lat0,PosLims);
    ppath = fullfile(path_raw_images,'track_map',['track_map_',DDD,'.png']);
    print('-dpng',ppath);

    close all
end
