rng('default'); % Reset the random number generator
simParameters = []; % Clear the simParameters variable

simParameters.NumFramesSim = 30; % Simulation time in terms of number of 10 ms frames

simParameters.NumCells = 3; % Number of cells
simParameters.CellRadius = 500; % Radius of each cell (in meters)

simParameters.GNBPosition = [1700 600;
    3000 600;
    2500 2000];

% Number of UEs in a cell. Each cell contains same number of UEs
simParameters.NumUEsCell = 4;

simParameters.GNBTxPower = 32; % Tx power for gNB (in dBm)

% Set the channel bandwidth to 10 MHz and subcarrier spacing (SCS) to 30 kHz
% as defined in 3GPP TS 38.104 Section 5.3.2.
simParameters.NumRBs = 24;
simParameters.SCS = 30; % kHz
simParameters.DLCarrierFreq = 2.635e9; % Hz
simParameters.DLBandwidth = 10e6; % Hz

% SINR to CQI mapping table for 0.1 BLER
simParameters.SINR90pc = [-5.46 -0.46 4.54 9.05 11.54 14.04 15.54 18.04 ...
    20.04 22.43 24.93 25.43 27.43 30.43 33.43];

% Maximum limit on the RBs allotted for PDSCH. Transmission limit is
% applicable for new PDSCH assignments and not for the retransmissions
simParameters.RBAllocationLimitDL = 15; % For PDSCH

% Specify the ID of cell of interest. All the visualizations and metrics are shown for this cell
simParameters.CellOfInterest = 2; % Set a value from 0 to NumCells-1

% The parameters CQIVisualization and RBVisualization control the display
% of these visualizations: (i) CQI visualization of RBs (ii) RB assignment
% visualization. By default, the 'RBVisualization' plot is disabled. You
% can enable it by setting to 'true'
simParameters.CQIVisualization = true;
simParameters.RBVisualization = false;
% The output metrics plots are updated periodically NumMetricsSteps times within the
% simulation duration
simParameters.NumMetricsSteps = 20;
% MAT-files to write the logs into. They are used for post simulation analysis and visualization
simParameters.ParametersLogFile = 'simParameters'; % For logging the simulation parameters
simParameters.SimulationLogFile = 'simulationLogs'; % For logging the simulation logs

% Set the periodic DL application traffic model for UEs. The following
% configuration applies for each cell
simParameters.DLPacketPeriodicityUEs = 10; % Periodicity (in ms) at which the DL packets are generated for UEs at gNB
simParameters.DLPacketSizesUEs = 2e4; % Size of the DL packets generated (in bytes) for UEs at gNB

% Validate the simulation configuration
hNRIntercellInterferenceValidateConfig(simParameters);

simParameters.DuplexMode = 0; % FDD
simParameters.ULBandwidth = 10e6; % Hz
simParameters.ULCarrierFreq = 2.515e9; % Hz
simParameters.NumUEs = simParameters.NumUEsCell; % Number of UEs in a cell
simParameters.NCellIDList = 0:simParameters.NumCells-1; % List of physical cell IDs
% To store waveforms at the receiver, set the value to twice the number of transmitters
simParameters.UERxBufferSize = 2 * simParameters.NumCells; % Reception buffer size

% CSI-RS resource configuration. All UEs are assumed to measure channel quality on same CSI-RS resource
simParameters.CSIRSRowNumber = 2; % Possible row numbers for single transmit antenna case are 1 and 2
simParameters.SubbandSize = 8; % Size of sub-band for CQI reporting in terms of number of RBs

% Set the BSRPeriodicity to 'inf' as there is no UL traffic
simParameters.BSRPeriodicity = inf; % In ms

% Slot duration for the selected SCS and number of slots in a 10 ms frame
slotDuration = 1/(simParameters.SCS/15); % In ms
numSlotsFrame = 10/slotDuration; % Number of slots in a 10 ms frame
numSlotsSim = simParameters.NumFramesSim * numSlotsFrame; % Number of slots in the simulation

