Is it possible to make realistic n-body solar system simulation in matter of size and mass?

I did program sol system simulation too so here are mine insights:

1. rendering

   I use OpenGL with 1:1 scaling. All units are in SI so [m,s,kg,…]. The problem gets started with Z-buffer. The usual Z-buffer bit wide is 16/24/32 bit which is nowhere near what you need. I am rendering from 0.1m up to 1000 AU so how to overcome this?

   I did manage it by rendering with 3 frustrums at once combining Z-sorting and Z-buffering (Z-sort is necessary because of transparent rings …and other effects). So first I render most distant parts up to zfar=1000AU. Sky dome is projected at z=750AU distance then clear the Z-buffer and render objects up to zfar=0.1AU. Then clear the Z-buffer again and render close objects up to zfar=100000 m.

   To get this work you have to have as precise projection matrix as you can. The gluPerspective has unprecise cotangens so it needs to repair concerned elements (get me a long time to spot that). Z near value is dependent on the Z-buffer bit width. When coded properly then this works good even with zoom 10000x. I use this program as navigation/searcher of objects for mine telescope 🙂 in real time from mine home view. I combine 3D stars, astro bodies, ships, real ground (via DTM and satellite texture). Its capable even of red-cyan anaglyph output :). Can render from surface,atmosphere,space … (not just locked to Earth). No other 3th party lib then OpenGL is used. Here is how it looks like:

   
   
   
   
   

   As you can see it works fine on any altitude or zoom the atmosphere is done like this atmosphere scattering shader

2. simulation

   I am not using n-body gravity simulation because for that you need a lot of data that is very very hard to get (and almost impossible in desired precision). The computations must be done very precisely.

   I use Kepler’s equation instead so see these:

   • Solving Kepler’s equation
   • C++ implementation.

   If you still want to use gravity model then use JPL horizons from NASA. I think they have also source codes in C/C++ there but they are using incompatible reference frame with mine maps so it is unusable for me.

   In general Kepler’s equation has bigger error but it is not increasing with time too much. The gravity model is more precise but its error is rising with time and you need to update the astro body data continuously to make it work …

[edit1] integration precision

your current implementation is like this:

// object variables
double  acc[3],vel[3],pos[3];
// timer iteration
double dt=timer.interval;
for (int i=0;i<3;i++)
 {
 vel[i]+=acc[i]*dt;
 pos[i]+=vel[i]*dt;
 }

The problem is when you are adding very small and very big value then they are
shifted to the same exponent before addition which will round off significant data
To avoid this just change it to this:

// object variables
double          vel0[3],pos0[3]; // low
double          vel1[3],pos1[3]; // high
double  acc [3],vel [3],pos [3]; // full
// timer iteration
double dt =timer.interval;
double max=10.0; // precision range constant
for (int i=0;i<3;i++)
 {
 vel0[i]+=acc[i]*dt; if (fabs(vel0[i]>=max)) { vel1[i]+=vel0[i]; vel0[i]=0.0; } vel[i]=vel0[i]+vel1[i];
 pos0[i]+=vel[i]*dt; if (fabs(pos0[i]>=max)) { pos1[i]+=pos0[i]; pos0[i]=0.0; } pos[i]=pos0[i]+pos1[i];
 }

Now xxx0 is integrated up to max and the whole thing is added to xxx1

The rounding is still there but it is not cumulative anymore. You have to select max value that the integration itself is safe and also the addition xxx0+xxx1 have to be safe. So if the numbers are too different for one spliting then split twice or more …

• like: xxx0+=yyy*dt; if (fabs(xxx0>max0))... if (fabs(xxx1>max1))...

This however does not fully use the dynamic range of added variables to their full potential. There are alternatives that does it like this:

• How to emulate double precision using two floats

[Edit2] Stars

• The Colors of the Stars The best star visualization I ever saw
• Star B-V color index to apparent RGB color all star catalogues uses B-V index
• Using Stelar catalogs also star names cross-reference link is there
• Skybox: combine different star data

[Edit3] Improving Newton D’ALembert integration precision even more

The basic problem with iterative integration is that applaying gravity based acceleration based on current body position will result in bigger orbits because durring integration step dt the position changes a bit which is not accounted in the naive integration. To remedy this take a look at this picture:

Let assume our body is at circular orbit and in the 0 deg position. Instead of using acceleration direction based on current position I used position after 0.5*dt. This augment the acceleration small bit resulting in much much higher precision (correspondence to Kepler orbits). With this tweak I was able to successfully convert from Kepler orbit into Newton D’Alembert for 2 body system. (doing this for n-body is the next step). Of coarse correlating with real data from our solar system is only possible for 2 body system not affected by tidal effects and or moons. To construct own fictional data you can use Kepler circular orbit and contripedal force equalizing gravity:

G = 6.67384e-11;
v = sqrt(G*M/a);                           // orbital speed
T = sqrt((4.0*M_PI*M_PI*a*a*a)/(G*(m+M))); // orbital period

where a is the circular orbit radius m is body mass , M is focal body mass (sun). To maintain precision in acceptable tolerance (for me) the integration step dt should be:

dt = 0.000001*T

So to put new body for testing just put it at:

pos = (a,0,0)
vel = (0,sqrt(G*M/a),0)

While main focal body (Sun) is at:

pos = (0,0,0)
vel = (0,0,0)

This will place your body in circular orbit so you can compare Kepler versus Newton D’Alembert to evaluate precision of your simulation.