function [CRLB_Jnt,CRLB_Ind,CRLB_NoOfst]=MSLocMOMTCRLB(u,t,vu,s,Qr,Qd,Qfr,Qfd)
% [CRLB_Jnt,CRLB_Ind,CRLB_NoOfst]=MSLocMOMTCRLB(u,t,vu,s,Qr,Qd,Qfr,Qfd)
%
% This function obtains the CRLBs of an object position estimate under the 
% MOMT 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.
% vu:        (Dim x 1), object velocity.
% 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.
% Qfr:       (M x M), covariance matrix of the indirect range-rate measurements.
% Qfd:       (M x M), covariance matrix of the direct range-rate measurements.
% 
% Output parameters:
% CRLB_Jnt:  (2 Dim x 2 Dim), CRLB of object position and velocity for 
%            joint estimation with transmitter location and offsets.
% CRLB_Ind:  (2 Dim x 2 Dim), CRLB of object position using indirect-path
%            measurements only.
% CRLB_NoOfst: (2 Dim x 2 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
%

[Dim,M]=size(s);
vt=zeros(Dim,1);

H=[-ones(M-1,1),eye(M-1)]';
Hq=[H,zeros(M,M-1);zeros(M,M-1),H];

rho_dt=[];
rho_rt=[];
rho_ru=[];
rho_fru=[];
rho_frt=[];
rho_fdt=[];

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'];
    rho_fdti=vt/norm(t-s(:,i))-((t-s(:,i))'*vt*(t-s(:,i)))/(norm(t-s(:,i))^3);
    rho_fdt=[rho_fdt;rho_fdti'];
    rho_frui=vu/norm(u-s(:,i))-((u-s(:,i))'*vu*(u-s(:,i)))/(norm(u-s(:,i))^3)+(vu-vt)/norm(u-t)-((u-t)'*(vu-vt)*(u-t))/(norm(u-t)^3);
    rho_fru=[rho_fru;rho_frui'];
    rho_frti=(vt-vu)/norm(u-t)-((t-u)'*(vt-vu)*(t-u))/(norm(u-t)^3);
    rho_frt=[rho_frt;rho_frti'];
end

Qrfr=blkdiag(Qr,Qfr);
Qdfd=blkdiag(Qd,Qfd);

rho_rfrtheta=[rho_ru,zeros(M,Dim);rho_fru,rho_ru];
rho_rfrphi=[rho_rt,ones(M,1),zeros(M,1);rho_frt,zeros(M,1),ones(M,1)];
rho_dfdphi=[rho_dt,ones(M,1),zeros(M,1);zeros(M,Dim),zeros(M,1),ones(M,1)];
rho_rfrt=[rho_rt;rho_frt];

K=Qrfr+rho_rfrphi*inv(rho_dfdphi'*inv(Qdfd)*rho_dfdphi)*rho_rfrphi';
Kn=Qrfr+rho_rfrt*inv(rho_dt'*inv(Qd)*rho_dt)*rho_rfrt';

Kq=Hq*inv(Hq'*Qrfr*Hq)*Hq';


FIM=rho_rfrtheta'*(inv(K))*rho_rfrtheta;
FIMq=rho_rfrtheta'*Kq*rho_rfrtheta;
FIMn=rho_rfrtheta'*inv(Kn)*rho_rfrtheta;

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

