propagateOrbit

Calculate position and velocity of spacecraft

Since R2024b

Syntax

[positions,velocities] = propagateOrbit(time,gpStruct)

[positions,velocities] = propagateOrbit(time,a,ecc,incl,RAAN,argp,nu)

[positions,velocities] = propagateOrbit(time,rEpoch,vEpoch)

[positions,velocities] = propagateOrbit(___,Name=Value)

Description

[positions,velocities] = propagateOrbit(time,gpStruct) calculates the positions and velocities in the International Celestial Reference Frame (ICRF) corresponding to the time specified by time using the two-line-element (TLE) or orbit mean-elements message (OMM) data.

[positions,velocities] = propagateOrbit(time,a,ecc,incl,RAAN,argp,nu) 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 time argument.

[positions,velocities] = propagateOrbit(time,rEpoch,vEpoch) calculates the positions and velocities at the epoch defined by rEpoch and vEpoch. The function assumes that the epoch is the first element in the time argument.

[positions,velocities] = propagateOrbit(___,Name=Value) calculates the positions and velocities using additional parameters specified by one or more name-value arguments.

Examples

Read Data from TLE File and Calculate Position and Velocity from Data

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);

Simulate a Lunar Free-Return Trajectory Using Numerical Propagator

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: [1x1 struct]
    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

time — Julian dates or datetime objects

1-by-m vector | m-by-1 vector

Julian dates or datetime objects, specified as a 1-by-m vector or m-by-1 vector of doubles or datetime objects. time must be strictly increasing.

When you specify a datetime object whose TimeZone property is empty, the propagateOrbit function assumes that the datetime object time zone is UTC.

Example: datetime(2023,6,16)

Example: juliandate(2023,6,16)

Data Types: double | datetime

a — Semimajor axes

numSc-element vector | scalar

Semimajor axes, specified as a numSc-element vector or a scalar, in meters. If all spacecraft have the same semimajor axis, you can specify a as a scalar. numSc is the number of spacecraft.

Example: 10000000

Example: [8000000 10000000]

Dependencies

If any element of a is negative, the corresponding element in ecc must be greater than or equal to 1. The only supported PropModel in this instance is "numerical".

Data Types: double

ecc — Orbit eccentricity

numSc-element vector | scalar

Orbit eccentricity, which is the deviation of the orbital curve from circular, specified as a numSc-element vector or scalar. If all spacecraft have the same eccentricity, you can specify ecc as a scalar. numSc is the number of spacecraft.

Example: 0.1

Example: [0.1 0.2]

Dependencies

If any element of ecc is greater than or equal to 1, the corresponding element in a must be negative. The only supported PropModel in this instance is numerical.

incl — Inclination

numSc-element vector | scalar

Inclination, which is the vertical tilt of the orbit with respect to the ICRF xy plane measured at the ascending node, specified as a numSc-element vector or a scalar in degrees. numSc is the number of spacecraft.

Example: 10

Example: [20 30]

RAAN — Right ascension of ascending node (RAAN)

numSc-element vector | scalar

Right ascension of ascending node (RAAN), specified as a numSc-element vector or a scalar, in degrees. If all spacecraft have the same right ascension of ascending node, you can specify RAAN as a scalar. numSc is the number of spacecraft.

Example: 10

Example: [20 30]

argp — Argument of periapsis

numSc-element vector | scalar

Argument of periapsis, which is the angle from the orbit ascending node to periapsis (closest point of orbit to the central body), specified as a numSc-element vector or a scalar in degrees. If all spacecraft have the same argument of periapsis, you can specify RAAN as a scalar. numSc is the number of spacecraft.

Example: 10

Example: [20 30], if there are two spacecraft

nu — True anomaly

numSc-element vector | scalar

True anomaly, which is the angle between periapsis and the initial position of the spacecraft, specified as a numSc-element vector or a scalar. If all spacecraft have the same true anomaly, you can specify nu as a scalar. numSc is the number of spacecraft.