% Interval at which metrics visualization updates in terms of number of
% slots. As one slot is the finest time-granularity of the simulation, make
% sure that MetricsStepSize is an integer
simParameters.MetricsStepSize = ceil(numSlotsSim / simParameters.NumMetricsSteps);
if mod(numSlotsSim, simParameters.NumMetricsSteps) ~= 0
    % Update the NumMetricsSteps parameter if NumSlotsSim is not
    % completely divisible by it
    simParameters.NumMetricsSteps = floor(numSlotsSim / simParameters.MetricsStepSize);
end

% DL packet periodicity for UEs in terms of number of slots
appPeriodicityUEsSlotsDL = simParameters.DLPacketPeriodicityUEs / slotDuration;

simParameters.NumLogicalChannels = 1; % Only 1 logical channel is assumed in each UE in this example

% Logical channel configuration applies for all the nodes (UEs and gNBs) in the simulation
% Mapping between logical channel and logical channel group ID
simParameters.LCHConfig.LCGID = 1;
% Priority of each logical channel
simParameters.LCHConfig.Priority = 1;
% Prioritized bitrate (PBR) of each logical channel (in kilo bytes per second)
simParameters.LCHConfig.PBR = 8;
% Bucket size duration (BSD) of each logical channel (in ms). However, the priority,
% PBR and BSD for logical channel is not relevant in this example as single logical
% channel is assumed
simParameters.LCHConfig.BSD = 10;
% Logical channel ID (logical channel ID of data radio bearers starts from 4)
simParameters.LCHConfig.LCID = 4;

% RLC entity direction. Value 0 represents DL only, 1
% represents UL only and 2 represents both UL and DL
% directions. Setting entity direction to have both UL and DL
simParameters.RLCConfig.EntityDir = 0;

% Maximum RLC service data unit (SDU) length (in bytes) as per 3GPP TS 38.323
simParameters.maxRLCSDULength = 9000;

% Generate the positions of UEs in each cell
simParameters.UEPosition = generateUEPositions(simParameters);

% Total number of UEs in the simulation
simParameters.MaxReceivers = simParameters.NumCells * simParameters.NumUEsCell;

gNB = cell(simParameters.NumCells, 1);
UEs = cell(simParameters.NumCells, simParameters.NumUEsCell);
% Create DL packet distribution object
dlPacketDistributionObj = hNRPacketDistribution(simParameters, 0); % 0 for DL
% Create UL packet distribution object
ulPacketDistributionObj = hNRPacketDistribution(simParameters, 1); % 1 for UL

