function [mot,dT_hist,dT_edges,missing] = read_motion_mru_FATIMA_2022(dfl,read_raw)
%{
    reads Intetial Labs MRU_P data file into 10 Hz matlab table structure
    adjusts axes/sign convention to PSL standard
    
    inputs: dfl: path to fixed-width format MRU_P data file
           read_raw: set true to force read of raw txt file

    raw file format:

1   MRU output log starting at 20220620000  (YYYYJJJHHMM)
2
3   Date  Time  Heading  Pitch  Roll  GyroX  GyroY GyroZ   AccX   AccY  AccZ  MagX  MagY  MagZV  Lat  Lon  Heave  Espeed  Nspeed  Vspeed 
4   
5   03/03/2022  00:00:00:.044  304.970000,-5.790000,2.000000,-1.700000,8.800000,0.400000,-0.014600,-0.039300,0.397100,3730.000000,18600.000000,-53920.000000,39.992902,-105.261563,0.000000,-0.020000,-0.010000,0.000000
6   03/03/2022  00:00:00:.094  304.960000,-5.810000,2.090000,-2.100000,8.100000,0.800000,-0.016300,-0.039900,0.396400,3790.000000,18590.000000,-53910.000000,39.992902,-105.261563,0.000000,-0.020000,-0.010000,0.000000
7   03/03/2022  00:00:00:.144  304.950000,-5.820000,2.180000,-1.400000,9.000000,0.500000,-0.015700,-0.041100,0.398100,3840.000000,18600.000000,-53900.000000,39.992902,-105.261563,0.000000,-0.020000,-0.010000,0.000000
8   03/03/2022  00:00:00:.195  304.950000,-5.840000,2.270000,-1.000000,8.900000,0.200000,-0.018000,-0.041800,0.398300,3910.000000,18610.000000,-53900.000000,39.992902,-105.261563,0.010000,-0.020000,-0.010000,0.000000

    mot format:
    1   dnum     matlab datenum timestamp
    2   yyyy     year
    3   mm       month
    4   dd       day of month
    5   hr       hour of day
    6   min      minute
    7   sec      decimal seconds
    8   hed      heading, deg
    9   pitch    pitch angle, deg
    10  roll     roll angle, deg
    11  ratex    x axis rotational rate (roll), rad/s
    12  ratey    y axis rotational rate (pitch), rad/s
    13  ratez    z axis rotational rate (yaw), rad/s
    14  accx     x axis acceleration (surge), m/s2
    15  accy     y axis acceleration (sway), m/s2
    16  accz     z axis acceleration, (heave) m/s2
    17  magx     x axis magnetic field strength, nT
    18  magy     y axis magnetic field strength, nT
    19  magz     z axis magnetic field strength, nT
    20  lat      latitude N, deg
    21  lon      longitude E, deg
    22  heave    vertical heave, m
    23  uplat    velocity to N, m/s, including ship speed
    24  vplat    velocity to W, m/s, including ship speed
    25  wplat    vertical velocity, m/s
%}

% define variable names for output table
vars = {'dnum','yyyy','mm','dd','hr','min','sec','hed','pitch','roll','ratex',...
        'ratey','ratez','accx','accy','accz','magx','magy','magz','lat',...
        'lon','heave','uplat','vplat','wplat'};
units = {'none','none','none','none','none','none','none','deg','deg','deg',...
         'rad/s','rad/s','rad/s','m/s2','m/s2','m/s2','nT','nT','nT','degN','degE',...
         'm','m/s','m/s','m/s'};

[fpath,fname,~] = fileparts(dfl);

dflm = [fname,'.mat'];
% load mru mat file if it exists, unless force read raw is indicated
if ~read_raw && exist(fullfile(fpath,dflm),'file')
    disp(['Reading MRU file ',fname]);
    load(fullfile(fpath,dflm),'mot','dT_hist','dT_edges','missing');
    return
end

