propagateOrbit
Syntax
Description
[
calculates the positions and velocities in the International Celestial Reference Frame
(ICRF) corresponding to the time specified by positions,velocities] = propagateOrbit(time,gpStruct)time using the
two-line-element (TLE) or orbit mean-elements message (OMM) data.
[
calculates the positions and velocities using specified orbital elements to define epoch
states. The function assumes that the epoch is the first element in the
positions,velocities] = propagateOrbit(time,a,ecc,incl,RAAN,argp,nu)time argument.
[
calculates the positions and velocities at the epoch defined by positions,velocities] = propagateOrbit(time,rEpoch,vEpoch)rEpoch
and vEpoch. The function assumes that the epoch is the first element in
the time argument.
[
calculates the positions and velocities using additional parameters specified by one or more
name-value arguments.positions,velocities] = propagateOrbit(___,Name=Value)
Examples
Read the data from the TLE (Two-Line Element) file named 'leoSatelliteConstellation.tle'. Calculate the position and velocity based on the TLE structure, 'tleStruct'. This file is located on the MATLAB® path and is provided with the Aerospace Toolbox.
tleStruct = tleread('leoSatelliteConstellation.tle')tleStruct=40×1 struct array with fields:
Name
SatelliteCatalogNumber
Epoch
BStar
RightAscensionOfAscendingNode
Eccentricity
Inclination
ArgumentOfPeriapsis
MeanAnomaly
MeanMotion
Calculate the position and velocity corresponding to the input time using the TLE data defined in tleStruct.
[r,v] = propagateOrbit(datetime(2022, 1, 3, 12, 0, 0),tleStruct);
This examples shows how to simulate a lunar free-return trajectory using numerical propagator. Define trajectory start and end times.
startTime = datetime(2022,11,17,14,40,0); endTime = datetime(2022,11,24,12,45,0);
Construct datetime vector representing the trajectory time stamps spaced at 60 second intervals.
sampleTime = 60; % s
time = startTime:seconds(sampleTime):endTime;Define the initial position and velocity of the spacecraft in ICRF. These represent the spacecraft states immediately after the translunar injection burn.
initialPosition = [5927630.386747557; ... 3087663.891097251; ... 1174446.969646237]; initialVelocity = [-5190.330300215130; ... 8212.486957313564; ... 4605.538019512981];
Define the options for numerical orbit propagation. Configure the ordinary differential equation (ODE) solver options, gravitational potential model and the third body gravity. Include third body gravity from all supported solar system bodies. Note that the numerical propagator also supports inclusion of perturbations due to drag from Earth atmosphere and solar radiation pressure. However, these effects are ignored in this example.
numericalPropOpts = Aero.spacecraft.NumericalPropagatorOptions( ... ODESet=odeset(RelTol=1e-8,AbsTol=1e-8,MaxStep=300), ... IncludeThirdBodyGravity=true, ... ThirdBodyGravitySource=[ ... "Sun" ... "Mercury" ... "Venus" ... "Moon" ... "Mars" ... "Jupiter" ... "Saturn" ... "Uranus" ... "Neptune" ... "Pluto"], ... GravitationalPotentialModel="point-mass")
numericalPropOpts =
NumericalPropagatorOptions with properties:
ODESolver: "ode45"
ODESet: [1×1 struct]
CentralBodyOptions: [1×1 Aero.spacecraft.CentralBodyOptions]
GravitationalPotentialModel: "point-mass"
IncludeAtmosDrag: 0
IncludeThirdBodyGravity: 1
ThirdBodyGravitySource: ["Sun" "Mercury" "Venus" "Moon" "Mars" "Jupiter" "Saturn" "Uranus" "Neptune" "Pluto"]
IncludeSRP: 0
PlanetaryEphemerisModel: "de405"
Define the physical properties (mass, reflectivity coefficient and solar radiation pressure area) of the spacecraft.
physicalProps = Aero.spacecraft.PhysicalProperties( ... Mass=10000, ... ReflectivityCoefficient=0.3, ... SRPArea=2)
physicalProps =
PhysicalProperties with properties:
Mass: 10000
DragCoefficient: 2.1790
DragArea: 1
ReflectivityCoefficient: 0.3000
SRPArea: 2
Propagate the spacecraft trajectory. While PropModle name-value argument defaults to "numerical" if you explicitly specify NumericalPropagatorOptions or PhysicalProperties name-value argument, PropModel has been explicitly specified for illustrative purposes.
position = propagateOrbit( ... time, ... initialPosition, ... initialVelocity, ... PropModel="numerical", ... NumericalPropagatorOptions=numericalPropOpts, ... PhysicalProperties=physicalProps);
Convert the position history to timetable.
positionTT = timetable(time',position');
Visualize the trajectory on a satellite scenario viewer. To do this, create a satellite scenario object.
sc = satelliteScenario(startTime,endTime,sampleTime);
Add the spacecraft to the scenario using the satellite function and the position timetable.
spacecraft = satellite(sc,positionTT,Name="Spacecraft");The satellite scenario viewer used for visualizing the scenario already renders visualization for the moon. However, you can gain better situational awareness if the lunar orbit is plotted as well. To do this, add a satellite to the scenario using the satellite function that is on the same trajectory as that of the Moon. To calculate the trajectory of the Moon, start by computing the barycentric dynamical times (TDB) for the simulation time samples as Julian dates.
leapSeconds = 37; ttMinusTAI = 32.184; terrestrialTime = time + seconds(leapSeconds + ttMinusTAI); tdbJD = tdbjuliandate([ ... terrestrialTime.Year' ... terrestrialTime.Month' ... terrestrialTime.Day' ... terrestrialTime.Hour' ... terrestrialTime.Minute' ... terrestrialTime.Second']);
Calculate the positions of the Moon for each scenario time sample using planetEphemeris and convert the positions into a timetable.
pMoonkm = planetEphemeris(tdbJD,"earth","moon"); % km pMoon = convlength(pMoonkm,'km','m'); % m pMoonTT = timetable(time',pMoon);
Add a satellite representing the Moon to the scenario using the satellite function. Set the orbit and marker color to red.
moon = satellite(sc,pMoonTT,Name="Moon"); moon.Orbit.LineColor="red"; moon.MarkerColor="red";
The scale of the distance involved in the scenario is large enough that the Earth may not be readily visible when viewing from the shadow side. To mitigate this, add a ground station at the North pole of the Earth and label it "Earth". Set the marker size of this ground station to 0.001. This way, the label "Earth" will always be visible near the position of the Earth.
earth = groundStation(sc, ... 90, ... % Latitude, deg 0, ... % Longitude, deg Name="Earth"); earth.MarkerSize = 0.001;
Run satelliteScenarioViewer to launch a satellite scenario viewer. Set the playback speed multiplier to 60,000. Set the camera reference frame to "inertial".
v = satelliteScenarioViewer(sc, ... CameraReferenceFrame="Inertial", ... PlaybackSpeedMultiplier=60000);
Set the camera position and orientation to visualize the free-return trajectory from a top-down view.
campos(v, ... 55.991361, ... 18.826434, ... 1163851259.541816); camheading(v, ... 359.7544952991605); campitch(v, ... -89.830968253450209); camroll(v, ... 0);
Play the scenario.
play(sc);
Input Arguments
Name-Value Arguments
Output Arguments
References
[1] Hoots, Felix R., and Ronald L. Roehrich. . Aerospace Defense Command, Peterson AFB CO Office of Astrodynamics, 1980.
[2] Vallado, David, et al. “Revisiting Spacetrack Report #3.” AIAA/AAS Astrodynamics Specialist Conference and Exhibit, American Institute of Aeronautics and Astronautics, 2006. , https://doi.org/10.2514/6.2006-6753.
