How to determine an algorithm for Low Earth Orbit propagation, considering perturbation from the moon, sun, etc?

Question

I'm trying to determine an algorithm for orbit propagation for a satellite on LEO, with high accuracy in a six month time interval. The propagation should consider the Sun and Moon, and any other effects worth accounting for.

I'm going to integrate the orbit numerically. Now I'm using the equation, that describes a Keplerian motion perturbed by main effects of an oblate Earth.

However, I have to consider also other effects to achieve realistic results.

Could anyone suggest literature, or provide some links & references? I'm trying to implement the methods and algorithms myself, please, don't suggest SW for this purpose.

Answer 1

Newtonian Gravity:

Your biggest effect after the monopole term and $J_2$ will be all of the other multipole terms representing the realistic gravitational field of the Earth expressed in spherical harmonics. $J_2$ is the largest of them by far (order 1E-03), but certainly not the only one. I believe the next three in line include $J_{22}$ which represents/includes the lowest azimuthal (longitudinal) term. See for example the Wikipedia subsection The deviations of Earth's gravitational field from that of a homogeneous sphere where you will see where $J_2$ fits in. You will find the coefficients of the multipole terms for the scalar potential field expressed in various places, out to different orders. You will then have to do the math to find the vector gradient of the field which will give you the acceleration terms as you've shown for $J_2$. Or, you might find tables with the expressions and coefficients already evaluated for you.

Some approximate values just to give you an idea (copied from 1, 2 and 3):

J2  =  1.1E-03
J3  = -2.5E-06
J4  = -1.6E-06
J22 = -1.8E-06

I've written more about using the Geopotential coefficients in the as-yet unanswered question For the mathematical relationship between J2 (km^5/s^2) and dimensionless J2 - which one is derived from the other?.

General Relativity:

I've extracted this bit about a commonly used expression for corrections due to General Relativity from this answer, which includes several sources/references where you can read more about this kind of calculation.

Remember that your coordinates for oblateness are for the Earth, and the ecliptic (where you find the Sun and Jupiter for example) is tilted by roughly 23 degrees, and the Moon's orbit has it's own plane and complex orbital motion.

You either need to treat at least the Sun (which also moves!), Earth and Moon (next in line would come Jupiter and Venus) as additional bodies whose motion needs to be integrated numerically (you can do that first, store it, then interpolate from it) or you can interpolate their positions from an existing ephemeris. You can search this site for references to spice or Skyfield for interpolators of the NASA JPL Development Ephemerides for example.

Although I'm not so familliar with GR, I'm going to recommend an equation that seems to work well and seems to be supported by several links. It is an approximate relativistic correction to Newtonian gravity that is used in orbital mechanics simulations. You will see various forms in the following links, most with additional terms not shown here:

• https://physics.stackexchange.com/q/313146/83380
• Eq. 1 in https://www.lpi.usra.edu/books/CometsII/7009.pdf
• Eq. 27 in https://ipnpr.jpl.nasa.gov/progress_report/42-196/196C.pdf
• Eq. 4-26 in https://descanso.jpl.nasa.gov/monograph/series2/Descanso2_all.pdf
• also see David Hammen's discussion/advice in this answer.

The following approximation should be added to your Newtonian gravity terms, where $\mathbf{r}$ and $\mathbf{v}$ represent the distance and velocity vectors between each object who's acceleration you are evaluating and each object for which the GR term you feel has a strong-enough gravitational field to consider. Or you can do it for all pairs if you don't want to bother with the housekeeping.

$$\mathbf{a_{GR}} = GM \frac{1}{c^2 |r|^3}\left(4 GM \frac{\mathbf{r}}{|r|} - (\mathbf{v} \cdot \mathbf{v}) \mathbf{r} + 4 (\mathbf{r} \cdot \mathbf{v}) \mathbf{v} \right).$$

Answer 2

You can use my simulator and add your J2 code. Start here: http://orbitsimulator.com/gravitySimulatorCloud/simulations/1517709441154_tarlan.html

This shows Earth with the ISS in orbit. As the ISS is propagated, perturbations from the Moon, Sun, and all the planets are taken into account. But Earth is treated as a sphere. As a result, low orbits are reasonably accurate for a day or two before Earth's bulge noticeably precesses the ISS's orbit. High orbits remain accurate for years to decades.

Your equations show the x,y,z accelerations due to the bulge. In the menu, open the Autopilot > Per Iteration. You can add code here that will be executed with every iteration.

As a pretty useless example, enter this code and press "Update". This causes the ISS to accelerate at 1 m/s/s along each component.

// replace the following 3 lines with your J2 code xDotDot = 1; // m/s^2

yDotDot = 1;

zDotDot = 1;

// The ISS is object 12. Update its velocity using the acceleration. (delta v = at)

objvx[12] += xDotDot * timeStep;

objvy[12] += yDotDot * timeStep;

objvz[12] += zDotDot * timeStep;

Answer 3

If you are propagating the orbit of an object in a low-Earth orbit and will numerically integrate it, then you should use a very precise gravity model. AFAIK, this will add more effects to the solution then the other planets gravity. My suggestion is that you use EGM2008 (or EGM96) to compute the Earth acceleration.

For the numerical propagation, basically you will have to add all the accelerations and use a very accurate differential equation solution. If you want to see an example of a high precision orbit propagator, you can look this solution in MATLAB:

https://www.mathworks.com/matlabcentral/fileexchange/55167-high-precision-orbit-propagator

Answer 4

As others have mentioned, you're talking about a n-body problem. As they have also mentioned, it's an unsolved problem. The best we've come up with is to integrate the system numerically. Tools exist for this purpose, but you're generally going to have to know a bit of programming to get them to work. For instance you can try Rebound. I've never worked with this library personally, but it seems like a very good bet. I would start with this example which models the solar system and add your hypothetical spacecraft to it.

Answer 5

If you substitute your EOMs for the first term here I think you'll get the N-body EOMs that you are after:

From Ocampo, Elements of a Software System for Spacecraft Trajectory Optimization in Conway, Spacecraft Trajectory Optimization page 85:

An example expression for $\vec g$ is the acceleration acting on a spacecraft due to a main gravitational body $c_b$ and possibly including the gravitational acceleration due to $n_b$ additional celestial bodies whose time dependent positions $\vec r_j(t)$ with respect to $c_b$ are obtained from a precomputed ephemeris

$$\vec g(\vec r,t) = - \frac{Gm_{c_b}}{r^3} \vec r- G\sum_{j=1}^{n_b} m_j \left(\frac{\vec r-\vec r_j(t)}{|\vec r-\vec r_j(t)|^3} + \frac{\vec r_j(t)}{r_j^3(t)}\right) \tag{4.13}$$

where $G$ is the universal constant of gravitation, $m_{c_b}$ is the mass of the main celestial body, and $m_j$ is the mass of celestial body $j$. Depending on the problem, additional terms accounting for other common accelerations such as atmospheric drag, solar radiation pressure, and non-spherical celestial bodies need to be added to the vector function $\vec g$.