fix physics

[FrMo]

Hi Gerald, here is my explaination why i think that the physics in gravit has some flaws currently.

Cheers, Frank.

min. Code Changes

Basicially, the bug can be corrected by changing this line in frame-pp.c:

        force = state.g * pd1->mass * pd2->mass / inverseSquareDistance;

to

        force = state.g  * pd2->mass / (inverseSquareDistance * sqrt(inverseSquareDistance));

While experimenting, i found that the following formula gives realistic simulations, too:

        force = state.g  * pd2->mass / inverseSquareDistance;

The effect is, that galaxy disks become a lot more stable. Not sure why...

Why wait ??

Because this little change has a lot of implications:

• the current value of G does not work any more, and the simulation will get much more sensitive to this value
• the current spawn scripts must be changed (involves more physics know-how...)

  Some Physics background:

http://en.wikipedia.org/wiki/Newton%27s_law_of_universal_gravitation

so the force between two particles is defined as: F= G * m1 * m2 / (distance^2)

How do we get the acceleration "a" from this ?

• F = m1 * a --> a = F / m1
• --> a = G * m2 / (distance^2)

For vectors, we need to solve another small problem; the acceleration vector has size a, and must point into the same direction as the distance vector. The new acceleration vector is:

• a_vector = a * distance_vector / length(distance_vector)
• --> a_vector = distance_vector * G * m2 / ( (distance^2) * distance)

After summing up all acceleration for each particle, we only need to to move each particle to it's new position:

• v_new = v_old + deltaT * a
• p_new = p_old + detaT * v_new

assuming detaT = 1, this part in gravit should be correct. But there is potential for improvements, as there are different possible "integration functions" that reduce the error here (leapfrog, etc). Also, the multiplication with G can be moved outside of the particle loop, and done when computing v_new.

To avoid "division by zero" and "stars shooting out of the galaxy" effects, is is a common trick to add a small number (less than the min. distance between any two particles) to the distance; sometimes this is called "smoothing length" or "softening length". The smoothing length has almost no effect for larger distances, and avoids problems with very small distances.

[gak]

Hey mate. I'm pretty busy right now, so I don't have time to check this out. I'll leave the ticket open for now.

The plan is to verify the info, implement it and create new spawn calculations.

[FrMo]

Hi,

i have found some physics papers describing our subject. Even if you skip the math, the code examples are relatively easy to read.

BTW, I am still surprised that my "a = G * m2 / distance" (instead of "a = G * m2 / (distance^2)" ) gives very stable, realisticly-looking simulations, too (similar to MOND, dark matter, or other effects that happen in reality ???). It is also faster, because we do not need a square root in the code. Maybe we could add this as an optional "almost-proper" mode of physics.

cheers, Frank.

http://arxiv.org/abs/astro-ph/0511062 "Performance Tuning of N-Body Codes on Modern Microprocessors" by Keigo Nitadori, Junichiro Makino, Piet Hut

http://academic.research.microsoft.com/Publication/6559173/direct-n-body-kernels-for-multicore-platforms "Direct N-body Kernels for Multicore Platforms" by Nitin Arora, Aashay Shringarpure, Richard W. Vuduc

http://arxiv.org/abs/0806.3950 "Gravitational N-body Simulations" (summary review of the field) by M. Trenti, P. Hut

http://www.artcompsci.org/msa/web/vol_1/v1_web/v1_web.html "Moving Stars around" - a step-by-step intro on building n-body simulations