Constellation Modeling with the Orbit Propagator Block

This example shows how to propagate the orbits of a constellation of satellites and compute and visualize access intervals between the individual satellites and several ground stations. It uses:

• Aerospace Blockset `Orbit Propagator` block

• Aerospace Toolbox `satelliteScenario` object

The Aerospace Toolbox `satelliteScenario` object lets you load previously generated, time-stamped ephemeris data into a scenario from a timeseries or timetable object. Data is interpolated in the scenario object to align with the scenario time steps, allowing you to incorporate data generated in a Simulink model into either a new or existing `satelliteScenario` object. This example shows how to propagate a constellation of satellites in Simulink with the Aerospace Blockset `Orbit Propagator` block, and load the logged ephemeris data into a `satelliteScenario` object for access analysis.

Define Mission Parameters and Constellation Initial Conditions

Specify a start date and duration for the mission. This example uses MATLAB structures to organize mission data. These structures make accessing data later in the example more intuitive. They also help declutter the global base workspace.

```mission.StartDate = datetime(2020, 11, 30, 22, 23, 24); mission.Duration = hours(24);```

The constellation in this example is a Walker-Delta constellation modeled similar to Galileo, the European GNSS (Global navigation satellite system) constellation. The constellation consists of 24 satellites in medium Earth orbit (MEO).

Walker-Delta constellations use the notation:

`$\mathit{i}:\mathit{T}/\mathit{P}/\mathit{F}$`

where

$\mathit{i}=$ inclination

$\mathit{T}=$ total number of satellites

$\mathit{P}=$ number of equally space geometric planes

$\mathit{F}=$ spacing between satellites in adjacent planes

Walker-Delta constellations are a common solution for maximizing geometric coverage over Earth while minimizing the number of satellites required to perform the mission. The Galileo navigation system is a Walker-Delta ${56}^{\circ }:24/3/1$ constellation (24 satellites in 3 planes inclined at 56 degrees) in a 29599.8 km orbit.

Specify Keplerian orbital elements for the constellation at `mission.StartDate`.

```mission.Satellites.SemiMajorAxis = 29599.8e3 * ones(24,1); % meters mission.Satellites.Eccentricity = 0.0005 * ones(24,1); mission.Satellites.Inclination = 56 * ones(24,1); % deg mission.Satellites.ArgOfPeriapsis = 350 * ones(24,1); % deg```

The ascending nodes of the orbital planes of a Walker-Delta constellation are uniformly distributed at intervals of $\frac{{360}^{\circ }}{\mathit{P}}$ around the equator. The number of satellites per plane, S, is given as $\mathit{S}=\frac{\mathit{T}}{\mathit{P}}$. With 24 satellites total, this results in 3 planes of 8 satellites at 120 degree intervals around the equator. The satellites in each orbital plane are distributed at intervals of $\frac{{360}^{\circ }}{\mathit{S}}$, or 45 degrees.

```mission.Satellites.RAAN = sort(repmat([0 120 240], 1,8))'; % deg mission.Satellites.TrueAnomaly = repmat(0:45:315, 1,3)'; % deg```

Lastly, account for the relative angular shift between adjacent orbital planes. The phase difference is given as $\Delta \varphi =\mathit{F}*\frac{360}{\mathit{T}}$, or 15 degrees in this case.

```mission.Satellites.TrueAnomaly(9:16) = mission.Satellites.TrueAnomaly(9:16) + 15; mission.Satellites.TrueAnomaly(17:24) = mission.Satellites.TrueAnomaly(17:24) + 30;```

Show the constellation nodes in a table.

