function [pl, phase] = raypl(ray, varargin)
%RAYPL Path loss and phase shift for RF propagation ray
%   [PL, PHASE] = raypl(RAY) returns the path loss, PL, in dB and phase
%   shift, PHASE, in radians for one comm.Ray object, RAY. The carrier
%   frequency to calculate the path loss and phase shift is specified by
%   the Frequency property of RAY. Phase uses exp(-iwt) convention.
%
%   The RAY object must have the PathSpecification property set to
%   'Locations' and the Type field of all structures assigned to the
%   Interactions property set to 'Reflection'. The output PL is a double
%   precision, non-negative scalar. The output PHASE is a double precision,
%   scalar between -pi and pi.
%
%   The path loss and phase shift calculations are for unpolarized antennas
%   by default. In that case, geometric coupling between horizontal and
%   vertical polarizations are not taken into account. To turn on
%   polarizations, use the Name-Value pairs to specify both 
%   TransmitterPolarization and ReceiverPolarization to a polarization type
%   other than 'none' or a Jones vector (see below).
%   
%   [PL, PHASE] = raypl(RAY,Name,Value) calculates the path loss and phase
%   shift with additional options specified by one or more Name-Value
%   pairs.
%
%   raypl Name-Value pairs:
%
%   ReflectionMaterials - Reflection materials for a non-line-of-sight
%       (NLOS) ray, specified as a string scalar, 1-by-NI string vector,
%       character vector, 1-by-NI cell array of character vectors, 2-by-1
%       numeric vector, or 2-by-NI numeric matrix, where NI equals
%       RAY.NumInteractions. When it is a string scalar or vector, each
%       reflection material must be one of "brick" | "ceiling-board" |
%       "chipboard" | "concrete" | "floorboard" | "glass" | "loam" |
%       "marble" | "metal" | "perfect-reflector" | "plasterboard" |
%       "plywood" | "vegetation" | "water" | "wood". When it is a character
%       vector or cell array of character vectors, the same set of
%       reflection materials as above applies. When it is a numeric vector
%       or matrix, the two rows represent the relative permittivity and
%       conductivity of the corresponding reflection material,
%       respectively. When it has only one column, the material applies to
%       all the reflections. The default value is "concrete".
%
%   TransmitterPolarization - Transmit antenna polarization type, specified 
%       as one of "H" | "V" | "LHCP" | "RHCP" | "none" or as a 2-by-1
%       (normalized) Jones vector in the form of [H; V]. The default value
%       is "none". Horizontal & vertical polarization are sometimes
%       referred to as E_phi & E_theta in literature or E_azimuth &
%       E_elevation in MATLAB, respectively.
%
%   ReceiverPolarization - Receive antenna polarization type, specified as 
%       one of "H" | "V" | "LHCP" | "RHCP" | "none" or as a 2-by-1
%       (normalized) Jones vector in the form of [H; V]. The default value
%       is "none". Horizontal & vertical polarization are sometimes
%       referred to as E_phi & E_theta in literature or E_azimuth &
%       E_elevation in MATLAB, respectively.
%
%   TransmitterAxes - Orientation axes of transmit antenna, specified
%       as a 3-by-3 unitary matrix indicating the rotation from the
%       transmitter local coordinate system (LCS) into the global
%       coordinate system (GCS), which is the local East-North-Up (ENU)
%       coordinate system at transmitter when CoordinateSystem is set to
%       "Geographic" in the ray object. The default value is a 3-by-3
%       identity matrix.
%
%   ReceiverAxes - Orientation axes of receiver antenna, specified
%       as a 3-by-3 unitary matrix indicating the rotation from the
%       receiver local coordinate system (LCS) into the global coordinate
%       system (GCS), which is the local East-North-Up (ENU) coordinate
%       system at receiver when CoordinateSystem is set to "Geographic" in
%       the ray object. The default value is a 3-by-3 identity matrix.
%
%   % Example: Re-evaluate path loss and phase shift for a ray with
%   % different frequency and reflection materials.
%   
%   % Launch Site Viewer with buildings in Hong Kong
%   viewer = siteviewer("Buildings", "hongkong.osm");
%
%   % Specify transmitter and receiver sites
%   tx = txsite( ...
%       "Latitude", 22.2789, ...
%       "Longitude", 114.1625, ...
%       "AntennaHeight", 10, ...
%       "TransmitterPower", 5, ...
%       "TransmitterFrequency", 28e9);
%   rx = rxsite( ...
%       "Latitude", 22.2799, ...
%       "Longitude", 114.1617, ...
%       "AntennaHeight", 1);
%
%   % Specify a ray tracing propagation model using SBR method with up to 2
%   % order reflections.
%   pm = propagationModel("raytracing", ...
%       "Method", "sbr", ...
%       "MaxNumReflections", 2);
%   
%   % Perform ray tracing between the sites
%   rays = raytrace(tx, rx, pm);
%   
%   % Find the first ray with 2 reflections
%   ray = rays{1}(find([rays{1}.NumInteractions] == 2, 1));
%   
%   % Plot the ray and observe it reflects on buildings twice
%   plot(ray);
%   
%   % The material is concrete by default for all the buildings. We can
%   % change the material to metal for the second reflection and
%   % re-evaluate path loss and phase shift, using the raypl function.
%   [ray.PathLoss, ray.PhaseShift] = raypl(ray, ...
%       "ReflectionMaterials", ["concrete", "metal"]);
%
%   % Plot the ray and observe the slight color change
%   plot(ray);
%
%   % We can also re-evaluate the path loss and phase shift at a different
%   % frequency using the raypl function. 
%   ray.Frequency = 2e9;
%   [ray.PathLoss, ray.PhaseShift] = raypl(ray, ...
%       "ReflectionMaterials", ["concrete", "metal"]);
%   
%   % Plot the ray again and observe the obvious color change
%   plot(ray);
%   
%   See also comm.Ray, raytrace, earthSurfacePermittivity,
%   buildingMaterialPermittivity.

