GPU-Accelerated N-Body Galaxy Collision Simulator

Galaxy mergers are one of the dominant drivers of galactic structure evolution, implicated in triggering starbursts, growing massive elliptical galaxies, and ripping long, dramatic tidal tails out of spiral disks.

While researching stellar dynamics and applied computer graphics, a compelling computational challenge emerged: Can we mathematically model and simulate these massive $O(N^2)$ gravitational interactions interactively, in real-time, on consumer hardware without relying on pre-rendered frames?

Today, I am detailing the architecture of a custom N-Body simulator built using OpenGL 4.3 and C++. Modeled after the Medusa galaxy (an ongoing spiral-elliptical merger), the simulation evolves over 110,000 particles in real-time.

🪐 The Physics: Generating Stable Initial Conditions

If initial velocities and positions aren't perfectly balanced against the system's potential energy, the galaxy will instantly fly apart or violently collapse.

To solve this, I engineered a pipeline sampling initial conditions mathematically derived from Per Bjerkeli's 2007 MSc thesis, treating the galaxies as collisionless systems.

The Elliptical Galaxy (Hernquist Model):

A spherical distribution that models the dense core and fading outer edges of an elliptical galaxy, using the Hernquist density profile:

1. Particle Distribution (Inverse Transform Sampling):

To generate the positions of particles, the following equation is used.

$\frac{M(r)}{M}=\left(\frac{r}{a}\right)^3\left[1+\left(\frac{r}{a}\right)^2\right]^{-\frac{3}{2}}$

where,

• $M(r)$ = Mass enclosed within radius $r$.
• $M$ = Total mass of the galaxy.
• $r$ = Radial distance of a particle from the center of the galaxy.
• $a$ = Scale radius of the Plummer sphere.

By drawing random numbers for the mass fraction, different radii can be calculated. Each particle is then randomly distributed on the corresponding sphere. This equation is based on the Plummer distribution; the same approach is now applied to the Hernquist distribution equations below.

The Density Distribution:

The Hernquist density distribution is given by

$\rho(r)=\frac{M}{2\pi}\frac{a}{r(r+a)^3}$

where,

• $\rho(r)$ = Mass density at radius $r$.
• $M$ = Total mass of the galaxy.
• $a$ = Scale radius.
• $r$ = Radial distance from the galactic center.

The Mass Distribution:

$\frac{M(r)}{M}=\frac{r^2}{(r+a)^2}$

where,

• $M(r)$ = Mass enclosed within radius $r$.
• $M$ = Total mass of the galaxy.
• $r$ = Radial distance.
• $a$ = Scale radius.

Solving the above equation by setting

$\frac{M(r)}{M}=u$

(where $u$ is a uniformly distributed random number in $[0,1]$) gives

$r=a\frac{\sqrt{u}}{1-\sqrt{u}}$

This is directly implemented in the code to generate particle radii.

Then, to generate the coordinates, we need $x$, $y$, and $z$, where

$z=r\cos(\theta).$

Therefore, we randomly generate a value of $z$ in the interval

$-r \le z \le r.$

Next, to generate $x$ and $y$, we use

$x=\sqrt{r^2-z^2}\cos(\phi)$

and

$y=\sqrt{r^2-z^2}\sin(\phi),$

where $\phi$ is generated randomly in the interval

$0 \le \phi < 2\pi.$

2. Velocity Distribution

For each particle, a random number between 0 and the escape velocity is drawn, and the total energy for that velocity is calculated. An acceptance-rejection method then keeps values falling under the distribution function and rejects those above it, evaluated over the possible energy interval for the current radius.

The Distribution Function

$f(E)=\frac{M}{8\sqrt{2}\pi^3a^3v_g^3}\;\frac{1}{(1-q^2)^{\frac{5}{2}}}\left[3\arcsin(q)+q\sqrt{1-q^2}\left(1-2q^2\right)\left(8q^4-8q^2-3\right)\right]$

