function [estParams] = jurman_2007(trajData, accel_MPU)
%JURMAN_2007_FCN Summary of this function goes here
%   Detailed explanation goes here
g = 9.81;

%% Get accelerations for each position
av_acc = computeAvAccels(trajData, accel_MPU);

U_apos = av_acc(1:3,:); U_aneg = av_acc(4:6,:);

U_as = U_apos + U_aneg; B_a = 0.5 * U_as; b_est = mean(B_a, 1); B_a = repmat(b_est, [3 1]);
U_ad = U_apos - U_aneg;

R = chol(1/(4*(g^2)) * (U_ad * U_ad'), 'lower'); % gives lower triangular matrix, R, such that

[T, S] = lu(R); % nonorthogonality and scale factor

Y = U_apos - B_a; % bias compensated measurements matrix.

M = inv(T) * inv(S) * Y * inv(9.81*eye(3)); % can estimate the scale factors, nonorthogonality and misalignment using the six position method.

S_est = S*T*M; S_est = S_est';

estParams.Bias = b_est; % estimated accelerometer biases
estParams.AxesMisalignment = S_est; % estimated accelerometer axis sensitivities

alpha = asin(T(2,1)); beta = asin(T(3,1)); gamma = asin(T(3,2)/cos(asin(T(3,1))));

estParams.ScaleFactors = diag(S);
estParams.NonorthogAngles = [alpha beta gamma]; % in rad

% Get misalignment angles back from S_est:
theta = asin(S_est(1,3) / S_est(1,1)); phi = asin(1/(cos(alpha)*cos(theta))*-(S_est(2,3)/S_est(2,2) - sin(alpha)*sin(theta))); psi = asin(-S_est(1,2) / (cos(theta)*S_est(1,1)));

estParams.MisalignAngles = [phi theta psi]; % in rad

end