```ConstellationDefinition = table(mission.Satellites.SemiMajorAxis, ... mission.Satellites.Eccentricity, ... mission.Satellites.Inclination, ... mission.Satellites.RAAN, ... mission.Satellites.ArgOfPeriapsis, ... mission.Satellites.TrueAnomaly, ... 'VariableNames', ["a (m)", "e", "i (deg)", "Ω (deg)", "ω (deg)", "ν (deg)"])```
```ConstellationDefinition=24×6 table a (m) e i (deg) Ω (deg) ω (deg) ν (deg) ________ ______ _______ _______ _______ _______ 2.96e+07 0.0005 56 0 350 0 2.96e+07 0.0005 56 0 350 45 2.96e+07 0.0005 56 0 350 90 2.96e+07 0.0005 56 0 350 135 2.96e+07 0.0005 56 0 350 180 2.96e+07 0.0005 56 0 350 225 2.96e+07 0.0005 56 0 350 270 2.96e+07 0.0005 56 0 350 315 2.96e+07 0.0005 56 120 350 15 2.96e+07 0.0005 56 120 350 60 2.96e+07 0.0005 56 120 350 105 2.96e+07 0.0005 56 120 350 150 2.96e+07 0.0005 56 120 350 195 2.96e+07 0.0005 56 120 350 240 2.96e+07 0.0005 56 120 350 285 2.96e+07 0.0005 56 120 350 330 ⋮ ```

Open and Configure the Orbit Propagation Model

Open the included Simulink model. This model contains an `Orbit Propagator` block connected to output ports. The` Orbit Propagator` block supports vectorization. This allows you to model multiple satellites in a single block by specifying arrays of initial conditions in the `Block Parameters` window or using `set_param`. The model also includes a "Mission Analysis and Visualization" section that contains a dashboard `Callback button`. When clicked, this button runs the model, creates a new satelliteScenario object in the global base workspace containing the satellite or constellation defined in the `Orbit Propagator` block, and opens a Satellite Scenario Viewer window for the new scenario. To view the source code for this action, double click the callback button. The "Mission Analysis and Visualization" section is a standalone workflow to create a new satelliteScenario object and is not used as part of this written example.

```mission.mdl = "OrbitPropagatorBlockExampleModel"; open_system(mission.mdl);```

Define the path to the `Orbit Propagator` block in the model.

`mission.Satellites.blk = mission.mdl + "/Orbit Propagator";`

Set satellite initial conditions. To assign the Keplerian orbital element set defined in the previous section, use `set_param`.

```set_param(mission.Satellites.blk, ... startDate = num2str(juliandate(mission.StartDate)), ... stateFormatNum = "Orbital elements", ... orbitType = "Keplerian", ... semiMajorAxis = "mission.Satellites.SemiMajorAxis", ... eccentricity = "mission.Satellites.Eccentricity", ... inclination = "mission.Satellites.Inclination", ... raan = "mission.Satellites.RAAN", ... argPeriapsis = "mission.Satellites.ArgOfPeriapsis", ... trueAnomaly = "mission.Satellites.TrueAnomaly");```

Set the position and velocity output ports of the block to use the Earth-centered Earth-fixed frame, which is the International Terrestrial Reference Frame (ITRF).

```set_param(mission.Satellites.blk, ... centralBody = "Earth", ... outportFrame = "Fixed-frame");```

Configure the propagator. This example uses the `Oblate ellipsoid (J2)` propagator which includes second order zonal harmonic perturbations in the satellite trajectory calculations, accounting for the oblateness of Earth.

```set_param(mission.Satellites.blk, ... propagator = "Numerical (high precision)", ... gravityModel = "Oblate ellipsoid (J2)", ... useEOPs = "off");```

Apply model-level solver setting using `set_param`. For best performance and accuracy when using a numerical propagator, use a variable-step solver.

```set_param(mission.mdl, ... SolverType = "Variable-step", ... SolverName = "VariableStepAuto", ... RelTol = "1e-6", ... AbsTol = "1e-7", ... StopTime = string(seconds(mission.Duration)));```

Save model output port data as a dataset of time series objects.

```set_param(mission.mdl, ... SaveOutput = "on", ... OutputSaveName = "yout", ... SaveFormat = "Dataset");```

Run the Model and Collect Satellite Ephemerides

Simulate the model.

`mission.SimOutput = sim(mission.mdl);`

Extract position and velocity data from the model output data structure.

```mission.Satellites.TimeseriesPosECEF = mission.SimOutput.yout{1}.Values; mission.Satellites.TimeseriesVelECEF = mission.SimOutput.yout{2}.Values;```

Set the start data from the mission in the timeseries object.

