Modeling Satellite Constellations using Ephemeris Data
This example demonstrates how to add time-stamped ephemeris data for a constellation of 24 satellites (similar to ESA Galileo GNSS constellation) to a satellite scenario for access analysis. The example uses data generated by the Aerospace Blockset Orbit Propagator block. For more information, see the Aerospace Blockset example Constellation Modeling with the Orbit Propagator Block.
The satelliteScenario object supports loading previously generated, time-stamped satellite ephemeris data into a scenario from a timeseries or timetable object. An ephemeris is a table containing position (and optionally velocity) state information of a satellite during a given period of time. Ephemeris data used to add satellites to the scenario object is interpolated via the makima interpolation method to align with the scenario time steps. This allows you to incorporate data generated by a Simulink model into either a new or existing satelliteScenario.
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 the ESA Galileo GNSS (Global navigation satellite system) constellation. The constellation consists of 24 satellites in medium Earth orbit (MEO). The satellites' Keplerian orbital elements at the mission start date epoch are:
mission.ConstellationDefinition = table( ... 23222e3 * ones(24,1), ... % Semi-major axis (m) 0.0005 * ones(24,1), ... % Eccentricity 56 * ones(24,1), ... % Inclination (deg) 350 * ones(24,1), ... % Right ascension of the ascending node (deg) sort(repmat([0 120 240], 1,8))', ... % Argument of periapsis (deg) [0:45:315, 15:45:330, 30:45:345]', ... % True anomoly (deg) 'VariableNames', ["a (m)", "e", "i (deg)", "Ω (deg)", "ω (deg)", "ν (deg)"]); mission.ConstellationDefinition
ans=24×6 table
a (m) e i (deg) Ω (deg) ω (deg) ν (deg)
__________ ______ _______ _______ _______ _______
2.3222e+07 0.0005 56 350 0 0
2.3222e+07 0.0005 56 350 0 45
2.3222e+07 0.0005 56 350 0 90
2.3222e+07 0.0005 56 350 0 135
2.3222e+07 0.0005 56 350 0 180
2.3222e+07 0.0005 56 350 0 225
2.3222e+07 0.0005 56 350 0 270
2.3222e+07 0.0005 56 350 0 315
2.3222e+07 0.0005 56 350 120 15
2.3222e+07 0.0005 56 350 120 60
2.3222e+07 0.0005 56 350 120 105
2.3222e+07 0.0005 56 350 120 150
2.3222e+07 0.0005 56 350 120 195
2.3222e+07 0.0005 56 350 120 240
2.3222e+07 0.0005 56 350 120 285
2.3222e+07 0.0005 56 350 120 330
⋮
Load Ephemeris Timeseries Data
The timeseries objects contain position and velocity data for all 24 satellites in the constellation. The data is referenced in the International Terrestrial Reference frame (ITRF), which is an Earth-centered Earth-fixed (ECEF) coordinate system. The data was generated using the Aerospace Blockset Orbit Propagator block. For more information, see the Aerospace Blockset example Constellation Modeling with the Orbit Propagator Block.
mission.Ephemeris = load("SatelliteScenarioEphemerisData.mat", "TimeseriesPosITRF", "TimeseriesVelITRF"); mission.Ephemeris.TimeseriesPosITRF
timeseries
Common Properties:
Name: ''
Time: [79x1 double]
TimeInfo: [1x1 tsdata.timemetadata]
Data: [24x3x79 double]
DataInfo: [1x1 tsdata.datametadata]
More properties, Methods
mission.Ephemeris.TimeseriesVelITRF
timeseries
Common Properties:
Name: ''
Time: [79x1 double]
TimeInfo: [1x1 tsdata.timemetadata]
Data: [24x3x79 double]
DataInfo: [1x1 tsdata.datametadata]
More properties, Methods
Load the Satellite Ephemerides into a satelliteScenario Object
Create a satellite scenario object for the analysis.
scenario = satelliteScenario(mission.StartDate, mission.StartDate + hours(24), 60);
Use the satellite method to add all 24 satellites to the satellite scenario from the ECEF position and velocity timeseries objects. This example uses position and velocity information; however satellites can also be added from position data only and velocity states are then estimated. Available coordinate frames for Name-Value pair CoordinateFrame are "ECEF", "Inertial", and "Geographic". If the timeseries object contains a value for ts.TimeInfo.StartDate, the method uses that value as the epoch for the timeseries object. If no StartDate is defined, the method uses the scenario start date by default.
sat = satellite(scenario, mission.Ephemeris.TimeseriesPosITRF, mission.Ephemeris.TimeseriesVelITRF, ... "CoordinateFrame", "ecef", "Name", "GALILEO " + (1:24))
sat =
1×24 Satellite array with properties:
Name
ID
ConicalSensors
Gimbals
Transmitters
Receivers
Accesses
GroundTrack
Orbit
OrbitPropagator
MarkerColor
MarkerSize
ShowLabel
LabelFontSize
LabelFontColor
disp(scenario)
satelliteScenario with properties:
StartTime: 30-Nov-2020 22:23:24
StopTime: 01-Dec-2020 22:23:24
SampleTime: 60
Viewers: [0×0 matlabshared.satellitescenario.Viewer]
Satellites: [1×24 matlabshared.satellitescenario.Satellite]
GroundStations: []
AutoShow: 1
Alternatively, satellites can also be added as ephemerides to the satellite scenario as a MATLAB timetable, table, or tscollection. For example, a timetable containing the first 3 satellites of the position timeseries object in the previous section, formatted for use with satelliteScenario objects is shown below.
Satellites are represented by variables (column headers).
Each row contains a position vector associated with the row's
Timeproperty.
timetable(... datetime(getabstime(mission.Ephemeris.TimeseriesPosITRF)), ... squeeze(mission.Ephemeris.TimeseriesPosITRF.Data(1,:,:))', ... squeeze(mission.Ephemeris.TimeseriesPosITRF.Data(2,:,:))', ... squeeze(mission.Ephemeris.TimeseriesPosITRF.Data(3,:,:))',... 'VariableNames', ["Satellite_1", "Satellite_2", "Satellite_3"])
ans=79×3 timetable
Time Satellite_1 Satellite_2 Satellite_3
____________________ ________________________________________ ________________________________________ _________________________________________
30-Nov-2020 22:23:24 1.4317e+07 -1.7969e+07 -3.2957e+06 1.8576e+07 -8.4335e+06 1.1077e+07 1.1955e+07 6.0469e+06 1.8967e+07
30-Nov-2020 22:23:34 1.4326e+07 -1.7968e+07 -3.2607e+06 1.8567e+07 -8.4162e+06 1.1106e+07 1.1934e+07 6.071e+06 1.8974e+07
30-Nov-2020 22:24:26 1.4367e+07 -1.7966e+07 -3.0856e+06 1.8519e+07 -8.3288e+06 1.125e+07 1.1824e+07 6.192e+06 1.9003e+07
30-Nov-2020 22:28:44 1.4562e+07 -1.7939e+07 -2.2062e+06 1.8271e+07 -7.878e+06 1.1958e+07 1.128e+07 6.8025e+06 1.9125e+07
30-Nov-2020 22:40:57 1.4983e+07 -1.7724e+07 3.0748e+05 1.7492e+07 -6.472e+06 1.3824e+07 9.757e+06 8.5761e+06 1.9251e+07
30-Nov-2020 22:54:44 1.5245e+07 -1.7219e+07 3.1374e+06 1.6506e+07 -4.6682e+06 1.5645e+07 8.1023e+06 1.0623e+07 1.8999e+07
30-Nov-2020 23:08:36 1.5308e+07 -1.6416e+07 5.9118e+06 1.5446e+07 -2.6463e+06 1.7132e+07 6.541e+06 1.2681e+07 1.8329e+07
30-Nov-2020 23:24:15 1.5183e+07 -1.5144e+07 8.8865e+06 1.4218e+07 -1.4887e+05 1.8358e+07 4.9309e+06 1.4942e+07 1.7089e+07
30-Nov-2020 23:40:33 1.4888e+07 -1.3412e+07 1.1718e+07 1.2964e+07 2.6318e+06 1.9086e+07 3.4517e+06 1.7146e+07 1.529e+07
30-Nov-2020 23:57:25 1.4477e+07 -1.1213e+07 1.4269e+07 1.1748e+07 5.6176e+06 1.923e+07 2.1441e+06 1.9174e+07 1.2944e+07
01-Dec-2020 00:15:24 1.3994e+07 -8.4611e+06 1.6479e+07 1.0588e+07 8.8233e+06 1.8695e+07 9.861e+05 2.096e+07 9.9781e+06
01-Dec-2020 00:33:57 1.3527e+07 -5.2723e+06 1.8119e+07 9.5707e+06 1.2021e+07 1.742e+07 14938 2.2298e+07 6.5344e+06
01-Dec-2020 00:53:11 1.3142e+07 -1.7143e+06 1.9066e+07 8.7131e+06 1.5075e+07 1.5378e+07 -8.1352e+05 2.3065e+07 2.6968e+06
01-Dec-2020 01:12:46 1.2901e+07 2.0218e+06 1.9203e+07 8.0192e+06 1.7771e+07 1.2634e+07 -1.5539e+06 2.3146e+07 -1.3246e+06
01-Dec-2020 01:32:20 1.2826e+07 5.7119e+06 1.8499e+07 7.449e+06 1.9927e+07 9.3368e+06 -2.2854e+06 2.251e+07 -5.2889e+06
01-Dec-2020 01:51:55 1.2899e+07 9.1955e+06 1.6986e+07 6.9227e+06 2.1452e+07 5.6297e+06 -3.1024e+06 2.1184e+07 -9.0243e+06
⋮
Set Graphical Properties on the Satellites
Viewer windows with many satellites can become crowded and difficult to read. To keep the window readable, manually control graphical properties of the scenario elements.
Hide the satellite labels and ground tracks.
set(sat, "ShowLabel", false);
hide([sat(:).GroundTrack]);Set satellite in the same orbital plane to have the same orbit color.
set(sat(1:8), "MarkerColor", "red"); set(sat(9:16), "MarkerColor", "blue"); set(sat(17:24), "MarkerColor", "green"); orbit = [sat(:).Orbit]; set(orbit(1:8), "LineColor", "red"); set(orbit(9:16), "LineColor", "blue"); set(orbit(17:24), "LineColor", "green");
Add Ground Stations to Scenario
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:
MathWorks headquarters in Natick, Massachusetts, USA (42.30048°, -71.34908°)
MathWorks München in Ismaning (near München), Germany (48.23206°, 11.68445°)
MathWorks Bangalore in Bangalore, India (12.94448°, 77.69256°)
gsUS = groundStation(scenario, 42.30048, -71.34908, ... "MinElevationAngle", 10, "Name", "MW Natick"); gsDE = groundStation(scenario, 48.23206, 11.68445, ... "MinElevationAngle", 10, "Name", "MW Munchen"); gsIN = groundStation(scenario, 12.94448, 77.69256, ... "MinElevationAngle", 10, "Name", "MW Bangalore"); figure geoscatter([gsUS.Latitude gsDE.Latitude gsIN.Latitude], ... [gsUS.Longitude gsDE.Longitude gsIN.Longitude], "red", "filled") geolimits([-75 75], [-180 180]) title("MathWorks 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.
for idx = 1:numel(sat) access(gsUS, sat(idx)); access(gsDE, sat(idx)); access(gsIN, sat(idx)); end accessUS = [gsUS(:).Accesses]; accessDE = [gsDE(:).Accesses]; accessIN = [gsIN(:).Accesses];
Set access colors to match orbital plane colors assigned earlier in the example.
set(accessUS(1:8), "LineColor", "red"); set(accessUS(9:16), "LineColor", "blue"); set(accessUS(17:24), "LineColor", "green"); set(accessDE(1:8), "LineColor", "red"); set(accessDE(9:16), "LineColor", "blue"); set(accessDE(17:24), "LineColor", "green"); set(accessIN(1:8), "LineColor", "red"); set(accessIN(9:16), "LineColor", "blue"); set(accessIN(17:24), "LineColor", "green");
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 EndOrbit.
intervalsUS = accessIntervals(accessUS); intervalsUS = sortrows(intervalsUS, "StartTime", "ascend")
intervalsUS=52×8 table
Source Target IntervalNumber StartTime EndTime Duration StartOrbit EndOrbit
___________ ____________ ______________ ____________________ ____________________ ________ __________ ________
"MW Natick" "GALILEO 1" 1 30-Nov-2020 22:23:24 01-Dec-2020 01:40:24 11820 NaN NaN
"MW Natick" "GALILEO 2" 1 30-Nov-2020 22:23:24 01-Dec-2020 00:08:24 6300 NaN NaN
"MW Natick" "GALILEO 3" 1 30-Nov-2020 22:23:24 30-Nov-2020 22:34:24 660 NaN NaN
"MW Natick" "GALILEO 12" 1 30-Nov-2020 22:23:24 30-Nov-2020 23:35:24 4320 NaN NaN
"MW Natick" "GALILEO 13" 1 30-Nov-2020 22:23:24 30-Nov-2020 22:47:24 1440 NaN NaN
"MW Natick" "GALILEO 18" 1 30-Nov-2020 22:23:24 01-Dec-2020 01:51:24 12480 NaN NaN
"MW Natick" "GALILEO 19" 1 30-Nov-2020 22:23:24 01-Dec-2020 00:13:24 6600 NaN NaN
"MW Natick" "GALILEO 20" 1 30-Nov-2020 22:23:24 30-Nov-2020 22:27:24 240 NaN NaN
"MW Natick" "GALILEO 11" 1 30-Nov-2020 22:38:24 01-Dec-2020 00:00:24 4920 NaN NaN
"MW Natick" "GALILEO 17" 1 30-Nov-2020 22:47:24 01-Dec-2020 03:19:24 16320 NaN NaN
"MW Natick" "GALILEO 8" 1 30-Nov-2020 23:09:24 01-Dec-2020 03:15:24 14760 NaN NaN
"MW Natick" "GALILEO 7" 1 01-Dec-2020 00:28:24 01-Dec-2020 04:55:24 16020 NaN NaN
"MW Natick" "GALILEO 24" 1 01-Dec-2020 00:29:24 01-Dec-2020 04:40:24 15060 NaN NaN
"MW Natick" "GALILEO 6" 1 01-Dec-2020 01:54:24 01-Dec-2020 06:41:24 17220 NaN NaN
"MW Natick" "GALILEO 23" 1 01-Dec-2020 02:04:24 01-Dec-2020 05:52:24 13680 NaN NaN
"MW Natick" "GALILEO 5" 1 01-Dec-2020 03:30:24 01-Dec-2020 08:25:24 17700 NaN NaN
⋮
intervalsDE = accessIntervals(accessDE); intervalsDE = sortrows(intervalsDE, "StartTime", "ascend")
intervalsDE=51×8 table
Source Target IntervalNumber StartTime EndTime Duration StartOrbit EndOrbit
____________ ____________ ______________ ____________________ ____________________ ________ __________ ________
"MW Munchen" "GALILEO 2" 1 30-Nov-2020 22:23:24 01-Dec-2020 02:01:24 13080 NaN NaN
"MW Munchen" "GALILEO 3" 1 30-Nov-2020 22:23:24 01-Dec-2020 00:20:24 7020 NaN NaN
"MW Munchen" "GALILEO 4" 1 30-Nov-2020 22:23:24 30-Nov-2020 22:40:24 1020 NaN NaN
"MW Munchen" "GALILEO 19" 1 30-Nov-2020 22:23:24 01-Dec-2020 00:42:24 8340 NaN NaN
"MW Munchen" "GALILEO 20" 1 30-Nov-2020 22:23:24 30-Nov-2020 23:35:24 4320 NaN NaN
"MW Munchen" "GALILEO 10" 1 30-Nov-2020 22:38:24 30-Nov-2020 22:49:24 660 NaN NaN
"MW Munchen" "GALILEO 18" 1 30-Nov-2020 22:42:24 01-Dec-2020 01:37:24 10500 NaN NaN
"MW Munchen" "GALILEO 9" 1 30-Nov-2020 22:44:24 01-Dec-2020 00:43:24 7140 NaN NaN
"MW Munchen" "GALILEO 1" 1 30-Nov-2020 22:55:24 01-Dec-2020 03:37:24 16920 NaN NaN
"MW Munchen" "GALILEO 16" 1 30-Nov-2020 23:24:24 01-Dec-2020 02:15:24 10260 NaN NaN
"MW Munchen" "GALILEO 17" 1 01-Dec-2020 00:13:24 01-Dec-2020 02:19:24 7560 NaN NaN
"MW Munchen" "GALILEO 15" 1 01-Dec-2020 00:18:24 01-Dec-2020 03:45:24 12420 NaN NaN
"MW Munchen" "GALILEO 8" 1 01-Dec-2020 00:36:24 01-Dec-2020 05:07:24 16260 NaN NaN
"MW Munchen" "GALILEO 14" 1 01-Dec-2020 01:24:24 01-Dec-2020 05:17:24 13980 NaN NaN
"MW Munchen" "GALILEO 24" 1 01-Dec-2020 01:57:24 01-Dec-2020 02:36:24 2340 NaN NaN
"MW Munchen" "GALILEO 7" 1 01-Dec-2020 02:16:24 01-Dec-2020 06:29:24 15180 NaN NaN
⋮
intervalsIN = accessIntervals(accessIN); intervalsIN = sortrows(intervalsIN, "StartTime", "ascend")
intervalsIN=45×8 table
Source Target IntervalNumber StartTime EndTime Duration StartOrbit EndOrbit
______________ ____________ ______________ ____________________ ____________________ ________ __________ ________
"MW Bangalore" "GALILEO 3" 1 30-Nov-2020 22:23:24 01-Dec-2020 02:38:24 15300 NaN NaN
"MW Bangalore" "GALILEO 4" 1 30-Nov-2020 22:23:24 01-Dec-2020 01:05:24 9720 NaN NaN
"MW Bangalore" "GALILEO 5" 1 30-Nov-2020 22:23:24 30-Nov-2020 23:17:24 3240 NaN NaN
"MW Bangalore" "GALILEO 9" 1 30-Nov-2020 22:23:24 01-Dec-2020 00:55:24 9120 NaN NaN
"MW Bangalore" "GALILEO 10" 1 30-Nov-2020 22:23:24 30-Nov-2020 23:19:24 3360 NaN NaN
"MW Bangalore" "GALILEO 16" 1 30-Nov-2020 22:23:24 01-Dec-2020 02:42:24 15540 NaN NaN
"MW Bangalore" "GALILEO 21" 1 30-Nov-2020 22:23:24 30-Nov-2020 23:02:24 2340 NaN NaN
"MW Bangalore" "GALILEO 22" 1 30-Nov-2020 22:23:24 30-Nov-2020 22:29:24 360 NaN NaN
"MW Bangalore" "GALILEO 15" 1 30-Nov-2020 23:36:24 01-Dec-2020 05:12:24 20160 NaN NaN
"MW Bangalore" "GALILEO 2" 1 30-Nov-2020 23:44:24 01-Dec-2020 04:05:24 15660 NaN NaN
"MW Bangalore" "GALILEO 21" 2 01-Dec-2020 01:15:24 01-Dec-2020 04:12:24 10620 NaN NaN
"MW Bangalore" "GALILEO 1" 1 01-Dec-2020 01:20:24 01-Dec-2020 05:26:24 14760 NaN NaN
"MW Bangalore" "GALILEO 14" 1 01-Dec-2020 01:29:24 01-Dec-2020 07:33:24 21840 NaN NaN
"MW Bangalore" "GALILEO 20" 1 01-Dec-2020 02:04:24 01-Dec-2020 05:56:24 13920 NaN NaN
"MW Bangalore" "GALILEO 8" 1 01-Dec-2020 02:56:24 01-Dec-2020 06:37:24 13260 NaN NaN
"MW Bangalore" "GALILEO 19" 1 01-Dec-2020 03:20:24 01-Dec-2020 07:32:24 15120 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.
viewer3D = satelliteScenarioViewer(scenario);
Compare Access Between Ground Stations
Calculate access status between each satellite and ground station using the accessStatus method. Plot cumulative access for each ground station over the one day analysis window.
% Initialize array with size equal to number of timesteps in scenario timeSteps = mission.StartDate:seconds(60):mission.StartDate+days(1); statusUS = zeros(1, numel(timeSteps)); statusDE = statusUS; statusIN = statusUS; % Sum cumulative access at each timestep for idx = 1:24 statusUS = statusUS + accessStatus(accessUS(idx)); statusDE = statusDE + accessStatus(accessDE(idx)); statusIN = statusIN + accessStatus(accessIN(idx)); end clear idx; subplot(3,1,1); plot(timeSteps, statusUS); title("MW Natick to GALILEO") ylabel("# of satellites") subplot(3,1,2); plot(timeSteps, statusDE); title("MW München to GALILEO") ylabel("# of satellites") subplot(3,1,3); plot(timeSteps, statusIN); title("MW 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 = ["MW Natick", "MW München", "MW 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
MW Natick MW München MW Bangalore
_________ __________ ____________
Total # of intervals 52 51 45
Total interval time (hrs) 153.92 153.55 162.73
Mean interval length (min) 177.6 180.65 216.98
Mean # of satellites in view 6.4448 6.4289 6.8071
Min # of satellites in view 4 4 4
Max # of satellites in view 8 9 9
Walker-Delta constellations like Galileo 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. 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 40 minutes). All locations always have a minimum of 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] The European Space Agency: Galileo Facts and Figures. https://www.esa.int/Applications/Navigation/Galileo/Facts_and_figures
See Also
Objects
access|conicalSensor|groundStation|receiver|satellite|satelliteScenario|satelliteScenarioViewer|transmitter
Functions
Related Examples
- Multi-Hop Satellite Communications Link Between Two Ground Stations
- Satellite Constellation Access to a Ground Station
- Comparison of Orbit Propagators
- Estimate GNSS Receiver Position with Simulated Satellite Constellations
- Model, Visualize, and Analyze Satellite Scenario