function box2Results = box2Simulation(PHYParameters,MACParameters,ScenarioParameters,SimParameters)
% box2Simulation Example helper function

%   Copyright 2019 The MathWorks, Inc.

%% Generate Transmitter Sites
% Before the main body of the simulation, the transmitter site objects
% |txsite| are generated and assigned room names of the form 'Room#-STA#'
% for ease of reference. One transmitter (AP) per room is assumed. Each
% transmitter is assumed to be isotropic.

% Total number of transmitters, assuming one transmitter (tx) per room
numTx = prod(ScenarioParameters.BuildingLayout);

% Create transmitter sites with an isotropic antenna element
roomNames = strings(1,numTx);
for siteInd = 1:numTx
    roomNames(siteInd) = "Room " + siteInd;
end
txs = txsite('Name',roomNames,...
    'TransmitterFrequency',PHYParameters.TransmitterFrequency, ...
    'TransmitterPower',10.^((PHYParameters.TxPower+PHYParameters.TxGain-30)/10));

%% Generate Receiver Sites
% The receive site objects |rxsite| are generated and assigned associated
% room names for reference (Room#-STA#). The Scenario parameter
% |NumRxPerRoom| is used to define how many receivers (STAs) are present in
% each room. Each receiver is assumed to be isotropic.

% Total number of receivers, assuming one transmitter (tx) per room
numRx = numTx*ScenarioParameters.NumRxPerRoom;

% Create receiver sites
roomNames = strings(1,numRx);
for siteInd = 1:numRx
    roomNames(siteInd) = "Room " + siteInd + "-" + (ceil(siteInd/numTx));
end
rxs = rxsite('Name',roomNames);

% Receiver noise power in dBm
T = 290; % Temperature (Kelvin)
k = physconst('Boltzmann'); % Boltzmann constant
% Sample rate (Hz)
fs = wlanSampleRate(wlanHESUConfig('ChannelBandwidth',PHYParameters.ChannelBandwidth));
rxNoisePower = 10*log10(k*T*fs)+30+PHYParameters.NoiseFigure;

%% Box 2 Simulation

Ntx = PHYParameters.NumTxAntennas;
Nrx = PHYParameters.NumRxAntennas;
Nsts = PHYParameters.NumSTS;

% Assume all links are HE
cfgHE = wlanHESUConfig;
cfgHE.ChannelBandwidth = PHYParameters.ChannelBandwidth;
cfgHE.NumTransmitAntennas = Ntx;
cfgHE.NumSpaceTimeStreams = Nsts;
cfgHE.SpatialMapping = 'Custom';
cfgHE.SpatialMappingMatrix = helperSpatialExpansionMatrix(cfgHE);
ofdmInfo = wlanHEOFDMInfo('HE-Data',cfgHE);
Nst = ofdmInfo.NumTones;

% Generate a channel
cfgChan = wlanTGaxChannel;
cfgChan.SampleRate = fs;
cfgChan.ChannelBandwidth = cfgHE.ChannelBandwidth;
cfgChan.DelayProfile = 'Model-D';
cfgChan.TransmitReceiveDistance = 15; % N-LOS
cfgChan.NumTransmitAntennas = Ntx;
cfgChan.NumReceiveAntennas = Nrx;
cfgChan.PathGainsOutputPort = true;
cfgChan.NormalizeChannelOutputs = false;
cfgChan.ChannelFiltering = false; % No need to filter input
cfgChan.EnvironmentalSpeed = 0; % Stationary channel
cfgChan.ChannelFiltering = false;
cfgChan.NumSamples = 1;

% Precoding includes spatial mapping matrix and cyclic shift
W = getPrecodingMatrix(cfgHE); % Nst-by-Nsts-by-Ntx

% Pre-allocate output
box2Results = struct;
box2Results.sinr = nan(Nst,Nsts,numTx,SimParameters.NumTxEventsPerDrop,SimParameters.NumDrops); % for storing SINR values
box2Results.sinrEff = nan(numTx,SimParameters.NumTxEventsPerDrop,SimParameters.NumDrops); % for storing effective SINR values

for drop = 1:SimParameters.NumDrops
    fprintf('  Running drop #%d/%d ...\n',drop,SimParameters.NumDrops);

    % Drop receivers in each room
    [association,txChannels,rxChannels] = tgaxDropNodes( ...
        txs,rxs,ScenarioParameters,MACParameters.NumChannels);
   
    % Generate propagation model
    propModel = TGaxResidential('roomSize',ScenarioParameters.RoomSize);

    % Calculate signal strength for all links
    signalStrength = sigstrength(rxs,txs,propModel,'Type','power',...
        'ReceiverGain',PHYParameters.RxGain,'Map','none'); % all signal strengths in dBm
       
    % Threshold signals that are not on same non-overlapping channel
    signalStrength(~(txChannels == rxChannels')) = -Inf;
    
    if SimParameters.Test == 1
        % Interference free
        signalStrength(~eye(numTx,numRx)) = -Inf;
    end
    
    % Mask the transmitter-receiver links that are non-negligible to simulate and get the linear indices
    nonnegligibleMask = signalStrength > -Inf; % Can set to skip very small channels
    
    Ptxrx = 10.^((signalStrength-30)/10); % Convert signal strength to Watts
    
    % Create channel manager object to get channel realizations
    fprintf('    Generating %d fading channel realizations ...\n',sum(nonnegligibleMask,'all'));
    chanManager = tgaxChannelManager(cfgChan,nonnegligibleMask);
    
    for txevent = 1:SimParameters.NumTxEventsPerDrop
        fprintf('    Running transmission event #%d/%d ...\n',txevent,SimParameters.NumTxEventsPerDrop);
    
        % Determine active transmitters and receivers with Clear Channel Assessment
        switch SimParameters.Test
            case 1
                % Interference free
                activeTx = true(size(signalStrength,1),1);
                activeRx = [true(size(signalStrength,1),1); false(size(signalStrength,2)-size(signalStrength,1))]; 
                ccaActiveMask = eye(numel(activeTx),numel(activeRx));
            case 2
                % DL only
                activeTx = true(size(signalStrength,1),1);
                activeRx = [true(size(signalStrength,1),1); false(size(signalStrength,2)-size(signalStrength,1))]; 
                ccaActiveMask = true(numel(activeTx),numel(activeRx));
            case 3
                % All nodes per channel access rules
                [activeTx,activeRx] = tgaxCCA(signalStrength,MACParameters.CCALevel);
                ccaActiveMask = logical(activeTx.*activeRx');
        end
        
        % Remove active links if transmissions not scheduled
        nonnegligibleMaskTxEvent = nonnegligibleMask&ccaActiveMask;

        % Calculate signal power (STXRX), inter-stream interference power
        % (ISTXRX), and interference power (IOTXRX)
        activeRxIdx = find(activeRx);
        numActiveRx = sum(activeRx);
        sinr = zeros(Nst,Nsts,numActiveRx);
        for irx = 1:numActiveRx
            % For each active receiver calculate the SINR
            rxIdx = activeRxIdx(irx); % Receiver of interest index
            txIdx = find(association(:,rxIdx)); % Transmitter of interest index
            
            % Find the index of active interferers
            tm = true(numTx,1); % Assume all transmitters are interferes
            tm(txIdx) = false; % Exclude transmitter of interest
            interferers = find(nonnegligibleMaskTxEvent(:,rxIdx)&tm);

            % Generate channel of interest
            Htxrx_soi = getChannelMatrix(chanManager,cfgHE,txIdx,rxIdx);
            Ptxrx_soi = Ptxrx(txIdx,rxIdx);
            Wtx_soi = W; % Assume same precoding used at all transmitters

            % For each interferer generate a channel and get the receive power and precoding used
            numInterf = numel(interferers);
            Htxrx_int = cell(numInterf,1);
            Wtx_int = cell(numInterf,1);
            Ptxrx_int = zeros(numInterf,1);
            for iint = 1:numInterf
                Htxrx_int{iint} = getChannelMatrix(chanManager,cfgHE,interferers(iint),rxIdx);
                Wtx_int{iint} = W; % Assume same precoding used at all transmitters
                Ptxrx_int(iint) = Ptxrx(interferers(iint),rxIdx); % Interference power
            end

            % Calculate SINR
            N0 = 10.^((rxNoisePower-30)/10);
            sinr(:,:,irx) = calculateSINR(Htxrx_soi,Ptxrx_soi,Wtx_soi,N0,Htxrx_int,Ptxrx_int,Wtx_int);
        end

        % Effective SINR calculation in TGax Evaluation Methodology calibration
        sinrEffective = 10*log10(2.^((1/(Nst*Nsts))*sum(sum(log2(1+permute(10.^(sinr/10),[3 1 2])),3),2))-1);
        box2Results.sinr(:,:,activeRx,txevent,drop) = sinr;
        box2Results.sinrEff(activeRx,txevent,drop) = sinrEffective;
    end
end

end