Analyze NTN Coverage and Capacity for LEO Mega-Constellation

Since R2024b

This example shows how to analyze the coverage and capacity of a very low population density region for a mega-constellation of low-Earth orbit (LEO) satellites for satellite-based non-terrestrial network (NTN) applications. The coverage analysis includes:

• Percentage of satellite visibility for all the users in a region

• Percentage of link availability for all the users in a region

• Percentage of users with capacity greater than a threshold

This example assumes:

• All satellite links are regenerative.

• Each satellite can communicate to all other visible satellites.

Simulation Setup

This example models a service link operating at an S-band frequency with a constellation of LEO satellites as transmitters and multiple user equipment (UEs) as receivers in the region of interest (ROI). This example uses the NTN simulation setup defined in the Hexa-X D5.3 Deliverable [1].

Create Satellite Scenario

Create a satellite scenario using the satelliteScenario object. Specify the start time, stop time, and sample time for the scenario. Choose the stop time and sample time based on the satellite altitude, to analyze coverage for one complete satellite revolution around the Earth. The example uses a stop time of 2 hours to cover the time period of revolution of an LEO satellite at an altitude of 1325 km.

sc = satelliteScenario;
sc.StartTime = datetime(2024,4,29,09,01,57); % In datetime
sc.StopTime = sc.StartTime + hours(2);       % In datetime
sc.SampleTime = 60;                          % In s

Specify LEO Satellite Constellation

Define the configurations for constellations of LEO satellites, where each constellation consists of multiple orbital shells. An orbital shell represents a specific layer surrounding the Earth, defined by its altitude and inclination, as well as the number of orbits and number of satellites within each of those orbits. If you have an Aerospace Toolbox license, you can also create the shell of satellites by using the walkerStar (Aerospace Toolbox) function. This example uses three satellite constellations, and treats each constellation as a separate entity, conducting an analysis for each one. The example uses the parameters specified in the Hexa-X D5.3 Deliverable Table 4-1 [1] and the shell parameters correspond to the Telesat system of LEO satellite constellation [2].

% Define the satellite constellations. 
% All these fields must have same number of entries: ShellNames,
% ShellAltitude, ShellInclination, NumOrbits, and NumSatellitesPerOrbit
constellation = struct;
constellation(1).Name = "Constellation 1";
constellation(1).ShellNames = ["Altitude 1", "Altitude 2"];
constellation(1).ShellAltitude = [1015 1325];              % In km
constellation(1).ShellInclination = [98.98 50.88];         % In degrees
constellation(1).NumOrbits = [3 6];                        % Number of orbital planes
constellation(1).NumSatellitesPerOrbit = [5 5];            % Number of satellites per orbital plane

constellation(2).Name = "Constellation 2";
constellation(2).ShellNames = ["Altitude 1", "Altitude 2"];
constellation(2).ShellAltitude = [1015 1325];              % In km
constellation(2).ShellInclination = [98.98 50.88];         % In degrees
constellation(2).NumOrbits = [6 20];                       % Number of orbital planes
constellation(2).NumSatellitesPerOrbit = [13 11];          % Number of satellites per orbital plane

constellation(3).Name = "Constellation 3";
constellation(3).ShellNames = ["Altitude 1", "Altitude 2"];
constellation(3).ShellAltitude = [1015 1325];              % In km
constellation(3).ShellInclination = [98.98 50.88];         % In degrees
constellation(3).NumOrbits = [12 20];                      % Number of orbital planes
constellation(3).NumSatellitesPerOrbit = [13 22];          % Number of satellites per orbital plane

% Validate dimensions of each constellation
validateConstellation(constellation)

Specify the transmitter parameters for all the satellites. The example operates all the satellites with the same transmitter configuration, using an isotropic antenna for the transmitter.

txConfig = struct;
txConfig.Frequency = 2e9;  % Carrier frequency in Hz
txConfig.Power = 20;       % Transmitter power in dBW
txConfig.BitRate = 10;     % Bit rate in Mbps
txConfig.SystemLoss = 0;   % Transmitter system loss in dB
txConfig.Bandwidth = 10e6; % Link bandwidth in Hz

Specify User Equipment in ROI

