function emissOut = radarChannel(emissIn, platsIn, varargin)
%radarChannel  Atmospheric propagation and reflections of radar signals
%   EMPROP = radarChannel(EMTX, PLATS) returns radar emissions, EMPROP,
%   which are the combination of the input radar emissions, EMTX, as well
%   as the reflections of each of the input emissions off each of the
%   platforms, PLATS. EMTX and EMPROP are arrays of radarEmission objects.
%   
%   If EMTX is a cell array of radarEmission objects, then EMPROP is a
%   cell array of radarEmission objects. If EMTX is an array of
%   radarEmission objects, then EMPROP is an array of radarEmission
%   objects.
%
%   PLATS is cell array of platform objects provided by the Platforms
%   property on the scenario or an array of platform structures with fields
%   as listed <a href="matlab:help('fusion.internal.interfaces.DataStructures/platformStruct')">here</a>.
%
%   If an array of platform structures is used, the PlatformID and Position
%   fields are required. All of the remaining fields will be assigned
%   default values if not specified.
%
%   radarChannel(... 'Name', Value, ...) specifies additional
%   name-value pair arguments described below:
%
%   'HasOcclusion'  A real, finite, scalar value. When true or non-zero,
%                   the channel models occlusion of the input radar
%                   emissions, EMTX, along the propagation path to each
%                   platform, PLATS. If another platform blocks the
%                   propagation to a reflecting platform, a reflection is
%                   not returned for that reflecting platform. When set to
%                   false or zero, occlusion is not modeled between the
%                   input radar emissions and each reflecting platform. The
%                   default is true.
%
%   EXAMPLE 1: Reflect a radar emission off a platform.
%
%   % Define a platform for the emitter.
%   plat = struct('PlatformID', 1, 'Position', [0 0 0], ...
%                 'Signatures', {rcsSignature});
%
%   % Create a radar emission.
%   emTx = radarEmission('EmitterIndex', 1, ...
%                        'PlatformID', plat.PlatformID, ...
%                        'OriginPosition', plat.Position);
%   
%   % Define a target to reflect the emission.
%   tgt = struct('PlatformID', 2, 'Position', [10 0 0], ...
%                'Signatures', {rcsSignature});
%
%   % Reflect the emission off the target.
%   emProp = radarChannel(emTx, [plat tgt]);
%
%   % The first propagated emission is the same as the emitter's emission.
%   emProp(1)
%
%   % The second emission is the reflected emission. Its PlatformID is set
%   % to the PlatformID of the target from which it was reflected.
%   emProp(2)
%
%   EXAMPLE 2: Reflect a radar emission using radarEmitter.
%
%   % Create a radarEmitter.
%   emitter = radarEmitter(1);
%
%   % Define a platform for the emitter.
%   plat = struct('PlatformID', 1, 'Position', [0 0 0], ...
%                 'Signatures', {rcsSignature});
%
%   % Add a target platform to reflect the emitted signal.
%   tgt = struct('PlatformID', 2, 'Position', [30 0 0], ...
%                'Signatures', {rcsSignature});
%
%   % Emit the signal.
%   simulationTime = 0;
%   emTx = step(emitter, plat, simulationTime)
%
%   % Reflect the signal off the target platform.
%   emProp = radarChannel(emTx, [plat tgt])
%
%   See also: radarEmitter, radarEmission.

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

% Parse name-value pairs
defaultPairs = struct('DefaultHasOcclusion',true,...
    'DefaultReferenceFrame','NED');
mcosParser = fusion.internal.mixin.MCOSConstructorParser;
value = parseConstructorInputs(mcosParser, defaultPairs, varargin{:});

% Check that the name-value pair values are valid
validateattributes(value{1}, {'double','single','logical'}, ...
    {'real','finite','scalar'}, ...
    mfilename, 'HasOcclusion');
hasOcclusion = logical(value{1});
reference = value{2};

% Parse emissions and platforms
emiss = fusion.internal.remotesensors.EmissionUtilities.parseEmissions(emissIn,reference,mfilename);
plats = fusion.internal.remotesensors.EmissionUtilities.parsePlatforms(platsIn,reference,mfilename);

% Are any emissions or platforms in scenario?
if isempty(plats) || fusion.internal.remotesensors.CommonUtilities.isIndexedEmpty(emiss)
    emissOut = emiss;
    return
end