where,

• $f(E)$ = Distribution function corresponding to the particle energy $E$.
• $E$ = Total specific energy (energy per unit mass) of the particle.
• $M$ = Total mass of the galaxy.
• $a$ = Scale radius of the Hernquist profile.
• $v_g$ = Characteristic velocity of the system.
• $q$ = Dimensionless energy parameter.

The parameters $q$ and $v_g$ are given by

$q=\sqrt{\frac{-aE}{GM}}$

and

$v_g=\sqrt{\frac{GM}{a}}$

where $G$ is the universal gravitational constant.

The escape velocity is given by

$v_e(r)=\sqrt{\frac{2GM}{r+a}}$

where $v_e(r)$ is the escape velocity at radius $r$.

The velocity components in $x$, $y$, and $z$ are then calculated the same way as the position on the sphere.

Disk Galaxy (followed Section 6.3 in the Thesis)

First, the Miyamoto model initializes position and velocity of particles in the disk; then the Hernquist distribution model spawns the dark matter halo that holds the galaxy together.

1. Miyamoto Particle Position Distribution (Acceptance Rejection Sampling)

It describes both the disk and the bulge of the galaxy.

The Density Function

$\rho(r,z)=\frac{b^2M}{4\pi}\frac{ar^2+(a+3\sqrt{z^2+b^2})(a+\sqrt{z^2+b^2})^2}{\left(r^2+\left(a+\sqrt{z^2+b^2}\right)^2\right)^{\frac{5}{2}}\left(z^2+b^2\right)^{\frac{3}{2}}}$

where,

• $\rho(r,z)$ = Mass density at cylindrical coordinates $(r,z)$.
• $M$ = Total mass of the disk galaxy.
• $r$ = Cylindrical radial distance from the galactic center.
• $z$ = Height above (or below) the galactic disk.
• $a$ = Radial scale length of the Miyamoto-Nagai disk.
• $b$ = Vertical scale height of the disk.

Using the cylindrical coordinate system, $z$ is generated randomly from

$(-z_{\max},z_{\max}),$

where

$z_{\max}=5b,$

then the random value of $r$ from

$(0,r_{\max}),$

where

$r_{\max}=10a,$

then the probability is found with the probability function

$P(r,z)=r*\rho(r,z)$

where $P(r,z)$ is the acceptance probability used for rejection sampling.

This probability is compared with a random value $u \in (0,P_{\max})$; if $u<P(r,z)$, the coordinates are kept and the particle spawned there. The angle $\theta$ is then generated from

$(0,2\pi)$

and

$x=r\cos(\theta), \quad y=r\sin(\theta).$

2. Miyamoto Particle Velocity Distribution

First let us define the disk potential.

$\phi_d(r,z)=-\frac{GM_d}{\sqrt{r^2+\left(a_d+\sqrt{z^2+b_d^2}\right)^2}}$

where,

• $\phi_d(r,z)$ = Gravitational potential due to the Miyamoto-Nagai disk.
• $G$ = Universal gravitational constant.
• $M_d$ = Mass of the disk.
• $a_d$ = Radial scale length of the disk.
• $b_d$ = Vertical scale height of the disk.

Now Halo potential.

$\phi_h(r,z)=-\frac{GM_h}{r+a_h}$

where $\phi_h$ is the Hernquist halo potential, $M_h$ the halo mass, and $a_h$ the halo scale radius.

Now calculate the $v_c$ by this formula and substitute $z=0$.

$v_c=\sqrt{r\left(\frac{\Delta\phi_d}{\Delta r}+\frac{\Delta\phi_h}{\Delta r}\right)}$

where $v_c$ is the circular velocity at radius $r$.

Now kappa.

$\kappa(r)=\sqrt{\frac{\Delta^2\phi_d}{\Delta r^2}+\frac{\Delta^2\phi_h}{\Delta r^2}+3\frac{v_c^2}{r^2}}$