ddd = str2double(fname(5:7));
hh = str2double(fname(8:9));
start = datenum(str2double(fname(1:4))-1,12,31) + ddd + hh/24;
% default 10 Hz timestamp
dnum = start:1/(86400*10):(start+1/24 - 1/(86400*10));  dnum = dnum(:);  % 10 Hz
dnum_bin = start:1/(86400*10):start+1/24; dnum_bin = dnum_bin(:);  % 10Hz bins
dnum_bin_10min = start:10/1440:start+1/24; dnum_bin_10min = dnum_bin_10min(:); % 10 min

missing = ones(6,1)*72000/6;
dT_hist = NaN(1,128);
dT_edges = NaN(1,129);

% read data if file exists, else return NaN array
if exist(dfl,'file')==2
    disp(['Reading MRU file ',fname]);
    
    flist = fopen(dfl,'r');
    temp = textscan(flist,['%f/%f/%f  %f:%f:%f:%f ',repmat('%f,',1,17),'%f %*[^\n]'],...
                           'headerlines',4,'emptyvalue',NaN);
    fclose(flist);
    mot0 = cell2mat(temp);
    mot0(:,6) = mot0(:,6) + mot0(:,7); % seconds
    dv = [mot0(:,3),mot0(:,1),mot0(:,2),mot0(:,4),mot0(:,5),mot0(:,6)]; % datevec
    dnum0 = datenum(dv);
    mot0 = [dnum0,dv,mot0(:,8:25)];

    % adust measurement units
    mot0(:,11:13) = mot0(:,11:13)*pi/180; % deg/s to rad/s
    mot0(:,14:16) = mot0(:,14:16)*9.80;   % g to m/s2

    % redefine axes to right handed convention:
    %{
        PSL standard axis/sign conventions:
         +x velocity is forward
         +y velocity is to port
         +z velocity is up
         +rotation about the x-axis is port up (phi = roll)
         +rotation about the y-axis is bow down (theta = pitch)
         +rotation about the z-axis is bow to port (psi = yaw)
         vplat is positive to W in earth coords, pitch/roll components only
         uplat is positive to N in earth coords, pitch/roll components only

        MRU convention, as mounted with connector facing bow (forward)
         +x acc/vel is to port
         +y acc/vel is to stern
         +z acc/vel is up
         +rotation about the x-axis is bow down (theta = pitch)
         +rotation about the y-axis is port down (phi = roll)
         +rotation about the z-axis is bow to port (psi = yaw)
         Espeed: includes ship speed and pitch/roll components
         Nspeed: includes ship speed and pitch/roll components
    %}
    mot1 = mot0; % duplicate array
    mot1(:,11) = -mot0(:,12);  % PSL ratex = -MRU ratey
    mot1(:,12) = mot0(:,11);   % PSL ratey = MRU ratex
    mot1(:,14) = -mot0(:,15);  % PSL accx = -MRU accy
    mot1(:,15) = mot0(:,14);   % PSL accy = MRU accx
    % not sure how to adjust sign of N/E speed since ship speed is included
    mot1(:,23) = -mot0(:,24);  % uplat = -MRU Nspeed???
    mot1(:,24) = mot0(:,23);   % vplat = MRU Espeed???
    mot1(:,17) = mot0(:,18);   % swap magx and magy ??? don't use these
    mot1(:,18) = mot0(:,17);
    
    % histogram of delta times
    [dT_hist,dT_edges] = histcounts(diff(mot0(:,1))*86400,-0.2:0.025:3);
    % count missing points in each 10-min interval
    % should be 12000 points per 10 min
    for ii=1:6
        missing(ii) = 12000 - numel(dnum0(dnum0>=dnum_bin_10min(ii) & dnum0<dnum_bin_10min(ii+1)));
    end

    % average to exact 10 Hz
    temp = interval_avg(mot0(:,1),mot0(:,8:25),dnum_bin);
    mot = [dnum,datevec(dnum),temp(:,2:end)];
    mot = array2table(mot,'VariableNames',vars);
    mot.Properties.VariableUnits = units;

    % save as .mat file for quicker loading next time
    save(fullfile(fpath,dflm),'mot','dT_hist','dT_edges','missing');
else
    mot = NaN(length(dnum),25);
    mot(:,1) = dnum;
    mot(:,2:7) = datevec(dnum);
    mot = array2table(mot,'VariableNames',vars);
    mot.Properties.VariableUnits = units;
end

end