numEmiss = numel(emiss);
numPlats = numel(plats);

% Identify platforms defined as point targets. If all targets are point
% targets, then occlusion can be ignored.
isPtTgts = matlabshared.tracking.internal.fusion.targetTruth.isPointTarget(plats);
hasOcclusion = hasOcclusion && ~all(isPtTgts);

% Generate binary space partition tree for occlusion calculations
fv = matlabshared.tracking.internal.fusion.targetTruth.cuboids(plats);
cubes = matlabshared.tracking.internal.fusion.compgeo.faceVertexToFaceMatrix(fv);
bsptree = matlabshared.tracking.internal.fusion.compgeo.BSPTree(cubes);

% The receiver's location is unknown, so all possible combinations are
% considered. This requires calculating the signal propagation on both the
% receive and transmit path:
%   1. Receive path:  transmitter -> platform's signature
%   2. Transmit path: platform's reflected signal -> a different platform
% Where, here, the "reflecting" platform is treated as an "isotropic
% antenna" that re-emits the reflected emission. The receive and transmit
% gain of the reflecting platform's "antenna" is determined by the
% reflecting platform's radar cross section.

% Receive path: signal to be reflected received at platform from a transmitter
isIlluminated = false(numPlats,numEmiss);
for iEmiss = 1:numEmiss
    % Figure out which platforms are illuminated by each transmitted
    % emission
    thisEmiss = fusion.internal.remotesensors.CommonUtilities.getIndexed(emiss,iEmiss);

    for iThisPlat = 1:numPlats
        thisPlat = fusion.internal.remotesensors.CommonUtilities.getIndexed(plats,iThisPlat);

        % Skip if emission originated from this platform
        if thisPlat.PlatformID == thisEmiss.PlatformID
            continue
        end

        % This method checks:
        % - For point targets:    if the point lies inside of this
        %                         emission's field of view
        % - For extended targets: if any portion of any face on the
        %                         extended target lies inside of the field
        %                         of view
        inFoV = fusion.internal.remotesensors.EmissionUtilities.isPlatformInFieldOfView(thisEmiss,thisPlat,isPtTgts(iThisPlat));

        if hasOcclusion && inFoV(1)
            % The target is inside of the emission's field of view, but is
            % it occluded by other targets?
            isOcc = isOccluded(bsptree,thisEmiss.OriginPosition,iThisPlat);
            isIlluminated(iThisPlat,iEmiss) = ~isOcc;
        else
            isIlluminated(iThisPlat,iEmiss) = inFoV(1);
        end
    end
end

% Are any platforms illuminated?
numRxEmiss = sum(isIlluminated(:));
if numRxEmiss<1
    emissOut = emiss;
    return % no
end

% Allocate output
maxNumRefl = numRxEmiss*numPlats;

% Output format follows input format
if iscell(emiss)
    thisEmiss = {emiss{1}};
else
    thisEmiss = emiss(1);
end
reflections = repmat(thisEmiss,maxNumRefl,1);

iRefl = 1; % Next index in allocated reflections array to add a reflection to
for iEmiss = 1:numel(emiss)
    thisEmiss = fusion.internal.remotesensors.CommonUtilities.getIndexed(emiss,iEmiss);
    
    % Find the platforms illuminated by this emission
    iIllumPlats = find(isIlluminated(:,iEmiss));
    for iPlat = 1:numel(iIllumPlats)
        iReflPlat = iIllumPlats(iPlat);

        % Reflect the emission by re-emitting the emission from the
        % reflecting platform to all of the other platforms visible from
        % the reflecting platform's location
        [reflections,iRefl] = reflectEmission(reflections,iRefl,thisEmiss,iReflPlat,plats,isPtTgts,cubes,hasOcclusion,bsptree);
    end
end

numRefl = iRefl-1;
emissOut = assembleOutput(emiss,reflections,numRefl);
end

function [reflections,iRefl] = reflectEmission(reflections,iRefl,thisEmiss,iReflPlat,plats,isPtTgts,cubes,hasOcclusion,bsptree)
% Adds reflections for this emission from every platform that is
% illuminated by this emission to the allocated reflections array. iRefl is
% the index where the next reflection should be added to the reflections
% array.

% Get the platform that is reflecting this emission
reflPlat = fusion.internal.remotesensors.CommonUtilities.getIndexed(plats,iReflPlat);
reflPlatCube = cubes(:,:,:,iReflPlat);

