rng('default'); % Reset the random number generator
simParameters = []; % Clear the simParameters variable
simParameters.NumFramesSim = 20; % Simulation time in terms of number of 10 ms frames
simParameters.SchedulingType = 0; % Set the value to 0 (slot based scheduling) or 1 (symbol based scheduling)

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

% 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

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
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;

% Create RLC channel configuration structure
rlcChannelConfigStruct.LCGID = 1; % Mapping between logical channel and logical channel group ID
rlcChannelConfigStruct.Priority = 1; % Priority of each logical channel
rlcChannelConfigStruct.PBR = 8; % Prioritized bitrate (PBR), in kilobytes per second, of each logical channel
rlcChannelConfigStruct.BSD = 10; % Bucket size duration (BSD), in ms, of each logical channel
rlcChannelConfigStruct.EntityType = simParameters.RLCConfig.EntityDir;
rlcChannelConfigStruct.LogicalChannelID = simParameters.LCHConfig.LCID;

% 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;

if ~isfield(simParameters, 'SchedulingType') || simParameters.SchedulingType == 0 % If no scheduling type is specified or slot based scheduling is specified
    rbAssignmentPlotPeriodicity = numSlotsFrame; % Update RB assignment visualization every frame (10 ms)
    tickGranularity = 14;
    simParameters.PUSCHMappingType = 'A';
    simParameters.PDSCHMappingType = 'A';
else % Symbol based scheduling
    rbAssignmentPlotPeriodicity = 1; % Update RB assignment visualization every slot
    tickGranularity = 1;
    simParameters.PUSCHMappingType = 'B';
    simParameters.PDSCHMappingType = 'B';
end

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
        
        % 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);

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 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);

% Create metrics visualization object for MAC and PHY metrics (DL direction is represented as 0)
visualizer = hNRMetricsVisualizer(simParameters, 'MACLogger', simSchedulingLogger{cellOfInterestIdx}, 'PhyLogger', simPhyLogger{cellOfInterestIdx}, 'VisualizationFlag', 0);

% Run processing loop
slotNum = 0;
numSymbolsSim = numSlotsSim * 14; % Simulation time in units of symbol duration
% Execute all the symbols in the simulation
for symbolNum = 1 : tickGranularity : numSymbolsSim
    if mod(symbolNum - 1, 14) == 0
        slotNum = slotNum + 1;
    end
    
    % 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});

        % Run MAC and PHY layers of UEs
        for ueIdx = 1:simParameters.NumUEsCell
            run(UEs{cellIdx, ueIdx});
        end

        % MAC logging
        logCellSchedulingStats(simSchedulingLogger{cellIdx}, symbolNum, gNB{cellIdx}, UEs(cellIdx, :), 0); % 0 for DL

        % PHY logging
        logCellPhyStats(simPhyLogger{cellIdx}, symbolNum, gNB{cellIdx}, UEs(cellIdx, :));
    end

    % Visualization
    % Check slot boundary
    if symbolNum > 1 && ((simParameters.SchedulingType == 1 && mod(symbolNum, 14) == 0) || (simParameters.SchedulingType == 0 && mod(symbolNum-1, 14) == 0))
        % RB assignment visualization (if enabled)
        if simParameters.RBVisualization
            if mod(slotNum, rbAssignmentPlotPeriodicity) == 0
                % Plot at slot boundary, if the update periodicity is reached
                plotRBGrids(simSchedulingLogger{cellOfInterestIdx});
            end
        end
        % CQI grid visualization (if enabled)
        if simParameters.CQIVisualization
            if mod(slotNum, numSlotsFrame) == 0 % Plot at frame boundary
                plotCQIRBGrids(simSchedulingLogger{cellOfInterestIdx});
            end
        end
        % If the update periodicity is reached, plot scheduler metrics and PHY metrics visualization
        % at slot boundary
        if mod(slotNum, simParameters.MetricsStepSize) == 0
            plotMetrics(visualizer, slotNum);
        end
    end

    for cellIdx = 1:simParameters.NumCells
        % Advance timer ticks for gNB and UEs by the number of symbols per slot
        advanceTimer(gNB{cellIdx}, tickGranularity);
        for ueIdx = 1:simParameters.NumUEsCell
            advanceTimer(UEs{cellIdx, ueIdx}, tickGranularity);
        end
    end
end

dlStats = getPerformanceIndicators(simSchedulingLogger{simParameters.CellOfInterest + 1});
[~, avgBLERLogs] = getBLERLogs(simPhyLogger{simParameters.CellOfInterest + 1});
fprintf('\nPeak DL throughput: %0.2f Mbps. Achieved average DL Throughput: %0.2f Mbps', dlStats(1, 1), dlStats(2, 1));

fprintf('\nAchieved average DL Goodput: %0.2f Mbps', dlStats(5, 1));

fprintf('\nPeak DL spectral efficiency: %0.2f bits/s/Hz. Achieved average DL spectral efficiency: %0.2f bits/s/Hz', dlStats(3, 1), dlStats(4, 1));

disp(['Block error rate for each UE in the downlink direction: [' num2str(round(avgBLERLogs(:, 1)', 2)) ']']);

% Get the logs
if(simParameters.DuplexMode == 0) % FDD
    logInfo = struct('NCellID', [], 'DLTimeStepLogs', [], 'SchedulingAssignmentLogs', [], 'BLERLogs', [], 'AvgBLERLogs', []);
else
    logInfo = struct('NCellID', [], 'TimeStepLogs', [], 'SchedulingAssignmentLogs', [], 'BLERLogs', [], 'AvgBLERLogs', []);
end
simulationLogs = cell(simParameters.NumCells, 1);
for cellIdx = 1:simParameters.NumCells
    logInfo.NCellID = simParameters.NCellIDList(cellIdx);
    if(simParameters.DuplexMode == 0) % FDD
        logInfo.DLTimeStepLogs = getSchedulingLogs(simSchedulingLogger{cellIdx});
    else % TDD
        logInfo.TimeStepLogs = getSchedulingLogs(simSchedulingLogger{cellIdx});
    end
    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