% Copyright 2019-2022 The MathWorks, Inc.
%#codegen

% Parse inputs
if isempty(coder.target)
    p = inputParser;
    % Use hidden ReflectionMaterials property if the N-V pair not specified
    p.addParameter("ReflectionMaterials", {'concrete'});
    p.addParameter('TransmitterPolarization', "none");
    p.addParameter('ReceiverPolarization', "none");
    p.addParameter('TransmitterAxes', eye(3));
    p.addParameter('ReceiverAxes', eye(3));
    p.addParameter('ValidateInputs', true);
    p.parse(varargin{:});
    
    mtlsChar = convertStringsToChars(p.Results.ReflectionMaterials);
    txPol = convertCharsToStrings(p.Results.TransmitterPolarization);
    rxPol = convertCharsToStrings(p.Results.ReceiverPolarization);
    txAxes = p.Results.TransmitterAxes;
    rxAxes = p.Results.ReceiverAxes;
    validateInputs = p.Results.ValidateInputs;
else
    pvPairs = struct(...
        'ReflectionMaterials', uint32(0), ...
        'TransmitterPolarization', uint32(0), ...
        'ReceiverPolarization', uint32(0), ...
        'TransmitterAxes', uint32(0), ...
        'ReceiverAxes', uint32(0), ...
        'ValidateInputs', uint32(0));
    popts = struct( ...
        'CaseSensitivity', false, ...
        'StructExpand', true, ...
        'PartialMatching', true);
    pStruct = coder.internal.parseParameterInputs( ...
        pvPairs, popts, varargin{:});
    mtlsChar = convertStringsToChars( ...
        coder.internal.getParameterValue( ...
        pStruct.ReflectionMaterials, {'concrete'}, varargin{:}));
    txPol = convertCharsToStrings( ...
        coder.internal.getParameterValue( ...
        pStruct.TransmitterPolarization, "none", varargin{:}));
    rxPol = convertCharsToStrings( ...
        coder.internal.getParameterValue( ...
        pStruct.ReceiverPolarization, "none", varargin{:}));
    txAxes = coder.internal.getParameterValue( ...
        pStruct.TransmitterAxes, eye(3), varargin{:});
    rxAxes = coder.internal.getParameterValue( ...
        pStruct.ReceiverAxes, eye(3), varargin{:});
    validateInputs = coder.internal.getParameterValue( ...
        pStruct.ValidateInputs, true, varargin{:});
end

% Validate ray object
if validateInputs
    validateattributes(ray, {'comm.Ray'}, {'scalar'}, ...
        'raypl', 'ray input');
    validateConfig(ray);
    coder.internal.errorIf( ...
        strcmp(ray.PathSpecification, 'Delay and angles'), ...
        'shared_channel:raypl:NotInLocSpec');
    coder.internal.errorIf( ...
        sum('R' == [ray.Interactions.Type], 2) ~= size(ray.Interactions, 2), ...
        'shared_channel:raypl:NotAllReflection');    
end

% Get carrier frequency
fc = ray.Frequency;

% Validate material input for NLOS
if ischar(mtlsChar)
    mtls = {mtlsChar};
else
    mtls = mtlsChar;
end

