function [snapshot, total_SNR_dB, R_nf, r_m_array, sv_mat, noise] = gen_sig_for_df(fc, target_rng, elev_deg, azi_deg, SNR_dB, txarray_loc, rxarray_loc, Nsnap, ch_error, noise_flag)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Generating raw data(snapshot) for direction finding 
% Start : 23.06.01
% End   : 23.08.25
% developed by Kwanggoo Yeo
% 
% Description
%   - Input
%       - 1) fc [vector], [Hz] : Center freqeuncy of radar Tx signal
%       - 2) target_rng [vector], [m] : Distance(range) from the radar to each target
%       - 3) elev_deg [vector], [deg] : Elevation angle for each taret
%       - 4) azi_deg [vector], [deg] : Azimuth angle for each taret
%       - 5) SNR_dB [vector], [dB] : SNR for each taret
%       - 6) txarray_loc [matrix], [m] : Locations of txarray elements in cartesian coordinate
%       - 7) rxarray_loc [matrix], [m] : Locations of rxarray elements in cartesian coordinate
%       - 8) Nsnap [scalar], [-] : The number of snapshots to generate
%       - 9) ch_error [vector], [-] : gain & phase error for each element in virtual array
%       - 10) noise_flag [scalar], [-] : Flag for noise which is included in snapshot (1) or not (0)
%
%   - Output
%       - 1) snapshot [matrix], [-] : Generated raw data(snapshot)
%       - 2) total_SNR_dB [scalar], [dB]   : SNR of snapshot, not SNR of signal in single element 
%       - 3) R_nf [scalar], [m]           : Distance of near-field (Frensel Region) 
%       - 4) r_m_array [matrix], [m]  : Distance from the targets to each element in virtual array
%       - 5) sv_amt [matrix], [-]  : Steering matrix
%       - 6) noise [matrix], [-]  : Generated noise in each element in virtual array
%
% History
% (23.08.25) Completed
%
% Referece
%   - 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Assumptions
%
% 1. Non-dispersion medium
% 2. Narrowband
% 3. No assumption on Near/Far-field
% 4. Geometry
%    4-1. Cartesian(x-y-z) 좌표계이며, 레이더 시스템의 local coordinate
%    4-2. Elevation angle : x-y 평면 기준으로 위 : + , 아래 : -
%    4-3. Azimuth angle : x-y 평면 상에서, x축의 왼쪽 : +, 오른쪽 : -
%    4-4. Virtual array (tx1, rx1) 채널의 위치를 원점으로 기준(coordinate origin)
%    4-5. x축 : 레이더 boresignt 방향
% 5. Monostatic radar system

%% Basic constants settings

c0 = physconst('LightSpeed');
lambda_c = c0/fc;
k_c = 2*pi / lambda_c;

Np = length(azi_deg); % Number of targets

% ch_error must be column vector
if size(ch_error, 1) == 1
    ch_error = ch_error.';
end

%% Array settings

Nt = size(txarray_loc, 2);
Nr = size(rxarray_loc, 2);

% Near field region calcuation
array_dist = zeros(Nt*Nr, 1);
for txidx = 1 : Nt
    for rxidx = 1 : Nr
        array_dist(Nr*(txidx-1) + rxidx) = sqrt(sum((txarray_loc(:,txidx) - rxarray_loc(:,rxidx)).^2));
    end
end
max_len_array = max(array_dist);
R_nf = 2 * max_len_array^2 / lambda_c;

% Target location calcuation in Cartesian
target_loc = zeros(3, Np);
for tidx = 1 : Np
    target_loc(:, tidx) = [target_rng(tidx)*cosd(elev_deg(tidx))*sind(azi_deg(tidx)) ; target_rng(tidx)*cosd(elev_deg(tidx))*cosd(azi_deg(tidx)); target_rng(tidx)*sind(elev_deg(tidx))];
end

% Generating steering vector(sv) matrix
r_m_array = complex(zeros(Nt*Nr, Np));
sv_mat = complex(zeros(Nt*Nr, Np));

for target_idx = 1 : Np
    for txarray_idx = 1 : Nt
        for rxarray_idx = 1 : Nr
            r_m_array(Nr*(txarray_idx-1)+rxarray_idx, target_idx) = sqrt(sum((txarray_loc(:, txarray_idx) - target_loc(:, target_idx)).^2)) + sqrt(sum((rxarray_loc(:, rxarray_idx) - target_loc(:, target_idx)).^2));
            sv_mat(Nr*(txarray_idx-1)+rxarray_idx, target_idx) = exp(1i * k_c * r_m_array(Nr*(txarray_idx-1)+rxarray_idx, target_idx));
        end
    end
    sv_mat(:, target_idx) = sv_mat(:, target_idx) .* ch_error;    
    sv_mat(:, target_idx) = sv_mat(:, target_idx) * conj(sv_mat(1, target_idx));
end


%% Signal Generation
% Calculating signal complex gain baed on SNR_dB
noise = (sqrt(0.5) * (randn(Nt*Nr, Nsnap) + 1i * randn(Nt*Nr, Nsnap)));
Pn = sum(diag(noise * noise')/Nsnap) / (Nt*Nr);

sig = zeros(Np, Nsnap);
for target_idx = 1 : Np
    Ps = 10^(SNR_dB(target_idx)/10) * Pn;
    sig(target_idx, :) = sqrt(Ps) * exp(1i * 2 * pi * rand(1, Nsnap));
end

% Signal model
if noise_flag == 1
    snapshot = sv_mat * sig + noise;
else
    snapshot = sv_mat * sig;
end

% Total SNR_dB
Ps_total = sum(diag((sv_mat * sig) * (sv_mat * sig)')/Nsnap) / (Nt*Nr);
total_SNR_dB = 10 * log10( Ps_total / Pn );

end