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" ^^

Freed

consider that physics engines for games aren't really trying to model the real world, they approximate. You can not expect to apply real world physics formula and get a precise result or even just something that comes very close. Once you understand and accept that, you will find that making it look believable by tweaking values and formulas as needed will be the better approach for making a game. If instead you want an actual simulation, you need software that actually tries to simulate real world physics rather than just emulating it to the point were the results seem familiar. – Farseeing

Sorry, deleted my comment on PhysX as I realized this question has absolutely no connection to using game engine physics. Cheers. – Monoplane

Maybe you might want to have a look at how they have done physics in KSP, here is a link : github.com/mhoram-kerbin/ksp-physics-documentation/releases (they are using unity if I remember correctly) – Gisela

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.

[Edit4]

Here link to win32 binaries of mine effemerids:

• Effemerids v1.118 + SVK 200m DTM

be sure to read the astro_ephemerides_new.txt for key shortcuts, note that efemerids are from ~1990 so the bodies positions will be slightly off and also moons are not working properly as they need different approach (keppler is more or less useless for them)

Project was written in win32 Borland/Embarcadero environments and needs proper OpenGL gfx driver to work. All the bodies are configured by ini file so you can tweak what you want and used GLSL shaders are also present as files...

Abaft

But doesnt Kepler's equations "assumes" that there is an orbit? What in my case there is no necessity to have an orbit. A body should be able to take sling shot from another mass object, which you can't achieve with kepler's equations, right? - I strictly want to make it in n-body. I am actually working in n-body, and to test it, i want to build solar system out of it. And the problem is not with randering but calculating with incredibly big values of mass, and incredibly small values of acceleration. Well, unity doesn't account the value of 0,00000000000000009f, it basicly returns 0. – Freed

Also, why gravity model increases in error with time? With the fake valued solar system simulation i make (link), it works quite consistent. There is a problem in the trajectory simulator but that is because the algorithm i wrote for it, not the n-body algortihm. I will fix it soon. – Freed

@BerkayDrsn 1. Yes Kepler's equation can not be used for slingshots and non periodic movements. 2. gravity model has increasing error due to continuous summing integration error (it could look realistic but when you compare it to the real thing then it does not match) 3. I strongly recommend not to use unity for such computation instead do it yourself but you need to use some kind of higher precision fixed/floating point numbers and small enough integration dt step. There are most likely some astronomical libs for that too. Unity was designed for games not for realistic simulation – Abaft

I don't use Unity's Physics, no rigidbody at all. I have my own physics code, which is similar to yours in the link. But the problem still exists. I have been logging recently, and i figured where the exact problem is. I have something like Velocity += Acceleration * Time.deltaTime; position += Velocity * Time.deltaTime; in the code, (vector3) And the acceleration is something like 0,0000000000000009.. right now. As the code flows, velocity increases as it should be. But in the beggining, the initial position of earth is (0,0,23500) You can see it in the chart i gave in the question. – Freed

Well now when i add speed * timedelta which is something like 0,00000000000000005 to 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. I don't know if it's unity problem or c#'s float problem. What should i do? – Freed

@BerkayDrsn you need to have at least 2 floating variables per each variable for example velocity=vel0+vel1; where vel0 are for example velocity part up to 10.000 and the vel1 is the rest. now you do vel0+=acc*dt; if (fabs(vel0>=10.0)) { vel1+=vel0; vel0=0.0; } you can add any number of variables and choose the range so it has as low error as it can and as long overflow as it can ... this way you eliminate the summing small and big numbers rounding error – Abaft

I couldn't understand the splitting method. Can you explain it a little more in detailed? ^^ – Freed

@BerkayDrsn added [edit3] with some precision tweaks I come up with during my attempts of Kepler <-> Newton d'Lambert orbit conversions finally I got it working :) (just for 2 rigid bodies for now) – Abaft

@Freed hi take a look at this How to emulate double precision using two floats its much more precise than just ranging floats ... (just added the link to answer too) – Abaft

@Abaft do you have a link to your project? Would love to see it in action – Pellet

@Pellet added link to edit4 ... – Abaft

Answer 2

Scaling things down will not necessarily help, as you have discovered. Here is some good reading on things to consider when using floating point numbers: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Basically doing a simulation from first principles (Newton's laws) is bad for numerical accuracy, because you don't imbue the numerical methods with an idea of the scale of important effects, so you end up throwing a whole bunch of different effects at different scales together and the result is low accuracy.

Usually things like ephemerides for planets, satellites etc don't start with Newton's law, They start by assuming that the orbits are keplerian, then and small perturbative corrections.

Here is an algo that calculate the position of the planets (semi-imperically). http://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf

If you want to do an N-body simulation you it seems like will need more precision. If unity prevents you from using double precision, then I suggest doing the calculations in plain C#, then convert to single precision when the job is done.

Thankless

Well as i said in other other answer, kepler equiations assumes the body has an orbit around the mass object. But in my case, there is no necessity to have one. For example, a body should be able to take sling shot from the mass object. I can't do that with kepler, right? I strictly want to use N-body. And about the last part, can you tell me how to do that? – Freed

To me it would appear that folks maybe should use a Kepler function for circular orbits (planets, moons), and for smaller objects that accelerate and have non-circular orbits, flag those objects and maybe switch to a high frequency Newtonian function with a high sample rate (maybe the NASA code) when they are nearby to other objects then back to Kepler in the deeper space? But I am taking a guess, thanks for the thoughts! – Andraandrade