function [h] = get_channel1(simPar, user)
% Return a channel vector of length N according to the pre-defined simulation parameters. 
% Ensure each element of h is of unit energy. 
h = zeros(simPar.N,simPar.T+1);

randPos = zeros(3, simPar.channel.numPaths); 
g       = zeros(1, simPar.channel.numPaths); 

for idx = 1:simPar.channel.numPaths
    randPos(:, idx)   = [20*rand()-5;0;10*rand()];
    g(idx)            = (randn() + 1i*randn())/sqrt(2); 
end

if simPar.channel.genMethod == "CDL"
    cdl_model = nrCDLChannel('ChannelFiltering', false);
    cdl_model.MaximumDopplerShift       = simPar.v/simPar.lambda;
    cdl_model.UTDirectionOfTravel       = [0;0];
    cdl_model.CarrierFrequency          = simPar.fc;
    cdl_model.SampleRate                = 1/simPar.Ts;
    cdl_model.SampleDensity             = Inf;
    cdl_model.NumTimeSamples            = simPar.T+simPar.T_pre;
    cdl_model.TransmitAntennaArray.Size = [simPar.N, 1, 1, 1, 1];   % [M N P Mg Ng]
    cdl_model.ReceiveAntennaArray.Size  = [1, 1, 1, 1, 1]; 
    cdl_model.Seed                      = randi(10000);
    [pathGains, ~] = cdl_model();
    for idx = 1:simPar.T+simPar.T_pre
         h1(:,idx) = reshape(sum(pathGains(idx, :, :, :), 2), [simPar.N, 1]);
    end
    h = h1;
elseif simPar.channel.genMethod == "Ray_tracing"
    mx_number = simPar.N;
% mz_number = 50;
% my_number = 50;
% threshold = 0.2;
% threhold = 1e-3;
% mulength = 40;
% mz_total = -mz_number:mz_number;
% my_total = -my_number:my_number;
lam = simPar.lambda;
dy = lam/8; dz =  dy;
% iter_test = 30;
% R = zeros(2*my_number+1,2*mz_number+1);
% N = (2*mz_number+1)*(2*my_number+1);

% bsPosition = [22.287495, 114.140706]; % lat, lon
% bsAntSize = [mx_number 1];                    % number of rows and columns in rectangular array (base station)
% bsArrayOrientation = [0 0].';       % azimuth (0 deg is East, 90 deg is North) and elevation (positive points upwards) in deg
% uePosition = [22.287551+rand(1)*0.000008,114.140702+rand(1)*8*0.000001];  % lat, lon
% ueAntSize = [1 1];                    % number of rows and columns in rectangular array (UE).
% ueArrayOrientation = [180 0].';      % azimuth (0 deg is East, 90 deg is North) and elevation (positive points upwards)  in deg
% reflectionsOrder = 2;                 % number of reflections for ray tracing analysis (0 for LOS)

bsPosition = [22.287+rand(1)*0.001, 114.140+rand(1)*0.001]; % lat, lon
bsAntSize = [mx_number 1];                    % number of rows and columns in rectangular array (base station)
bsArrayOrientation = [90 0].';       % azimuth (0 deg is East, 90 deg is North) and elevation (positive points upwards) in deg
uePosition = [bsPosition(1)-0.000001*4+rand(1)*0.000008,bsPosition(2)+0.000001*20+rand(1)*8*0.000001-0.000001*4];  % lat, lon
ueAntSize = [1 1];                    % number of rows and columns in rectangular array (UE).
ueArrayOrientation = [270 0].';      % azimuth (0 deg is East, 90 deg is North) and elevation (positive points upwards)  in deg
reflectionsOrder = 2; 

bsSite = txsite("Name","Base station", ...
    "Latitude",bsPosition(1),"Longitude",bsPosition(2),...
    "AntennaAngle",bsArrayOrientation(1:2),...
    "AntennaHeight",4,...  % in m
    "TransmitterFrequency",3e8/lam);

ueSite = rxsite("Name","UE", ...
    "Latitude",uePosition(1),"Longitude",uePosition(2),...
    "AntennaHeight",4,... % in m
    "AntennaAngle",ueArrayOrientation(1:2));