% Find the faces that are visible on the reflecting platform from the
% emitter
if isPtTgts(iReflPlat)
    isFacing = ~isnan(reflPlatCube(1,:));
    iReflFaces = find(isFacing(:),1);
else
    thisCube = bsxfun(@minus,reflPlatCube,thisEmiss.OriginPosition(:));
    isFacing = matlabshared.tracking.internal.fusion.compgeo.isFacing(thisCube);
    iReflFaces = find(isFacing(:)); % convert to column to eliminate run-time warning in codegen about row/column conversion for find
end

% Remove faces on the reflecting platform that are not visible from the
% position of the emitter
if hasOcclusion
    isOcc = false(size(iReflFaces));
    for m = 1:numel(iReflFaces)
        isOcc(m) = isOccluded(bsptree,thisEmiss.OriginPosition,iReflPlat,iReflFaces(m));
    end
    if all(isOcc)
        return
    end
    iReflFaces = iReflFaces(~isOcc);
end

sol = fusion.internal.Constants.SpeedOfLight;
for iNextPlat = 1:numel(plats)
    % Next platform to which the emission is reflected
    nextPlat = fusion.internal.remotesensors.CommonUtilities.getIndexed(plats,iNextPlat);
    
    % Is platform towards which the reflection is transmitted the same
    % platform that received the emission from the transmitter?
    if nextPlat.PlatformID==reflPlat.PlatformID
        continue % Yes, skip to next platform
    end
    
    % Receive reflected emission at this next platform
    iReflFacesNext = iReflFaces;
    if hasOcclusion
        % Remove faces on the reflecting platform that are not visible from
        % the position of the next platform
        nextPos = nextPlat.Position;
        isOcc = false(size(iReflFacesNext));
        for m = 1:numel(iReflFacesNext)
            isOcc(m) = isOccluded(bsptree,nextPos,iReflPlat,iReflFacesNext(m));
        end
        if all(isOcc)
            continue
        end
        iReflFacesNext = iReflFacesNext(~isOcc);
    end
    
    % Use a point on the surface of the visible faces on the reflecting
    % platform as the point of reflection
    [reflPos,reflVel] = reflectedPosAndVel(reflPlat,isPtTgts(iReflPlat),reflPlatCube,iReflFacesNext);
    [rgIn,rrIn,urIn] = receiveEmission(thisEmiss,reflPos,reflVel);

    % Free space propagation loss from transmitter to receiver
    % - For extended objects, the loss is computed as the distance from
    %   the transmitter to the platform's center of rotation (e.g.,
    %   origin offset).
    lambda = sol/thisEmiss.CenterFrequency;
    propLoss = (4*pi.*rgIn/lambda).^2;
    propLossdB = fusion.internal.units.pow2db(propLoss);

    % Re-emit received energy towards next platform
    nextPos = nextPlat.Position(:);
    [rcsdB,orient] = transmitEmission(thisEmiss,urIn,reflPlat,reflPos,nextPos);

    % Update fields for the reflection to the receiving platform
    thisRefl = thisEmiss;
    thisRefl.PlatformID = reflPlat.PlatformID;
    thisRefl.OriginPosition(:) = reflPos;
    thisRefl.OriginVelocity(:) = reflVel;
    thisRefl.Orientation(:) = orient;
    thisRefl.FieldOfView(:) = [180, 180];
    thisRefl.EIRP = thisRefl.EIRP-propLossdB;
    thisRefl.RCS = thisRefl.RCS+rcsdB;
    thisRefl.PropagationRange = thisRefl.PropagationRange+rgIn;
    thisRefl.PropagationRangeRate = thisRefl.PropagationRangeRate+rrIn;
        
    if ~isPtTgts(iReflPlat)
        % Save off the faces on the reflecting platform that were visible
        % from this next platform
        thisRefl.IlluminatedSurface(:,1:size(reflPlatCube,2),1:numel(iReflFacesNext)) = reflPlatCube(:,:,iReflFacesNext);
    end

    reflections = fusion.internal.remotesensors.CommonUtilities.setIndexed(reflections,iRefl,thisRefl);
    iRefl = iRefl+1;
end
end

