%timeVaryingCDLChannel Add fading (Doppler) to TR 38.901 and TR 36.873 CDL channel
%
%   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 [H,sampleTimes] = timeVaryingCDLChannel(model,staticChannelInfo,alph,scattererSpeeds,theTime,insize,outputtype)  

    % get dimensionality information for a single cluster
    [D,t] = getClusterDimensionality(model,insize,theTime);

    % get number of clusters N and set up the corresponding dimension of
    % the channel matrix output 'H' ('H' will have dimension sizes 'D')
    staticPathGains = staticChannelInfo.StaticPathGains;
    N = size(staticPathGains,2);
    D(2) = N;

    % create empty channel matrix 
    H = complex(zeros(D,outputtype),zeros(D,outputtype));

    % get Doppler term according to TR 38.901 Step 11 Eqs. (7.5-28/29)
    departureUnitVectors = staticChannelInfo.DepartureUnitVectors;
    arrivalUnitVectors = staticChannelInfo.ArrivalUnitVectors;
    dopplerTerm = wireless.internal.channelmodels.getDopplerTerm(model,arrivalUnitVectors,departureUnitVectors,alph,scattererSpeeds,t); % size T_H-by-N-by-M

    % Introduce the Doppler effect of each ray in the static (cluster) path
    % gains by adding the time-varying ray complex amplitudes within each
    % cluster. The size of this cluster response H is D. Note the implicit
    % cast to the type of H, defined by outputtype.

    % for each receive antenna:
    U = D(4);
    for u = 1:U
        H(:,:,:,u) = permute(sum(dopplerTerm .*staticPathGains(:,:,:,:,u),3),[1 2 4 3]);
    end

    % time instants when the channel response is sampled
    sampleTimes = getSampleTimes(model,insize,theTime);

end

function [D,t] = getClusterDimensionality(model,insize,theTime)

    t = getSampleTimes(model,insize,theTime); % sample times
    
    T = numel(t);               % number of time-domain channel samples
    P = insize(2);              % number of input signals
    R = model.NumOutputSignals; % number of output signals
    
    D = [T 1 P R];
    
end

function [t,F_cg] = getSampleTimes(model,insize,theTime)

    persistent c0;
    if isempty(c0)
        c0 = physconst('LightSpeed');
    end

    if (isfield(model,'SampleTimes'))
        
        % if filtering is off, or filtering is on and SampleTimes are
        % present, use SampleTimes to define time points
        t = model.SampleTimes(:);
        F_cg = [];
        
    else
        
        if (model.SampleDensity == Inf)
            
            % set coefficient generation sampling rate F_cg equal to input
            % sampling rate (SampleRate is mandatory if filtering is on)
            % and set the number of samples to the input waveform length
            F_cg = model.SampleRate;
            T = insize(1);
            
        else

            % set coefficient generation sampling rate F_cg according to
            % Doppler frequency and SampleDensity parameter
            lambda_0 = c0 / model.CarrierFrequency;    % carrier wavelength, m

            F_cg = (sum(model.MaximumDopplerShift) + (2 * (model.MaximumScattererSpeed/lambda_0))) * 2 * model.SampleDensity;

            % calculate the number of time points at which to take channel
            % samples, including at least an extra half sample period
            % beyond the end of the waveform to allow for zero order hold
            % "balancing" (see zohChannelMatrixTimes in getPathGainIndices)
            T = ceil(insize(1) * F_cg / model.SampleRate + 0.5);
            
        end

        % Limit F_cg to sample rate and T to input waveform length
        F_cg = min(F_cg,model.SampleRate);
        T = min(T,insize(1));
        
        % calculate vector of time points at sampling rate F_cg
        % (with special case for F_cg=0 for zero Doppler (T=1,t=0))
        if (T==1)
            t = 0;
        else
            t = (0:(T-1)).'/F_cg;
        end
        
        % add coefficients generation time offset
        t = t + theTime;

    end
    
end