High Precision Orbit Propagator

Precise modeling of satellite's perturbed motion (special perturbations approach)


Updated Sat, 17 Dec 2022 12:01:10 +0000

View License

The motion of a near-Earth satellite is affected by various forces. One of these forces is the Earth's central gravitation and the others are known as perturbations. These perturbations are classified into gravitational and non-gravitational forces. In this case, the equation of motion can be written as:
r ̈=-(GM/r^3)*r+γ_p
γ_p is the vector of additional accelerations induced by the disturbing forces.
γ_p=r ̈_E+r ̈_S+r ̈_M+r ̈_P+r ̈_e+r ̈_o+r ̈_D+r ̈_SP+r ̈_A+r ̈_emp
r ̈_E = Accelerations due to the non-spherically and inhomogeneous mass distribution within Earth (central body)
r ̈_S, r ̈_M, r ̈_P = Accelerations due to other celestial bodies (Sun, Moon, and planets)
r ̈_e, r ̈_o = Accelerations due to Earth and oceanic tides
r ̈_D = Accelerations due to atmospheric drag
r ̈_SP, r ̈_A = Accelerations due to direct and Earth-reflected solar radiation pressure
r ̈_emp = Accelerations due to unmodeled forces
Here I have used the following integrator and force model to simulate the satellite's perturbed motion:
Integrator: Variable-order Radau IIA integrator with step-size control
Force Model:
- gravity field of the Earth (GGM03C model)
- gravity of the solar system planets (positions of the planets are computed by JPLDE440)
- drag effect using NRLMSISE-00, Jacchia-Bowman 2006 & 2008, MSIS-86, Jacchia 70, or modified Harris-Priester atmospheric density model (in Accel.m you can uncomment your favorite model)
- solar radiation pressure using conical, geometrical, or cylindrical shadow model
- solid Earth tides (IERS Conventions 2010)
- ocean tides
- general relativity
- ECEF2ECI and ECI2ECEF transformations using IAU 2006 Resolutions
The Simulation starts by running the test_HPOP.m. In the InitialState.txt set initial values for your favorite satellite; lines 2-7 are the state vector of satellite/spacecraft in the International Terrestrial Reference Frame (ITRF). Lines 8 to 12 are the satellite parameters related to the forces from atmospheric drag and solar radiation pressure. Lines 8-10 are in units m^2 and kg. Line 11: Cr is the radiation pressure coefficient (Cr = 1 + reflectivity of satellite). Line 12: Cd is the atmospheric drag coefficient of the satellite. In the test_HPOP.m you can consider different perturbations by setting them 1 as follows:
AuxParam.n = 70; % maximum degree of central body's gravitation field
AuxParam.m = 70; % minimum order of central body's gravitation field
AuxParam.sun = 1; % perturbation of the Sun
AuxParam.moon = 1; % perturbation of Moon
AuxParam.planets = 1; % perturbations of planets
AuxParam.sRad = 1; % solar radiation pressure
AuxParam.drag = 1; % atmospheric drag
AuxParam.SolidEarthTides = 1; % solid Earth tides
AuxParam.OceanTides = 1; % ocean tides
AuxParam.Relativity = 1; % general relativity
Montenbruck O., Gill E.; Satellite Orbits: Models, Methods and Applications; Springer Verlag, Heidelberg; Corrected 3rd Printing (2005).
Montenbruck O., Pfleger T.; Astronomy on the Personal Computer; Springer Verlag, Heidelberg; 4th edition (2000).
Seeber G.; Satellite Geodesy; Walter de Gruyter, Berlin, New York; 2nd completely revised and extended edition (2003).
Vallado D. A; Fundamentals of Astrodynamics and Applications; McGraw-Hill, New York; 4th edition (2013).
NIMA. 2000. Department of Defense World Geodetic System 1984. NIMA-TR 8350.2, 3rd ed, Amendment 1. Washington, DC: Headquarters, National Imagery, and Mapping Agency.

Cite As

Meysam Mahooti (2023). High Precision Orbit Propagator (https://www.mathworks.com/matlabcentral/fileexchange/55167-high-precision-orbit-propagator), MATLAB Central File Exchange. Retrieved .

MATLAB Release Compatibility
Created with R2021b
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
Version Published Release Notes

test_HPOP, Density_Jacchia70.m, msis86.m, and nrlmsise00.m were modified.


IAU SOFA 2021-05-12 release was used.


Accel.m, AccelHarmonic_AnelasticEarth.m, AccelHarmonic_ElasticEarth.m, ECEF2ECI.m, ECI2ECEF.m, and nrlmsise00.m were modified to decrease CPU time.


Conical shadow model was added.


It was revised on 2022-10-12.


JB2006.m is added.


The AccelHarmonic_AnelasticEarth.m and the AccelHarmonic_ElasticEarth.m are modified.

AccelHarmonic_ElasticEarth.m and Accel.m are modified.

The AccelHarmonic_AnelasticEarth.m, the AccelHarmonic_ElasticEarth.m, and the Accel.m are modified.

The AccelHarmonic_AnelasticEarth.m and the AccelHarmonic_ElasticEarth.m are modified.

The AccelHarmonic_AnelasticEarth.m, the AccelHarmonic_ElasticEarth.m, and the Accel.m are modified.

Accel.m and JB2008.m are modified.

Mjday_TDB, test_HPOP, and nrlmsise00 are modified.

The DE436 full matrix is added.

JB2008 is modified.

AccelHarmonic_ElasticEarth.m, AccelHarmonic_AnelasticEarth.m, and HPOP.m are modified.

AccelHarmonic_ElasticEarth.m and AccelHarmonic_AnelasticEarth.m are modified.

HPOP.m, JPL_Eph_DE436.m and SAT_Const.m are modified.

Ephemeris.m is modified to decrease CPU time.
DE430 (planetary and lunar ephemerides) is replaced by DE436.
Satellite's ground track is plotted.
Effect of Solid Earth Tides is computed for elastic and anelastic Earth based on IERS Conventions 2010 (AccelHarmonic_ElasticEarth.m, AccelHarmonic_AnelasticEarth.m).
Revised on 2017-12-22.
IERS.m, JB2008.m, nrlmsise00.m, msis86.m, Density_Jacchia70.m, and HPOP.m are modified to decrease CPU time.
SAT_Const.m and constants.m are called once to decrease CPU time.

Revised on 2016-11-30.
Description is updated.
Integrator, shadow model and atmospheric density model are changed.
The image is added.