Define the region of interest (ROI) by providing the limits of latitude and longitude. Specify the number of UEs deployed in this ROI, as these UEs act as notional reference to compute the coverage metrics in the ROI. This example uses ground stations as UEs, and operates all the ground stations with the same receiver configuration. Assume each ground station is equipped with an isotropic antenna. This example deploys 50 users in a region of the North Atlantic Ocean that contains a portion of North America and Africa.

lat = [20 70];     % Limits of latitude in degrees ([min max])
lon = [-90 20];    % Limits of longitude in degrees ([min max])
numUEs = 50;       % Number of UE
minElevAngle = 30; % Minimum elevation angle in degrees

Specify the receiver parameters for all the ground stations.

rxConfig = struct;
rxConfig.MaxGByT = -5;        % Maximum gain-to-noise-temperature in dB/K
rxConfig.SystemLoss = 0;      % Receiver system loss in dB
rxConfig.PreReceiverLoss = 0; % Pre-receiver loss in dB
rxConfig.RequiredEbNo = 11;   % Required bit energy to noise power spectral density ratio in dB

Show ROI in Satellite Scenario

Use meshgrid to find the region covering the specified latitude and longitude limits. Find the uniformly spaced coordinate points (in degrees) within the region, and place the ground stations. View the ground stations in the ROI using the satellite scenario.

% Find the coordinates of uniformly spaced UEs within the ROI
[latSpacing,lonSpacing] = findClosestFactors(numUEs);
latPts = linspace(lat(1),lat(2),latSpacing);
lonPts = linspace(lon(1),lon(2),lonSpacing);
[latMesh,lonMesh] = meshgrid(latPts,lonPts);
latCoord = latMesh(:);
lonCoord = lonMesh(:);

% Add the ground station to the scenario
ue = groundStation(sc,latCoord,lonCoord);
[ue.MinElevationAngle] = minElevAngle;

% Add the receive antenna to the ground station
isotropic = arrayConfig(Size=[1 1]);
rx = receiver(ue,Antenna=isotropic, ...
    SystemLoss=rxConfig.SystemLoss, ...
    PreReceiverLoss=rxConfig.PreReceiverLoss, ...
    GainToNoiseTemperatureRatio=rxConfig.MaxGByT);

% Show the ROI
s = satelliteScenarioViewer(sc,ShowDetails=false);

Calculate Coverage Metrics

Calculate the coverage metrics for the three specified constellations. In each constellation, add the satellites in the specified orbital shells, add the transmitter with the specified configuration to the satellites, and for each UE follow these steps:

• Find the satellites visible to the UE using elevation angle.

• Create a temporary link to the UEs for all visible satellites in the entire simulation duration.

• Use ebno to calculate the received bit energy to noise power spectral density ratio (Eb/No) at the UE for all the links.

• Store the maximum Eb/No at the UE for each time instance.

• Compute the maximum carrier-to-noise (CNR) ratio at the UE for each time instance.

• Compute the maximum possible link capacity for the maximum CNR and bandwidth using the Shannon theorem.

% Delete viewer to avoid animating all the constellation and temporary
% links
delete(s)

% For each timestamp and each UE, find the maximum signal power
timeSteps = sc.StartTime:seconds(sc.SampleTime):sc.StopTime;
numTimeSteps = numel(timeSteps);
% Store the results in an output structure
resultPerConstellation = struct("SatelliteVisibility",[], ...
    "LinkAvailability",[], ...
    "MinLinkAvailability",[], ...
    "MinSatelliteVisibility",[], ...
    "MaxEbNo",[], ...
    "MaxCNR",[], ...
    "Capacity",[], ...
    "NumTotalSatellites",[]);
results = repmat(resultPerConstellation,numel(constellation),1);