where $\kappa(r)$ is the epicyclic frequency.

Radial Velocity Distribution

$\sigma_r=q\;\frac{3.36\,G}{\kappa}\;\frac{M_d}{2\pi a^2}\;\exp\left(-\frac{r}{a}\right)$

where $\sigma_r$ is the radial velocity dispersion and $q$ the Toomre stability parameter.

Vertical Dispersion

$\sigma_z=\sqrt{\pi G b\;\frac{M_d}{2\pi a^2}\;\exp\left(-\frac{r}{a}\right)}$

where $\sigma_z$ is the vertical velocity dispersion.

Azimuthal Dispersion

$\sigma_\phi=\sigma_r\frac{\kappa}{2\Omega}$

where $\sigma_\phi$ is the azimuthal velocity dispersion and $\Omega$ the angular frequency.

Mean Stream Velocity

$v_\phi=\sqrt{v_c^2+\sigma_r^2\left(1-\frac{\kappa^2}{4\Omega^2}-\frac{2r}{a}\right)}$

where $v_\phi$ is the mean azimuthal streaming velocity.

$v_r$ (radial velocity) is drawn randomly between $(0,\sigma_r)$, $v_z$ (vertical velocity) between $(0,\sigma_z)$, and the azimuthal velocity is randomly generated from the mean stream velocity.

The velocity components in $x$ and $y$ direction are then

$v_x=v_r\cos(\theta)-v_\phi\sin(\theta)$

$v_y=v_r\sin(\theta)+v_\phi\cos(\theta)$

where $\theta$ is the azimuthal angle in the interval

$0\le\theta<2\pi.$

The halo in the disk galaxy is spawned using the same Hernquist distribution as the elliptical galaxy above; since these are dark matter particles, the model doesn't render them in the simulation.

💻 The Engineering Bottleneck: Barnes-Hut & Compute Shaders

The Barnes-Hut Treecode:

The greatest computational bottleneck in an N-body simulation lies in evaluating the gravitational interaction between every pair of particles. Under Newtonian gravity, every particle exerts a force on every other particle, requiring

pairwise force evaluations each simulation step, a computational complexity of $O(N^2)$. For a system containing approximately 110,000 particles, a brute-force implementation would require over 12 billion gravitational interactions every frame, making real-time execution on consumer hardware impossible.

To overcome this limitation, the simulation employs the Barnes–Hut Treecode, originally proposed by Barnes and Hut. Rather than evaluating every interaction individually, the algorithm recursively subdivides the domain into an adaptive octree, where each node stores the total mass and center of mass for that region. During force evaluation, distant clusters can then be approximated as a single massive body instead of computed particle-by-particle.

This approximation is controlled through the Multipole Acceptance Criterion (MAC),

where $s$ is the width of the current octree cell, $d$ is the distance between the particle and the cell's center of mass, and $\theta$ is the opening angle.

If satisfied, the entire node is treated as a single body at its center of mass; otherwise the node is opened and its children examined recursively. Nearby particles retain full accuracy, while distant regions are efficiently approximated, reducing complexity from $O(N^2)$ to approximately $O(N\log N)$ and allowing systems of over one hundred thousand particles to be simulated interactively.

Unlike uniform spatial grids, the octree is adaptive — dense regions near galactic cores subdivide into many small cells, while sparse regions occupy only a few large nodes, reducing memory usage while keeping force accuracy where it matters most.

Octree Data Structure

Each node of the Barnes–Hut tree stores the information required for hierarchical force approximation.

struct alignas(16) OctreeNode
{
    glm::vec4 centerOfMassMass;
    glm::vec4 boundsMin;
    glm::vec4 boundsMax;
 
    int children[8];
 
    int isLeaf;
    int particleIndex;
};

Each node stores its total mass, center of mass, axis-aligned bounding box, pointers to its eight child octants, and whether it's a leaf or internal node.

