function [uvw,Ts,accplat,uvwplat,xyzplat,euler,plat_rate,lagPnts,mu,uvw_raw,...
         uvw_motcorr,f1,f2,f3] = motcorr4(son,mot,heading,sens_disp,fsam,decorr)
%{
    motcorr2.m - computes motion-corrected wind from raw sonic and motionPak
    data files.  Edson-98 method.
    
    motcorr3.m - decorrelation uses multiple regression approach over 10-min segments
    
    motcorr4.m - use with Inertial Labs MRU motion sensor.
    
    Required inputs are arrays son (36000x5), mot(36000x7), heading (36000x1) and
    sens_disp (1x3).  Additional input is sampling rate, typically 10Hz
    and decorr (boolean).
    
    Outputs are 36000x3 arrays uvw (winds), accplat (platform accelerations),
    uvwplat(platform velocities), xyzplat (platform displacements), lagPnts
    and mu, a 6x33 array of before/after corr coefs from the decorrelation.
    
    Requires functions angles.m, accels.m Filtcoef.m interval_decorr.m and neaqs_trans.m
    
    Plots comparison between E98 plat velocities and MRU-derived velocities.
%}

% define 10 Hz variables - transpose to row vectors to match legacy code
son_vel = son(:,2:4)';    % sonic wind velocities
plat_acc = mot(:,2:4)';   % platform accelerations, m/s
plat_rate = mot(:,5:7)';  % platform rotational rates, rad/s
uvw_mru = mot(:,8:10)';   % mru derived velocities
uvw_mru(2:3,:);
Ts = son(:,5)';

% unwrapped slow rotation about z axis
heading = unwrap(heading)';

% filter coefficients
[bhi,ahi,~,~,~,~] = Filtcoef(fsam,1/30);

% euler angles
[euler,~] = angles(ahi,bhi,fsam,plat_acc,plat_rate,heading);

% platform velocities, displacements and accelerations
[uvwplat,xyzplat,accplat] = accels2(bhi,ahi,fsam,plat_acc,euler);

% rotate uvw into earth frame
R = [sens_disp(1); sens_disp(2); sens_disp(3)] * ones(1,length(son_vel));
uvw = neaqs_trans(son_vel+cross(plat_rate,R),euler,0);
uvw_raw = uvw;

% compute lag between W and platform vertical velocity
% shift uvw accordingly and sum with uvwplat to correct winds
% uvw are then in earth coordinates, corrected for platform pitch/roll/yaw
lagPnts = lagcorrection(uvw(3,:),uvwplat(3,:),16);
if lagPnts > 0
    uvw(1,1:length(uvw(3,:))-lagPnts) = uvw(1,1:length(uvw(3,:))-lagPnts) + ...
                                        uvwplat(1,lagPnts+1:length(uvw(3,:)));
    uvw(2,1:length(uvw(3,:))-lagPnts) = uvw(2,1:length(uvw(3,:))-lagPnts) + ...
                                        uvwplat(2,lagPnts+1:length(uvw(3,:)));
    uvw(3,1:length(uvw(3,:))-lagPnts) = uvw(3,1:length(uvw(3,:))-lagPnts) + ...
                                        uvwplat(3,lagPnts+1:length(uvw(3,:)));
else
    uvw(1,-lagPnts+1:length(uvw(3,:))) = uvw(1,-lagPnts +1:length(uvw(3,:))) + ...
                                         uvwplat(1,1:length(uvw(3,:))+lagPnts);
    uvw(2,-lagPnts+1:length(uvw(3,:))) = uvw(2,-lagPnts +1:length(uvw(3,:))) + ...
                                         uvwplat(2,1:length(uvw(3,:))+lagPnts);
    uvw(3,-lagPnts+1:length(uvw(3,:))) = uvw(3,-lagPnts +1:length(uvw(3,:))) + ...
                                         uvwplat(3,1:length(uvw(3,:))+lagPnts);
end
uvw_motcorr = uvw;

Ts_raw = Ts;

