Programming an Orbital Simulator

[Davidoux]

Hello,

I plan to write a (modest) orbital simulator in C++, I am an experienced developer with an average knowledge of classic physics;
I would like to have some recommendations (recommended book reading, eventual existing C++ libraries that can help) before starting knowing that I need near-real time performance for 3D display.
The scope of this simulator is to predict orbit of one or more orbiting spacecraft s
(two would be enough) around the earth, the spacecraft should be able to perform some manoeuvers (attitude correction, injection, orbit circularization, ...) and eventually rendezvousing.

Thanks again for your help, I am new to this forum and I hope I posted in the right topic

David

 

[Sane]

It might help to tell us what programming knowledge you already have.

If you plan on using C++, you will want to look at OpenGL. OpenGL can meet all of your three-dimensional needs. But if you've never programmed anything in 3D, doing an orbital simulator is going to be very difficult.

You'll want to look back over vectors in R3. But this will probably come back to you as you need it, since most of it's common sense. There are many resources for OpenGL and programming in three dimensions. I'm not sure what else to suggest.

I wouldn't worry about it right now, but you'll need to find 3D models for the planets, or else design them yourself. It would be straight-forward, depending on what level of detail you require, if you have the resources.

Depending on what you've done before, you might not want to start with a project this deep. What is your experience?

 

[chroot]

This simulator need not actually use any 3D graphics, though, admittedly, it'd be nice.

Simulating the physics is nothing more than dividing time up into small steps. At each step, you need the position and velocity vector for all of your objects. At each step, you calculate the forces on the objects due to gravity, sum them all up into a single force vector, and then solve where the object will be in the next time step, given its current position, velocity, and experienced force.

- Warren

 

[Sane]

Oh, that is true. I assumed that camera movements were required. But since this is never mentioned, I guess its only ideal.

Then I'd suggest you should look into the Simple Direct Media Layer for C++. It's reasonable to learn and use. More than OpenGL at least. But you still might want to do some easier graphical projects before tackling this.

 

[Davidoux]

Hi Sane, chroot,

actually I don' t want to focus on the 3D and graphical stuff, my concern is about developing first an orbital physics library, which would support one ore more spacecraft orbiting Earth and being able to exert thrust forces, so my concern is about choosing an integrator algorithm (RK, simplectic, euler, ..) simulate the spacecraft movement (rotations, translations,) .. and morevover in the future to support gravitational perturbation. So it' s all about choosing the correct math, equations and integration method.

David

 

[Sane]

For the method of integration, it would all depend on how you plan on instructing your planets. Since this seems to be more of a question directed towards Physics, than it is a question about Programming, you would get more help asking for the easiest equations to calculate all of this on one of the Physics forums. Making two posts, in different forums, of the same question, is against the rules, so you may have to clear this with a moderator first, or rephrase your question.

 

[D H]

In your first post you said you wanted a modest orbit simulator. Here are some suggestions to make it simple:

• Make your orbital maneuvers impulsive (i.e., apply the desired delta-v instantaneously).
• Forgo rotational dynamics. Simply apply your delta-v in the desired direction. Don't worry about how a real spacecraft would have to orient itself to achieve this delta-v.
• Model the Earth as a spherical body. You can add quite a bit more realism by modeling the Earth as an oblate spheroid (i.e., account for J2). Anything beyond J2 will require a considerably more complex geopotential model.
• Forgo atmospheric drag. This is a significant effect, but it is a bear to model.
• Forgo everything but the Earth. Pretend the Earth-centered J2000 frame is an inertial frame. There is, afterall, no reason to simulate third-body effects when you are treating the Earth as a spherical body.
• Use a simple integrator, but not too simple. Euler integration is unstable at all time steps in a central force field. I suggest using a fixed time step RK-4 integrator as a good start. A more advanced integrator, can make much larger time steps than RK-4, but those big time steps aren't of much use when you move on to modeled a real spacecraft (non-impulsive burns, rotational dynamics, jet errors, navigation errors, ...) or when you want to display the vehicles graphically.

Is this too simple?

Most NASA centers and most university aerospace engineering departments use much more elaborate simulators. For example, look into the Trick simulation environment at NASA's Johnson Space Center or SPICE at the Jet Propulsion Laboratory.

 

[Davidoux]

(to the moderators : forgive me if I posted in the wrong topic, but maybe there should be an 'numerical analysis' topic in the programming sub-forum, and dedicated to physics simulation programming ?)

Hi D H,

Thanks this is exactly the kind of suggestions I sought for :)
I agree with most of your suggestions as a good starting point and RK4 seems to be indeed a reasonable choice.

