%getFieldTerm Field term of antenna element for spherical angles
%
%   Note: This is an internal undocumented function and its API and/or
%   functionality may change in subsequent releases.
% 
%   TR 38.901 Equation 7.5-22/7.5-28/7.5-29 and TR 36.873 Equation
%   7.3-22/7.5-27. Field term (F) which represents the field pattern of a
%   transmit or receive antenna element in the GCS, including any
%   transformation from the LCS due to the antenna orientation.

%   Copyright 2020 The MathWorks, Inc.

%#codegen
function fieldTerm_ant = getFieldTerm(array,theta,phi,ant)

    % reshape theta and phi to combine cluster and ray dimensions into a
    % single row
    theta = theta(:).';
    phi = phi(:).';
    
    % Angles of the orientation of the array
    % alpha: bearing angle of LCS w.r.t. GCS
    % beta: downtilt angle of LCS w.r.t. GCS
    % gamma: slant angle of LCS w.r.t GCS
    % note that the bearing (alpha) and downtilt (beta) angles in the LCS
    % are formed from the combination of the overall array orientation and
    % the orientation of the current element
    alpha = array.Orientation(1) + array.ElementOrientations(1,ant);
    beta  = array.Orientation(2) + array.ElementOrientations(2,ant);
    gamma = array.Orientation(3);
    
    % Angles of the initial orientation of the array. The final orientation
    % of the array is the result of combining the InitialArrayOrientation
    % and Orientation angles.
    alphaI = array.InitialArrayOrientation(1);
    betaI  = array.InitialArrayOrientation(2);
    gammaI = array.InitialArrayOrientation(3);
            
    % TR 38.901 Section 7.1.3 and TR 36.873 Section 5.1.3. R = Rz*Ry*Rx is
    % a 3-by-3 ZYX rotation matrix, which is a sequential set of intrinsic
    % rotations around Z, Y, and X axis used for the transformation from
    % the local coordinate system (LCS) to the global coordinate system
    % (GCS). When the array's initial orientation is not aligned with
    % x-axis, the transformation is corrected using the initial orientation
    % angles.
    R0 = wireless.internal.channelmodels.rotationMatrix([alphaI; betaI; gammaI])';
    R = R0'*wireless.internal.channelmodels.rotationMatrix([alpha; beta; gamma])*R0;
    Rp = R';
    
    sinTheta = sind(theta);
    cosTheta = cosd(theta);
    sinPhi = sind(phi);
    cosPhi = cosd(phi);
    
    rho = [sinTheta.*cosPhi;
           sinTheta.*sinPhi;
           cosTheta];
    
    % TR 38.901 Section 7.1 Equations 7.1-7 and 7.1-8
    % LCS zenith (theta) and azimuth (phi) angles for a given LCS
    % orientation alpha/beta/gamma and Cartesian representation (rhat) of
    % GCS zenith and azimuth angles theta and phi
    theta_p = real(acosd(Rp(3,:)*rho));
    phi_p = atan2d(Rp(2,:)*rho, Rp(1,:)*rho);
    phi_p(theta_p==0) = 0; % Set ambiguous phi to 0 

    cosThetaP = cosd(theta_p);
    % theta vector in the LCS (Equation 7.1-13)
    theta_p_vec = [cosThetaP.*cosd(phi_p);
                   cosThetaP.*sind(phi_p);
                  -sind(theta_p)];            
              
    % theta and phi vectors in the GCS (Equation 7.1-13 and -14)
    phi_vec   = [-sinPhi; 
                  cosPhi; 
                  zeros(size(phi))];
    theta_vec = [cosTheta.*cosPhi;
                 cosTheta.*sinPhi;
                -sinTheta];
        
    % TR 38.901 Section 7.1 Equation 7.1-12
    % Rotation of the spherical basis vectors theta_hat and phi_hat by an
    % angle psi due to the orientation of the LCS w.r.t. the GCS
    psi = atan2d(sum(phi_vec.*(R*theta_p_vec)),sum(theta_vec.*(R*theta_p_vec)));
    
    % TR 38.901 Section 7.3.2
    % Polarization slant angle
    zeta = array.ElementOrientations(3,ant);
    
    % F = [F_theta_prime; F_phi_prime], where F_theta_prime and F_phi_prime
    % are the field components in the direction of the LCS spherical unit
    % vectors theta_hat_prime and phi_hat_prime respectively. For
    % heterogenous arrays, the fields depend on the antenna element.
    F = array.ElementPattern(theta_p,phi_p,zeta,ant);
    
    % TR 38.901 Section 7.1 Equation 7.1-11
    % Rotation of element field pattern from LCS to GCS
    fieldTerm_ant = [F(1,:).*cosd(psi) - F(2,:).*sind(psi); F(1,:).*sind(psi) + F(2,:).*cosd(psi)];

end
