🛰 osculating2mean

🎯 Features

• Convert osculating orbital elements to/from mean orbital elements using spherical harmonics Earth gravity potential in MATLAB.

---

🚀 Index

• 💡 Description
• ✍🏼 Authors
• 📞 Contact
• 💿 Installation
• 📖 Documentation
• 🦆 Example
• ✨ Contributing to osculating2mean
• 📄 Lincense
• 💥 References

---

💡 Description

The osculating2mean toolbox is developed for MATLAB. The backbone computations of the spherical harmonics Earth gravity potential perturbations are written in FORTRAN and called from MATLAB.

The EGM96 NASA GSFC and NIMA Joint Geopotential Model is used.

This package was developed for non-singular orbital elements, for near-circular orbits, but it may be expanded to general orbits easily.

---

This repository is part of the work

L. Pedroso and P. Batista, Distributed decentralized receding horizon control for very large-scale networks with application to satellite mega-constellations, Control Engineering Practice, vol. 141, pp. 105728, 2023. doi: 10.1016/j.conengprac.2023.105728.

and based on

Hwang, C. and Hwang, L.S., 2002. Satellite orbit error due to geopotential model error using perturbation theory: applications to ROCSAT-2 and COSMIC missions. Computers & geosciences, 28(3), pp.357-367.

If you use this repository, please reference the publications above.

---

✍🏼 Authors

Leonardo Pedroso1
Pedro Batista1
Cheinway Hwang2
1Institute for Systems and Robotics, Instituto Superior Técnico, Universidade de Lisboa, Portugal
2Department of Civil Engineering, National Yang Ming Chiao Tung University, Taiwan

---

📞 Contact

osculating2mean toolbox is currently maintained by Leonardo Pedroso (leonardo.pedroso@tecnico.ulisboa.pt).

---

💿 Installation

To install osculating2mean:

• Clone or download the repository to the desired installation directory;
• Open the osculating2mean installation directory in MATLAB instance and run

  make_osculating2mean
  

• Add osculating2mean to the MATLAB path

  addpath('[installation directory]/osculating2mean');
  

Note

• The command make_osculating2mean must be run in the installation directory. If it is run (either when intalling osculating2mean or later) in other directory, the EGM96 model data will not be available.
• Ignore the massive dump of warnings when running make_osculating2mean :)

---

📖 Documentation

The documentation is divided into the following categories:

• Osculating elements from/to position-vector
• Eckstein-Ustinov J2 Perturbations
• Kaula Spherical Harmonics Geopotential Perturbations

Osculating elements from/to position-vector

To compute the osculating orbital elements from a position-velocity vector use

OE = rv2OEOsc(x)

and to compute the position-velocity vector from osculating orbital elements use

x = OEOsc2rv(OE)

where x is a $6\times 1$ position-velocity vector in SI units and OE is a $6\times 1$ vector of non-singular orbital elements, for near-circular orbits, i.e.,

• $a$: semi-major axis [m]
• $u$: $\omega + M$ ($\omega$: argument of perigee; $M$: mean anomaly) [rad]
• $e_x$: $e\cos(\omega)$ ($e$: excentricity)
• $e_y$: $e\sin(\omega)$
• $i$: inclination [rad]
• $\Omega$: longitude of ascending node [rad]

Implementation of these functions was adapted from (Vallado, 1997).

To adjust the maximum number of iterations and convergence criteria of the numerical solution of the Kepler equation in OEOsc2rv, it is also possible to set additional arguments

x = OEOsc2rv(OE, MaxIt, epsl)

For more details, see the thorough comments in the source code of this function.

Example: compute the osculating orbital elements from/to a position-velocity vector

>> x = 1.0e+06 * ...
  [6.4329;
  -1.5777;
  -2.0041;
   0.0028;
   0.0042;
   0.0056];
>> OE = rv2OEOsc(x);
>> OE(1)
ans =
  6.9195e+06
>> OE(2:6)
ans =
   5.9111
  -0.0003
  -0.0004
   0.9249
   6.2728