function [reflPos,reflVel] = reflectedPosAndVel(reflPlat,isPtTgt,illumSurf,iVisTxRx)
% Returns the origin position and velocity of the reflected emission. For
% point objects, this is simply the position and velocity of the point
% object reflecting the emission. For extended objects, the origin position
% is the centroid of the object's illuminated surface.
if isPtTgt
    reflPos = reflPlat.Position(:);
else
    % Compute geometric center of illuminated surfaces
    if isempty(iVisTxRx)
        reflPos = NaN(3,1,'like',reflPlat.Position);
    else
        reflPos = zeros(3,1,'like',reflPlat.Position);
        totArea = 0;
        for m = 1:numel(iVisTxRx)
            thisFace = illumSurf(:,:,iVisTxRx(m));
            isGd = ~isnan(thisFace(1,:));
            if sum(isGd)>=3 % Need at least 3 points to have a face
                [thisRot, thisOffset] = matlabshared.tracking.internal.fusion.compgeo.rotateIntoFace(thisFace(:,isGd));
                thisFaceRot = thisRot'*bsxfun(@minus,thisFace,thisOffset);
                thisArea = polyarea(thisFaceRot(1,isGd),thisFaceRot(2,isGd));
                reflPos = reflPos+thisArea*mean(thisFace(:,isGd),2);
                totArea = totArea+thisArea;
            end
        end
        reflPos = reflPos/totArea;
    end
end
reflVel = reflPlat.Velocity(:);
end

function emissOut = assembleOutput(emiss,reflections,numRefl)
% Concatenate input and reflected emissions
% - If input was a cell array then output is also a cell array, otherwise
%   it is an array of emission objects.

numIn = numel(emiss);
usesCell = iscell(emiss);
if usesCell
    % Allocate space
    emissOut = repmat({emiss{1}},numIn+numRefl,1);
    
    % Add inputs to the output
    for m = 1:numIn
        emissOut{m} = emiss{m};
    end
    
    % Add reflections to the output
    for m = numIn+1:numel(emissOut)
        emissOut{m} = reflections{m-numIn};
    end
else
    % Allocate space
    emissOut = repmat(emiss(1),numIn+numRefl,1);

    % Add inputs to the output
    for m = 1:numIn
        emissOut(m) = emiss(m);
    end
    
    % Add reflections to the output
    for m = numIn+1:numel(emissOut)
        emissOut(m) = reflections(m-numIn);
    end
end
end

function [rcsdB, orientOut] = transmitEmission(emiss,urIn,reflPlat,reflPos,nextPos)
% Compute the bistatic RCS to account for the receive gain and the transmit
% gain of the platform reflecting the emission and the direction of the
% reflected emission's propagation to the next platform.

vecOut = nextPos-reflPos;

% Convert to spherical coordinates (az,el)
[azOut,elOut,rgOut] = cart2sph(vecOut(1),vecOut(2),vecOut(3));
azOut = rad2deg(azOut);
elOut = rad2deg(elOut);
azOut = fusion.internal.units.interval(azOut,[-180 180]);
elOut = fusion.internal.units.interval(elOut,[-180 180]);

% Add in atmospheric refraction
if rgOut > 0
    rgTrue = rgOut;
    elTrue = elOut;
    elOut = fusion.internal.remotesensors. ...
        SingleExponentialRefraction.modifiedDepressionModel(rgTrue,elTrue);
end
[x,y,z] = sph2cart(deg2rad(azOut),deg2rad(elOut),1);
urOut = [x(:) y(:) z(:)]';

% Compute RCS including refraction effects
rcsdB = getBistaticRCS(reflPlat,urIn,urOut,emiss.CenterFrequency);

% Direction of emission propagation to the next platform as an orientation
% matrix.
orient = lookOrientation(urOut);

% Use same format for orientation in output signals as was used
% in input signals
if isa(emiss.Orientation,'quaternion')
    orientOut = quaternion(orient,'rotmat','frame');
else
    orientOut = orient;
end
end

function orient = lookOrientation(vLook)
% Returns the look direction as an orientation matrix.
u1 = vLook/norm(vLook);
if abs(u1(1))>eps || abs(u1(2))>eps
    v2 = [-u1(2);u1(1);0];
    u2 = v2/norm(v2);
    u3 = cross(u1,u2);
else
    v3 = [-u1(3);0;u1(1)];
    u3 = v3/norm(v3);
    u2 = cross(u3,u1);
end
orient = [u1 u2 u3]';
end