```mission.Satellites.TimeseriesPosECEF.TimeInfo.StartDate = mission.StartDate; mission.Satellites.TimeseriesVelECEF.TimeInfo.StartDate = mission.StartDate;```

The timeseries objects contain position and velocity data for all 24 satellites.

`mission.Satellites.TimeseriesPosECEF`
``` timeseries Common Properties: Name: '' Time: [57x1 double] TimeInfo: [1x1 tsdata.timemetadata] Data: [24x3x57 double] DataInfo: [1x1 tsdata.datametadata] More properties, Methods ```

Load the Satellite Ephemerides into a `satelliteScenario O`bject

Create a satellite scenario object for the analysis.

`scenario = satelliteScenario(mission.StartDate, mission.StartDate + hours(24), 60);`

Add all 24 satellites to the satellite scenario from the ECEF position and velocity timeseries objects using the `satellite` method.

```sat = satellite(scenario, mission.Satellites.TimeseriesPosECEF, mission.Satellites.TimeseriesVelECEF, ... CoordinateFrame="ecef", Name="GALILEO " + (1:24))```
```sat = 1x24 Satellite array with properties: Name ID ConicalSensors Gimbals Transmitters Receivers Accesses GroundTrack Orbit OrbitPropagator MarkerColor MarkerSize ShowLabel LabelFontColor LabelFontSize ```
`disp(scenario)`
``` satelliteScenario with properties: StartTime: 30-Nov-2020 22:23:24 StopTime: 01-Dec-2020 22:23:24 SampleTime: 60 AutoSimulate: 1 Satellites: [1×24 matlabshared.satellitescenario.Satellite] GroundStations: [1×0 matlabshared.satellitescenario.GroundStation] Viewers: [0×0 matlabshared.satellitescenario.Viewer] AutoShow: 1 ```

Set Graphical Properties on the Satellites

Set satellites in each orbital plane to have the same orbit color.

```set(sat(1:8), MarkerColor="#FF6929"); set(sat(9:16), MarkerColor="#139FFF"); set(sat(17:24), MarkerColor="#64D413"); orbit = [sat(:).Orbit]; set(orbit(1:8), LineColor="#FF6929"); set(orbit(9:16), LineColor="#139FFF"); set(orbit(17:24), LineColor="#64D413");```

To provide accurate positioning data, a location on Earth must have access to at least 4 satellites in the constellation at any given time. In this example, use three MathWorks locations to compare total constellation access over the 1 day analysis window to different regions of Earth:

• Natick, Massachusetts, USA (42.30048°, -71.34908°)

• München, Germany (48.23206°, 11.68445°)

• Bangalore, India (12.94448°, 77.69256°)

```gsUS = groundStation(scenario, 42.30048, -71.34908, ... MinElevationAngle=10, Name="Natick"); gsUS.MarkerColor = "red"; gsDE = groundStation(scenario, 48.23206, 11.68445, ... MinElevationAngle=10, Name="Munchen"); gsDE.MarkerColor = "red"; gsIN = groundStation(scenario, 12.94448, 77.69256, ... MinElevationAngle=10, Name="Bangalore"); gsIN.MarkerColor = "red"; figure geoscatter([gsUS.Latitude gsDE.Latitude gsIN.Latitude], ... [gsUS.Longitude gsDE.Longitude gsIN.Longitude], "red", "filled") geolimits([-75 75], [-180 180]) title("Ground Stations")```

Compute Ground Station to Satellite Access (Line-of-Sight Visibility)

Calculate line-of-sight access between the ground stations and each individual satellite using the `access` method.

```accessUS = access(gsUS, sat); accessDE = access(gsDE, sat); accessIN = access(gsIN, sat);```

Set access colors to match orbital plane colors assigned earlier in the example.

```set(accessUS, LineWidth="1"); set(accessUS(1:8), LineColor="#FF6929"); set(accessUS(9:16), LineColor="#139FFF"); set(accessUS(17:24), LineColor="#64D413"); set(accessDE, LineWidth="1"); set(accessDE(1:8), LineColor="#FF6929"); set(accessDE(9:16), LineColor="#139FFF"); set(accessDE(17:24), LineColor="#64D413"); set(accessIN, LineWidth="1"); set(accessIN(1:8), LineColor="#FF6929"); set(accessIN(9:16), LineColor="#139FFF"); set(accessIN(17:24), LineColor="#64D413");```

