%getCDLInitialPhases TR 38.900 / 38.901 Section 7.5 initial phase setup (step 10)
%
%   Note: This is an internal undocumented function and its API and/or
%   functionality may change in subsequent releases.
%
% Phi is an array of initial phases of size L-by-M-by-4, where L is the
% number of clusters and M is the number of rays. The 3rd dimension
% contains polarization combinations (4 elements) as follows:
% theta/theta, theta/phi, phi/theta, phi/phi

%   Copyright 2016-2021 The MathWorks, Inc.

%#codegen

function Phi = getCDLInitialPhases(model,randomStream,clusterTypes)

    persistent c0;
    if isempty(c0)
        c0 = physconst('LightSpeed');
    end
    
    % If InitPhase is 'Random' replace with '38.900' to give backwards
    % compatibility prior to adding 38.901 support
    if (strcmpi(model.InitPhase,'Random'))
        model.InitPhase = '38.900';
    end
    
    M = wireless.internal.channelmodels.getNumberOfRays(model);
    L = numel(clusterTypes);
    LOS = find(strcmpi(clusterTypes,'LOS'));
    splitNLOS = find(strcmpi(clusterTypes,'SubclusteredNLOS'));
    siz = [L M 4];
    
    if (any(strcmpi(model.InitPhase,{'38.900','38.901','36.873'})))
        % Note that in the case of 38.901, the phases drawn here for the
        % LOS cluster will be ignored, but are still drawn from the random
        % number generator for consistency of the NLOS cluster phases
        % between 38.900 and 38.901
        if (isempty(randomStream))
            Phi = rand(siz);
        else
            Phi = rand(randomStream,siz);
        end
        Phi = (Phi*2*pi) - pi;
    elseif (isscalar(model.InitPhase)) % model.InitPhase is a numeric scalar
        Phi = repmat(double(model.InitPhase),siz);
    else % model.InitPhase is a numeric matrix
        initPhase = double(model.InitPhase);
        % verify phase array size
        if (isempty(coder.target))
            wireless.internal.channelmodels.verify3DArraySize(initPhase,L - numel(splitNLOS)/3*2 - numel(LOS),M,4,'Phase Initialization');
        end

        % For configurations with a LOS path (HasLOSCluster=true), the
        % first row of InitialPhases has no effect as it corresponds to the
        % LOS (specular) path for which no phase initialization angles are
        % applied. Because the function skips the first row in case of LOS,
        % add here an additional value corresponding to the LOS cluster, if
        % any. Specifically, the function manipulates the array to deal
        % with the LOS component and to configure the subsequent channel
        % generation. The array is manipulated as follows.
        %  * Values of -Inf are inserted in phase initialization
        % locations not used for the LOS cluster (rays 2...20, and
        % polarization combinations theta/phi, phi/theta and phi/phi for
        % the first ray).
        %  * Values of -Inf are inserted in rays not used by a particular
        %  subcluster to capture the mapping of sub-clusters to rays in
        %  Table 7.5-5.
        PhiTmp = [zeros(numel(LOS),M,4); initPhase];

        if (~isempty(splitNLOS))
            % duplicate 1st dimension elements in Phi that belong to
            % clusters with sub-clustering
            Phi = wireless.internal.channelmodels.duplicateSubclusterRows(PhiTmp,LOS,splitNLOS);
        else
            Phi = PhiTmp;
        end
    end
    
    if (~isempty(LOS))
        Phi(LOS,2:20,:) = -Inf;
        Phi(LOS,1,2:4) = -Inf;
        if (strcmpi(model.InitPhase,'38.901')) || isnumeric(model.InitPhase)
            % TR 38.901 Equation 7.5-29, phase of exponential term involving d_3D.
            % The phase here is incorporated into the polarization matrix
            % calculated in wireless.internal.channelmodels.CDLChannel.getLOSPolarizationMatrix
            lambda_0 = c0 / model.CarrierFrequency;  % carrier wavelength, m
            d_3D = sqrt(sum((model.TransmitAntennaArray.Position - model.ReceiveAntennaArray.Position).^2));
            Phi(LOS,1) = -2 * pi * d_3D / lambda_0;
        end 
    end
    
    if (~isempty(splitNLOS))
        splitNLOS = reshape(splitNLOS,3,numel(splitNLOS)/3);
        [R_1, R_2, R_3] = wireless.internal.channelmodels.getRayMapping(M);
        Phi(splitNLOS(1,:),[R_2 R_3],:) = -Inf; % [R_2 R_3] = not R_1
        Phi(splitNLOS(2,:),[R_1 R_3],:) = -Inf; % [R_1 R_3] = not R_2
        Phi(splitNLOS(3,:),[R_1 R_2],:) = -Inf; % [R_1 R_2] = not R_3
    end

end