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

Question

Important note: this question has utterly no relation to "PhysX", which is a computer-game-physics system (useful for the physics in arcade games such as ball games, etc); PhysX is a system built-in to Unity3D and other game engines; PhysX is totally irrelevant here.

//////////////////// UPDATE (read bottom first) /////////////////////

I have been logging the values and searching where the exact problem is, and i think i found it. I have something like this in my code

Velocity += Acceleration * Time.deltaTime;
position += Velocity * Time.deltaTime;

The acceleration is something like 0,0000000000000009.. right now. As the code flows, velocity increases as it should be, no problem with the float. But in the beggining, the initial position of earth is (0,0,23500f) You can see this in the chart i gave at the end.

Well now when i add speed * timedelta (which is something like 0,00000000000000005 at this point) to position which is 23500, it basicly doesn't adds it. position is still (0, 0, 23500) not something like (0,0, 23500.00000000000005), thus the earth doesn't move, thus the acceleration doesn't change.

If i set the initial position of earth to 0,0,0 and still, set the acceleration to 0.0000000000000000009 to assuume it's position is (0,0,23500) It then "ADDS" the velocity * timedelta. It becomes something like (0,0,000000000000000000005) and keep increases. When the float is 0, there is no problem with adding such small value. But if the float is something like 23500, then it doesn't adds up the small values.

I don't know if it's exactly unity's problem or c#'s float.

And that's why i can't make it work with small values. If i can overcome this, my problem will be solved.

///////////////////////////////////////////////////////////////////////////////

i have been develeping n-body phyics to simulate our solar system, so i have been gathering data around to make it as realistic as possible. But there is a problem with the data size. I searched every tiny bit of internet and i couldn't find a single explanation how people overcomes this. (If they so) So i trying my shot here.

So, to keep the ratio of distance, radius and "mass" between planets fixed, i created an excel file to calculate all the datas. (Because why the heck would someone put "what would be the earth's mass if it had "that" radius chart" on the internet?) I will give the ss as attachment. It basicly "normalizes" or in other words "scales" every property of a planet to a given reference. In this case, i took the reference as "earth's radius."

I work in unity, and you know, you can't work with "too big" or "too small" values in unity. So i had to scale the solar system down, "a lot!"

So i use the Newton's law of universal gravitation, which is F = GMm/r^2, to make it simple, i am directly calculating the a = GM/r^2, for a given body from all other bodies.

So, the real value of earth's gravitational acceleration "towards sun" is roughly 0,000006 km/s^2, which is even incredibly small value to work with in unity, but it could work. Yet, to get this value,1 i need to set earth's radius (scale) to 6371 unit, and sun to scale of 696,342!, which is TOO big to render it in unity.

So i said, let the earth's radius be 1, in unity units. So, when the radius changes, everything changes, the mass, the distance... I kept the density of the planet and calculate the mass from the new volume with the new radius. All the calculations are in the attachment.

So the thing is, when i take the earth's radius as 1, the gravitational accelaration towards sun becomes is something like 0,0000000000009 which is ridiculously small. And of course Unity doesn't work with that value.

So, if i increase the earth's radius instead, then the mass and radius of the Sun gets ridiculously big and then again, i can't work with it.

I don't know how other people fixed this, what they did to overcome this problem but as i see from here, it looks impossible to make realistic n-body simulation of solar system. (in unity atleast)

So i need to have 10 rep to post images -_-, i will give link instead. http://berkaydursun.com/solar_system_simulator/data.PNG Also one directory up is the working experimental solar system simulation with n-body calculations but with UNREALISTIC values. It works quite well, and it even looks somehow close to real, but no, it doesn't have the right ratios ^^ You can test it here if you wish http://berkaydursun.com/solar_system_simulator/

Edit: WoW I almost started every paragraph with "So" ^^

Answer 1

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.