function [h] = get_channel(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. 

    UE = struct();
%     UE.theta_UE     = -pi/2 + pi*rand();
%     UE.R_UE         = 10;
    v = [0;0;20];
    UE.pos=user+v.*simPar.xb(1, simPar.N)/simPar.c;
    UE.R_UE         = sqrt(UE.pos(1)^2+UE.pos(3)^2);
    UE.theta_UE     = atan(UE.pos(1)/UE.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                = 10e2;
    cdl_model.SampleDensity             = Inf;
    cdl_model.NumTimeSamples            = 1;
    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"
        h = exp(1j*simPar.k0*vecnorm(simPar.BS.antennaPos-UE.pos).');

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

        for idx = 1:simPar.channel.numPaths
            g = (randn() + 1i*randn())/sqrt(2); 
            rand_pos = [20*rand()-10;0;10*rand()];
            h = h + sqrt(1/(1+K)/simPar.channel.numPaths)*g*exp(1j*simPar.k0*vecnorm(simPar.BS.antennaPos-rand_pos).');
        end

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

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

    elseif simPar.channel.genMethod == "EM_Correlation"
        mu = 5*[-sin(UE.theta_UE); 0; cos(UE.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);
        h = V*sqrt(D)*w;

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

        for idx = 1:simPar.channel.numPaths
            g = (randn() + 1i*randn())/sqrt(2); 
            h = h + 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


end