>> x = OEOsc2rv(OE)
x =
  1.0e+06 *

   6.4329
  -1.5777
  -2.0041
   0.0028
   0.0042
   0.0056

---

Eckstein-Ustinov J2 Perturbations

To compute the Eckstein-Ustinov first order perturbations due to $J_2$ use

EUPerturbations =  EcksteinUstinovPerturbations(OEMean)

where OEMean are the non-singular mean orbital elements for which the perturbation is computed and EUPerturbations are the perturbations to the mean orbital elements. This function is implemented according to (Eckstein and Hechler, 1970).

To convert from mean orbital elements to osculating orbital elements (or, equivalently, position-velocity) taking into account the Eckstein-Ustinov $J_2$ perturbations use

OEOsc = OEMeanEU2OEOsc(OEMean)

and to iteratively convert from osculating orbital elements (or, equivalently, position-velocity) to mean orbital elements, taking into account the Eckstein-Ustinov $J_2$ perturbations use

OEMean = OEOsc2OEMeanEU(OEOsc)

where OEOsc is a $6\times 1$ vector of non-singular osculating orbital elements for near-circular orbits and OEMean is a $6\times 1$ vector of non-singular mean orbital elements for near-circular orbits, both in SI units.

To adjust the maximum number of iterations and convergence criteria, it is also possible to set additional arguments

OEMean = OEOsc2OEMeanEU(OEOsc, MaxIt, epslPos, epslVel)

For more details, see the thorough comments in the source code of this function.

Example: compute the mean orbital elements taking into account the Eckstein-Ustinov J2 perturbations

>> x = 1.0e+06 * ...
  [6.4329;
  -1.5777;
  -2.0041;
   0.0028;
   0.0042;
   0.0056];
>> OEOsc = rv2OEOsc(x);
>> OEMean = OEOsc2OEMeanEU(OEOsc);
>> OEMean(1)
ans =
  6.9150e+06
>> OEMean(2:6)
ans =
   5.9114
  -0.0007
  -0.0000
   0.9247
   6.2730
>> OEOsc = OEMeanEU2OEOsc(OEMean);
>> OEOsc(1)
ans =
  6.9195e+06
>> OEMean(2:6)
ans =
   5.9111
  -0.0003
  -0.0004
   0.9249
   6.2728

---

Kaula Spherical Harmonics Geopotential Perturbations

To compute the spherical harmonics geopotential perturbations to the mean orbital elements according to (Kaula, 2013) use

dOE = KaulaGeopotentialPerturbations(t_tdb,OEmean,degree)

where t_tdb is the dynamic baricentric time since J2000 in seconds, OEMean are the non-singular mean orbital elements for which the perturbation is computed, degree is the maximum degree of the spherical harmonics geopotential model, and dOE are the perturbations to the following ordered set of orbital elements:

• $a$ (semi-major axis) [m]
• $e$ (excentricity)
• $i$ (inclination) [rad]
• $\Omega$ (longitude of ascending node) [rad]
• $\omega$ (argument of perigee) [rad]
• $M$ (mean anomaly) [rad]

The EGM96 NASA GSFC and NIMA Joint Geopotential Model is used.

The backbone of this function is implemented in FORTRAN, whose source code is based on the framework of (Hwang, 2001) and the programs published in (Hwang and Hwang, 2002). The FORTRAN code is called from MATLAB using a MEX function.

To convert from mean orbital elements to osculating orbital elements (or, equivalently, position-velocity ) taking into account the spherical harmonics geopotential perturbations use

OEOsc = OEMeanEUK2OEOsc(t_tdb,OEMean,degree) 

and to convert from position-velocity (or, equivalently, osculating orbital elements) to mean orbital elements, taking into account the spherical harmonics geopotential perturbations use

OEMean = OEOsc2OEMeanEUK(t_tdb,x,degree) 