Example: 10

Example: [20 30]

rEpoch — Position at epoch

3-by-1 vector | 3-by-numSc matrix

Position at epoch, specified as a 3-by-1 vector or 3-by-numSc matrix. If all spacecraft have the same initial position at epoch, you can specify rEpoch as a 3-by-1 vector. numSc is the number of spacecraft.

Example: [10000000;1000;2000]

Example: [10000000 12000000;1000 2000;2000 4000]

Dependencies

rEpoch and the corresponding vEpoch vectors cannot be parallel when PropModel is set to 'sgp4', 'sdp4', 'general-perturbations' or 'two-body-keplerian'.

vEpoch — Velocity at epoch

3-by-1 vector | 3-by-numSc matrix

Velocity at epoch, specified as a 3-by-1 vector or 3-by-numSc matrix. If all spacecraft have the same initial velocity at epoch, you can specify vEpoch as a 3-by-1 vector. numSc is the number of spacecraft.

Example: [0;8000;0]

Example: [0 100;8000 9000;200 300]

Dependencies

vEpoch and the corresponding rEpoch vectors cannot be parallel when PropModel is set to 'sgp4', 'sdp4', 'general-perturbations' or 'two-body-keplerian'.

gpStruct — General perturbations data structures

numSc-element vector

General perturbations structures representing TLE or OMM data, specified as a numSc-element vector. numSc is the number of spacecraft. To create the TLE or OMM data structure, consider using tleread or ommread respectively. For more information on TLE and OMM file structures, see Two Line Element (TLE) Files.

Example: tleread("leoSatelliteConstellation.tle")

Data Types: struct

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: PropModel="two-body-keplerian" specifies two-body-keplerian as the propagation method.

Epoch — Epoch

first value of time argument (default) | scalar

Epoch, specified as a scalar datetimeobject. The data type of 'Epoch' is double.

When you specify a datetime object whose TimeZone property is empty, the propagateOrbit function assumes that the datetime object time zone is UTC.

Example: datetime(2023,6,19)

Default value

• Epoch defined in gpStruct when using gpStruct input and PropModel is set to 'sgp4', 'sdp4' or 'general-perturbations'.

• First element of time input in all other cases.

Dependencies

When you specify the epoch states using gpStruct, do not specify this name-value argument.

Data Types: double

NumericalPropagatorOptions — Options used by the numerical orbit propagator

NumericalPropagatorOptions object

Options used by the numerical orbit propagator, specified as a scalar Aero.spacecraft.NumericalPropagatorOptions object.

Default Value

The default value is the object returned when calling the Aero.spacecraft.NumericalPropagatorOptions object without any inputs.

Dependencies

If you specify PropModel to 'sgp4', 'sdp4', 'general-perturbations' or 'two-body-keplerian', do not specify this name-value argument.

PhysicalProperties — Physical properties of spacecraft used by numerical orbit propagator

Aero.spacecraft.PhysicalProperties object | array of Aero.spacecraft.PhysicalProperties objects

Physical properties of spacecraft used by the numerical orbit propagator, specified as a scalar or array of Aero.spacecraft.PhysicalProperties objects. If all spacecraft use the same physical properties, you can specify 'PhysicalProperties' as a scalar.

Default Value

The default value is the object returned when calling the Aero.spacecraft.PhysicalProperties object without any inputs.

Dependencies

If you specify PropModel to 'sgp4', 'sdp4', 'general-perturbations' or 'two-body-keplerian', do not specify this name-value argument.

PropModel — Orbital state propagation method

"general-perturbations" | "two-body-keplerian" | "sgp4" | "sdp4" | "numerical"

Orbital state propagation method, specified as one of these values:

• "general-perturbations" — "sgp4" or "sdp4", depending on orbital period. When the orbital period corresponding to the epoch states is less than 225 minutes, the propagation model is "sgp4". Otherwise, it is "sdp4".