% Loop each constellation
for constIdx = 1:numel(constellation)
    % Display status message at the start of processing loop
    fprintf("Computing coverage statistics for constellation " + constIdx)

    % Initialize the variables
    [linkAvailability,satelliteVisibility, ...
       maxEbNo] = deal(zeros(numUEs,numTimeSteps));

    % Add satellites to the scenario
    sat = addSatellites(sc,constellation(constIdx));
    % Add transmitter to the satellites
    tx = transmitter(sat, ...
        Frequency=txConfig.Frequency, ...
        Power=txConfig.Power, ...
        SystemLoss=txConfig.SystemLoss, ...
        BitRate=txConfig.BitRate, ...
        Antenna=arrayConfig(Size=[1 1]));

    for ueIdx = 1:numUEs
        % Display progress through dots in intervals of 20%
        if mod(ueIdx,ceil(0.2*numUEs)) == 0
            fprintf(".")
        end

        % Find the elevation angle of all the satellites with respect to
        % UE
        [~,el] = aer(ue(ueIdx),sat);
        % Find the satellites that have elevation angles greater than or equal
        % to minimum elevation angle
        elIdx = el >= ue(ueIdx).MinElevationAngle;
        % Satellite visibility: Set to 1 when a satellite is visible to the
        % UE
        satelliteVisibility(ueIdx,:) = any(elIdx,1);
        % Create temporary links for all the satellites that are visible over
        % the simulation duration
        validTx = tx(any(elIdx,2));
        if ~isempty(validTx)
            links = link(validTx,rx(ueIdx));
            % Calculate the received Eb/No
            ebByNo = ebno(links);
            % Get the maximum Eb/No when link is closed
            maxEbNo(ueIdx,:) = max(ebByNo);
            % Link Availability: Set to 1 when the link is closed
            linkAvailability(ueIdx,:) = any(linkStatus(links));
            % Delete links
            delete(links)
        end
    end

    % Compute the max CNR based on max Eb/No
    cno = maxEbNo + pow2db(tx(1).BitRate) + 60;  % Bit rate is in Mbps
    cnrdB = cno - pow2db(txConfig.Bandwidth);
    cnr = db2pow(cnrdB);
    % Compute capacity using Shannon theorem, only for available links
    capacity = txConfig.Bandwidth*log2(1+cnr); % bps

    % Capture the results to the output structure
    results(constIdx).SatelliteVisibility = satelliteVisibility;
    results(constIdx).MinSatelliteVisibility = min(mean(satelliteVisibility,2));
    results(constIdx).LinkAvailability = linkAvailability;
    results(constIdx).MinLinkAvailability = min(mean(linkAvailability,2));
    results(constIdx).MaxEbNo = maxEbNo;
    results(constIdx).MaxCNR = cnr;
    results(constIdx).Capacity = capacity;
    results(constIdx).NumTotalSatellites = numel(sat);

    % Delete the satellites
    delete(sat)

    % Display status message at end of processing loop
    fprintf(newline + "Computed coverage statistics for constellation " + constIdx + newline)
end

Computing coverage statistics for constellation 1

.....

Computed coverage statistics for constellation 1

Computing coverage statistics for constellation 2

.....

Computed coverage statistics for constellation 2

Computing coverage statistics for constellation 3

.....

Computed coverage statistics for constellation 3

Calculate these three coverage metrics: satellite visibility, link availability, and the probability of UEs with capacity greater than the threshold.

Satellite Visibility

The satellite visibility for a UE is defined as the ratio of the number of time instances with visible satellites to the total number of time instances. The visibility of a satellite over a region is determined by the lowest visibility level among all the UEs in that area. If the value of satellite visibility over region is x, it implies that all the UEs in a region have satellite visibility for at least x percent of the time.

satelliteVisibility = [results.MinSatelliteVisibility]'*100; % Value in percentage

Link Availability

The link availability of a UE is defined as the ratio of the number of time instances with nonnegative link margin to the total number of time instances. The availability of a link over a region is determined by the lowest availability level among all the UEs in that area. If the value of link availability over region is x, it implies that all the UEs in a region have valid communication links for at least x percent of the time.

linkAvailability = [results.MinLinkAvailability]'*100; % Value in percentage

Probability of UEs with Capacity Greater Than Threshold

The probability of UEs with capacity greater than a threshold at a particular time is defined as the ratio of the number of UEs with capacity greater than that threshold to the total number of UEs. The probability of UEs with capacity greater than a threshold over the simulation time is the mean of probability of UEs with capacity greater than a threshold across all the time instances. If the value of the probability is x, it implies that an average of x percent of users have capacity greater than the threshold over all the entries. This example uses a threshold of 1 Mbps.

