MultiHopPathSelection.m

% In a satellite constellation, a multi-hop path is a communication route
% where a data signal travels through multiple intermediate satellites 
% (nodes) to reach its destination, rather than directly from the source 
% to the target. This approach allows for greater coverage than a 
% single-hop path, connecting distant points or overcoming limitations in 
% a satellite's direct radio range by using other satellites as relays.

% This example shows how to determine the path through a large 
% constellation consisting of 1,000 low-Earth orbit (LEO) satellites to 
% gain access between two ground stations. Following this, it demonstrates 
% how to calculate the intervals during the next 3 hour period when this 
% path can be used.

% Here, the two ground stations that we'll be using are Tohoku Univ and 
% IEM Kolkata :P

startTime = datetime(2025,1,7,14,32,0); % 10 December 2021, 6:27:57 PM UTC
stopTime = startTime + hours(3);        % 10 December 2021, 9:27:57 PM UTC
sampleTime = 60;                        % Seconds
sc = satelliteScenario(startTime,stopTime,sampleTime,"AutoSimulate",false);

% Since the path must be determined only for the first time step, we'll set 
% AutoSimulate of the scenario to false to prevent it from automatically 
% advancing through the time steps until StopTime. When AutoSimulate is 
% false, SimulationStatus becomes available.

% Adding a large constellation of satellites; this consists of 1000 LEO
% satellites
sat = satellite(sc,"largeConstellation.tle");
numSatellites = numel(sat);

% Adding ground stations :)

gsSource = groundStation(sc,22.56263 ,88.36304 ... 
    ); % Latitude and longitude 
gsTarget = groundStation(sc,38.2538,140.8741 ...  
    ); % Latitude and longitude

% Determining elevation angles of ground station wrt ground stations

% Calculating the scenario state corresponding to StartTime.
advance(sc);

% Retrieving the elevation angle of each satellite with respect to the 
% ground stations.
[~,elSourceToSat] = aer(gsSource,sat);
[~,elTargetToSat] = aer(gsTarget,sat);

% Determine the elevation angles that are greater than or equal to 30
% degrees.
elSourceToSatGreaterThanOrEqual30 = (elSourceToSat >= 30)';
elTargetToSatGreaterThanOrEqual30 = (elTargetToSat >= 30)';

% Determining the best satellite for initial access to constellation

% The best satellite to be used for the initial access to the large 
% constellation is assumed to be the one that satisfies the following 
% simultaneously:
% 1.Has the closest range to "Target Ground Station".
% 2.Has an elevation angle of at least 30 degrees with respect to 
% "Source Ground Station".

% Finding the indices of the elements of elSourceToSatGreaterThanOrEqual30
% whose value is true.
trueID = find(elSourceToSatGreaterThanOrEqual30 == true);

% These indices are essentially the indices of satellites in sat whose
% elevation angle with respect to "Source Ground Station" is at least 30
% degrees.
% Determining the range of these satellites to "Target Ground
% Station".
[~,~,r] = aer(sat(trueID), gsTarget);

% Determining the index of the element in r bearing the minimum value.
[~,minRangeID] = min(r);

% Determining the element in trueID at the index minRangeID.
id = trueID(minRangeID);

% This is the index of the best satellite for initial access to the
% constellation. This will be the first hop in the path. We'll initialize a
% variable 'node' that stores the first two nodes of the routing - namely,
% "Source Ground Station" and the best satellite for initial constellation
% access.
nodes = {gsSource sat(id)};

% Determining remaining nodes in path to the ground station


% The remaining nodes in the path are determined using a similar logic as
% what was used for determining the first satellite. If the elevation angle
% of the satellite in the current node is already at least 30 degrees wrt 
% "Target Ground Station", the next node is "Target Ground Station", 
% thereby completing the path. Otherwise, the next node in the 
% constellation is chosen using a logic that is similar to what was used 
% for determining the first satellite in the path.
% The next node is the one that simultaneously satisfies the following:
%  1.Has the closest range to "Target Ground Station".
%  2.Elevation angle with respect to the satellite in the current node is 
%    greater than or equal to -15 degrees.

% The elevation value of -15 degrees is chosen because the horizon wrt
% each satellite in the constellation is about -21.9813 degrees. This value
% can be derived by assuming a spherical Earth geometry and the fact that 
% these satellites are in near-circular orbits at an altitude of roughly 
% 500 kilometers.

% Note:the spherical Earth assumption is used only for 
% computing the elevation angle of the horizon below. The satellite 
% scenario simulation itself assumes a WGS84 ellipsoid model for the Earth.