function [rgIn, rrIn, urIn] = receiveEmission(emiss,posIn,velIn)
% Compute the propagation range, range-rate, and direction at which the
% emission is received at the platform from which it will be reflected.

emissPos = emiss.OriginPosition(:);
emissVel = emiss.OriginVelocity(:);
emissOrient = fusion.internal.remotesensors. ...
    CommonUtilities.getRotmat(emiss.Orientation);

vecIn = posIn-emissPos;

if emiss.PropagationRange==0
    % Rotate receiving platform into the transmitter's frame
    vecTx = emissOrient*vecIn;
    
    % Convert to spherical coordinates (az,el)
    [azTx,elTx,rgTx] = cart2sph(vecTx(1),vecTx(2),vecTx(3));
    azTx = rad2deg(azTx);
    elTx = rad2deg(elTx);
    azTx = fusion.internal.units.interval(azTx,[-180 180]);
    elTx = fusion.internal.units.interval(elTx,[-180 180]);
    
    % Add in atmospheric refraction
    rgIn = rgTx;
    if rgIn > 0
        rgTrue = rgTx;
        elTrue = elTx;
        rgIn = fusion.internal.remotesensors. ...
            SingleExponentialRefraction.modifiedRangeModel(rgTrue,elTrue);
        elTx = fusion.internal.remotesensors. ...
            SingleExponentialRefraction.modifiedDepressionModel(rgTrue,elTrue);
    end
    
    % Compute range-rate and RCS including refraction effects
    [x,y,z] = sph2cart(deg2rad(azTx),deg2rad(elTx),1);
    urIn = [x(:) y(:) z(:)]';
    urIn = emissOrient'*urIn;
else
    rgIn = norm(vecIn);
    urIn = vecIn/rgIn;
end

rrIn = urIn'*(velIn-emissVel);
end

function rcsdB = getBistaticRCS(platIn,vecIn,vecOut,freq)
% Returns the RCS of the platform using the monostatic-bistatic
% equivalence theorem (MBET) in decibels.
%
% plat:     Pose for platform whose RCS is being calculated
% vecIn:    Vector defining direction of input signal to
%           platform in scenario coordinates
% vecOut:   Vector defining direction of output signal from
%           platform in scenario coordinates
% freq:     Frequency in hertz of input signal

if iscell(platIn)
    plat = platIn{1};
else
    plat = platIn;
end

% Rotate vectors into platform's body frame
orient = fusion.internal.remotesensors. ...
    CommonUtilities.getRotmat(plat.Orientation);
vecIn = orient*vecIn;
vecOut = orient*vecOut;

% Make both vecIn and vecOut point away from platform
vecIn = -vecIn;

% Convert vectors into spherical coordinates (az,el) in platform's body frame
[azIn,elIn] = cart2sph(vecIn(1),vecIn(2),vecIn(3));
azIn = rad2deg(azIn);
elIn = rad2deg(elIn);

[azOut,elOut] = cart2sph(vecOut(1),vecOut(2),vecOut(3));
azOut = rad2deg(azOut);
elOut = rad2deg(elOut);

% Compute lookup angle
% - MBET says that equivalent monostatic angle is 1/2 of the
%   bistatic angle when bistatic angle is small (e.g. < 90 deg)
% - For monostatic case, bsAng == 0
% - Note: MBET will be used to estimate the bistatic RCS for angles > 90
%   degrees as well
bsAng = azOut-azIn;
bsAng = fusion.internal.units.interval(bsAng,[-180 180]);
azAng = azIn+bsAng/2;
azAng = fusion.internal.units.interval(azAng,[-180 180]);

bsAng = elOut-elIn;
bsAng = fusion.internal.units.interval(bsAng,[-180 180]);
elAng = elIn+bsAng/2;
elAng = fusion.internal.units.interval(elAng,[-180 180]);

% Lookup RCS at computed angle
signatures = plat.Signatures;
for m = 1:numel(signatures)
    thisSig = signatures{m};
    cond = isa(thisSig,'rcsSignature');
    if cond
        break
    end
end
coder.internal.errorIf(~cond,'shared_radarfusion:RemoteSensors:SignatureNotFound','RCS',plat.PlatformID);
tmp = value(thisSig, azAng, elAng, freq,'Parent',mfilename);
rcsdB = real(tmp(1)); 
end

% LocalWords:  emprop emtx tgt refl az el rotmat mbet vec bs ang