if validateInputs && ~ray.LineOfSight
    if iscell(mtls)
        % Force mtls to be a homogeneous cell array for code generation
        coder.varsize('mtls', [1 length(mtls)], [0 0]);

        validateattributes(mtls, {'cell'}, {'row'}, ...
            'raypl', 'material input');
        supportedMtls = {'plasterboard','ceiling-board','chipboard', ...
            'floorboard','concrete','brick','wood','glass', 'marble', 'plywood', ...
            'metal','water','vegetation','loam','perfect-reflector'};
        for i = 1:length(mtls)
            coder.internal.errorIf(~any(strcmpi(mtls{i}, supportedMtls)), ...
                'shared_channel:raypl:MaterialNameNotSupported');
        end
    else
        validateattributes(mtls, {'numeric'}, ...
            {'real','nrows',2}, ...
            'raypl', 'material input');
        mtls = double(mtls);
    end
    
    % Check number of materials matches number of reflections
    coder.internal.errorIf((size(mtls, 2) > 1) && ...
        (size(mtls, 2) ~= ray.NumInteractions), ...
        'shared_channel:raypl:NumRefNumMtlMismatch');
end

% Validate Tx orientation axes
if validateInputs 
    validateAxes(txAxes, 'transmitter');
end

% Validate Rx orientation axes
if validateInputs 
    validateAxes(rxAxes, 'receiver');
end

% Validate Tx polarization input and derive Jones vector from it
[txJV,isTxPol] = getJonesVector(txPol, 'transmitter', validateInputs);

% Validate Rx polarization input and derive Jones vector from it
[rxJV,isRxPol] = getJonesVector(rxPol, 'receiver', validateInputs);

% For unpolarized antennas (e.g. "isotropic"), ignore txAxes' and rxAxes'
% rotation effect on local polarization
if ~isTxPol
    txAxes = eye(3);
end
if ~isRxPol
    rxAxes = eye(3);
end

% Get polarization matrix
polMtx = getPolMatrix(ray, mtls, txAxes, rxAxes);

% Calculate total loss due to interactions, depending on polarization
if ~isTxPol || ~isRxPol
    % To handle non-/un-polarized antenna(s), which by definition have
    % random/unknown phase coherency between polarizations, treat both
    % antennas as unpolarized with no polarization mismatch or phase
    % coherency.
    intactLoss = sum(abs(polMtx),'all')/2; % = sqrt(2)/2*[1;1]' * abs(polMtx) * sqrt(2)/2*[1;1]
else
    intactLoss = rxJV' * polMtx * txJV;
end

% Calculate path loss in dB from interactions and free-space
lambda = 299792458/ray.Frequency; % wavelength in free-space
pl = -20*log10(abs((intactLoss))) + ...     % Reflection loss
      fspl(ray.PropagationDistance,lambda); % Free space loss = fspl
pl = max(0, pl);

% Calculate phase change in radians from reflections and free space
numCycles = ray.PropagationDelay*fc;
phase = mod(angle(intactLoss) ...   % Reflection phase change, exp(-iwt) convention
       + 2*pi*numCycles, 2*pi);     % Free space phase change, exp(-iwt) convention
end

function validateAxes(antennaAxes, antennaStr)
% Validate tx/rx orientation axes input

validateattributes(antennaAxes, {'double'}, ...
    {'real','finite','size',[3 3]}, ...
    'raypl', [antennaStr, ' orientation axes input']);
coder.internal.errorIf(...
    max(max(abs(antennaAxes'*antennaAxes-eye(3)))) > sqrt(eps), ...
    'shared_channel:raypl:AxesNotUnitary', antennaStr);

end

function [JV,isPol] = getJonesVector(pol, antennaStr, validateInputs)
% Validate tx/rx polarization input and derive Jones vector in the form of
% [H; V]

if validateInputs
    validateattributes(pol, {'string','numeric'}, {}, ...
        'raypl', [antennaStr, ' polarization input']);
end

if isstring(pol)
    if validateInputs
        validateattributes(pol, {'string'}, {'scalar'}, ...
            'raypl', [antennaStr, ' polarization input']);
        coder.internal.errorIf( ...
            ~any(strcmpi(pol, {'H','V','LHCP','RHCP','none'})), ...
            'shared_channel:raypl:PolTypeNotSupported', antennaStr);
    end
    switch pol  % This matches the top-of-file & online documentation
        case 'H'
            JV = [1; 0];
            isPol = true;
        case 'V'
            JV = [0; 1];
            isPol = true;
        case 'LHCP'
            JV = 1/sqrt(2)*[1; 1i]; % exp(-iwt) transmit convention
            isPol = true;
        case 'RHCP'
            JV = 1/sqrt(2)*[1; -1i]; % exp(-iwt) transmit convention
            isPol = true;
        otherwise % 'none' = randomly polarized
            JV = 1/sqrt(2)*[1; 1];
            isPol = false;
    end
else % Numeric
    if validateInputs
        validateattributes(pol, {'double'}, {'finite','size',[2,1]}, ...
            'raypl', [antennaStr, ' polarization input']);
        coder.internal.errorIf(abs(norm(pol) - 1) > sqrt(eps), ...
            'shared_channel:raypl:JonesVectorNotNormalized', antennaStr);
    end
    JV = pol;
    isPol = true;
end

end

% [EOF]