View the full access table between each ground station and all satellites in the constellation as tables. Sort the access intervals by interval start time. Satellites added from ephemeris data do not display values for StartOrbit and Stop orbit.

```intervalsUS = accessIntervals(accessUS); intervalsUS = sortrows(intervalsUS, "StartTime", "ascend")```
```intervalsUS=40×8 table Source Target IntervalNumber StartTime EndTime Duration StartOrbit EndOrbit ________ ____________ ______________ ____________________ ____________________ ________ __________ ________ "Natick" "GALILEO 1" 1 30-Nov-2020 22:23:24 01-Dec-2020 04:04:24 20460 NaN NaN "Natick" "GALILEO 2" 1 30-Nov-2020 22:23:24 01-Dec-2020 01:24:24 10860 NaN NaN "Natick" "GALILEO 3" 1 30-Nov-2020 22:23:24 30-Nov-2020 22:57:24 2040 NaN NaN "Natick" "GALILEO 12" 1 30-Nov-2020 22:23:24 01-Dec-2020 00:00:24 5820 NaN NaN "Natick" "GALILEO 13" 1 30-Nov-2020 22:23:24 30-Nov-2020 23:05:24 2520 NaN NaN "Natick" "GALILEO 18" 1 30-Nov-2020 22:23:24 01-Dec-2020 04:00:24 20220 NaN NaN "Natick" "GALILEO 19" 1 30-Nov-2020 22:23:24 01-Dec-2020 01:42:24 11940 NaN NaN "Natick" "GALILEO 20" 1 30-Nov-2020 22:23:24 30-Nov-2020 22:46:24 1380 NaN NaN "Natick" "GALILEO 11" 1 30-Nov-2020 22:25:24 01-Dec-2020 00:18:24 6780 NaN NaN "Natick" "GALILEO 17" 1 30-Nov-2020 22:50:24 01-Dec-2020 05:50:24 25200 NaN NaN "Natick" "GALILEO 8" 1 30-Nov-2020 23:20:24 01-Dec-2020 07:09:24 28140 NaN NaN "Natick" "GALILEO 7" 1 01-Dec-2020 01:26:24 01-Dec-2020 10:00:24 30840 NaN NaN "Natick" "GALILEO 24" 1 01-Dec-2020 01:40:24 01-Dec-2020 07:12:24 19920 NaN NaN "Natick" "GALILEO 14" 1 01-Dec-2020 03:56:24 01-Dec-2020 07:15:24 11940 NaN NaN "Natick" "GALILEO 6" 1 01-Dec-2020 04:05:24 01-Dec-2020 12:14:24 29340 NaN NaN "Natick" "GALILEO 23" 1 01-Dec-2020 04:10:24 01-Dec-2020 08:03:24 13980 NaN NaN ⋮ ```
```intervalsDE = accessIntervals(accessDE); intervalsDE = sortrows(intervalsDE, "StartTime", "ascend")```
```intervalsDE=40×8 table Source Target IntervalNumber StartTime EndTime Duration StartOrbit EndOrbit _________ ____________ ______________ ____________________ ____________________ ________ __________ ________ "Munchen" "GALILEO 2" 1 30-Nov-2020 22:23:24 01-Dec-2020 04:34:24 22260 NaN NaN "Munchen" "GALILEO 3" 1 30-Nov-2020 22:23:24 01-Dec-2020 01:58:24 12900 NaN NaN "Munchen" "GALILEO 4" 1 30-Nov-2020 22:23:24 30-Nov-2020 23:05:24 2520 NaN NaN "Munchen" "GALILEO 10" 1 30-Nov-2020 22:23:24 30-Nov-2020 23:58:24 5700 NaN NaN "Munchen" "GALILEO 19" 1 30-Nov-2020 22:23:24 01-Dec-2020 01:36:24 11580 NaN NaN "Munchen" "GALILEO 20" 1 30-Nov-2020 22:23:24 01-Dec-2020 00:15:24 6720 NaN NaN "Munchen" "GALILEO 21" 1 30-Nov-2020 22:23:24 30-Nov-2020 22:28:24 300 NaN NaN "Munchen" "GALILEO 9" 1 30-Nov-2020 22:34:24 01-Dec-2020 02:22:24 13680 NaN NaN "Munchen" "GALILEO 18" 1 30-Nov-2020 22:41:24 01-Dec-2020 02:31:24 13800 NaN NaN "Munchen" "GALILEO 1" 1 30-Nov-2020 23:05:24 01-Dec-2020 06:42:24 27420 NaN NaN "Munchen" "GALILEO 16" 1 30-Nov-2020 23:29:24 01-Dec-2020 04:47:24 19080 NaN NaN "Munchen" "GALILEO 15" 1 01-Dec-2020 00:50:24 01-Dec-2020 07:27:24 23820 NaN NaN "Munchen" "GALILEO 17" 1 01-Dec-2020 01:05:24 01-Dec-2020 03:00:24 6900 NaN NaN "Munchen" "GALILEO 8" 1 01-Dec-2020 01:57:24 01-Dec-2020 08:25:24 23280 NaN NaN "Munchen" "GALILEO 14" 1 01-Dec-2020 02:36:24 01-Dec-2020 10:19:24 27780 NaN NaN "Munchen" "GALILEO 7" 1 01-Dec-2020 04:35:24 01-Dec-2020 09:43:24 18480 NaN NaN ⋮ ```
```intervalsIN = accessIntervals(accessIN); intervalsIN = sortrows(intervalsIN, "StartTime", "ascend")```
```intervalsIN=31×8 table Source Target IntervalNumber StartTime EndTime Duration StartOrbit EndOrbit ___________ ____________ ______________ ____________________ ____________________ ________ __________ ________ "Bangalore" "GALILEO 3" 1 30-Nov-2020 22:23:24 01-Dec-2020 05:12:24 24540 NaN NaN "Bangalore" "GALILEO 4" 1 30-Nov-2020 22:23:24 01-Dec-2020 02:59:24 16560 NaN NaN "Bangalore" "GALILEO 5" 1 30-Nov-2020 22:23:24 01-Dec-2020 00:22:24 7140 NaN NaN "Bangalore" "GALILEO 9" 1 30-Nov-2020 22:23:24 01-Dec-2020 03:37:24 18840 NaN NaN "Bangalore" "GALILEO 10" 1 30-Nov-2020 22:23:24 01-Dec-2020 00:09:24 6360 NaN NaN "Bangalore" "GALILEO 16" 1 30-Nov-2020 22:23:24 01-Dec-2020 08:44:24 37260 NaN NaN "Bangalore" "GALILEO 21" 1 30-Nov-2020 22:23:24 30-Nov-2020 23:25:24 3720 NaN NaN "Bangalore" "GALILEO 22" 1 30-Nov-2020 22:23:24 30-Nov-2020 22:58:24 2100 NaN NaN "Bangalore" "GALILEO 15" 1 01-Dec-2020 00:17:24 01-Dec-2020 11:16:24 39540 NaN NaN "Bangalore" "GALILEO 2" 1 01-Dec-2020 00:25:24 01-Dec-2020 07:10:24 24300 NaN NaN "Bangalore" "GALILEO 22" 2 01-Dec-2020 00:48:24 01-Dec-2020 05:50:24 18120 NaN NaN "Bangalore" "GALILEO 21" 2 01-Dec-2020 01:32:24 01-Dec-2020 08:29:24 25020 NaN NaN "Bangalore" "GALILEO 1" 1 01-Dec-2020 03:06:24 01-Dec-2020 07:17:24 15060 NaN NaN "Bangalore" "GALILEO 20" 1 01-Dec-2020 03:36:24 01-Dec-2020 12:38:24 32520 NaN NaN "Bangalore" "GALILEO 14" 1 01-Dec-2020 05:48:24 01-Dec-2020 13:29:24 27660 NaN NaN "Bangalore" "GALILEO 19" 1 01-Dec-2020 05:53:24 01-Dec-2020 17:06:24 40380 NaN NaN ⋮ ```

