How to Make Realistic N-Body Solar System Simulation in Matter of Size and Mass

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.

N-Body Gravity / Solar System Javascript Simulation

Using NASA values is not a problem when using separate coordinates for drawing. Using an appropriate linear transfomration from real coordinates to screen coordinatees for displaying does not influence the physical values and computations.

For simulating the motion of a planet with iterative updates one can assume that the gravitational force and the velocity are constant for a small portion of time dt. This factor dt is missing in your conversions from accelration to velocity and from velocity to distance. Choosing an appropriate value for dt may need some experiments. If the value is too big the approximation will be too far off from reality. If the value is too small you may not see any movement or rounding errors may influence the result.

For the beginning let us assume that the sun is always at (0,0). Also for a start let us ignore the forces between the planets. Then here are the necessary formulas for a first not too bad approximation:

• scalar acceleration of a planet at position (x,y) by the gravitational force of the sun (with mass M): a = G*M/(d*d) where d=sqrt(x*x+y*y). Note that this is indepent of the planet's mass.
• acceleration vector: ax = -a*x/d, ay = -a*y/d (the vector (-x,-y) is pointing towards the sun and must be brought the length a)
• change of the planet's velocity (vx,vy): vx += ax*dt, vy += ay*dt
• change of the planet's position: x += vx*dt, y += vy*dt

Unable to get a proper orbit on my solar system model

the first problem is usually caused by too big time step of the simulation in conjuction with absent collision handling. when your objects get close the forces get big and the incremental step in simulation became too big so the next iterated position is after the collision and usually further away then before it so the breaking force is smaller leading to bigger and bigger orbits in time...

To remedy that you can do:

• add collision handling
• limit max speed or force to some sane limit

I never encountered your second problem and without code I can not even guess ... other than rounding errors

Take a look at this:

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

an also all sublinks in there ...

[Edit1] after I saw your code

The architecture of the code looks OK the problem is that your equations are slightly off. it should be:

vel+=acc*dt;
pos+=vel*dt;
acc=0.0;

instead of yours:

pos+=vel;
vel+=acc;
acc=0.0;

so you got wrong order and missing the *dt where dt is the time step. Because of this no matter how you change your timer interval the result is the same (just slower/faster) and also the direction of acceleration is applied one step latter than should causing the spinning of orbit (because acceleration was computed from different position than it was applied to the final position so its direction is always off).

---

Related Topics