for cellIdx = 1:simParameters.NumCells % For each cell
    
    simParameters.NCellID = simParameters.NCellIDList(cellIdx); % Cell ID
    simParameters.Position = [simParameters.GNBPosition(cellIdx, :) 0]; % gNB position in (x,y,z) coordinates
    gNB{cellIdx} = hNRGNB(simParameters); % Create gNB node
    scheduler = hNRSchedulerProportionalFair(simParameters); % Create proportional fair scheduler
    addScheduler(gNB{cellIdx}, scheduler); % Add scheduler to gNB
    
    gNB{cellIdx}.PhyEntity = hNRGNBPhy(simParameters); % Create the PHY layer instance
    configurePhy(gNB{cellIdx}, simParameters); % Configure the PHY layer
    setPhyInterface(gNB{cellIdx}); % Set up the interface to PHY layer
    
    % For each cell, create the set of UE nodes and place them randomly within the cell radius
    for ueIdx = 1:simParameters.NumUEsCell
        simParameters.Position = [simParameters.UEPosition{cellIdx}(ueIdx, :) 0]; % Position of UE in (x,y,z) coordinates
        UEs{cellIdx, ueIdx} = hNRUE(simParameters, ueIdx);
        UEs{cellIdx, ueIdx}.PhyEntity = hNRUEPhy(simParameters, ueIdx); % Create the PHY layer instance
        configurePhy(UEs{cellIdx, ueIdx}, simParameters); % Configure the PHY layer
        setPhyInterface(UEs{cellIdx, ueIdx}); % Set up the interface to PHY
        
        % Create RLC channel configuration structure
        rlcChannelConfigStruct.EntityType = simParameters.RLCConfig.EntityDir;
        rlcChannelConfigStruct.LogicalChannelID = simParameters.LCHConfig.LCID;
        rlcChannelConfigStruct.LCGID = simParameters.LCHConfig.LCGID;
        rlcChannelConfigStruct.Priority = simParameters.LCHConfig.Priority;
        rlcChannelConfigStruct.PBR = simParameters.LCHConfig.PBR;
        rlcChannelConfigStruct.BSD = simParameters.LCHConfig.BSD;
        % Setup logical channel at gNB for the UE
        configureLogicalChannel(gNB{cellIdx}, ueIdx, rlcChannelConfigStruct);
        % Setup logical channel at UE
        configureLogicalChannel(UEs{cellIdx, ueIdx}, ueIdx, rlcChannelConfigStruct);
        
        % Add data traffic pattern generators to gNB node
        packetSize = simParameters.DLPacketSizesUEs;
        % Calculate the data rate (in kbps) of On-Off traffic pattern using
        % packet size (in bytes) and packet interval (in ms)
        dataRate = ceil(1000/simParameters.DLPacketPeriodicityUEs) * packetSize * 8e-3;
        % Limit the size of the generated application packet to the maximum
        % RLC SDU size. The maximum supported RLC SDU size is 9000 bytes
        if packetSize > simParameters.maxRLCSDULength
            packetSize = simParameters.maxRLCSDULength;
        end
        % Create an object for On-Off network traffic pattern for the specified
        % UE and add it to the gNB. This object generates the downlink data
        % traffic on the gNB for the UE
        app = networkTrafficOnOff('PacketSize', packetSize, 'GeneratePacket', true, ...
            'OnTime', simParameters.NumFramesSim/100, 'OffTime', 0, 'DataRate', dataRate);
        gNB{cellIdx}.addApplication(ueIdx, simParameters.LCHConfig.LCID, app);
    end
    
    % Setup the UL and DL packet distribution mechanism
    hNRSetUpPacketDistribution(simParameters, gNB{cellIdx}, UEs(cellIdx, :), dlPacketDistributionObj, ulPacketDistributionObj);
end

% Display network topology
plotNetwork(simParameters);

% To store these UE metrics for each slot: throughput bytes
% transmitted, goodput bytes transmitted, and pending buffer amount bytes.
% The number of goodput bytes is calculated by excluding the
% retransmissions from the total transmissions
UESlotMetricsDL = zeros(simParameters.NumUEsCell, 3);

% To store last received new data indicator (NDI) values for DL HARQ process
HARQProcessStatusDL = zeros(simParameters.NumUEsCell, 16); % Max 16 HARQ process

% To store current DL CQI values on the RBs for different UEs
downlinkChannelQuality = zeros(simParameters.NumUEsCell, simParameters.NumRBs);

% Store the DL BLER for each UE
dlBLERStats = zeros(simParameters.NumUEsCell, 2);

simSchedulingLogger = cell(simParameters.NumCells, 1);
simPhyLogger = cell(simParameters.NumCells, 1);
for cellIdx = 1:simParameters.NumCells
    % Create an object for MAC DL scheduling information visualization and logging
    simParameters.NCellID = simParameters.NCellIDList(cellIdx);
    simSchedulingLogger{cellIdx} = hNRSchedulingLogger(simParameters, 0); % 0 for DL
    % Create an object for BLER visualization and PHY layer metrics logging
    simPhyLogger{cellIdx} = hNRPhyLogger(simParameters, 0); % 0 for DL
end

% Store the index of cell ID of interest
cellOfInterestIdx = find(simParameters.CellOfInterest == simParameters.NCellIDList);