coveragePercent = zeros(numel(results),1);
threshold = 1;                                                 % In Mbps
for i = 1:numel(results)
    cTemp = results(i).Capacity;
    coveragePercent(i) = mean(cTemp > threshold*1e6,[1 2])*100;
end

Visualize Results

Tabulate the coverage metrics.

numTotalSatellites = [results.NumTotalSatellites];
constellationNames = [constellation.Name];
resultTable = table(constellationNames',numTotalSatellites',satelliteVisibility,linkAvailability,coveragePercent, ...
    VariableNames={'Constellation Name','Num. Satellites','Satellite Visibility (%)','Link Availability (%)', ...
    ['Capacity > ',num2str(threshold),' Mbps (%x)']})

resultTable=3×5 table
    Constellation Name    Num. Satellites    Satellite Visibility (%)    Link Availability (%)    Capacity > 1 Mbps (%x)
    __________________    _______________    ________________________    _____________________    ______________________

    "Constellation 1"            45                   9.0909                    3.3058                    36.893        
    "Constellation 2"           298                   73.554                    63.636                    94.893        
    "Constellation 3"           596                   99.174                    97.521                    99.983        

Plot the satellite visibility and link availability.

figure(1)
% Plot the percentage of satellite visibility of the region for the whole
% simulation duration
plot(numTotalSatellites,satelliteVisibility,"-*")
grid on
box on
xlabel("Number of satellites")
ylabel("Coverage Metric (%)")
hold on
% Plot the percentage of link availability of the region for the whole
% simulation duration
plot(numTotalSatellites,linkAvailability,"-*")
% Add legend to the plot
legend("Satellite Visibility","Link Availability", ...
    Location="southeast")
hold off
title("Coverage Metric As a Function of Number of Satellites")

Based on the average number of UEs in the region, plot the percentage of UEs with capacity greater than the threshold over the complete simulation time, whenever a link is available.

figure(2)
plot(numTotalSatellites,coveragePercent,"-*")
grid on
box on
xlabel("Number of satellites")
ylabel("Capacity > " + num2str(threshold) + " Mbps (% Number of users)")
title("Percentage of Users Exceeding Capacity Threshold vs. Number of Satellites")

Plot the link availability chart for each UE across all the timestamps for the first specified constellation.

figure(3)
% Scale the status with the UE number for better visualization
colors = colororder;
firstLineColor = colors(1,:);
ueIndex = (1:numUEs)';
linkAvailabilityFirstConst = results(1).LinkAvailability;
linkAvailabilityFirstConst(linkAvailabilityFirstConst == 0) = nan;
plot(timeSteps,linkAvailabilityFirstConst.*ueIndex, ...
    Color=firstLineColor,LineWidth=1)
xlim([timeSteps(1) timeSteps(end)])
ylim([0 numUEs+1])
xlabel("Time")
ylabel("UE Index")
title("Link Availability Chart")

Plot the regional coverage for the third specified constellation at the stop time on siteviewer, which supports RF propagation analysis and visualization.

capacityLastInstance = results(end).Capacity(:,end)./1e6;                % In Mbps
% Launch siteviewer for each example run to plot the capacity
siteviewer;
pd = propagationData(latCoord,lonCoord,"Capacity",capacityLastInstance);
legendTitle = "Capacity" + newline + "(Mbps)";
contour(pd,LegendTitle=legendTitle)

Compare your plotted regional coverage to this figure, which shows the regional coverage for the third specified constellation at the stop time, obtained by running 1000 UEs in the region with the default parameter settings. With a large number of UEs in the region, the coverage results show better resolution across the area.

Further Exploration

Using the coverage analysis over a particular region, you can extend this example in these ways:

• Compare satellite constellations with different orbital parameters.

• Map CNR with throughput for a standard-based system, such as narrowband Internet of Things (NB-IoT) and 5G New Radio.

• Analyze the impact of large user density in a region by increasing the number of UEs in the ROI.

• Model the impact of antenna pattern by using different antennas from Antenna Toolbox™ and Phased Array System Toolbox™ at the transmitter and receiver.

