function [CRLB_Jnt,CRLB_Ind,CRLB_NoOfst]=MSLocSOSTCRLB(u,t,s,Qr,Qd)
% [CRLB_Jnt,CRLB_Ind,CRLB_NoOfst]=MSLocSOSTCRLB(u,t,s,Qr,Qd)
%
% This function obtains the CRLBs of an object position estimate under the 
% SOST scenario for the cases of (i) using both indirect- and direct-path 
% measurements, (ii) using indirect-path measurements only, and 
% (iii) without time and frequency offsets.
%
% Input parameters:
% u:         (Dim x 1), object position.
% t:         (Dim x 1), transmitter position.
% s:         (Dim x M), receiver position matrix, M is the number of receivers.         
% Qr:        (M x M), covariance matrix of the indirect range measurements.
% Qd:        (M x M), covariance matrix of the direct range measurements.
% 
% Output parameters:
% CRLB_Jnt:  (Dim x Dim), CRLB of object position for joint estimation 
%            with transmitter location and offsets.
% CRLB_Ind:  (Dim x Dim), CRLB of object position using indirect-path
%            measurements only without estimation of transmitter location
%            and offsets.
% CRLB_NoOfst: (Dim x Dim), CRLB of object position for joint estimation 
%              with transmitter location without time and frequency
%              offsets.
%
% The program can be used for 2D(Dim=2) or 3D(Dim=3) localization.
%
% Reference:
% Y. Zhang and K. C. Ho, "Multistatic localization in partially dynamic 
% scenario with only sensor positions svailable," IEEE Trans. Aerosp. 
% Electron. Syst., vol. 57, no. 5, pp. 3416-3432, Oct. 2021.
% 
% Yang Zhang and K. C. Ho   10-28-2021
%
%       Copyright (C) 2021
%       Computational Intelligence Signal Processing Laboratory
%       University of Missouri
%       Columbia, MO 65211, USA.
%       hod@missouri.edu
%

[~,M]=size(s);
H=[-ones(M-1,1),eye(M-1)]';

rho_dt=[];
rho_rt=[];
rho_ru=[];

for i=1:M
    rho_dti=(t-s(:,i))/norm(t-s(:,i));
    rho_dt=[rho_dt;rho_dti'];
    rho_rti=(t-u)/norm(t-u);
    rho_rt=[rho_rt;rho_rti'];
    rho_rui=(u-s(:,i))/norm(u-s(:,i))+(u-t)/norm(u-t);
    rho_ru=[rho_ru;rho_rui'];
end


rho_mrphi=[rho_rt,ones(M,1)];
rho_mdphi=[rho_dt,ones(M,1)];

K_MrMd=Qr+rho_mrphi*inv(rho_mdphi'*inv(Qd)*rho_mdphi)*rho_mrphi';
K_MrMdn=Qr+rho_rt*inv(rho_dt'*inv(Qd)*rho_dt)*rho_rt';
K_Mr=H*inv(H'*Qr*H)*H';

FIM=rho_ru'*inv(K_MrMd)*rho_ru;
FIMq=rho_ru'*K_Mr*rho_ru;
FIMn=rho_ru'*inv(K_MrMdn)*rho_ru;

CRLB_Jnt=inv(FIM);
CRLB_Ind=inv(FIMq);
CRLB_NoOfst=inv(FIMn);