dnum_bins = son(1,1):10/1440:son(1,1)+60/1440; dnum_bins = dnum_bins';
if decorr
    % Decorrelate wind-motion and Ts-motion.  Decorrelate in 10-min
    % segments with multiple regression
    X = [ones(length(uvw),1),accplat(1,:)',uvwplat(1,:)'];   % U
    [u_decorr,mu_u] = interval_decorr(son(:,1),uvw(1,:)',dnum_bins,X);
    X = [ones(length(uvw),1),accplat(2,:)',uvwplat(2,:)'];   % V
    [v_decorr,mu_v] = interval_decorr(son(:,1),uvw(2,:)',dnum_bins,X);
    X = [ones(length(uvw),1),accplat(3,:)',uvwplat(3,:)'];   % W
    [w_decorr,mu_w] = interval_decorr(son(:,1),uvw(3,:)',dnum_bins,X);
    X = [ones(length(Ts),1),accplat',uvwplat'];              % Tsonic
    [Ts_decorr,mu_Ts] = interval_decorr(son(:,1),Ts',dnum_bins,X);
    uvw = [u_decorr';v_decorr';w_decorr'];
    Ts = Ts_decorr';
    mu = [dnum_bins(1:end-1),mu_u,mu_v,mu_w,mu_Ts]; % concat all mu coefs
else
    mu = [dnum_bins(1:end-1),NaN(6,25)];
end

N = floor(length(son)/2);
dstr = datestr(son(1:1));
dt = datetime(son(:,1),'ConvertFrom','datenum');

seg = N-150:N+150;
f1 = figure('position',[200,400,600,1000]);
subtightplot(3,1,1,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
    plot(dt(seg),uvwplat(1,seg),'r',dt(seg),uvw_raw(1,seg),'b',...
         dt(seg),uvw_motcorr(1,seg),'g',dt(seg),uvw(1,seg),'k',...
         dt(seg),uvw_mru(1,seg),'--r','LineWidth',1); grid on; grid minor;
    ylabel('U velocity'); xlim([dt(N-150),dt(N+150)]);
    title([dstr,' Winds (Earth Frame): Raw, Motcorr, Decorr']);
    legend('uplat','raw','motcorr','decorr','mru uplat',...
           'location','northwest');
subtightplot(3,1,2,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
    plot(dt(seg),uvwplat(2,seg),'r',dt(seg),uvw_raw(2,seg),'b',...
         dt(seg),uvw_motcorr(2,seg),'g',dt(seg),uvw(2,seg),'k',...
         dt(seg),uvw_mru(2,seg),'--r','LineWidth',1); grid on; grid minor;
    ylabel('V velocity'); xlim([dt(N-150),dt(N+150)]);
    legend('vplat','raw','motcorr','decorr','mru vplat',...
           'location','northwest');
subtightplot(3,1,3,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
    plot(dt(seg),uvwplat(3,seg),'r',dt(seg),uvw_raw(3,seg),'b',...
         dt(seg),uvw_motcorr(3,seg),'g',dt(seg),uvw(3,seg),'k',...
         dt(seg),uvw_mru(3,seg),'--r','LineWidth',1); hold on;
    plot([dt(N-150),dt(N+150)],[0,0],'k');  hold off; grid on; grid minor;
    ylabel('W velocity'); xlim([dt(N-150),dt(N+150)]);
    legend('wplat','raw','motcorr','decorr','mru wplat',...
           'location','northwest');

f2 = figure('position',[200,400,600,1000]);
subtightplot(3,1,1,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
    plot(dt,uvw_raw(1,:),'b',dt,uvw(1,:),'g'); grid on; grid minor;
    ylabel('U velocity'); title([dstr,' Winds (Earth Frame): Raw, Motcorr, Decorr']);
    xlim([dt(1),dt(1)+hours(1)]); legend('raw','motcorr/decorr','location','northwest');
subtightplot(3,1,2,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
    plot(dt,uvw_raw(2,:),'b',dt,uvw(2,:),'g'); grid on; grid minor;
    ylabel('V velocity'); xlim([dt(1),dt(1)+hours(1)]);
subtightplot(3,1,3,[0.05,0.05],[0.05,0.05],[0.08,0.05]);
    plot(dt,uvw_raw(3,:),'b',dt,uvw(3,:),'g'); hold on;
    plot([dt(1),dt(end)],[0,0],'k'); grid on; grid minor;
    ylabel('W velocity'); xlim([dt(1),dt(1)+hours(1)]);

seg = N-300:N+300;
f3 = figure('position',[200,400,800,600]);
plot(dt(seg),Ts_raw(seg),'b',dt(seg),Ts(seg),'k'); grid on; grid minor;
ylabel('Tson, degC'); title([dstr,' Sonic Temperature: Raw, Decorr']);
xlim([dt(N-300),dt(N+300)]);legend('raw','decorr','location','northwest');

end