The structure is explicitly aligned to 16 bytes, ensuring efficient memory transfers between the CPU and GPU using Shader Storage Buffer Objects (SSBOs).

Building the Octree

At the start of every simulation step, the tree is reconstructed using the latest particle positions.

glm::vec3 minB(1e9f), maxB(-1e9f);
 
for (const auto &p : particles)
{
    minB = glm::min(minB, glm::vec3(p.posMass));
    maxB = glm::max(maxB, glm::vec3(p.posMass));
}
 
createNode(minB, maxB);
 
for (int i = 0; i < particles.size(); ++i)
{
    insert(0, i, particles);
}

The simulation first determines the global bounding box enclosing every particle, which becomes the root node, then inserts each particle individually. Whenever two particles occupy the same leaf node, it's subdivided into eight child octants and both particles are recursively redistributed into their corresponding regions.

Recursive Particle Insertion

The insertion routine updates the center of mass during traversal.

glm::vec3 newCoM =(glm::vec3(currentCoM)*currentMass+pPos*pMass)/(currentMass+pMass);
 
nodes[currNode].centerOfMassMass =glm::vec4(newCoM,currentMass+pMass);

Rather than computing centers of mass after the tree is built, the implementation updates both mass and center of mass incrementally during insertion, eliminating an extra traversal over the completed tree.

If two particles occupy identical coordinates, a small numerical tolerance prevents infinite subdivision caused by floating-point precision.

Hybrid CPU-GPU Pipeline

The simulation combines CPU spatial data structures with GPU parallel force computation. Each frame follows the sequence

1. Read updated particle positions back from the GPU.
2. Construct a new Barnes–Hut octree on the CPU.
3. Upload the octree into a Shader Storage Buffer Object.
4. Launch the OpenGL Compute Shader.
5. Compute gravitational accelerations in parallel.
6. Synchronize GPU memory before rendering.

The synchronization stage is implemented as

Particle* ptr =(Particle*)glMapBufferRange(GL_SHADER_STORAGE_BUFFER,0,particles.size()*sizeof(Particle),GL_MAP_READ_BIT);
 
std::copy(ptr,ptr+particles.size(),particles.begin());
 
tree.build(particles);
 
glBufferData(GL_SHADER_STORAGE_BUFFER,tree.nodes.size()*sizeof(OctreeNode),tree.nodes.data(),GL_DYNAMIC_DRAW);

Once the tree has been transferred to the GPU, each compute shader invocation evaluates the gravitational acceleration for a single particle using the hierarchical octree representation.

Compute Shader Dispatch

Finally, the simulation launches one compute thread for every particle.

glUniform1f(
glGetUniformLocation(computeProgram,"theta"),
0.5f);
 
int numGroups =
(NUM_PARTICLES+255)/256;
 
glDispatchCompute(
numGroups,
1,
1);
 
glMemoryBarrier(
GL_SHADER_STORAGE_BARRIER_BIT |
GL_VERTEX_ATTRIB_ARRAY_BARRIER_BIT);

The opening angle $\theta = 0.5$ balances accuracy against speed: smaller values open more nodes for greater accuracy, larger values approximate more aggressively for speed.

⚡ Parallelizing the N-Body Solver with OpenGL Compute Shaders

Although the Barnes–Hut tree reduces complexity from $O(N^2)$ to approximately $O(N\log N)$, evaluating the force acting on every particle remains the most computationally demanding stage of the simulation.

Modern GPUs expose thousands of programmable execution units capable of arbitrary numerical computation, not just rendering. Since each particle's acceleration can be evaluated independently, the N-body problem exhibits a high degree of data parallelism, making it an ideal candidate for GPU execution.

To exploit this parallelism, the force evaluation and time integration stage was implemented entirely with OpenGL 4.3 Compute Shaders. Rather than iterating over particles sequentially on the CPU, thousands of GPU threads execute simultaneously, each computing acceleration, velocity, and position for a single particle.