where the arguments t_tdband degree are the same as those of function KaulaGeopotentialPerturbations, OEOsc is a $6\times 1$ vector of non-singular osculating orbital elements for near-circular orbits and OEMean is a $6\times 1$ vector of non-singular mean orbital elements for near-circular orbits, both in SI units.

The conversion from osculating to mean orbital elements is performed in 3 steps as proposed in (Spiridonova, Kirschner and Hugentobler, 2014):

• Iteratively compute the mean orbital elements taking into account the Eckstein-Ustinov J2 perturbations
• Compute the Kaula geopotential perturbations corresponding to the Eckstein-Ustinov mean orbital elements
• Compute the mean orbital elements by subtracting the perturbations to the osculating orbital elements

Example: compute the mean orbital elements taking into account the Kaula Spherical Harmonics Geopotential Perturbations

>> x = 1.0e+06 * ...
  [6.4329;
  -1.5777;
  -2.0041;
   0.0028;
   0.0042;
   0.0056];
>> OEOsc = rv2OEOsc(x);
>> OEMean = OEOsc2OEMeanEUK(11100,OEOsc,10);
>> OEMean(1)
ans =
  6.9150e+06
>> OEMean(2:6)
ans =
   5.9114
  -0.0007
  -0.0000
   0.9247
   6.2730
>> OEOsc = OEMeanEUK2OEOsc(11100,OEMean,10);
>> OEOsc(1)
ans =
  6.9195e+06
>> OEMean(2:6)
ans =
   5.9111
  -0.0003
  -0.0004
   0.9249
   6.2728

---

🦆 Example

In this example, which can be run with

example_osculating2mean

the osculating, mean Eckstein-Ustinov, and mean Eckstein-Ustinov-Kaula orbital elements are compared for a time-series of roughly 10 orbits of a satellite in LEO. The time-series were obtain with a simulation in TUDAT considering:

• atmospheric drag
• cannon ball solar radiation pressure
• third body perturbations from the Sun, Moon, Mars, Venus
• spherical harmonic gravity up to degree and order 12

The evolution of the non-singular orbital elements is depicted below

---

✨ Contributing to osculating2mean

The community is encouraged to contribute with

• Suggestions
• Addition of tools

To contribute to osculating2mean

• Open an issue (tutorial on how to create an issue)
• Make a pull request (tutorial on how to contribute to GitHub projects)
• Or, if you are not familiar with GitHub, contact the authors

---

📄 License

MIT License

---

💥 References

L. Pedroso and P. Batista, Distributed decentralized receding horizon control for very large-scale networks with application to satellite mega-constellations, Control Engineering Practice, vol. 141, pp. 105728, 2023. doi: 10.1016/j.conengprac.2023.105728.

@article{PedrosoBatista2023DistributedRHC,
	author = {Leonardo Pedroso and Pedro Batista},
	title = {Distributed decentralized receding horizon control for very large-scale networks with application to satellite mega-constellations},
	journal = {Control Engineering Practice},
	year = {2023},
	volume = {141},
	pages = {105728},
	doi = {10.1016/j.conengprac.2023.105728}
}

Eckstein, M.C., Hechler, H., 1970. A reliable derivation of the perturbations due to any zonal and tesseral harmonics of the geopotential for nearly-circular satellite orbits, ESOC, ESRO SR-13.

Hwang, C., 2001. Gravity recovery using COSMIC GPS data: application of orbital perturbation theory. Journal of Geodesy, 75(2), pp.117-136.

Hwang, C. and Hwang, L.S., 2002. Satellite orbit error due to geopotential model error using perturbation theory: applications to ROCSAT-2 and COSMIC missions. Computers & geosciences, 28(3), pp.357-367.

Kaula, W.M., 2013. Theory of satellite geodesy: applications of satellites to geodesy. Courier Corporation.

Spiridonova, S., Kirschner, M. and Hugentobler, U., 2014. Precise mean orbital elements determination for LEO monitoring and maintenance.

Vallado, D.A., 1997. Fundamentals of astrodynamics and applications. McGraw-Hill.