Define Mission Parameters and Satellite 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(2019, 1, 4, 12, 0, 0); mission.Duration = hours(6);

Specify Keplerian orbital elements for the satellite(s) at themission.StartDate.

mission.Satellite.SemiMajorAxis = 6786233.13;% metersmission.Satellite.Eccentricity = 0.0010537; mission.Satellite.Inclination = 51.7519;% degmission.Satellite.RAAN = 95.2562;% degmission.Satellite.ArgOfPeriapsis = 93.4872;% degmission.Satellite.TrueAnomaly = 202.9234;% deg

Specify the latitude and longitude of a ground station to use in access analysis below.

mission.GroundStation.Latitude =42;% degmission.GroundStation.Longitude =-71;% deg

Open and Configure the Orbit Propagation Model

Open the included Simulink model. This model contains anOrbit Propagatorblock connected to output ports. TheOrbit Propagatorblock supports vectorization. This allows you to model multiple satellites in a single block by specifying arrays of initial conditions in theBlock Parameterswindow or usingset_param. The model also includes a "Mission Analysis and Visualization" section that contains a dashboard回调按钮. When clicked, this button runs the model, creates a new satelliteScenario object in the global base workspace containing the satellite or constellation defined in theOrbit Propagatorblock, 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 example.

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

Define the path to theOrbit Propagatorblock in the model.

mission.Satellite.blk = mission.mdl +"/Orbit Propagator";

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