As most of the spacecraft maneuvers will be at close range from a station (final approach of docking), I plan not to simulate J2 as its influence is negligible during final approach; However drag is most important, but it can be 'loosely' simulated at first (if it is possible).

Finally, I know that there is an Euler variant called the 'Explicit' (the one you are talking is the 'Implicit' and is indeed unstable) but I don't really know the advantages of using such method against RK4

D

 

[D H]

Whoa!

You said modest in your first post, now you want to simulate rendezvous? Rendezvous is an active phase of flight. This makes for a very elaborate simulation.

A spacecraft inherently has imperfect sensors, imperfect effectors, and imperfect navigation. These errors are significant contributors to the success (or lack thereof) of a docking.

Proper attitude is obviously a critical factor in rendezvous. While orbital maneuvers can be simulated with a 3-DOF simulation, rendezvous cannot.

Drag is important during final approach. The station is an ungainly butterfly compared to the bullet-like spacecraft that dock with the station.

 

[Jheriko]

One thing to bear in mind when writing real-time simulations is that the frame rate is not constant, which means that the \Delta ts in some of the numerical integrators will need to be carefully dealt with to get the correct result. In addition to this, any changes you make to create framerate independence will change the order of accuracy of the method. If you aren't relying on more than one frame of previous data (Euler, RK4 etc...) then you are probably safe from this.

If you want speed though I would choose an algorithm like Verlet integration, which may add some complications, but gives similar convergence to RK4 yet runs much faster (not quite as good once the algorithm is made framerate independent, but it "tends towards" fourth order accuracy if the timesteps are kept constant) and allows computationally cheap addition of distance constraints (useful for rigid body, rather than particle, models). Be carefull though, since Verlet uses the expansion of f(x + \Delta t) and f(x - \Delta t), there are actually two time differences that need to be kept track of, this frame's and the last frame's.