% bsSite = txsite("Name","Base station", ...
%     "Latitude",bsPosition(1),"Longitude",bsPosition(2),...
%     "AntennaHeight",0,...  % in m
%     "TransmitterFrequency",3e8/lam);
% 
% ueSite = rxsite("Name","UE", ...
%     "Latitude",uePosition(1),"Longitude",uePosition(2),...
%     "AntennaHeight",0);

% show(bsSite);
% show(ueSite);

pm = propagationModel("raytracing","Method","sbr","MaxNumReflections",reflectionsOrder);
rays = raytrace(bsSite,ueSite,pm,"Type","pathloss");

% plot(rays{1})

pathToAs = [rays{1}.PropagationDelay]-min([rays{1}.PropagationDelay]);  % Time of arrival of each ray (normalized to 0 sec)
avgPathGains  = -[rays{1}.PathLoss];                                    % Average path gains of each ray
pathAoDs = [rays{1}.AngleOfDeparture];                                  % AoD of each ray
pathAoAs = [rays{1}.AngleOfArrival];                                    % AoA of each ray
isLOS = any([rays{1}.LineOfSight]);                                     % Line of sight flag

% for i = 1:40
cdl_model = nrCDLChannel('ChannelFiltering', false);
cdl_model.DelayProfile = 'Custom';
cdl_model.PathDelays = pathToAs;
cdl_model.AveragePathGains = avgPathGains;%!
cdl_model.AnglesAoD = pathAoDs(1,:);       % azimuth of departure
cdl_model.AnglesZoD = 90-pathAoDs(2,:);    % channel uses zenith angle, rays use elevation
cdl_model.AnglesAoA = pathAoAs(1,:);       % azimuth of arrival
cdl_model.AnglesZoA = 90-pathAoAs(2,:);    % channel uses zenith angle, rays use elevation
cdl_model.HasLOSCluster = isLOS;
% cdl_model.MaximumDopplerShift = 300;
% cdl_model.CarrierFrequency = 3e8/lam;
% cdl_model.SampleRate = 10e2;
%cdl_model.XPR = 0;
% cdl_model.SampleDensity = Inf;
% cdl_model.NumTimeSamples = 1;
% cdl_model.TransmitAntennaArray.Size = [1, 1, 1, 1, 1]; % [M N P Mg Ng]
% cdl_model.ReceiveAntennaArray.Size = [2*my_number+1, 2*mz_number+1, 1, 1, 1];
% cdl_model.ReceiveAntennaArray.PolarizationAngles = [45,0];%?
% cdl_model.ReceiveAntennaArray.ElementSpacing = [dy/lam,dz/lam,1,1];%?
%cdl_model.ChannelFiltering = false;
% cdl_model.Seed = randi(100*i);
% [pathGains, ~] = cdl_model();
% h = reshape(sum(pathGains(1, :, :, :), 2), [2*my_number+1, 2*mz_number+1]);
    cdl_model.MaximumDopplerShift       = simPar.v/simPar.lambda;
    cdl_model.UTDirectionOfTravel       = [0;0];
    cdl_model.CarrierFrequency          = simPar.fc;
    cdl_model.SampleRate                = 1/simPar.Ts;
    cdl_model.SampleDensity             = Inf;
    cdl_model.NumTimeSamples            = simPar.T+simPar.T_pre;
    cdl_model.TransmitAntennaArray.Size = [simPar.N, 1, 1, 1, 1];   % [M N P Mg Ng]
    cdl_model.ReceiveAntennaArray.Size  = [1, 1, 1, 1, 1]; 
    cdl_model.Seed                      = randi(10000);
    [pathGains, ~] = cdl_model();
    for idx = 1:simPar.T+simPar.T_pre
         h1(:,idx) = reshape(sum(pathGains(idx, :, :, :), 2), [simPar.N, 1]);
    end
    h = h1;
else
for i = 1:simPar.T
    UE1 = struct();
%     UE.theta_UE     = -pi/2 + pi*rand();
%     UE.R_UE         = 10;
    v = [0; 0; 20];
    t = simPar.x(1, (i-1)*simPar.N+1)/simPar.c; 
    UE1.pos=user+v*t; 

    UE1.R_UE         = sqrt(UE1.pos(1)^2+UE1.pos(3)^2);
    UE1.theta_UE     = atan(UE1.pos(1)/UE1.pos(3));
