function [crlb] = computeCRLB_MPM(sigma, trueParams, us)
% for method of Tee et al. 
%   sigma - std of the accelerometer measurements

M = length(us); % number of positions

N = 100; % samples in the averaging?

s_x = trueParams.AxesMisalignment(1, 1); s_y = trueParams.AxesMisalignment(2, 2); s_z = trueParams.AxesMisalignment(3, 3);
alpha = trueParams.NonorthogAngles(1); beta = trueParams.NonorthogAngles(2); gamma = trueParams.NonorthogAngles(3);
b_x = trueParams.Bias(1); b_y = trueParams.Bias(2); b_z = trueParams.Bias(3);

As = zeros(9, 9, M); % to store A matrices below

for m = 1:M % for each position in the data

    u = us(m, :); u_x = u(1); u_y = u(2); u_z = u(3); % true specific force values in platform frame.

    As(:, :, m) = [u_x^2 0 0 0 0 0 u_x 0 0 % [d(mu)/d(s_x)...
                   0 (alpha*u_z + u_y)^2 0 s_y*u_z*(alpha*u_z + u_y) 0 0 0 alpha*u_z + u_y 0 % [d(mu)/d(s_y)...
                   0 0 (beta*u_x + gamma*u_y + u_z)^2 0 s_z*u_x*(beta*u_x + gamma*u_y + u_z) s_z*u_y*(beta*u_x + gamma*u_y + u_z) 0 0 beta*u_x + gamma*u_y + u_z % [d(mu)/d(s_z)...
                   0 s_y*u_z*(u_y + alpha*u_z) 0 (s_y*u_z)^2 0 0 0 s_y*u_z 0 % [d(mu)/d(alpha)...
                   0 0 s_z*u_x*(beta*u_x + gamma*u_y + u_z) 0 (s_z*u_x)^2 u_y*u_x*(s_z^2) 0 0 s_z*u_x % [d(mu)/d(beta)...
                   0 0 (s_y*u_y)*(beta*u_x + gamma*u_y + u_z) 0 u_x*u_y*(s_z^2) (s_z*u_y)^2 0 0 s_z*u_y % [d(mu)/d(gamma)...
                   u_x 0 0 0 0 0 1 0 0 % [d(mu)/d(b_x)...
                   0 alpha*u_z + u_y 0 s_y*u_z 0 0 0 1 0 % [d(mu)/d(b_y)...
                   0 0 beta*u_x + gamma*u_y + u_z 0 s_z*u_x s_z*u_y 0 0 1]; % [d(mu)/d(b_z)...
end

A = sum(As, 3);
J = (N / sigma^2) .* A;

crlb = diag(inv(J)); % returns a variance

end