The simulation stores both particle data and the Barnes–Hut tree inside Shader Storage Buffer Objects (SSBOs), giving arbitrary read/write access directly from GPU memory without transferring intermediate data back to the CPU during force evaluation.

The particle buffer and tree buffer are declared as

layout(std430, binding = 0) buffer ParticleBuffer
{
    Particle particles[];
};
 
layout(std430, binding = 1) buffer TreeBuffer
{
    OctreeNode tree[];
};

Using the std430 memory layout guarantees the GLSL structures match their C++ counterparts' memory layout, so both processors access exactly the same data representation.

Massive Parallel Execution

Each invocation of the compute shader is responsible for evolving exactly one particle.

layout(local_size_x = 256) in;
 
void main()
{
    uint index = gl_GlobalInvocationID.x;
 
    if(index >= numParticles)
        return;
 
    vec3 myPos = particles[index].posMass.xyz;

This creates workgroups containing 256 parallel threads. When dispatched, the GPU automatically launches enough workgroups to process the complete particle system, calculated on the CPU as

int numGroups =(NUM_PARTICLES + 255) / 256;
 
glDispatchCompute(numGroups,1,1);

This guarantees that every particle is assigned exactly one compute shader invocation, even when the particle count is not an exact multiple of the workgroup size.

Hierarchical Force Evaluation

Rather than iterating over every particle in the system, each shader invocation traverses the Barnes–Hut octree stored inside GPU memory.

int stack[128];
int stackPtr = 0;
 
stack[stackPtr++] = 0;
 
while(stackPtr > 0)
{
    int nodeIdx = stack[--stackPtr];
    OctreeNode node = tree[nodeIdx];

Instead of recursion, the traversal uses an explicit stack allocated in local shader memory, avoiding recursive calls, which are unsupported or inefficient on many GPU architectures.

Whenever an internal node is encountered, the shader evaluates the Barnes–Hut Multipole Acceptance Criterion

float width =max(max(boundsSize.x,boundsSize.y),boundsSize.z);
 
if((width/dist) < theta)
{
    accel_nbody +=(G*mass/(distSq*dist))*distVec;
}

If the node is sufficiently distant, the subtree is approximated using its center of mass; otherwise its children are pushed onto the traversal stack and evaluated individually.

Time Integration on the GPU

After the total gravitational acceleration has been accumulated, the particle state is immediately updated inside the same compute shader.

particles[index].velocity.xyz += accel * dt;
particles[index].posMass.xyz += particles[index].velocity.xyz * dt;

The simulation employs a Semi-Implicit Euler Integrator, updating velocities with the newly computed acceleration before advancing positions — significantly better numerical stability than explicit Euler for long-running simulations, at very low computational cost.

Softening Length & Symplectic Time Integration:

Although stars in a galaxy interact through Newtonian gravity, a real galaxy is fundamentally a collisionless system rather than a collection of particles undergoing frequent close encounters — each simulation particle represents a large ensemble of stars, not an individual star. The objective is therefore to reproduce the smooth, large-scale gravitational potential governing galactic evolution, rather than resolve every microscopic stellar interaction.

A direct implementation of Newton's inverse-square law introduces two major numerical challenges. First, when two particles pass extremely close to one another, the gravitational force increases rapidly as

$F \propto \frac{1}{r^2},$

causing unrealistically large accelerations that can dominate the timestep and destabilize the entire simulation. These close encounters are numerical artifacts arising from the finite number of particles used to approximate a continuous stellar distribution and do not accurately represent the collisionless dynamics of real galaxies.

Secondly, advancing particle trajectories over thousands of timesteps requires a stable integration scheme — conventional explicit methods gradually accumulate truncation errors, causing artificial gains or losses of total energy that manifest as unrealistic orbital decay, expanding trajectories, or distorted galactic structures.

To address these issues, the simulator incorporates two complementary numerical techniques: Plummer softening, which regularizes the gravitational force during close encounters, and a Semi-Implicit (Symplectic) Euler Integrator, which provides significantly improved long-term energy behaviour while maintaining a computationally efficient integration scheme suitable for real-time execution.

Direct evaluation of the Newtonian gravitational force,

$\mathbf{F}=G\frac{m_1m_2}{r^2},$

becomes numerically unstable as the inter-particle separation approaches zero, diverging as $r \rightarrow 0$ and forcing tiny timesteps with unrealistically large accelerations due to finite-precision arithmetic. To regularize this singularity, the simulation employs Plummer Softening, replacing the squared separation with

$r^2+\epsilon^2,$

where $\epsilon$ is the softening length. Consequently, the gravitational acceleration is evaluated as

$\mathbf{a}=Gm\frac{\mathbf{r}}{\left(r^2+\epsilon^2\right)^{3/2}},$

which removes the divergence at short distances while preserving inverse-square behaviour at larger scales, and suppresses artificial two-body scattering so the particle system more accurately represents a collisionless stellar distribution.

The softening is implemented directly inside the compute shader during force evaluation.

const float softeningSq = 0.05f * 0.05f;
 
vec3 distVec = com - myPos;
float distSq = dot(distVec, distVec) + softeningSq;
float dist   = sqrt(distSq);
 
accel_nbody +=(G * mass / (distSq * dist))* distVec;

Once the total gravitational acceleration has been accumulated through the Barnes-Hut traversal, particle trajectories are advanced using a Semi-Implicit Euler (Symplectic Euler) Integrator. Unlike standard explicit Euler, which exhibits continuous numerical energy drift in Hamiltonian systems, the symplectic formulation preserves the geometric structure of phase space considerably better, giving significantly improved long-term orbital stability.

The integration scheme first updates the particle velocity using the newly computed acceleration,

$\mathbf{v}_{t+\Delta t}=\mathbf{v}_t+\mathbf{a}\Delta t,$

followed by advancing the particle position using the updated velocity,

$\mathbf{x}_{t+\Delta t}=\mathbf{x}_t+\mathbf{v}_{t+\Delta t}\Delta t.$

This ordering yields substantially better energy conservation over long runs while keeping the simplicity and efficiency needed for real-time simulation, executed entirely on the GPU:

particles[index].velocity.xyz += accel * dt;
particles[index].posMass.xyz  += particles[index].velocity.xyz * dt;

By combining Plummer softening with a symplectic integration scheme, the simulator remains numerically stable during close particle encounters while preserving physically plausible orbital evolution over thousands of simulation timesteps.

⚙️ Code Architecture Pipeline

The C++ and GLSL codebase is approximately 1,200 lines, structured specifically to minimize CPU-GPU memory bottlenecks. Here is the operational pipeline:

1. GLFW & Application Initialization

GLFW handles the creation of the OpenGL 4.3 core profile context, window management, and input polling. GLAD resolves all OpenGL function pointers, and camera controls (mouse look and WASD movement) are set up. It serves as the application's backbone, capturing user inputs and calculating delta time for frame-rate-independent physics integration, and also seeds the random number generators and invokes the inverse transform sampling that generates the initial position and velocity vectors for the 110,000 particles.

2. Shader Storage Buffer Objects (SSBOs)

To avoid the performance penalty of sending 110,000 particles from CPU to GPU every frame, data is stored in SSBOs. The C++ particle array is uploaded once at initialization; every frame, the CPU maps the particle buffer to read positions, builds the Barnes-Hut Octree locally, and uploads the newly constructed tree to a second SSBO. This gives the Compute Shader instantaneous, high-bandwidth access to both particle states and spatial acceleration structures without stalling the graphics pipeline.

3. Compute Shaders (The Physics Engine)

The Compute Shader is the heavy lifter of this architecture. Dispatched in localized workgroups (256 threads per group), it assigns exactly one GPU thread to one particle. Each thread reads the Octree SSBO and executes the explicit stack-based Barnes-Hut traversal, accumulating gravitational acceleration, then applies the semi-implicit Euler integration, updating the particle's velocity and position directly in the SSBO. A memory barrier ensures all physics calculations finish before rendering.

4. Vertex Shaders

The Vertex Shader operates directly on the same SSBO updated by the Compute Shader, so zero particle data is sent back to the CPU for rendering. It reads the updated vec4 posMass attributes and applies the Model-View-Projection (MVP) matrix transformation to map 3D coordinates to screen space, and uses gl_VertexID to assign a color per galactic population: cool blue/white for the spiral disk, reddish-yellow for the elliptical galaxy, and full transparency (discarding) for the dark matter halo.

5. Fragment Shaders

The Fragment Shader turns the mathematical points into radiating stars. Instead of harsh square pixels, it uses gl_PointCoord to measure the fragment's distance from the center of the point sprite, discarding anything past the radius to form a perfect circle. Remaining fragments get a radial intensity gradient — fading alpha towards the edges — and additive blending accumulates light where particles overlap, producing the glowing dense cores of the galaxies.

🔬 Research Significance

Simulators of this nature are not just visual toys; galaxy mergers are fundamental to modern astrophysics and cosmological structure formation. This tool specifically replicates the dynamics of systems like the Medusa galaxy, allowing us to observe the non-linear growth of tidal structures over millions of years.

Because every initial condition is exposed to the user, the simulation functions as a lightweight numerical experiment platform. Varying the mass ratio, disk geometry, or approach velocity lets researchers explore how each factor shapes the resulting tidal tail, and whether the encounter ends as a clean fly-by or a bound merger. Beyond astrophysics, implementing Barnes-Hut treecodes inside GPU compute shaders is a pattern directly transferable to resolving $O(N^2)$ bottlenecks in plasma physics and molecular dynamics.

🎛️ Interactive UI & Conclusion

Galaxy mergers exhibit an extraordinary sensitivity to their initial conditions. Small variations in the mass ratio, impact parameter, relative velocity, or disk geometry can produce entirely different tidal tails, bridges, and merger remnants. Consequently, recompiling the application after every parameter adjustment becomes an inefficient workflow during experimentation.

To address this, the simulator integrates Dear ImGui directly into the OpenGL rendering pipeline, providing an interactive graphical interface for configuring the physical parameters of both galaxies before the simulation begins. Rather than modifying hard-coded constants within the source code, users can dynamically alter the structural properties of the galaxies and immediately launch a new numerical experiment.

The interface exposes several physically meaningful parameters, including

• Spiral disk radius
• Spiral disk thickness
• Total disk mass
• Elliptical galaxy mass
• Side-to-side collision offset (impact parameter)
• Initial approach velocity

Additionally, several predefined collision configurations are available, allowing common merger scenarios such as head-on collisions and glancing encounters to be reproduced with a single button press.

The Dear ImGui context is initialized directly on top of the existing GLFW/OpenGL rendering context.

IMGUI_CHECKVERSION();
 
ImGui::CreateContext();
 
ImGui::StyleColorsDark();
 
ImGui_ImplGlfw_InitForOpenGL(window, true);
 
ImGui_ImplOpenGL3_Init("#version 430");

Once initialized, the setup window is constructed entirely using Dear ImGui widgets, allowing the simulation parameters to be modified interactively before particle initialization.

ImGui::SliderFloat("Size (Radius)", &a, 0.3f, 3.0f);
 
ImGui::SliderFloat("Thickness", &b, 0.02f, 0.5f);
 
ImGui::SliderFloat("Total Mass", &M_disk, 1.0f, 20.0f);
 
ImGui::SliderFloat("Elliptical Mass", &M_ell, 50.0f, 800.0f);
 
ImGui::SliderFloat("Side-to-Side Offset", &v_x_special, -3.0f, 3.0f);
 
ImGui::SliderFloat("Approach Speed", &v_y_special, 1.0f, 5.0f);

To further simplify experimentation, predefined presets are provided for commonly studied collision geometries.

if (ImGui::Button("Preset: Head-On Collision"))
{
    v_x_special = 3.0f;
    v_y_special = 2.0f;
}
 
if (ImGui::Button("Preset: Glancing Blow"))
{
    v_x_special = 2.0f;
    v_y_special = 2.0f;
}

Once the desired configuration has been selected, the simulation is initialized only after the Start Simulation button is pressed, ensuring that all particle distributions are generated using the updated physical parameters.

if (ImGui::Button("Start Simulation", ImVec2(120,0)))
{
    simulationStarted = true;
}

Unlike conventional N-body simulators that rely on configuration files or source-code modifications, this interface enables rapid exploration of parameter space, making it easier to investigate how different collision geometries influence the resulting galactic morphology.

Free-Fly Camera System

Since galaxy mergers evolve in three spatial dimensions, a single fixed viewpoint gives only limited physical insight. To enable detailed inspection of tidal bridges, stellar streams, and disk deformation, a fully interactive first-person camera system was implemented using GLFW input callbacks.

The camera orientation is controlled through yaw and pitch Euler angles, while mouse movement continuously updates the viewing direction.

direction.x =cos(glm::radians(yaw))*cos(glm::radians(pitch));
 
direction.y =sin(glm::radians(pitch));
 
direction.z =sin(glm::radians(yaw))*cos(glm::radians(pitch));
 
cameraFront =glm::normalize(direction);

Zooming is implemented by dynamically modifying the camera's field of view using the mouse scroll wheel.

void scroll_callback(...)
{
    fov -= (float)yoffset;
 
    if(fov < 1.0f)
        fov = 1.0f;
 
    if(fov > 45.0f)
        fov = 45.0f;
}

Real-time navigation throughout the simulation domain is performed using the standard WASD control scheme.

if(glfwGetKey(window, GLFW_KEY_W) == GLFW_PRESS)
    cameraPos += cameraSpeed * cameraFront;
 
if(glfwGetKey(window, GLFW_KEY_S) == GLFW_PRESS)
    cameraPos -= cameraSpeed * cameraFront;
 
if(glfwGetKey(window, GLFW_KEY_A) == GLFW_PRESS)
    cameraPos -= glm::normalize(glm::cross(cameraFront,cameraUp))* cameraSpeed;
 
if(glfwGetKey(window, GLFW_KEY_D) == GLFW_PRESS)
    cameraPos += glm::normalize(glm::cross(cameraFront,cameraUp)) * cameraSpeed;

The resulting free-fly camera allows the evolving merger to be inspected from arbitrary viewpoints, enabling close examination of dynamically evolving tidal tails, stellar bridges, disk warping, and the eventual relaxation of the merged galaxy.

Conclusion

Developing this simulator required combining concepts from computational astrophysics, numerical methods, modern OpenGL, and GPU programming into a unified real-time application. Stable initial conditions were generated using physically motivated galaxy models, gravitational interactions were accelerated through the Barnes-Hut Treecode, and the computational workload was further parallelized using OpenGL Compute Shaders executing thousands of particle updates simultaneously.

Beyond the numerical solver itself, the interactive Dear ImGui interface and free-fly visualization system transformed the simulator into an experimental platform rather than a static demonstration, letting users immediately observe how changes in galaxy mass, disk structure, impact parameter, and encounter velocity influence the resulting merger dynamics.

By combining physically motivated galaxy models with hardware-accelerated numerical algorithms, the simulator demonstrates that research-inspired N-body simulations involving over 110,000 particles can be performed interactively on commodity GPU hardware while maintaining both numerical stability and real-time visualization.