• "two-body-keplerian" — Two-body Keplerian orbit propagator.

• "sgp4" — Simplified General Perturbations-4 orbit propagator.

• "sdp4" — Simplified Deep-Space Perturbations-4 orbit propagator.

• "numerical" — Numerical propagation model.

Example: PropModel = "two-body-keplerian"

Default Value

When the epoch states are specified using one of the orbital elements such as a, ecc, incl, RAAN, argp and nu or position or velocity states such as rEpoch or vEpoch, the default value is "general-perturbations" if the orbital energy corresponding to every epoch state is negative. Otherwise, the default value is "numerical". To calculate the orbital energy, the epoch % states are converted to inertial position and velocity vectors r_eci and v_eci. Assuming μ is the standard gravitational parameter of Earth, the orbital energy is $$E=\frac{|v_{eci}{|}^2}{2}-\frac{\mu }{\left|r_{eci}\right|}$$

Data Types: char | string

BStar — Aerodynamic drag term

0 (default) | numSc-element vector | scalar

Aerodynamic drag term, specified as a numSc-element vector or scalar. If all spacecraft use the same Bstar value, you can specify this argument as a scalar. numSc is the number of spacecraft.

Example: BStar = 0

Example: BStar = [0.0001 0.0002], if there are two spacecraft

Dependencies

propagateOrbit uses this value when:

• PropModel is set to "general-perturbations", "sgp4", or "sdp4".

• Epoch states are not specified using gpStruct.

Data Types: double

InputCoordinateFrame — Input coordinate frame of rEpoch and vEpoch

"icrf" (default) | "fixed-frame" | "geographic"

Input coordinate frame of rEpoch and vEpoch, specified as one of these values:

• "icrf" — rEpoch and vEpoch are defined in the ICRF in m and m/s, respectively.

• "fixed-frame" — rEpoch and vEpoch are defined as the Earth-centered Earth-fixed (ECEF) frame in m and m/s, respectively.

• "geographic" — rEpoch is defined as latitude (deg), longitude (deg), and altitude (m). vEpoch is defined in the fixed-frame velocity in m/s represented in the north-east-down (NED) frame defined by rEpoch.

Example: InputCoordinateFrame="fixed-frame"

Data Types: char | string

OutputCoordinateFrame — Output coordinate frame of posititions and velocities

"icrf" (default) | "fixed-frame" | "geographic"

Output coordinate frame of output positions and velocities, specified as one of these values:

• "icrf" — positions and velocities are defined in the ICRF in m and m/s respectively.

• "fixed-frame" — positions and velocities are defined in the Earth-centered Earth-fixed (ECEF) frame in m and m/s respectively.

• "geographic" — positions is defined as latitude (deg), longitude (deg), and altitude (m).. velocities is defined as the fixed-frame velocity in m/s represented in the north-east-down (NED) frame defined by positions.

Example: OutputCoordinateFrame = "geographic"

Data Types: char | string

Output Arguments

positions — Positions

3-by-m-by-numSc array

Positions, returned as a 3-by-m-by-numSc array. m is the number of time samples. numSc is the number of spacecraft.

velocities — Velocities

3-by-m-by-numSc array

Velocities, returned as a 3-by-m-by-numSc array. m is the number of time samples. numSc is the number of spacecraft.

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.

Version History

Introduced in R2024b

R2024b: Version History

The Aerospace Toolbox now incorporates the latest and most accurate equations from https://celestrak.org/publications/AIAA/2006-6753/ in propagateOrbit. Both the SGP4 and SDP4 orbit propagators use these equations and return identical results. The algorithms are detailed in the work of Vallado, David, et al, “Revisiting Spacetrack Report #3.”, AIAA/AAS Astrodynamics Specialist Conference and Exhibit, American Institute of Aeronautics and Astronautics, 2006, https://celestrak.org/publications/AIAA/2006-6753/.