Wow, it was still running, but it did converge to form a line distribution:

Like real galaxies often form disc shaped spirals.

I often wonder why this is?

At some point, i came up with some answer, but i forgot it.

Do all planets of a solar system tend to orbit on the same plane? I guess it's some net effect of all forces, which create some order from chaos. And the simplest kind of order from chaos is to reduce one dimension?

There is a book called Galactic Dynamics. You might like it!

P.S. I googled for “why are galaxies flat”, and it comes up with an answer. ? I assume the same happens for protoplanetary discs, although there is a big difference — galaxies are bound by gravitation, protoplanetary discs are bound by gravitation and electromagnetism.

@joej Hahah, I just saw this comment. Nice.

Also, I think I see what you are doing now. I missed the addition of your calcPairForce . Now I think I understand how this reduces calculations. 

taby said:
“why are galaxies flat”

Google:

Galaxies are flat because of their rotation. All of the stars, planets and other objects in a galaxy are rotating around the core of the said galaxy, and the conservation of angular momentum allows these objects to spread outward, but not any other direction, which is why they are flat.

Wolfram Alpha: No answer? I expected more…

I'd need to check more of the search results. It's clear that if points form a plane, the plane acts like an averaged attractor, and moving away from the plane costs extra energy.
But why does chaos converge to planes in the first place? ‘Rotation’ is no answer, but itself a result of such processes: Natural linear motion plus some distance constraint causes circular trajectories.

I guess it's this: Imagine we have stars uniformly distributed on the surface of a sphere. The randomness is so perfect, any trial to calculate an average plane direction gives us a zero vector.
Then we perturb the position of just one star by a small random offset. So the uniformity is broken and our average plane direction has a tiny magnitude.
This will now cause all stars to move a tiny little bit towards this given plane. And over time, the universe becomes more and more planar. Ofc. it would just collapse as we had no initial velocities, but before that, it would show some flattened elliptical shapes. that's somehow how i imagine it.

If i would understand better, i would maybe also understand how to calculate the average plane from a set of random points. I rarely needed this, so i don't understand any of the methods a had found.
If i had time to read a book, i would prefer a math book over galactic dynamics, as this would give me more practical knowledge.
But to read a math book, they assume you already know math. <:\

jonny-b said:
Now I think I understand how this reduces calculations.

You are done when you realize that this does not reduce calculations, but you had found a way to increase calculations for no benefit.

Ofc. everybody starts with your approach, because it requires the least thinking about it. So that's naturally human. But every experienced programmer will immediately think you're a beginner.
So it's not just about efficiency, but i was indeed a bit serious with my comment.
My impression is you're not a beginner (no shame if so), but for some strange reason n-body problems never occurred to you before, although it's quite generic and appears in many non simulation problems too.
I think you really want to adopt this faster loop for the future, never ever using the naive loop again.

We both found some enlightenment here. I work quite a lot on simulation stuff, and it's embarrassing to me i did not get this ‘zero net acceleration inside the planet thing’ on my own long before. So we're in the same boat.

@JoeJ I'm a web developer. These kind of problems are not something that I encounter or solve on a daily basis. I deliberately took on this little project as a way to push myself into a space that I don't normally get to go (pun intended?).

@joej Also, nice open faced compliment sandwich ( ; . I know I've said it before but I'll say it again. I really appreciate you sticking with this thread for so long. It has been a joy. I plan on continuing working on this little project for a while but now that the simulation is complete (more or less) I will being doing more web design and UX improvements. However, I will keep stopping back from time to time.

@joej Are you familiar at all with linear regression?

jonny-b said:
I'm a web developer.

I worked as web dev for one year. Iirc, such n-body problem appeared at least once, when the problem was to place some set of banners, available at two image sizes each, fitting a variable page size but picking larger images for more important banners.

I guess you'll often think back on me in the future, no matter what you work on.

But i don't want to be pedantic. It's just for good i'm sure. ; )

taby said:
@joej Are you familiar at all with linear regression?

Yes, a bit. I used it to calculate the average line from a point set. It worked using some least squares method, and i did not understand how exactly.
But then i found out it does not work well at all. It works best if the resulting line is aligned to the x or y axis. So i made an iterative solution, by rotating the point set to the line result of the former iteration. This kind of worked, but it was too slow. I saw this method often on the internet, but nobody pointed out it's accuracy depends on the orientation of input. It seemed most people were just unaware.

So i came up with a better solution: Calculate the averaged center of the points, then convert vectors from there to the points to angles, multiply all angles by two, calculate their averaged angle, finally divide the averaged angle by two, which gives the perfect average line i want in one step. We can get rid of trig by using complex numbers or vectors instead angles, and we can also have a weight for each point.
This is awesome. I use it to calculate curvature directions on meshes, crossfields, and all that remeshing stuff.

But it does not work in 3D. To do the same idea in 3D, we can use spherical harmonics. The third band of SH represents lines. But i found no way to extract the local maxima of the spherical function analytically. So i have to search for it with iterations, which is again slow. There must be a better way.

There was a thread here, about orienting bounding boxes to a good fit. Very similar problem, and @vilem otte gave some code to a solution based on matrices. I tried this, but could not get it to work.
Similar matrix math happens in my MPM fluid simulator. I do not fully understand the used math yet, but it somehow accumulates the adjacent velocity field to 3x3 matrix rows and columns, and then, doing a singular value decomposition we get a matrix which represents an ‘ellipsoid’ to describe the deformation of surrounding space due to pressure. That's awesome stuff, because SVD is pretty fast so we can do this for massive data.
I need to learn how this works in detail, so i can use it for other things… when i need it.

That's where i am.
@alvaro pointed out there exist regression methods which give accurate results in one step, iirc.
And David Eberly's codebase has many examples to fit certain shapes from points. Whenever i was stuck at math, his code did help me the most. I could not thank him enough, and would be still stuck at PacMan otherwise.

So i make my way, but i still lack many basics, like eigen vectors and values. And i lack experience with expressing problems in matrix form. I reinvent much more wheels than necessary.

Ok, then I’m sure that there is a way to do planar regression.