earthRadius = 6378137; % meters
altitude = 500000; % meters
horizonElevationAngle = asind(earthRadius/(earthRadius + altitude)) - 90;

% Any satellite whose elevation angle with respect to another satellite is
% greater than -21.9813 degrees is guaranteed to be visible to the latter. 
% However, choosing -15 degrees provides an adequate margin.

% Minimum elevation angle of satellite nodes with respect to the prior
% node.
minSatElevation = -15; % degrees

% Flag to specify if the complete multi-hop path has been found.
pathFound = false;

% Determine nodes of the path in a loop. Exit the loop once the complete
% multi-hop path has been found.
while ~pathFound
    % Index of the satellite in sat corresponding to current node is
    % updated to the value calculated as index for the next node in the
    % prior loop iteration. Essentially, the satellite in the next node in
    % prior iteration becomes the satellite in the current node in this
    % iteration.
    idCurrent = id;

    % This is the index of the element in elTargetToSatGreaterThanOrEqual30
    % tells if the elevation angle of this satellite is at least 30 degrees
    % with respect to "Target Ground Station". If this element is true, the
    % routing is complete, and the next node is the target ground station.
    if elTargetToSatGreaterThanOrEqual30(idCurrent)
        nodes = {nodes{:} gsTarget}; %#ok<CCAT> 
        pathFound = true;
        continue
    end

    % If the element is false, the path is not complete yet. The next node
    % in the path must be determined from the constellation. Determine
    % which satellites have elevation angle that is greater than or equal
    % to -15 degrees with respect to the current node. To do this, first
    % we'll determine the elevation angle of each satellite wrt the
    % current node.
    [~,els] = aer(sat(idCurrent),sat); 

    % Overwrite the elevation angle of the satellite with respect to itself
    % to be -90 degrees to ensure it does not get re-selected as the next
    % node.
    els(idCurrent) = -90; 

    % Determine the elevation angles that are greater than or equal to -15
    % degrees.
    s = els >= minSatElevation;

    % Find the indices of the elements in s whose value is true.
    trueID = find(s == true);

    % These indices are essentially the indices of satellites in sat whose
    % elevation angle with respect to the current node is greater than or
    % equal to -15 degrees. Determine the range of these satellites to
    % "Target Ground Station".
    [~,~,r] = aer(sat(trueID), gsTarget);

    % Determine the index of the element in r bearing the minimum value.
    [~,minRangeID] = min(r);

    % Determine the element in trueID at the index minRangeID.
    id = trueID(minRangeID);

    % This is the index of the best satellite among those in sat to be used
    % for the next node in the path. Append this satellite to the 'nodes'
    % variable.
    nodes = {nodes{:} sat(id)}; %#ok<CCAT>
end

% Determining intervals when calculated path can be used

% We must now determine the intervals over the next 3 hours during which 
% calculated path can be used. To accomplish this, manually stepping 
% through each time step of the scenario is not necessary. Instead, we can 
% make the scenario auto-simulate from StartTime to StopTime. Therefore, 
% we'll set AutoSimulate of the scenario to true.

sc.AutoSimulate = true;

% We'll add access analysis with calculated nodes in the path and set the 
% line colour to green
ac = access(nodes{:});
ac.LineColor = "green";

% Determining the access intervals using accessIntervals fn. Since 
% AutoSimulate has been set to true, the scenario will automatically 
% simulate from StartTime to StopTime at the specified SampleTime before 
% calculating the access intervals.These intervals are the times when the 
% calculated multi-hop path can be used.

intvls = accessIntervals(ac);

% Note that this does not necessarily result in the shortest overall 
% distance along the multi-hop path. One way to find the shortest path is 
% to find all possible paths using a graph search algorithm and then choose
% the shortest path among them. To find the distance between the different 
% satellites and ground stations, use the third output of the aer function.

% Additionally, this example computed the path only once for the scenario 
% StartTime, which was then used for the duration of the scenario. As a 
% result, the ground stations had access only for 5 minutes out of the 3 
% hours spanning the scenario duration. To increase the access duration 
% between the ground stations, you can modify the above code to re-compute 
% the path at each scenario time step using advance after setting 
% AutoSimulate to false. Calling advance advances SimulationTime by one 
% SampleTime. Once you compute all the paths, set AutoSimulate back to 
% true, create access objects corresponding to each unique path, and 
% re-compute the access intervals by calling accessIntervals on all the 
% access objects.