function dopplerTerm = getDopplerTerm(model,rhat_rx,rhat_tx,alph,D,t)
%getDopplerTerm get Doppler term
%
%   Note: This is an internal undocumented function and its API and/or
%   functionality may change in subsequent releases.
%

%   Copyright 2023 The MathWorks, Inc.

%#codegen

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

        % Reshape rhat_rx, rhat_tx, alph and D to combine cluster and ray dimensions into a single
        % row
        s = size(rhat_rx);
        rhat_rx = reshape(rhat_rx,[3 numel(rhat_rx)/3]);
        rhat_tx = reshape(rhat_tx,[3 numel(rhat_tx)/3]);
        alph = alph(:);
        D = D(:);
        
        lambda_0 = c0 / model.CarrierFrequency;    % Carrier wavelength, m
        
        v_rx = model.MaximumDopplerShift(1) * lambda_0; % Rx speed, m/s
        v_tx = model.MaximumDopplerShift(2) * lambda_0; % Tx speed, m/s
        
        theta_v_rx = model.UTDirectionOfTravel(2,1);   % Rx zenith angle of travel, degrees
        phi_v_rx = model.UTDirectionOfTravel(1,1);     % Rx azimuth angle of travel, degrees
        theta_v_tx = model.UTDirectionOfTravel(2,2);   % Tx zenith angle of travel, degrees
        phi_v_tx = model.UTDirectionOfTravel(1,2);     % Tx azimuth angle of travel, degrees
        
        vbar_rx = v_rx * wireless.internal.channelmodels.getRhoHat(phi_v_rx,theta_v_rx);    % Rx velocity
        vbar_tx = v_tx * wireless.internal.channelmodels.getRhoHat(phi_v_tx,theta_v_tx);    % Tx velocity

        % Calculate Doppler term
        dopplerTerm = exp(kron(t,(1i*2*pi*((rhat_rx.'*vbar_rx)+(rhat_tx.'*vbar_tx)+(2*(alph.*D)))/lambda_0).'));

        dopplerTerm = reshape(dopplerTerm,[length(t) s(2:end)]);
        
end