symbolNum = 0;
% Run processing loop
for slotNum = 1:numSlotsSim
    
    % All the cells operating on same SCS hence slot durations are same
    for cellIdx = 1:simParameters.NumCells % For each cell
        
        % Run MAC and PHY layers of gNB
        run(gNB{cellIdx}.MACEntity);
        run(gNB{cellIdx}.PhyEntity);
        
        % Run MAC and PHY layers of UEs
        for ueIdx = 1:simParameters.NumUEsCell
            % Read the last received NDI flags for DL HARQ processes for
            % logging (Reading it before it gets overwritten by run function of MAC)
            HARQProcessStatusDL(ueIdx, :) = getLastNDIFlagHarq(UEs{cellIdx, ueIdx}.MACEntity, 0); % 0 for DL
            run(UEs{cellIdx, ueIdx}.MACEntity);
            run(UEs{cellIdx, ueIdx}.PhyEntity);
        end
        
        % MAC logging
        % Read DL assignments done by gNB MAC scheduler at current time.
        % Resource assignments returned by a scheduler is empty, if
        % scheduler was not scheduled to run at the current time or no
        % resources got scheduled
        [~, resourceAssignmentsDL] = getCurrentSchedulingAssignments(gNB{cellIdx}.MACEntity);
        % Read throughput and goodput bytes sent for each UE
        [UESlotMetricsDL(:, 1), UESlotMetricsDL(:, 2)] = getTTIBytes(gNB{cellIdx});
        UESlotMetricsDL(:, 3) = getBufferStatus(gNB{cellIdx}); % Read pending buffer (in bytes) on gNB, for all the UEs
        
        for ueIdx = 1:simParameters.NumUEsCell
            % Read the DL channel quality at gNB for each of the UEs for logging
            downlinkChannelQuality(ueIdx,:) = getChannelQuality(gNB{cellIdx}, 0, ueIdx); % 0 for DL
        end
        
        % Update DL scheduling logs based on the current slot run of UEs
        % and gNB. Logs are updated in each slot, RB grid visualizations
        % are updated every frame, and metrics plots are updated every
        % metricsStepSize slots
        logScheduling(simSchedulingLogger{cellIdx}, symbolNum + 1, resourceAssignmentsDL, UESlotMetricsDL, downlinkChannelQuality, HARQProcessStatusDL, 0); % 0 for DL
        
        % PHY logging
        for ueIdx = 1:simParameters.NumUEsCell
            dlBLERStats(ueIdx, :) = getDLBLER(UEs{cellIdx, ueIdx}.PhyEntity);
        end
        % Log the DL BLER statistics
        logBLERStats(simPhyLogger{cellIdx}, dlBLERStats, []);
    end
    
    for cellIdx = 1:simParameters.NumCells
        % Advance timer ticks for gNB and UEs by the number of symbols per slot
        advanceTimer(gNB{cellIdx}, 14);
        for ueIdx = 1:simParameters.NumUEsCell
            advanceTimer(UEs{cellIdx, ueIdx}, 14);
        end
    end
    
    % RB assignment visualization (if enabled)
    if simParameters.RBVisualization
        if mod(slotNum, numSlotsFrame) == 0
            plotRBGrids(simSchedulingLogger{cellOfInterestIdx});
        end
    end
    
    % CQI grid visualization (if enabled)
    if simParameters.CQIVisualization
        if mod(slotNum, numSlotsFrame) == 0
            plotCQIRBGrids(simSchedulingLogger{cellOfInterestIdx});
        end
    end
    
    % Plot scheduler metrics and PHY metrics visualization at slot
    % boundary, if the update periodicity is reached
    if mod(slotNum, simParameters.MetricsStepSize) == 0
        plotMetrics(simSchedulingLogger{cellOfInterestIdx});
        plotMetrics(simPhyLogger{cellOfInterestIdx});
    end
    
    % Symbol number in the simulation
    symbolNum = symbolNum + 14;
end

% Get the logs
logInfo = struct('NCellID',[], 'DLTimeStepLogs',[], 'SchedulingAssignmentLogs',[] ,'BLERLogs',[]);
simulationLogs = cell(simParameters.NumCells, 1);
for cellIdx = 1:simParameters.NumCells
    logInfo.NCellID = simParameters.NCellIDList(cellIdx);
    [logInfo.DLTimeStepLogs, ~] = getSchedulingLogs(simSchedulingLogger{cellIdx});
    logInfo.SchedulingAssignmentLogs = getGrantLogs(simSchedulingLogger{cellIdx}); % Scheduling assignments log
    logInfo.BLERLogs = getBLERLogs(simPhyLogger{cellIdx}); % Block error rate logs
    simulationLogs{cellIdx, 1} = logInfo;
