Reasonably Fast N-Body Gravity Simulation

Simulate 2D N-Body systems using the O(N log N) algorithm of Barnes & Hut to compute forces between particles efficiently. MEX accelerated.

This code enables, depending on your hardware, the simulation of tens of thousands to millions of particles in real time. For this the Barnes-Hut algorithm (https://doi.org/10.48550/arXiv.1301.3342) was implemented using MEX functions. For optimal performance recompile (to profit from -march=native) the C++ functions to MEX using the prepared file compileMEX.mlx.

Features:

- Barnes-Hut gravity simulator using optimized MEX functions written in C++ (also included)

- Velocity-Verlet integration scheme

- Convienient initialization of particle clusters as shapes like rotating discs or squares, triangles etc.

- Plot kinetic and potential energy to check for energy conservation

- Make a video using the MP4 writing function.

- Matlab functions used as template for quadtree building and traversal are still inculded to better understand the C++ files.

The main idea is to construct at each time step a quadtree to store the particles in nodes down to single leafs in which a single particle resides:

To compute the forces we traverse the tree for each particle to check at each level if the criterion is fulfuilled. It checks basically if the node/cluster is far away enough to treat it as a single big particle with the sum of mass of all particles in that node and located at the node's center of mass. Depending on the opening angle theta one can trade accuracy vs speed. The following pcitures illustrate the effect of varing theta from 0.3 to 0.5 to 1 from left to right with the yellow dot being the currently considered particle.

The code implements softened gravity (Plummer potential) where the force is: . Where epsilon is the softening parameter to avoid exploding forces at very close encounters. In the limit when epsilon goes to zero we obtain the correct Netwonian gravity equation.

See demo:

https://www.youtube.com/watch?v=BHtQJeBWOt8&lc=UgyomFsZFVz3DpRVvUR4AaABAg

Optimization notes:

Tree builder is now about 90x faster than the first version!

Force computation was pretty hard to further optimze. And credit where credit is due: Claude that crazy mofo came up with the Morton sort. BTW you might want to recompile using the radix sort instead of std::sort in the C++ file for the force. Just uncomment to test the radix sort which can be faster depending on hardware.

Ah and the colormap infinity.m was taken from https://github.com/ThomasBury/scicomap