I'm currently considering a project to simulate the (Earth-seeing) coverage of 1 or more satellites over the course of a year. The point would be to compare how different orbits result in different amounts of coverage, both of the Earth's surface and objects placed there.

The code for this would be written using Golang, in an attempt to increase efficiency/speed of calculation over other languages like Python. As there aren't many orbital propagators for Golang, this would mean writing my own orbit propagation algorithm (both kinetic and kinematic), or using one of the few available ones, such as the SGP4-based go-satellite package. There would also be multi-processing involved.

As such, I'm trying to wrap my head around which orbit propagation methodology would be useful for this kind of simulation. Given that the simulation would be, well, simulating an entire year, errors would likely accumulate for any propagation method I choose. However, Fidelity should be high for short-term analysis, and reasonable for long-term. I'm also interested in what kind of CPU resource requirements are necessary for the different propagator models/methods/algorithms, too - I'd ideally not require a supercomputing cluster.

As such, my question is as the title stated it: how important is the choice I make in orbit propagation model for this kind of simulation what would really be the importance/magnitude of error, and which propagation method is best suited for this kind of thing?


2 Answers 2


I'm not an expert and this is not an expert answer, but these points may be helpful.

Don't even think about using SGP4, per my comments below this question and Wikipedia it's a circa 1980's clever approximation to get approximate state vectors within a few weeks of any given TLE's epoch. It's an approximator based on evolution of orbital elements, not a real propagator.

Orbital integration for a year is numerically pretty easy. It's not at all a stiff problem, most things are slowly varying unless you use a really high order gravity model. You could even implement a simple variable step size RK45, some higher order RK's and some canned numerical integrators available in Golang and compare them and I think you'll find that a few minutes on a laptop is all you need to run a year's worth of orbit propagation unless you have a complicated physical model for all of the small accelerations besides Earth's central field.

considering that laptops are gigaflops and Golang is C-like in speed; e.g. 10,000 flops per step with a 10 second step size, a year of propagation would only need 30 seconds at 1 gigaflop. However a gigaflop per laptop may be optimistic because it involves parallelized calculations (e.g. arrays) and to take advantage of that you might want to run several cases in parallel if you're really pinched for time.

But be sure to check your numerical integration technique against a symplectic integrator. For more on that see all the wonderful answers to What does “symplectic” mean in reference to numerical integrators, and does SciPy's odeint use them?

If you want to include general relativistic effects which are small but should be checked, see the answers to How to calculate the planets and moons beyond Newtons's gravitational force?

The real problem you have is in the model for all of those smaller accelerations.

For a longer discussion on these as well as a great discussion on the varieties of numerical integrators and how you might implement them I entheusiastically recommend anyone to look at Satellite Orbits; Models, Methods, Applications by Oliver Montenbruck and Eberhard Gill, Springer, 2000 which can also be found in google books and is discussed here and here. Also see answers to When calculating the future orbit of an artificial Earth satellite, is the Moon's gravity significant or insignificant?

  1. Gravity model; how lumpy of a gravity field are you going to use? Just J2, or several low order terms, or a larger set of spherical harmonics for the geopotential?
  2. Drag model; atmospheric drag is really hard to model accurately. You can get some simple drag coefficients from TLEs of similarly shaped objects, but remember that the Sun's activity heats the atmosphere, raising the density at LEO altitudes, and so drag is fundamentally unpredictable. No matter how good your model is, you will have to run several cases with different patterns of solar activity to see how much they change the orbit's altitude and therefore phase and therefore ground track.
  3. Sun and Moon: their gravitational effects will be small and you can probably add simple models for those by including their motion in your simulation.
  4. Sunlight's photon pressure: this will be small but not non-existent.

I think others may post answers or comments that recommend existing software you can run to just get a feel for the problem and also to check your calculations against.


Propagating for an entire year is extremely difficult to get right. There are ways to do it, but you should not attempt to rederive or reimplement them yourself. Before I describe them, however, I have a different suggestion: don't propagate at all.

If you are designing a constellation of satellites that have a mission to perform, then the operators of those satellites are going to make station-keeping maneuvers, spending fuel to maintain their assigned orbits. Trying to model those maneuvers in addition to the perturbing forces would be even harder, but you don't have to, because the whole point of making the maneuvers is to cancel out the effect of the forces.

For example, due to the lumpiness of the Earth, satellites in geostationary orbits gradually drift east or west in subpoint longitude, but nobody lets that add up for a whole year unless the satellite is dead. In order to stay in their assigned boxes, all active, so-called stationary satellites regularly (roughly weekly) execute short burns (just a few seconds at a time is plenty) in order to correct that drift and move from the "downhill" edge of their box back to the "uphill" edge, and then start sliding down again until it comes time for the next burn back up the potential slope.

If you're trying to calculate what this looks like for a year, then your most accurate answer will be obtained by figuring out how many days of drift you have between burns, and just replaying those few days of motion enough times to fill out a year. Similar tricks are useful for low-orbiting satellites: find how many days it takes to nearly repeat a ground track, and simulate only that many days. This may not be an integer, and is unlikely to evenly divide 365, so you will actually get a better answer for coverage if you simulate only one repeat period than exactly a year.

Doing that to high precision requires some propagation, but not for a year, which you should work hard to avoid. Under many circumstances, though, even that much isn't necessary. For basic constellation design, figuring out the inclinations and altitudes of the orbit planes and how many satellites to put in each is the main thing, and it can be done reasonably well just from nominal orbits, as in a Walker circular layout.

If, however, you are stuck simulating objects which aren't actively controlled (whether you are tracking debris, or have a customer who for some reason is willing to pay for east-west burns but not for north-south burns to control inclination), then you have to do what you actually asked about. In that case, the core idea is to avoid position and velocity. They change much too fast, so your integrator will need a very small step size. This will force it to take so many steps that it will both take far too long and also produce a very incorrect answer. You need to make your integrator's life easier by enabling it to take many fewer, and thus much larger, steps. To maintain accuracy, you need to switch to different coordinates, that change much more slowly.

This is pretty much what Kepler's orbital elements already are. The reason we always refer to Kepler is that even with every imaginable perturbation turned on, most orbits don't vary much from an ellipse, at least for a while. I apologize for the way this becomes a link-only answer, but I am not writing an entire textbook in orbital mechanics. Even if I were, I still wouldn't try to shoehorn it all into one post. :) At some point, you just have to hit the library and read the back issues. The way to proceed is to look at Lagrange's planetary equations showing how the inclination, eccentricity, and so on change over time, as described here and here and here.

Once you've gone that far, stop to wonder whether there might be some other set of coordinates, in which the variation from long-term periodic motion is even easier to analyze. The answer is yes, and there are actually quite a few different kinds, but they get progressively more difficult even to describe, much less understand, derive, or implement well. Now take a look at Delaunay variables, and Brouwer's and Kozai's mean elements, and try to follow the path that led from there to the SGP4 we all love to hate, through the Goddard Trajectory Determination System (here, here, and here) to the Draper Semi-Analytic Satellite Theory (here and here), including developments published just two weeks ago.

Finally, note that SST, which has been under development for decades and published in dozens of journal articles and doctoral theses, is now part of the open-source orbit toolset called OreKit. Download that, and have a go at wrapping a Java library to leverage their already tested and peer-reviewed code, rather than trying to type in literally hundreds of pages of equations and somehow finding enough time or assistants to test them at all adequately.

  • $\begingroup$ Nailed it - great answer $\endgroup$
    – Chris
    Commented Aug 28, 2021 at 0:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.