View the Satellite Scenario

Open a 3-D viewer window of the scenario. The viewer window contains all 24 satellites and the three ground stations defined earlier in this example. A line is drawn between each ground station and satellite during their corresponding access intervals. Hide the details of the satellites and ground stations by setting the `ShowDetails` name-value pair to `false`. Show satellite orbits and labels for the ground station locations.

```viewer3D = satelliteScenarioViewer(scenario, ShowDetails=false); show(sat.Orbit); gsUS.ShowLabel = true; gsUS.LabelFontSize = 11; gsDE.ShowLabel = true; gsDE.LabelFontSize = 11; gsIN.ShowLabel = true; gsIN.LabelFontSize = 11;```

Compare Access Between Ground Stations

Calculate access status between each satellite and ground station using the accessStatus method. Each row of the output array corresponds with a satellite in the constellation. Each column corresponds with time steps in the scenario. A value of `True` indicates that the satellite can access the aircraft at that specific time sample. The second output of `accessStatus` contains the time steps of the scenario. Plot cumulative access for each ground station over the one day analysis window.

```[statusUS, timeSteps] = accessStatus(accessUS); statusDE = accessStatus(accessDE); statusIN = accessStatus(accessIN); % Sum cumulative access at each timestep statusUS = sum(statusUS, 1); statusDE = sum(statusDE, 1); statusIN = sum(statusIN, 1); subplot(3,1,1); stairs(timeSteps, statusUS); title("Natick to GALILEO") ylabel("# of satellites") subplot(3,1,2); stairs(timeSteps, statusDE); title("München to GALILEO") ylabel("# of satellites") subplot(3,1,3); stairs(timeSteps, statusIN); title("Bangalore to GALILEO") ylabel("# of satellites")```

