Satellite Orbit Modeling

Modeling of satellite's perturbed motion (special perturbations approach)

Orbit modeling is the process of creating mathematical models to simulate motion of a massive body as it moves in orbit around another massive body due to gravity. Other forces such as gravitational attraction from tertiary bodies, air resistance, solar pressure, or thrust from a propulsion system are typically modeled as secondary effects. Directly modeling an orbit can push the limits of machine precision due to the need to model small perturbations to very large orbits. Because of this, perturbation methods are often used to model the orbit in order to achieve better accuracy. Orbit models are typically propagated in time and space using special perturbation methods. This is performed by first modeling the orbit as a Keplerian orbit. Then perturbations are added to the model to account for the various perturbations that affect the orbit. Special perturbations can be applied to any problem in celestial mechanics, as it is not limited to cases where the perturbing forces are small. Special perturbation methods are the basis of the most accurate machine-generated planetary ephemerides such as Jet Propulsion Laboratory Development Ephemerides.

Here I have used the following integrator and force model to simulate satellite's perturbed motion:

Integrator: Runge-Kutta-Fehlberg 7(8).

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 NRLMSISE00 atmospheric density model, solar radiation pressure using cylindrical model, solid Earth tides (IERS Conventions 2010), ocean tides and general relativity.

The Simulation starts by running the test_Envisat.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_Envisat.m you can consider different perturbations by setting them 1 as follows:

AuxParam.n = 40; % maximum degree of central body's gravitation field

AuxParam.m = 40; % maximum 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

References:

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.

https://ssd.jpl.nasa.gov/planets/eph_export.html

http://celestrak.org/SpaceData/EOP-All.txt

http://celestrak.org/SpaceData/SW-All.txt