If you have no stringent limits on data storage I would encourage precalculating the gravitational field(s) too, this way you can include all the ellipticity and lumpyness you want with no overhead at run-time. You could also include other information, such as a coefficient at each point for atmospheric drag, relativistic effects... etc. Basically any field quantity generated by the Earth in someway could be thrown in there for little to no extra slowdown since the Earth is essentially static from the perspective of your model. This may sound very static, but the only dynamic thing (other than your spacecraft /stations which aren't going to produce all of these effects) to model here really is the rotation of the Earth, which can be achieved by rotating the position to sample data from... If you want to include the effects of the Moon (an example, not necessarily justified) then you could have two sets of data which are sampled differently to account for the Earth's rotation... etc... although having said this, a Newtonian particle model will run happily in real-time using RK4 with several hundred particles at well over 60fps on an old computer.

I would also recommend (if you don't have them already) writing a vector and matrix class using the SSE instruction set and inline asm (i don't like the intel intrinsics but there is nothing wrong with them either). The latest version of compiler X will probably generate code with SSE/2/3 instructions if you want it to but the performance gain will not be as good as if you write a vector class to use them properly, the optimisers in compilers will spot the simple cases only, such as vectored addition, whilst missing less obvious candidates, like matrix inversion. Anyway... this may sound like some effort but if you plan to do physics simulations in C++ you will probably find yourself using them over and over again... I don't know of any particular library that does this (I wrote my own code as a learning experience), but there should be one out there considering how useful it is for graphics applications, specifically games.

Good luck with your research for this.

 

[D H]

One thing to bear in mind when writing real-time simulations is that the frame rate is not constant

This is just wrong. NASA, for one, writes many real-time simulations with a constant frame rate. The flight software runs at a constant rate. The simulation has to reflect that.

If you have no stringent limits on data storage I would encourage precalculating the gravitational field(s) too, this way you can include all the ellipticity and lumpyness you want with no overhead at run-time. You could also include other information, such as a coefficient at each point for atmospheric drag, relativistic effects... etc. Basically any field quantity generated by the Earth in someway could be thrown in there for little to no extra slowdown since the Earth is essentially static from the perspective of your model.

This is a very bad suggestion for an orbital simulation. The Earth is rotating. A non-spherical gravity field must be written in terms of an Earth-fixed frame; vehicle state must be represented in an inertial (or pseudo-inertial frame). A realistic atmospheric model must account for the diurnal atmospheric bulge, which approximate has the same rotation rate as the Earth. The Earth is anything but static.

I would also recommend (if you don't have them already) writing a vector and matrix class using the SSE instruction set and inline asm (i don't like the intel intrinsics but there is nothing wrong with them either).

Another very bad suggestion. Firstly, it ties the implementation to a specific architecture. Secondly, the only reason to go after performance is when performance matters. Tripling the efficiency of an infrequently-called function doesn't accomplish much. Thirdly, there are usually much better ways to improve the efficiency of an algorithm than write it in assembler.

 

[Jheriko]

...

You make several valid points but I think your criticism is a little undeserved.

The constant frame rate thing is interesting. Personally I use variable framerate for almost everything simply because in my experience I don't have enough resources to ensure a constant high framerate, and although 30 or 40 fps is acceptable to watch I would rather have it use 60 or 75 if it is possible. Having a constant framerate does simplify the calculations required and decreases the errors introduced by making an integration approximation use variable step sizes, so I can see the benefit of doing it, I am just unsure if it is a practical constraint.

As for the complaints about the fixed data I suggest you re-read my comments. Accounting for the rotation when sampling from a 3d grid of data is trivial and I mentioned this in my post. Recalculating gravitational potentials dynamically is more expensive than multiplication of a vector by a rotation matrix, especially if you wish to include 'lumpyness'.

The point on the SSE instruction set is interesting, although it doesn't necessarily tie anything to anyone architecture (#define, #ifdef etc..), and there are other sets for other CPUs (AltiVec for example). Its still worth researching since most of any simulation is vector maths and matrix manipulation... yes there are better ways to optimise but they require thought and analysis etc... another option which I forgot to mention before is openmp "#pragma omp parallel for" is quite handy and easy to use, although i would recommend it keeping out of the deeper parts of a program (e.g. vector addition).

 

[D H]

Personally I use variable framerate for almost everything simply because in my experience I don't have enough resources to ensure a constant high framerate, and although 30 or 40 fps is acceptable to watch I would rather have it use 60 or 75 if it is possible.

What are you simulating? A complex on-orbit scene (sun and backscattered illumination of multiple vehicles with moving solar arrays, semi-realistic Earth underneath) is incredibly expensive to generate. 10-15 Hz is a fast graphics rate for even a powerful computer dedicated to scene generation. The simulation itself runs at a much higher rate.

As for the complaints about the fixed data I suggest you re-read my comments. Accounting for the rotation when sampling from a 3d grid of data is trivial and I mentioned this in my post. Recalculating gravitational potentials dynamically is more expensive than multiplication of a vector by a rotation matrix, especially if you wish to include 'lumpyness'.

Unless the grid is very fine, the errors in gridding will dominate unless the grid is extremely fine. Have you examined the errors that result from using this approach?

The point on the SSE instruction set is interesting, although it doesn't necessarily tie anything to anyone architecture (#define, #ifdef etc..), and there are other sets for other CPUs (AltiVec for example). Its still worth researching since most of any simulation is vector maths and matrix manipulation... yes there are better ways to optimise but they require thought and analysis etc... another option which I forgot to mention before is openmp "#pragma omp parallel for" is quite handy and easy to use, although i would recommend it keeping out of the deeper parts of a program (e.g. vector addition).

Thought and analysis shows that most vector and matrix operations in an orbital simulation involve 3-vectors and 3x3 matrices. The most obvious optimization is to write inline functions or macros for 3x3 matrix operations with fully unrolled loop (i.e., compute each element explicitly). The pipelined architectures used by modern CPUs zip through such linear sequences of code quickly. Enable the compiler's optimizer will make this even faster.

Compare this to a set of generic matrix and vector hand-optimized in assembler. Function calls are very expensive. The stack needs to be loaded before the call and unloaded after the call. Moreover, the function call, the looping in the functions, and the return from the function stall the pipeline. The simple, non-assembly optimization will beat the pants off a set of assembler functions that use loops.

 

[Jheriko]

What are you simulating? A complex on-orbit scene (sun and backscattered illumination of multiple vehicles with moving solar arrays, semi-realistic Earth underneath) is incredibly expensive to generate. 10-15 Hz is a fast graphics rate for even a powerful computer dedicated to scene generation. The simulation itself runs at a much higher rate.

In "normal" graphics work 30fps is considered to be about the minimum that is acceptable, so I don't see how 10-15fps is fast. Maybe for a complex simulation it is fast, but it barely qualifies for the term "real-time". 60fps is often quoted as an ideal framerate for graphics applications.

Unless the grid is very fine, the errors in gridding will dominate unless the grid is extremely fine. Have you examined the errors that result from using this approach?

Yes, interpolating over a grid isn't as good as having the exact value but it is still a good optimisation in the sense that it removes the most expensive parts of the program from the render loop. Errors from using the grid will accumulate too you are correct, however with some good interpolation you can reduce the problems substantially (not totally of course).

Thought and analysis shows that most vector and matrix operations in an orbital simulation involve 3-vectors and 3x3 matrices. The most obvious optimization is to write inline functions or macros for 3x3 matrix operations with fully unrolled loop (i.e., compute each element explicitly). The pipelined architectures used by modern CPUs zip through such linear sequences of code quickly. Enable the compiler's optimizer will make this even faster.

Compare this to a set of generic matrix and vector hand-optimized in assembler. Function calls are very expensive. The stack needs to be loaded before the call and unloaded after the call. Moreover, the function call, the looping in the functions, and the return from the function stall the pipeline. The simple, non-assembly optimization will beat the pants off a set of assembler functions that use loops.

I'm sorry to disagree with this but I have written several programs where I need to do trivial operations with 4-vectors of 8-bit integers. Even for the simplest of cases like vector addition there is a noticable difference between a clean C++ implementation and an assembly routine using an MMX opcode. This is exactly the sort of operation that these opcodes were designed for in the first place. Ask anyone who has written blending code for software rendered graphics of some kind (even video players and such), MMX makes it substantially faster even for simple cases like additive or subtractive blending.

Where did I recommend loops in assembly? I think maybe you may have misinterpreted something that I said... I would also like to point out that your "most obvious optimisation" is how I would have written it in the first place... I don't consider that an optimisation, its just the obvious way to do it. Using a loop is obviously heavy handed, and more to the point, I never recommended it.

Also, what makes you think that I don't turn on the compiler's optimiser?

Function calls are very expensive

Function calls are very expensive? I used to have this same misconception, it took several discussions with a good friend of mine over the course of a year or so for me to shake it off.

Function calls are cheap and stack operations are relatively fast. A good compiler will inline a function call for one thing. Secondly there are __fastcall and similar non MS-specific specifiers which hint at the compiler to use the registers (although infact it only uses 2 because the compiler likes to use the others as it sees fit).

If we ignore these then we are left with the stack issues, which you can't avoid these by "letting the compiler optimise" or not using assembler... even a = b + c requires stack operations. Operators are after all just functions, and they use the same calling conventions as any other function (you can test this by manually examining the stack). The only way around that is to inline, which I (mistakingly) assumed was too obvious to bother mentioning explicitly. I just assume that people will inline or not in order to optimise their program, most compilers will even inline things for you where it is appropriate.

The fact is that in the case of a vector class with inline operators this all becomes a non-issue anyway, regardless as to whether or not function calls are slow, since there are no function call overheads and no stack operations for the parameters.

Maybe I jumped the gun a bit by suggesting asm code/libraries first, but at least it is one optimisation which you can apply from the beginning without too much thought, and it always makes the program faster... even if it is by some tiny amount.

Anyway... this has grown quite long and rantish... sorry if I seem stubborn, but I have been using these opcodes for quite a while and they do provide a substantial performance increase where you have to do lots of vector maths. If you don't believe me download the source code for any relatively recent 3d game, even a mod SDK will do since they often contain the math library. Media players are similar again, and I know that at least Winamp and WMP use MMX very heavily.

I will accept that perhaps my advice was bad though...

 

[whatta]

check out some open source stuff.

 

[TrickSimNewbe]

I'm learning how to use the Trick Simulation Environment for a simulation at work, and I am having trouble getting the stripchart feature to work. There is a call from

/trick/bin/tcl/stripchart.tcl
to
/trick/lib/blt-2.4/libBLT24-Linux_3.4_234.so

but it generates the following message:
"Error in startup script: couldn't load file: cannot open shared object file: No such file or directory"

I can see the file, but it may be missing some contents. Other things work ok.

Any clue where to find a fix?

 

[D H]

I'm learning how to use the Trick Simulation Environment for a simulation at work, and I am having trouble getting the stripchart feature to work. ...

Your error is something along the lines of

Error in startup script: couldn't load file "/trick/lib/blt-2.4/libBLT24-3.4_234.so" ... No such file or directory
    while executing
"load $library BLT"

I suspect you are working in an environment where different machines have different versions of Linux and/or gcc, and that you are trying to use Trick on a machine for which Trick has not been properly configured.

Some preliminary questions:

• Are you sure the specified .so file exists, or just one that has a very similar name?
• What is your TRICK_HOST_CPU environment variable?
• Does the file ${TRICK_HOME}/lib/blt-2.4/libBLT24-${TRICK_HOST_CPU}.so exist?