Collect access interval metrics for each ground station in a table for comparison.

```statusTable = [table(height(intervalsUS), height(intervalsDE), height(intervalsIN)); ... table(sum(intervalsUS.Duration)/3600, sum(intervalsDE.Duration)/3600, sum(intervalsIN.Duration)/3600); ... table(mean(intervalsUS.Duration/60), mean(intervalsDE.Duration/60), mean(intervalsIN.Duration/60)); ... table(mean(statusUS, 2), mean(statusDE, 2), mean(statusIN, 2)); ... table(min(statusUS), min(statusDE), min(statusIN)); ... table(max(statusUS), max(statusDE), max(statusIN))]; statusTable.Properties.VariableNames = ["Natick", "München", "Bangalore"]; statusTable.Properties.RowNames = ["Total # of intervals", "Total interval time (hrs)",... "Mean interval length (min)", "Mean # of satellites in view", ... "Min # of satellites in view", "Max # of satellites in view"]; statusTable```
```statusTable=6×3 table Natick München Bangalore ______ _______ _________ Total # of intervals 40 40 31 Total interval time (hrs) 167.88 169.95 180.42 Mean interval length (min) 251.82 254.93 349.19 Mean # of satellites in view 7.018 7.1041 7.5337 Min # of satellites in view 5 5 5 Max # of satellites in view 9 10 9 ```

Walker-Delta constellations are evenly distributed across longitudes. Natick and München are located at similar latitudes, and therefore have very similar access characteristics with respect to the constellation. Bangalore is at a latitude closer to the equator, and despite having a lower number of individual access intervals, it has the highest average number of satellites in view, the highest overall interval time, and the longest average interval duration (by about 95 minutes). All locations always have at least 4 satellites in view, as is required for GNSS trilateration.

References

[1] Wertz, James R, David F. Everett, and Jeffery J. Puschell. Space Mission Engineering: The New Smad. Hawthorne, CA: Microcosm Press, 2011. Print.

[2] Beech, Theresa W., Sefania Cornana, Miguel B. Mora. A Study of Three Satellite Constellation Design Algorithms. 14th International Symposium on Space Flight Dynamics, Foz do Iguaçu, Brazil 1999.

[3] The European Space Agency: Galileo Facts and Figures. https://www.esa.int/Applications/Navigation/Galileo/Facts_and_figures