%CDLChannel TR 38.901 and TR 36.873 CDL channel model without Doppler
%
%   Note: This is an internal undocumented function and its API and/or
%   functionality may change in subsequent releases.

%   Copyright 2023 The MathWorks, Inc.

%#codegen

function staticChannelInfo = staticCDLChannel(model,info,pdp,coupling,Phi)

    %----------------------------------------------------------------------
    % dimensionality information
    
    insize = [1 model.NumInputSignals];
    D = getClusterDimensionality(insize,model.NumOutputSignals);
    L = numel(info.ClusterTypes);   % number of clusters
    M = D(3);
    P = D(4);                       % number of transmit antennas
    R = D(5);                       % number of receive antennas
    
    % get XPR matrix
    XPR = getXPR(model,info.ClusterTypes);

    %----------------------------------------------------------------------
    % static channel matrix (no Doppler shifts), departure and arrival unit
    % vectors for LOS clusters, sub-clustered (split) NLOS clusters and
    % other NLOS clusters (step 11 and ray offset angles of step 7). These
    % calculations do not account for Doppler shifts of each ray, which
    % are added in wireless.internal.channelmodels.timeVaryingCDLChannel.

    [HLOS,rhoTxLOS,rhoRxLOS] = generateClusterGains(model,insize,pdp,coupling,XPR,Phi,info.ClusterTypes,'LOS');                                 % size 1-by-LOS-M-by-P-by-R
    [HSplitNLOS,rhoTxSplitNLOS,rhoRxSplitNLOS] = generateClusterGains(model,insize,pdp,coupling,XPR,Phi,info.ClusterTypes,'SubclusteredNLOS');  % size 1-by(NSC*3)-by-M-by-P-by-R
    [HNLOS,rhoTxNLOS,rhoRxNLOS] = generateClusterGains(model,insize,pdp,coupling,XPR,Phi,info.ClusterTypes,'NLOS');                             % size 1-by-(Lin-LOS-NSC)-by-M-by-P-by-R

    %----------------------------------------------------------------------
    % assembly of overall static channel matrix and unit vectors

    LOS = strcmpi(info.ClusterTypes,'LOS');
    SubclusterNLOS = strcmpi(info.ClusterTypes,'SubclusteredNLOS');
    NLOS = strcmpi(info.ClusterTypes,'NLOS');
    
    pathGains = complex(zeros(1,L,M,P,R),zeros(1,L,M,P,R));
    pathGains(:,LOS,:,:,:) = HLOS;
    pathGains(:,SubclusterNLOS,:,:,:) = HSplitNLOS;
    pathGains(:,NLOS,:,:,:) = HNLOS;

    % Unit vector of each departure ray
    departureUnitVectors = zeros(3,L,M); % rho_hat_tx
    departureUnitVectors(:,LOS,:) = rhoTxLOS;
    departureUnitVectors(:,SubclusterNLOS,:) = rhoTxSplitNLOS;
    departureUnitVectors(:,NLOS,:) = rhoTxNLOS;

    % Unit vector of each arrival ray
    arrivalUnitVectors = zeros(3,L,M); % rho_hat_rx
    arrivalUnitVectors(:,LOS,:) = rhoRxLOS;
    arrivalUnitVectors(:,SubclusterNLOS,:) = rhoRxSplitNLOS;
    arrivalUnitVectors(:,NLOS,:) = rhoRxNLOS;

    %----------------------------------------------------------------------
    %  normalizations
    
    if (model.NormalizeChannelOutputs)
        pathGains = pathGains / sqrt(size(pathGains,5));
    end

    if (model.NormalizePathGains)
        [~,power] = wireless.internal.channelmodels.namedDelayProfileColumns();
        p = sum(10.^(pdp(:,power)/10));
        pathGains = pathGains / sqrt(p);
    end

    %----------------------------------------------------------------------
    %  Collect path gains and unit vectors in output struct

    staticChannelInfo = struct();
    staticChannelInfo.StaticPathGains = pathGains;
    staticChannelInfo.DepartureUnitVectors = departureUnitVectors;
    staticChannelInfo.ArrivalUnitVectors = arrivalUnitVectors;

