function mot = read_motion_SeaState_2015(dfl,ddd,hhh)
%{
reads mot file specified by dfl and returns array with the
following columns:

1   jd_ref      10 Hz timestamp
2   accx        x axis acceleration, m/s
3   accy        y axis acceleration, m/s
4   accz        z axis acceleration, m/s
5   ratex       x axis rotation rate, rad/s
6   ratey       y axis rotation rate, rad/s
7   ratez       z axis rotation rate, rad/s

input parameters: dfl = string path to scs file
                  ddd = julian date
                  hhh = hour
%}

%% open mot file and read
disp(['Reading mot file for hour ',int2str(hhh)]);
flist = fopen(dfl,'r');
temp = {};    % cell array will hold one line of data, as read from file
z = [];       % destination matrix for data read from file
while ~feof(flist)
    try
        % temp will be 1x9 cell array, cells will be ~36000x1 arrays
        temp = textscan(flist,'%2f%2f%3f %f %f %f %f %f %f %*[^\n]','delimiter', ', ','headerlines', 1,'emptyvalue',NaN,'treatAsEmpty','Ship Motion');
        z = [z cell2mat(temp)'];  % 9x~36000 array
    catch
        for ii=1:9      % length of last cell reflects missing values from any/all fields
            if length(temp{1,ii})~=length(temp{1,9})
                temp{1,ii}(length(temp{1,9})+1) = [];  %truncate length of all cells to length(temp{1,9})
            end;
        end;
        if ~isempty(temp{1,1});
            z = [z cell2mat(temp)'];  % dim 2 will be < 36000 here
        end;
    end;
end;
fclose(flist);
tmotion =  ddd + (hhh + (z(1,:)+(z(2,:)+z(3,:)/1000)/60)/60)/24;

%% apply calibration constants & swap axes/change sign as needed
%  Change signs and swap axes such that:
%  +X is toward bow
%  +Y is toward port
%  +Z 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)

% MIRAI-2014 STRATUS-2015 values.
% calibration for Sundstrand rate meters and Systron accelerometers:
% Note arrays are ~36000x1
% accx = (-z(7,:))./3.9969293*9.80;         % forward is positive
% accy = (-z(8,:))./3.9604125.*9.80;        % to port is positive
% accz = (-z(9,:))./3.58409.*9.80;          % up is positive
% ratex = (+z(5,:))/0.02502/180*pi;         % port up is positive phi
% ratey = (-z(4,:))/0.02494/180*pi;         % bow down is positive theta
% ratez = (-z(6,:))/0.02503/180*pi;         % bow to port is positive psi (right-handed)
% Note: rate calib. is ~1/(100deg/s / 2.5Vdc) = 1/40 = 0.025

%MotionPak values as used on WHOTS-2014 & CALWATER-2015.
accx = (+z(8,:)+0)./3.7756.*9.80;         % forward is positive 
accy = (+z(7,:)+0)./3.8215.*9.80;         % to port is positive
accz = (-z(9,:)+0)./3.7005.*9.80;         % up is positive

ratex = (+z(5,:)+0)/0.025/180*pi;         % port up is positive phi
ratey = (+z(4,:)+0)/0.025/180*pi;         % bow down is positive theta
ratez = (-z(6,:)+0)/0.025023/180*pi;      % bow to port is positive psi (right-handed)

accx = replace_NaN_nearest_neighbor(accx);  % replace nans with nearest neighbor
accy = replace_NaN_nearest_neighbor(accy);
accz = replace_NaN_nearest_neighbor(accz);
ratex = replace_NaN_nearest_neighbor(ratex);
ratey = replace_NaN_nearest_neighbor(ratey);
ratez = replace_NaN_nearest_neighbor(ratez);

% figure;
% subplot(3,2,1);plot(accx);title('accx');
% subplot(3,2,2);plot(accy);title('accy');
% subplot(3,2,3);plot(accz);title('accz');
% subplot(3,2,4);plot(ratex);title('ratex');
% subplot(3,2,5);plot(ratey);title('ratey');
% subplot(3,2,6);plot(ratez);title('ratez');

%% compute gravity and format arrays
G = sqrt(mean(accx)^2+mean(accy)^2+mean(accz)^2);
disp(['*** Gravity = ', num2str(G,'%6.5f')]);

%% Interpolate
start = (ddd+hhh/24);
delta = double(1.0/864000);
last = start + 35999*delta;
tref = [start:delta:last];
[aa,ii,bb] = unique(tmotion);
saccx = interp1(tmotion(ii)',accx(ii)',tref','linear','extrap');    
saccy = interp1(tmotion(ii)',accy(ii)',tref','linear','extrap');
saccz = interp1(tmotion(ii)',accz(ii)',tref','linear','extrap');
sratex = interp1(tmotion(ii)',ratex(ii)',tref','linear','extrap');
sratey = interp1(tmotion(ii)',ratey(ii)',tref','linear','extrap');
sratez = interp1(tmotion(ii)',ratez(ii)',tref','linear','extrap');

mot = [tref' saccx saccy saccz sratex sratey sratez];