%     UE.pos          = UE.R_UE * [-sin(UE.theta_UE); 0; cos(UE.theta_UE)];
 
    
%     cdl_model = nrCDLChannel('ChannelFiltering', false);
%     cdl_model.MaximumDopplerShift       = 2*20/simPar.lambda;
%     cdl_model.CarrierFrequency          = simPar.fc;
%     cdl_model.SampleRate                = 1/t;
%     cdl_model.SampleDensity             = Inf;
%     cdl_model.NumTimeSamples            = 6;
%     cdl_model.TransmitAntennaArray.Size = [simPar.N, 1, 1, 1, 1];   % [M N P Mg Ng]
%     cdl_model.ReceiveAntennaArray.Size  = [1, 1, 1, 1, 1]; 
%     cdl_model.Seed                      = randi(10000);
%     

    if simPar.channel.genMethod == "NearField"
        h1 = exp(1j*simPar.k0*vecnorm(simPar.BS.antennaPos-UE1.pos).');

    elseif simPar.channel.genMethod == "NearField_SV"
        K = db2pow(simPar.channel.RicianFactor_dB); 
        h1 = sqrt(K/(1+K))*exp(1j*simPar.k0*vecnorm(simPar.BS.antennaPos-UE1.pos).'); 

        for idx = 1:simPar.channel.numPaths
            h1 = h1 + sqrt(1/(1+K)/simPar.channel.numPaths)*g(idx)*exp(1j*simPar.k0*vecnorm(simPar.BS.antennaPos-(randPos(:, idx) + v*t)).');
        end

    elseif simPar.channel.genMethod == "SV"
%         h1 = exp(1i*(0:simPar.N-1).'*simPar.k0*simPar.BS.spacing*sin(UE1.theta_UE));
        h1 = exp(1i*2*pi*(i-1)*simPar.Ts*20/simPar.lambda)*exp(-1i*(0:simPar.N-1).'*simPar.k0*simPar.BS.spacing*sin(atan(user(1)/user(3))));
%         h1 = h1/norm(h1);

%     elseif simPar.channel.genMethod == "CDL"
%         [pathGains, ~] = cdl_model();
%         h1 = reshape(sum(pathGains(1, :, :, :), 2), [simPar.N, 1]);

    elseif simPar.channel.genMethod == "EM_Correlation"
        mu = 5*[-sin(UE1.theta_UE); 0; cos(UE1.theta_UE)];

        theta = struct('rho', 2, 'v', [0; 0; 0], 'mu', mu);

        [Kh, ~] = get_K(simPar.x, simPar.p, theta, simPar.k0, simPar.c);
        assert(ishermitian(Kh), 'Hermitian condition not satisfied for GPR kernel'); 

        [V, D] = eig(Kh);
        w = (randn([simPar.N, 1]) + 1i*randn([simPar.N, 1]))/sqrt(2);
        h1 = V*sqrt(D)*w;

    elseif simPar.channel.genMethod == "VirtualChannel"
        % Convention: Phase Lag <-> exp(ikx) 
        K = db2pow(simPar.channel.RicianFactor_dB); 
        h1 = sqrt(K/(1+K))*exp(1i*(0:simPar.N-1).'*simPar.k0*simPar.BS.spacing*sin(UE1.theta_UE)); 

        for idx = 1:simPar.channel.numPaths
            g = (randn() + 1i*randn())/sqrt(2); 
            h1 = h1 + sqrt(1/(1+K)/simPar.channel.numPaths)*g*exp(1i*(0:simPar.N-1).'*simPar.k0*simPar.BS.spacing*sin(-pi/2+pi*rand()));
        end

    else
        error('Not Implemented.');

    end
    
    h(:,i) = h1;

end
h(:,simPar.T+1) = get_channel(simPar,user);
end
% check
% plot(angle(h(:,1)))
% hold on;
% plot(angle(h(:,2)))
% plot(angle(h(:,3)))
% legend; 

end

