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

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.

• See this answer space.stackexchange.com/a/23409/12102 and also try searching this site for tags or terms like jpl-horizons, spice, pyephem and skyfield for starters.
– uhoh
Feb 2, 2018 at 0:20
• You're basically asking for a solution the n-body problem. It hasn't been solved. All you can do is simulate. Feb 2, 2018 at 16:24
• @LorenPechtel I'm going to integrate numerically. Feb 3, 2018 at 19:26

## 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:

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:

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).$$

• corrections, improvements, edits welcome!
– uhoh
Feb 4, 2018 at 12:28
• Thanks for the answer! What about the atmospheric drag, solar radiation, etc? Feb 4, 2018 at 15:39
• @TarlanMammadzada I left both of those out because they are highly dependent on the shape and construction of each individual spacecraft. You can get simplifications from Wikipedia articles on those particular subjects, but detailed calculations are quite involved. At some altitudes, drag is also is highly dependent on space weather, which leads to big changes in the profile of the upper atmosphere. Why don't you try simple models for those first, and once your calculation is working, compare to actual spacecraft data?
– uhoh
Feb 4, 2018 at 17:08
• @TarlanMammadzada Oh, there is a really helpful list of resources in the meta section of this site, take a look: space.meta.stackexchange.com/questions/249/…
– uhoh
Feb 4, 2018 at 17:08
• So, you mean, I have to simulate the orbits, and make corrections regarding to the acceleration? The question is, how to simulate the orbit of, let's say, Jupiter? The equation I described in the problem description is related to the Earth. Feb 7, 2018 at 8:56

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 += xDotDot * timeStep;

objvy += yDotDot * timeStep;

objvz += zDotDot * timeStep;

• Thanks, but İ'm trying to implement the code myself. If you developed the simulator, could you please explain how to deal with perturbations? Feb 4, 2018 at 13:56
• I used Runge-Kutta 4. The state of the solar system at a given date is obtained from JPL Horizons. Load it into arrays objx, objy, objz, objvx, objvy, objvz and set a time step. Every time you call RK4(), it propagates the system forwards or backwards in time depending on the time step. Feb 4, 2018 at 17:53
• Ok, I can propagate the orbits of the planets, moon and the satellite if I know the effects of perturbations. Which equations do you solve by RK4? Feb 4, 2018 at 20:52
• deltaV = acceleration * time, position = v * time. If you do this for every planet on every other planet, you will get the perturbations automatically. Even if you do not know the perturbations, this will reveal them. Feb 7, 2018 at 16:06

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

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.

• Yes, I'm going to integrate numerically. However, I'm trying to understand the principles, not to use already developed SW. Feb 3, 2018 at 19:25
• To numerically integrate: Feb 5, 2018 at 15:57
• To numerically integrate set up a double loop so that each object can pull on all the other objects. Inside the 2nd loop, compute acceleration due to gravity on inner object from outer object: a=GM/d^2. Decompose this acceleration into x,y,z components. Update inner object's velocity, component by component by deltaV = a * timestep. Then update inner object's position, component by component, and update its position: deltaR = v * time step. Feb 5, 2018 at 16:03
• This is called Euler-Cromer method and is good enough for many purposes if the time step is not too high. Higher order methods such as RK4 work basically the same way, but they use half steps and quarter steps to improve accuracy. Feb 5, 2018 at 16:04
• What you described looks simple, however what about the equation I mentioned in the problem description? Feb 7, 2018 at 8:49

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$.