end
save(simParameters.ParametersLogFile, 'simParameters'); % Save simulation parameters in a MAT-file
save(simParameters.SimulationLogFile, 'simulationLogs'); % Save simulation logs in a MAT-file

function plotNetwork(simParameters)
% Create the figure
figure('Name', 'Network Topology Visualization', 'units', 'normalized', 'outerposition', [0 0 1 1], 'Visible', "on");
title('Network Topology Visualization');
hold on;

for cellIdx = 1:simParameters.NumCells
    
    % Plot the circle
    th = 0:pi/60:2*pi;
    xunit = simParameters.CellRadius * cos(th) + simParameters.GNBPosition(cellIdx, 1);
    yunit = simParameters.CellRadius * sin(th) + simParameters.GNBPosition(cellIdx, 2);
    if simParameters.CellOfInterest == simParameters.NCellIDList(cellIdx)
        h1 =  plot(xunit, yunit, 'Color', 'green'); % Cell of interest
    else
        h2 =  plot(xunit, yunit, 'Color', 'red');
    end
    xlabel('X-Position (meters)')
    ylabel('Y-Position (meters)')
    % Add tool tip data for gNBs
    s1 = scatter(simParameters.GNBPosition(cellIdx, 1), simParameters.GNBPosition(cellIdx, 2), '^','MarkerEdgeColor', 'magenta');
    cellIdRow = dataTipTextRow('Cell - ',{num2str(simParameters.NCellIDList(cellIdx))});
    s1.DataTipTemplate.DataTipRows(1) = cellIdRow;
    posRow = dataTipTextRow('Position[X, Y]: ',{['[' num2str(simParameters.GNBPosition(cellIdx, :)) ']']});
    s1.DataTipTemplate.DataTipRows(2) = posRow;
    
    % Add tool tip data for UEs
    uePosition = simParameters.UEPosition{cellIdx};
    for ueIdx = 1:size(uePosition, 1)
        s2 = scatter(uePosition(ueIdx, 1), uePosition(ueIdx, 2), '.','MarkerEdgeColor', 'blue');
        ueIdRow = dataTipTextRow('UE - ',{num2str(ueIdx)});
        s2.DataTipTemplate.DataTipRows(1) = ueIdRow;
        posRow = dataTipTextRow('Position[X, Y]: ',{['[' num2str(uePosition(ueIdx, :)) ']']});
        s2.DataTipTemplate.DataTipRows(2) = posRow;
    end
end
% Create the legend
if simParameters.NumCells > 1
    legend([h1 h2 s1 s2], 'Cell of interest', 'Interfering cells', 'gNodeB', 'UE', 'Location', 'northeastoutside')
else
    legend([h1 s1 s2], 'Cell of interest', 'gNodeB', 'UE', 'Location', 'northeastoutside')
end
axis([0 4000 0 3000]); % Set axis limits
hold off;
daspect([1000,1000,1]); % Set data aspect ratio
end

function uePositions = generateUEPositions(simParameters)
% Return the position of UEs in each cell

uePositions = cell(simParameters.NumCells, 1);
for cellIdx=1:simParameters.NumCells
    gnbXCo = simParameters.GNBPosition(cellIdx, 1); % gNB X-coordinate
    gnbYCo = simParameters.GNBPosition(cellIdx, 2); % gNB Y-coordinate
    theta = rand(simParameters.NumUEsCell, 1)*(2*pi);
    % Expression to calculate position of UEs with in the cell. By default,
    % it will place the UEs randomly with in the cell
    r = sqrt(rand(simParameters.NumUEsCell, 1))*simParameters.CellRadius;
    x = round(gnbXCo + r.*cos(theta));
    y = round(gnbYCo + r.*sin(theta));
    uePositions{cellIdx} = [x y];
end
end