function [motm,dT_hist_ddd,dT_bins_ddd,missing_ddd] = read_motion_day_FATIMA_2022(mru0_path_ddd,read_raw,...
                prtit,path_raw_images)
%{
    Reads hourly mru files for one day and plots diagnostics
    Loops through 24 hours, calling read_motion_mru_FATIMA_2022 for each hour
    10Hz data is averaged to 1Hz output
    
    Inputs: mru0_path_ddd: path to daily folder containing hourly mru files
            read_raw: set true to force read of raw text file
            prtit: set true to make plots
            path_raw_images: string path to root daily plot folder
    
    Output: motm: 86400 x 19 array of 1-Hz avg motion data
            Returns NaN for data gaps.
       
        1   dnum      matlab datenum timestamp
        2   hed       heading, deg
        3   pitch     pitch, deg
        4   roll      roll, deg
        5   ratex     x axis rotational rate, deg/s
        6   ratey     y axis rotational rate, deg/s
        7   ratezn    z axis rotational rate, deg/s
        8   accx      x axis acceleration, m/s2
        9   accy      y axis acceleration, m/s2
        10  accz      z axis acceleration, m/s2
        11  magx      x axis mag field strength, nT
        12  magy      y axis mag field strength, nT
        13  magz      z axis mag field strength, nT
        14  lat       latitude N, deg
        15  lon       longitude E, deg
        16  heave     vertial heave, m
        17  Espeed    easward speed, m/s
        18  Nspeed    northward speed, m/s
        19  Vspeed    vertical speed, m/s
%}

fclose all;
warning ('off','MATLAB:MKDIR:DirectoryExists');
mkdir(fullfile(path_raw_images,'rot_rate'));
mkdir(fullfile(path_raw_images,'accels'));
mkdir(fullfile(path_raw_images,'pitch_roll'));
mkdir(fullfile(path_raw_images,'heave_&_v_vel'));
mkdir(fullfile(path_raw_images,'n_e_velocity'));
mkdir(fullfile(path_raw_images,'QC_motion'));

[~,fldr_name,~] = fileparts(mru0_path_ddd);
YY = fldr_name(1:2); yyyy = 2000 + str2double(YY);
YYYY = sprintf('%4i',yyyy);
DDD = fldr_name(3:5); ddd = str2double(DDD);
disp(['MRU_P for day ',DDD]);


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

for hh = 0:23              % cycle thru 24 hourly motm files
    HH = sprintf('%02i',hh);
    dfl = fullfile(mru0_path_ddd,[YYYY,DDD,HH,'00MRU_P.txt']);
    [mot,dT_hist,dT_edges,missing] = read_motion_mru_FATIMA_2022(dfl,read_raw);
    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
    % average this hour's data into 1Hz motm array
    begin = start + hh/24.0; % start time this hour
    dnum_bin = begin:1/86400:begin + 1/24; % 3600 averaging bins this hour
    dnum_bin = dnum_bin(:);
    ii = hh*3600 + 1; % index in hourly array for start of this hour
    temp = interval_avg(mot.dnum, mot{:,8:25},dnum_bin);
    motm(ii:ii+3599,2:19) = temp(:,2:19);
end

vars = {'dnum','hed','pitch','roll','ratex','ratey','ratez','accx',...
        'accy','accz','magx','magy','magz','lat','lon','heave','Espeed',...
        'Nspeed','Vspeed'};
units = {'none','deg','deg','deg','deg/s','deg/s','deg/s','m/s2','m/s2',...
         'm/s2','nT','nT','nT','degN','degE','m','m/s','m/s','m/s'};
motm = array2table(motm,'VariableNames',vars);
motm.Properties.VariableUnits = units;

dt1 = datetime(dnum,'ConvertFrom','datenum');
dnum_10min = start:10/1440:start + 1 - 10/1440;
dnum_10min = dnum_10min(:); % 10 min timestamp
if prtit
    figure; plot(dt1,motm.ratex,dt1,motm.ratey,dt1,motm.ratez,'LineWidth',1);
    ylabel('Rotational Rate, deg/s');
    title(['FATIMA Day ',DDD,' MRU Rotational Rate']);
    legend('x axis','y axis','z axis','location','north'); grid on; grid minor;
    ppath = fullfile(path_raw_images,'rot_rate',['rot_rate_',DDD,'.png']);
    print('-dpng',ppath);

    figure; plot(dt1,motm.accx,dt1,motm.accy,dt1,motm.accz-9.8,'LineWidth',1);
    ylabel('Acceleratiuon, g');
    title(['FATIMA Day ',DDD,' MRU Accelerations']);
    legend('x axis','y axis','z axis - 9.8','location','north'); grid on; grid minor;
    ppath = fullfile(path_raw_images,'accels',['accels_',DDD,'.png']);
    print('-dpng',ppath);
    
    figure; plot(dt1,motm.pitch,dt1,motm.roll,'LineWidth',1);
    ylabel('Pitch & Roll Angles, deg');
    title(['FATIMA Day ',DDD,' MRU Pitch & Roll']);
    legend('pitch','roll','location','north'); grid on; grid minor;
    ppath = fullfile(path_raw_images,'pitch_roll',['pitch_roll_',DDD,'.png']);
    print('-dpng',ppath);

    figure; plot(dt1,motm.heave,dt1,motm.Vspeed,'LineWidth',1);
    ylabel('Heave & Vert. Vel, m & m/s');
    title(['FATIMA Day ',DDD,' MRU Heave & Vertical Velocity']);
    legend('heave','velocity','location','north'); grid on; grid minor;
    ppath = fullfile(path_raw_images,'heave_&_v_vel',['heave_v_vel_',DDD,'.png']);
    print('-dpng',ppath);

    figure; plot(dt1,motm.Nspeed,dt1,motm.Espeed,'LineWidth',1);
    ylabel('North Vel. & East Vel., m/s');
    title(['FATIMA Day ',DDD,' MRU North and East Velocity']);
    legend('north','east','location','north'); grid on; grid minor;
    ppath = fullfile(path_raw_images,'n_e_velocity',['n_e_velocity',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,' MRU Missing Data']); xlabel('hour of day');
    ppath = fullfile(path_raw_images,'QC_motion',['missing_mru_',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,7200]);
    title(['FATIMA Day ',DDD,' Daily Mean MRU Delta-T Distribution']);
    ppath = fullfile(path_raw_images,'QC_motion',['dT_hist_mru_',DDD,'.png']);
    print('-dpng',ppath);

    close all;
end
end