end

% Generation of cluster gains, covering:
%    * Section 7.5 Step 7: ray offset angles
%    * Section 7.5 Step 11: Generate channel coefficients (except Doppler)
% Note that the phases produced by "Step 10: Draw initial random phases"
% are provided as an input 'Phi' here and used to created the polarization
% term (also, 'Phi' can be externally specified rather than randomly
% generated depending on MODEL.InitPhase). Ray coupling indices produced by
% "Step 8: coupling of rays" are provided as an input 'coupling' here. Note
% that steps 5 (Generate delays), 6 (Generate cluster powers), 7 (Generate
% arrival and departure angles) and 9 (Generate XPRs) are not performed and
% the values of these small-scale parameters for a single link are provided
% either from the tabulated CDL profiles or from the custom CDL option
% (MODEL.DelayProfile='Custom'). The Doppler effect is added to the channel
% coefficients in wireless.internal.channelmodels.timeVaryingCDLChannel.
function [Hstatic,rhat_tx,rhat_rx] = generateClusterGains(model,insize,pdp,coupling,XPR,Phi,allClusterTypes,thisClusterType)
    
    persistent c0;
    if isempty(c0)
        c0 = physconst('LightSpeed');
    end
        
    % named column subscripts
    [~,power,AoD,AoA,ZoD,ZoA] = wireless.internal.channelmodels.namedDelayProfileColumns();
    
    % get dimensionality information for a single cluster
    D = getClusterDimensionality(insize,model.NumOutputSignals);
    M = D(3);   % number of rays per cluster (see R_i below)
    S = D(4);   % number of input signals
    U = D(5);   % number of output signals
    
    % get part of delay profile corresponding to this cluster type
    pdp = pdp(strcmpi(allClusterTypes,thisClusterType),:);

    % get number of clusters and set up the corresponding dimension of the
    % channel matrix output 'H' ('H' will have dimension sizes 'D')
    N = size(pdp,1);
    D(2) = N;
    
    % Initialize channel matrix and unit vectors
    Hstatic = zeros(D,'like',1i);
    rhat_tx = zeros([3 D(2:3)]); % 3-by-N-by-M
    rhat_rx = zeros([3 D(2:3)]); % 3-by-N-by-M

    % if no clusters are active for this cluster type, return empty channel
    % matrix
    if (N==0)
        return;
    end
    
    % convert cluster powers to a linear scale
    P = 10.^(pdp(:,power).'/10);    

    % get per-cluster parameters
    if isfield(model,'DelayProfile') && ~strcmpi(model.DelayProfile,'Custom') % nrCDLChannel model with DelayProfile different from Custom
        C_ASD = model.AngleSpreads(1);
        C_ASA = model.AngleSpreads(2);
        C_ZSD = model.AngleSpreads(3);
        C_ZSA = model.AngleSpreads(4);
    else % lte3DChannel model or custom DelayProfile of nrCDLChannel
        C_ASA = model.DesiredASA;
        C_ASD = model.DesiredASD;
        C_ZSA = model.DesiredZSA;
        C_ZSD = model.DesiredZSD;
    end
    
    % get part of ray coupling matrix corresponding to this cluster type
    coupling = coupling(strcmpi(allClusterTypes,thisClusterType),:,:);

    % get part of XPR matrix corresponding to this cluster type
    XPR = XPR(strcmpi(allClusterTypes,thisClusterType),:,:);

    % get part of Phi matrix corresponding to this cluster type
    Phi = Phi(strcmpi(allClusterTypes,thisClusterType),:,:);

    % apply ray offset angles to AOA and AOD angles according to Section
    % 7.5 Equation 7.5-13
    ray_offset_alpha = kron([0.0447 0.1413 0.2492 0.3715 0.5129 0.6797 0.8844 1.1481 1.5195 2.1551],[1 -1]); % size 1-by-M
    if (strcmpi(thisClusterType,'LOS'))
        ray_offset_alpha(:) = 0; % no ray offseting for LOS path (only first ray will be used anyway)
    end
    phi_AOA = pdp(:,AoA) + C_ASA * ray_offset_alpha;   % size N-by-M
    phi_AOD = pdp(:,AoD) + C_ASD * ray_offset_alpha;   % size N-by-M
    
    % apply ray offset angles to ZOA and ZOD angles according to Section
    % 7.5 Equation 7.5-18
    theta_ZOA = pdp(:,ZoA) + C_ZSA * ray_offset_alpha; % size N-by-M
    theta_ZOD = pdp(:,ZoD) + C_ZSD * ray_offset_alpha; % size N-by-M
    
    % scale angles according to Section 7.7.5.1
    if(model.AngleScaling)
        phi_AOA = wireless.internal.channelmodels.scaleCDLAngles(model,pdp,phi_AOA,'DesiredASA','DesiredMeanAoA');
        phi_AOD = wireless.internal.channelmodels.scaleCDLAngles(model,pdp,phi_AOD,'DesiredASD','DesiredMeanAoD');
        theta_ZOA = wireless.internal.channelmodels.scaleCDLAngles(model,pdp,theta_ZOA,'DesiredZSA','DesiredMeanZoA');
        theta_ZOD = wireless.internal.channelmodels.scaleCDLAngles(model,pdp,theta_ZOD,'DesiredZSD','DesiredMeanZoD');
    end
    
    % wrap azimuth angles to [-180,180]
    phi_AOA = wrapAzimuthAngles(phi_AOA);
    phi_AOD = wrapAzimuthAngles(phi_AOD);
    
    % wrap zenith angles to [0,360] and map [180,360] to [180,0]
    theta_ZOA = wrapZenithAngles(theta_ZOA);
    theta_ZOD = wrapZenithAngles(theta_ZOD);
    
    % perform coupling of rays (step 8)
    % * AoD remains in the original order
    % * AoA rearranged w.r.t. AoD by page 1
    % * ZoD rearranged w.r.t. AoD by page 3
    % * ZoA rearranged w.r.t. ZoD by page 2. Therefore, ZoA is rearranged 
    %   w.r.t. AoD firstly by page 2 and then by page 3
    %
    %           1
    %        .-----> AoA
    %       /
    %     AoD
    %       \   3          2
    %        '-----> ZoD -----> ZoA
    %
    % AoD to AoA coupling
    phi_AOA = phi_AOA(coupling(:,:,1));
    % AoD to ZoD coupling
    theta_ZOD = theta_ZOD(coupling(:,:,3));
    % ZoD to ZoA coupling
    theta_ZOA = theta_ZOA(coupling(:,:,2));
    % rearrange ZoA in AoD order
    theta_ZOA = theta_ZOA(coupling(:,:,3));
    
    % spherical unit vectors of departure for each cluster and each ray
    rhat_tx = wireless.internal.channelmodels.getRhoHat(phi_AOD,theta_ZOD); % size 3-by-N-by-M
    
    % transmit antenna/subarray location vector
    lambda_0 = c0 / model.CarrierFrequency;  % carrier wavelength, m
    if ~isempty(model.TransmitAntennaArray.SubarrayPositions) % Partitioned/Replicated subarrays
        txRadiatorPositions = model.TransmitAntennaArray.SubarrayPositions;
    else
        txRadiatorPositions = model.TransmitAntennaArray.ElementPositions;
    end
    dbar_tx = reshape(txRadiatorPositions,3,[]) * lambda_0 + model.TransmitAntennaArray.Position; % size 3-by-S
    
    transmitFieldTerm = complex(zeros(2,N*M,S));
    transmitLocationTerm = complex(zeros(N*M,S),zeros(N*M,S));
    % for each transmit antenna:
    for s = 1:S
        
        % get transmit field term
        transmitFieldTerm(:,:,s) = wireless.internal.channelmodels.getFieldTerm(model.TransmitAntennaArray,theta_ZOD,phi_AOD,s);
        
        % get transmit location term
        transmitLocationTerm(:,s) = getLocationTerm(rhat_tx,dbar_tx,s,model.CarrierFrequency);

    end
    
    % spherical unit vectors of arrival for each cluster and each ray
    rhat_rx = wireless.internal.channelmodels.getRhoHat(phi_AOA,theta_ZOA); % size 3-by-N-by-M
    
    % receive antenna/subarray location vectors
    if ~isempty(model.ReceiveAntennaArray.SubarrayPositions) % Partitioned/Replicated subarrays
        rxRadiatorPositions = model.ReceiveAntennaArray.SubarrayPositions;
    else
        rxRadiatorPositions = model.ReceiveAntennaArray.ElementPositions;
    end
    dbar_rx = reshape(rxRadiatorPositions,3,[]) * lambda_0 + model.ReceiveAntennaArray.Position; % size 3-by-U
    
    receiveFieldTerm = complex(zeros(2,N*M,U));
    receiveLocationTerm = complex(zeros(N*M,U),zeros(N*M,U));
    % for each receive antenna:
    for u = 1:U
        
        % get receive field term
        receiveFieldTerm(:,:,u) = wireless.internal.channelmodels.getFieldTerm(model.ReceiveAntennaArray,theta_ZOA,phi_AOA,u);
        
        % get receive location term
        receiveLocationTerm(:,u) = getLocationTerm(rhat_rx,dbar_rx,u,model.CarrierFrequency);
        
    end
    
    % get polarization term, size 4-by-(N*M)
    if (strcmpi(thisClusterType,'LOS'))
        polarizationTerm = getLOSPolarizationMatrix(Phi);
    else
        polarizationTerm = getNLOSPolarizationMatrix(XPR,Phi);
    end

    % establish set of active rays for each cluster
    R = getActiveRays(Phi);
    
    % for each transmit antenna:
    for s = 1:S

        % for each receive antenna:
        for u = 1:U

            % calculate the term of the coefficients combining the transmit
            % and receive antenna fields and the polarization matrix;
            % output is of size (N*M)-by-1
            fieldTerms = ((receiveFieldTerm(1,:,u).*polarizationTerm(1,:) + receiveFieldTerm(2,:,u).*polarizationTerm(2,:)) .* transmitFieldTerm(1,:,s) + ...
                          (receiveFieldTerm(1,:,u).*polarizationTerm(3,:) + receiveFieldTerm(2,:,u).*polarizationTerm(4,:)) .* transmitFieldTerm(2,:,s)).';

            % combine the field terms with the location terms for the
            % current transmit and receive antenna; output is of size
            % (N*M)-by-1. This does not account for Doppler shifts, added
            % in wireless.internal.channelmodels.timeVaryingCDLChannel.
            allTerms = fieldTerms .* receiveLocationTerm(:,u) .* transmitLocationTerm(:,s);

            % zero out any entries corresponding to unused rays for the
            % current set of clusters; output is of size (N*M)-by-1
            allTerms(R(:)==0) = 0;
            Hstatic(1,:,:,s,u) = reshape(allTerms,N,M);

        end
    end

    % calculate number of active rays for each cluster, |R_i|
    modR_i = sum(R,2).';

    scaling = sqrt(P./modR_i);

    % apply final scaling term to H
    % the scaling here includes |R_i| because the desired cluster power
    % |R_i|/M for sub-clusters has already been applied to 'pdp', during
    % its initial construction
    Hstatic = Hstatic .* scaling;

end

% establish set of active rays for each cluster in matrix 'R'
% each row of R corresponds to one cluster, with the columns corresponding
% to the ray number - the 'i'th element of a row is set to 1 if ray R_i is
% active as follows:
% * R_i = 1 for LOS
% * R_i according to Table 7.5-5 for sub-clustered NLOS
% * R_i = 1...20 for other NLOS
% (R_i is only defined in TR 38.901 for sub-clustered NLOS clusters, but
% the variable is given the values above to allow uniform processing for
% the three cases)
function R = getActiveRays(Phi)

    [N,M,~] = size(Phi);
    R = ones(N,M);
    R(Phi(:,:,1)==-Inf) = 0;

end

% TR 38.901 Equation 7.5-22/7.5-28/7.5-29, exponential term which is a
% function of the spherical unit vector 'rhat' for arrival (receive) or
% departure (transmit) azimuth and zenith angles and the location of the
% receive or transmit antenna element 'dbar'
function locationTerm_ant = getLocationTerm(rhat,dbar,ant,freq)

    persistent c0;
    if isempty(c0)
        c0 = physconst('LightSpeed');
    end
    
    % reshape rhat to combine cluster and ray dimensions into a single row
    rhat = reshape(rhat,[3 numel(rhat)/3]);
    
    lambda_0 = c0 / freq;  % carrier wavelength, m
    locationTerm_ant = exp(1i*2*pi*rhat.'*dbar(:,ant)/lambda_0);

end

% TR 38.901 Section 7.5 step 11
% polarization matrix in Equations 7.5-22 and 7.5-28
function polarizationMatrix = getNLOSPolarizationMatrix(XPR,Phi)

    % reshape Phi to combine cluster and ray dimensions into a single row
    Phi = reshape(permute(Phi,[3 1 2]),[4 numel(Phi)/4]);
    
    % get cross polarization power ratio in linear scale, kappa, and
    % reshape it to combine cluster and ray dimensions into a single row
    kappa = 10.^(XPR/10);
    kappa = reshape(permute(kappa,[3 1 2]),[1 numel(kappa)]);
    
    % get polarization matrix in Equations 7.5-22 and 7.5-28
    % note that the 1st dimension of Phi, when interpreted as a 2-by-2
    % matrix, is transposed during this operation (1st dimension element
    % order [1 3 2 4]), as the order of the polarization terms in Phi
    % described by the standard and the MODEL.InitPhase parameter
    % (theta/theta, theta/phi, phi/theta, phi/phi) is effectively row major
    polarizationMatrix = [ones(size(kappa)); sqrt(1./kappa); sqrt(1./kappa); ones(size(kappa))] .* exp(1i*Phi([1 3 2 4],:,:));

end

% TR 38.901 Section 7.5 step 11
% polarization matrix in Equation 7.5-29
function polarizationMatrix = getLOSPolarizationMatrix(Phi)

    polarizationMatrix = [1; 0; 0; -1] .* exp(1i*Phi(:,:,1));
    
end

function theta = wrapZenithAngles(theta)

    theta = mod(theta,360);
    if (~isempty(coder.target) && isvector(theta))
        % force column indexing as code generation 
        % assumes that vectors are columns
        theta(:) = wrap(theta(:));
    else
        theta = wrap(theta); 
    end
    
    function theta = wrap(theta)
        theta(theta>180) = 360 - theta(theta>180);
    end

end

function phi = wrapAzimuthAngles(phi)

    phi =  (mod(phi + 180,360) - 180);

end

function D = getClusterDimensionality(insize,numOutputSignals)

    T = 1;                      % 1 time-domain channel sample
    P = insize(2);              % number of input signals
    R = numOutputSignals;       % number of output signals
    M = wireless.internal.channelmodels.getNumberOfRays(); % Number of rays per cluster    

    D = [T 1 M P R];
    
end

% Get XPR to be a matrix of size numel(clusterTypes)-by-M
function XPR = getXPR(model, clusterTypes)

    % get number of rays
    M = wireless.internal.channelmodels.getNumberOfRays(model);

    % get temporary XPR value
    xpr = model.XPR;

    % Get XPR matrix
    if isscalar(xpr)
        XPR = repmat(xpr,numel(clusterTypes),M);
    else
        LOS = find(strcmpi(clusterTypes,'LOS'));
        splitNLOS = find(strcmpi(clusterTypes,'SubclusteredNLOS'));
        % For configurations with a LOS path (HasLOSCluster=true), the
        % first row of XPR has no effect as it corresponds to the LOS
        % (specular) path for which no cross polarization power ratio is
        % applied. Because the function skips the first row in case of LOS,
        % add here an additional value corresponding to the LOS cluster, if
        % any.
        xprTmp = [xpr(LOS,:); xpr]; % Split the first path into a LOS part and Rayleigh part
        if ~isempty(splitNLOS)
            % duplicate 1st dimension elements in XPR that belong to
            % clusters with sub-clustering
            XPR = wireless.internal.channelmodels.duplicateSubclusterRows(xprTmp,LOS,splitNLOS);
        else
            XPR = xprTmp;
        end
    end
end