set_param(mission.Satellite.blk,..."startDate", num2str(juliandate(mission.StartDate)),..."stateFormatNum","Orbital elements",..."orbitType","Keplerian",..."semiMajorAxis","mission.Satellite.SemiMajorAxis",..."eccentricity","mission.Satellite.Eccentricity",..."inclination","mission.Satellite.Inclination",..."raan","mission.Satellite.RAAN",..."argPeriapsis","mission.Satellite.ArgOfPeriapsis",..."trueAnomaly","mission.Satellite.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.Satellite.blk,..."centralBody","Earth",..."outportFrame","Fixed-frame");

Configure the propagator. This example uses a numerical propagator for higher position accuracy. Use numerical propagators to model Earth gravitational potential using the equation for universal gravitation ("Pt-mass"), a second order zonal harmonic model ("Oblate Ellipsoid (J2)"),或一个球面谐波模型("Spherical Harmonics"). Spherical harmonics are the most accurate, but trade accuracy for speed. For increased accuracy, you can also specify whether to use Earth orientation parameters (EOP's) in the internal transformations between inertial (ICRF) and fixed (ITRF) coordinate systems.

set_param(mission.Satellite.blk,..."propagator","Numerical (high precision)",..."gravityModel","Spherical Harmonics",..."earthSH","EGM2008",...% Earth spherical harmonic potential model"shDegree","120",...% Spherical harmonic model degree and order"useEOPs","on",...% Use EOP's in ECI to ECEF transformations"eopFile","aeroiersdata.mat");% EOP data file

Apply model-level solver setting usingset_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. In this example, theOrbit Propagatorblock is set to output position and velocity states in the ECEF (ITRF) coordinate frame.

mission.SimOutput = sim(mission.mdl);

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

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

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

mission.Satellite.TimeseriesPosECEF.TimeInfo.StartDate = mission.StartDate; mission.Satellite.TimeseriesVelECEF.TimeInfo.StartDate = mission.StartDate;

Load the Satellite Ephemerides into asatelliteScenarioObject

Create a satellite scenario object to use during the analysis portion of this example.

scenario = satelliteScenario;

Add the satellites to the satellite scenario as ECEF position and velocity timeseries using thesatellitemethod.

sat = satellite(scenario, mission.Satellite.TimeseriesPosECEF, mission.Satellite.TimeseriesVelECEF,..."CoordinateFrame","ecef")

sat = Satellite with properties: Name: "Satellite" ID: 1 ConicalSensors: [] Gimbals: [] Transmitters: [] Receivers: [] Accesses: [] GroundTrack: [1×1 matlabshared.satellitescenario.GroundTrack] Orbit: [1×1 matlabshared.satellitescenario.Orbit] OrbitPropagator: "ephemeris" MarkerColor: [1 0 0] MarkerSize: 10 ShowLabel: 1 LabelFontSize: 15 LabelFontColor: [1 0 0]

disp(scenario)

satelliteScenario with properties: StartTime: 04-Jan-2019 12:00:00 StopTime: 04-Jan-2019 18:00:00 SampleTime: 60 Viewers: [0×0 matlabshared.satellitescenario.Viewer] Satellites: [1×1 matlabshared.satellitescenario.Satellite] GroundStations: [] AutoShow: 1

Preview latitude (deg), longitude (deg), and altitude (m) for each satellite. Use thestatesmethod to query satellite states at each scenario time step.

foridx = numel(sat):-1:1% Retrieve states in geographic coordinates[llaData, ~, llaTimeStamps] = states(sat(idx),"CoordinateFrame","geographic");% Organize state data for each satellite in a seperate timetablemission.Satellite.LLATable{idx} = timetable(llaTimeStamps', llaData(1,:)', llaData(2,:)', llaData(3,:)',...'VariableNames', {'Lat_deg','Lon_deg',“Alt_m”}); mission.Satellite.LLATable{idx}end

ans=361×3 timetableTime Lat_deg Lon_deg Alt_m ____________________ _______ _______ __________ 04-Jan-2019 12:00:00 -44.804 120.35 4.2526e+05 04-Jan-2019 12:01:00 -42.797 124.73 4.2229e+05 04-Jan-2019 12:02:00 -40.626 128.77 4.2393e+05 04-Jan-2019 12:03:00 -38.322 132.53 4.2005e+05 04-Jan-2019 12:04:00 -35.848 136.07 4.2004e+05 04-Jan-2019 12:05:00 -33.289 139.35 4.203e+05 04-Jan-2019 12:06:00 -30.655 142.41 4.187e+05 04-Jan-2019 12:07:00 -27.884 145.34 4.1982e+05 04-Jan-2019 12:08:00 -25.069 148.09 4.1831e+05 04-Jan-2019 12:09:00 -22.234 150.68 4.1404e+05 04-Jan-2019 12:10:00 -19.297 153.19 4.1829e+05 04-Jan-2019 12:11:00 -16.343 155.58 4.1713e+05 04-Jan-2019 12:12:00 -13.388 157.89 4.07e+05 04-Jan-2019 12:13:00 -10.354 160.15 4.104e+05 04-Jan-2019 12:14:00 -7.3077 162.37 4.1291e+05 04-Jan-2019 12:15:00 -4.2622 164.55 4.0487e+05 ⋮

clearllaDatallaTimeStamps;

Display Satellite Trajectories Over the 3D Globe

显示卫星在地球轨道(WGS84 ellipsoid), use helper functionplot3DTrajectory.

mission.ColorMap = lines(256);% Define colormap for satellite trajectoriesmission.ColorMap(1:3,:) = []; plot3DTrajectories(mission.Satellite, mission.ColorMap);

Display Global and Regional 2D Ground Traces

View the global ground trace as a 2D projection using helper functionplot2DTrajectories:

plot2DTrajectories(mission.Satellite, mission.GroundStation, mission.ColorMap);

View regional ground trace. Select the region of interest from the dropdown menu:

plot2DTrajectories(mission.Satellite, mission.GroundStation, mission.ColorMap,"usa");

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

Add the ground station to the satelliteScenario object using thegroundStationmethod.

gs = groundStation(scenario, mission.GroundStation.Latitude, mission.GroundStation.Longitude,..."MinElevationAngle", 10,"Name","Ground Station")

gs = GroundStation with properties: Name: "Ground Station" ID: 2 Latitude: 42 Longitude: -71 Altitude: 0 MinElevationAngle: 10 ConicalSensors: [] Gimbals: [] Transmitters: [] Receivers: [] Accesses: [] MarkerColor: [0 1 1] MarkerSize: 10 ShowLabel: 1 LabelFontSize: 15 LabelFontColor: [0 1 1]

Attach line-of-sight access analyses between all individual satellites and the ground station using theaccessmethod.

foridx = 1:numel(sat) access(sat(idx), gs);endac = [sat(:).Accesses]; [ac(:).LineColor] = deal("green");

Display Access Intervals

Display access intervals for each satellite as atimetable. UseaccessStatusandaccessIntervalssatellite methods to interact with access analysis results.

foridx = numel(ac):-1:1 mission.Satellite.AccessStatus{idx} = accessStatus(ac(idx)); mission.Satellite.AccessTable{idx} = accessIntervals(ac(idx));% Use local function addLLAToTimetable to add geographic positions and% closest approach range to the Access Intervals timetablemission.Satellite.AccessTable{idx} = addLLAToTimetable(...mission.Satellite.AccessTable{idx}, mission.Satellite.LLATable{idx}, mission.GroundStation);endclearidx;

Display access intervals overlaying 2D ground traces of the satellite trajectories using helper functionplotAccessIntervals.

plotAccessIntervals(mission.Satellite, mission.GroundStation, mission.ColorMap);

mission.Satellite.AccessTable{:}

ans=2×8 tableSource Target IntervalNumber StartTime EndTime Duration LLA (deg, deg, m) ClosestApproach (m) ___________ ________________ ______________ ____________________ ____________________ ________ _________________ ___________________ "Satellite" "Ground Station" 1 04-Jan-2019 12:44:00 04-Jan-2019 12:50:00 360 {6×3 double} 5.0087e+05 "Satellite" "Ground Station" 2 04-Jan-2019 14:21:00 04-Jan-2019 14:25:00 240 {4×3 double} 1.102e+06

Further Analysis

Play thesatelliteScenarioobject to open and animate the scenario in asatelliteScenarioViewerwindow.

play(scenario) disp(scenario.Viewers(1))

Viewer with properties: Name: 'Satellite Scenario Viewer' Position: [560 240 800 600] Basemap: 'satellite' PlaybackSpeedMultiplier: 50 CameraReferenceFrame: 'ECEF' CurrentTime: 04-Jan-2019 12:00:25 Dimension: '3D'

References

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