• Analyze coverage and capacity for other satellite-based applications not limited to NTN.

References

[1] Hexa-X. Final 6G Architectural Enablers and Technological Solutions. Deliverable D5.3. Edited by Mårten Ericson, Hannu Flinck, Panagiotis Vlacheas, and Stefan Wänstedt. Hexa-X, April 30, 2023.

[2] Pachler, Nils, Inigo del Portillo, Edward F. Crawley, and Bruce G. Cameron. “An Updated Comparison of Four Low Earth Orbit Satellite Constellation Systems to Provide Global Broadband.” In 2021 IEEE International Conference on Communications Workshops (ICC Workshops), 1–7, 2021. https://doi.org/10.1109/ICCWorkshops50388.2021.9473799.

Local Functions

function [factor1,factor2] = findClosestFactors(n)
% This function finds the two factors of n that are closest to each other

% Initialize the factors from 1 and n
factor1 = 1;
factor2 = n;

% Loop through possible factors up to the square root of n
for i = round(sqrt(n)):-1:1
    if mod(n,i) == 0  % If i is a factor
        factor1 = i;
        factor2 = n/i;
        break % Exit the loop as you find the closest factors
    end
end

end

function sat = addSatellites(sc,const)
% Add the satellites to the scenario with the input constellation

sat = [];
for i = 1:numel(const.ShellAltitude)
    % Get orbital elements for each satellite in the shell
    [semimajoraxis,eccentricity, ...
        inclination,RAAN,argofperiapsis, ...
        trueanomaly,numTotSat] = getOrbitalElements(const.ShellAltitude(i),const.ShellInclination(i), ...
        const.NumOrbits(i),const.NumSatellitesPerOrbit(i));

    % Add satellites to the scenario
    satellitesOrbit = satellite(sc, ...
        semimajoraxis,eccentricity,inclination,RAAN,argofperiapsis,trueanomaly, ...
        Name=(const.Name + " " + const.ShellNames(i) + " " + string(1:numTotSat)'));

    % Output the satellite objects
    sat = [sat(:); satellitesOrbit(:)];
end

end

function [semimajoraxis,eccentricity, ...
    inclination,RAAN,argofperiapsis, ...
    trueanomaly,numTotSat] = getOrbitalElements(shellAltitude,shellInclination, ...
    numOrbits,numSatellitesPerPlane)
% Get the orbital elements based on the shell altitude, shell inclination,
% number of orbits and number of satellites per orbital plane

% Number of satellites in the constellation
numTotSat = numOrbits*numSatellitesPerPlane;

% Get orbital elements
orbitIdx = repelem(1:numOrbits,1,numSatellitesPerPlane);
planeIdx = repmat(1:numSatellitesPerPlane,1,numOrbits);
RAAN = 180*(orbitIdx-1)/numOrbits;
trueanomaly = 360*(planeIdx-1 + 0.5*(mod(orbitIdx,2)-1)) ...
    /numSatellitesPerPlane;
semimajoraxis = repmat((6371 + shellAltitude)*1e3,1,numTotSat); % meters
inclination = repmat(shellInclination,1,numTotSat);             % degrees
eccentricity = zeros(1,numTotSat);                              % degrees
argofperiapsis = zeros(1,numTotSat);                            % degrees

end

function validateConstellation(constellation)
% Validates the dimensions of all fields in each constellation

for constIdx = 1:numel(constellation)
    validateattributes(constellation(constIdx).Name,["char","string"],"scalartext")
    numShellName = numel(constellation(constIdx).ShellNames);
    validateattributes(constellation(constIdx).ShellNames,"string","nonempty")
    validateattributes(constellation(constIdx).ShellAltitude,"double", ...
        {'positive','integer','numel',numShellName})
    validateattributes(constellation(constIdx).ShellInclination,"double", ...
        {'numel',numShellName})
    validateattributes(constellation(constIdx).NumOrbits,"double", ...
        {'positive','integer','numel',numShellName})
    validateattributes(constellation(constIdx).NumSatellitesPerOrbit,"double", ...
        {'positive','integer','numel',numShellName})
end

end