function cost_MLE2 = TOARBLAncErrCalSrc_CostEval_MLE2(r, ANsy, d, pNsy, C, Q_r, Q_a, Q_d, Q_p, qHat, tHat, AHat)
%
%  cost_MLE2 = TOARBLAncErrCalSrc_CostEval_MLE2(r, ANsy, d, pNsy, C, Q_r, Q_a, Q_d, Q_p, qHat, tHat, AHat)
%
% This program evaluates the residual cost function of MLE2.  The code 
% is applicable for 2-D and 3-D localizations.
% Details of the algorithm development can be found in
% Benjian Hao, K. C. Ho and Zan Li,  "Range based rigid body localization 
% with a calibration emitter for mitigating anchor position uncertainties," 
% IEEE Trans. Wireless Commun., vol. 18, no. 12, pp. 5734-5748, Dec. 2019.
%
% Input Parameter List:
% r:        MxN, range (TOA) measurements, M is the num. of anchors, N is the number of sensors.
% ANsy:     KxM, noisy anchor positions, K is localization dimension.
% d:        Mx1, range (TOA) measurements for the calibration emitter.
% pNsy:     Kx1, noisy calibration emitter position.
% C:        KxN, sensor positions local to the rigid body.
% Q_r:      MNxMN, covariance matrix of range measurements.
% Q_a:      KMxKM, covariance matrix of anchor positions.
% Q_d:      MxM, covariance matrix of calibration emitter range measurements.
% Q_p:      KxK, covariance matrix of calibration emitter position.
% qHat:     Kx1, Euler angle vector estimate from MLE2, in degrees.
% tHat:     Kx1, position estimate from MLE2.
% AHat:     KxM, anchor position estimates from MLE2.
%
% Output Parameter List:
% cost_MLE2:    ML residual cost from MLE2.
%
% Benjian Hao, K. C. Ho    12-31-2019
%
%       Copyright (C) 2019
%       Computational Intelligence Signal Processing Laboratory
%       University of Missouri
%       Columbia, MO 65211, USA.
%       hod@missouri.edu
%

M = size(ANsy,2);   % number of anchors
N = size(C,2);      % number of sensors
K = size(ANsy,1);   % number of anchors

if (K==2)
    theta = qHat(1);
    QHat = [cosd(theta) -sind(theta); sind(theta) cosd(theta)];
else
    alpha = qHat(1); Q_ap = [1 0 0; 0 cosd(alpha) -sind(alpha); 0 sind(alpha) cosd(alpha)];
    beta = qHat(2);  Q_bt = [cosd(beta) 0  sind(beta); 0  1 0; -sind(beta)  0 cosd(beta)];
    gamma = qHat(3); Q_gm = [cosd(gamma) -sind(gamma) 0; sind(gamma) cosd(gamma) 0; 0 0 1];
    QHat = Q_gm * Q_bt * Q_ap;
end;

for n = 1:N
    for m = 1:M
        r_calculated(m,n)  = norm(QHat * C(:,n) + tHat - AHat(:,m));
    end
end

dd_dp = zeros(K,M);
for m = 1:M
    d_calculated(m)  = norm(pNsy - AHat(:,m));
    dd_dp(:,m) = -1/d_calculated(m) * (AHat(:,m) - pNsy);
end

R = [Q_r, zeros(M*N,K*M), zeros(M*N,M); zeros(K*M,M*N), Q_a, zeros(K*M,M); zeros(M,M*N), zeros(M,K*M), Q_d + dd_dp' * Q_p * dd_dp];

Theta = [r(:); ANsy(:); d];
Theta_Calculated = [r_calculated(:); AHat(:); d_calculated'];

cost_MLE2 = (Theta-Theta_Calculated)' * inv(R) * (Theta-Theta_Calculated);
