A high-performance N-body gravitational simulator written in C++, featuring a pluggable gravity model system, a symplectic leapfrog integrator, and OpenMP parallelisation.
During my Master's thesis at Ludwig-Maximilians-Universität München (2025), I investigated whether late-time modifications to gravity could resolve the Hubble tension — the ~5σ discrepancy between early- and late-universe measurements of the Hubble constant. That work was analytical and statistical, using Bayesian inference on cosmological datasets.
This simulator was built in March 2026 as a personal C++ project, inspired by that research. It brings the same class of gravity modifications down to the N-body scale: rather than constraining models statistically, here their orbital effects can be observed directly in simulation.
- Leapfrog (symplectic) integrator — time-reversible and energy-conserving; energy drift < 1e-12% over one simulated year
- Pluggable gravity models via abstract base class — swap physics with a
single line in
main.cpp - Yukawa-screened modified gravity — parameterised by strength
αand screening scaleλ, directly motivated by modified gravity research - OpenMP parallelisation of the O(N²) force loop — tested across 10 threads on Apple M-series hardware
- CSV trajectory output with Python/matplotlib visualisation
The Yukawa correction adds an extra attractive force at solar-system scales, pulling Earth into a tighter orbit. Both simulations conserve energy well.
| Model | Energy drift (1 year) |
|---|---|
| Newtonian gravity | 9.58 × 10⁻¹³ % |
| Modified gravity (α=0.1, λ=1 AU) | 0.014 % |
.
├── Vec3.h # 3D vector type with operator overloading
├── Particle.h # Particle class (pos, vel, force, mass)
├── GravityModel.h # Abstract base class for gravity models
├── NewtonianGravity.h # Standard 1/r² gravity
├── ModifiedGravity.h # Yukawa-screened modified gravity
├── Simulation.h # Leapfrog integrator + OpenMP force loop
├── FileWriter.h # CSV trajectory output
├── main.cpp # Entry point — configure and run simulations
├── plot.py # Python visualisation (matplotlib)
├── Makefile # Build system
└── images/ # Output plots
- C++17 compiler (g++ or Apple Clang)
- Homebrew (macOS only, for OpenMP)
- Python 3 with
pandasandmatplotlib
brew install libomp
make
./maing++ -O2 -std=c++17 -fopenmp main.cpp -o main
./mainpython3 plot.pyStandard inverse-square law with softening to prevent force singularities at close range:
F = G · m₁ · m₂ / (r² + ε²)
Adds a scale-dependent correction motivated by scalar-tensor theories of gravity:
F = F_Newton · (1 + α · exp(−r / λ))
| Parameter | Meaning | Default |
|---|---|---|
α (alpha) |
Modification strength | 0.1 |
λ (lambda) |
Screening length (metres) | 1.496e11 (1 AU) |
Setting α = 0 recovers standard Newtonian gravity exactly.
Subclass GravityModel and implement one method:
class MyGravity : public GravityModel {
public:
Vec3 force(const Vec3& r, double m1, double m2) const override {
// your implementation
}
};Pass it to Simulation in main.cpp:
Simulation sim(std::make_shared<MyGravity>());The Hubble constant H₀ — the rate at which the universe is expanding — can be measured two ways: from the early universe (CMB, giving ~67 km/s/Mpc) and from the late universe (distance ladder, giving ~73 km/s/Mpc). The ~5σ discrepancy between these values is one of the most significant open problems in cosmology.
One proposed resolution is that gravity itself behaves differently at late cosmological times. My thesis used the MGCAMB package and Bayesian inference on cosmological datasets to constrain such models. This simulator explores the same class of modifications at the N-body scale.
Ege Özmeral — MSc Physics, Ludwig-Maximilians-Universität München
linkedin.com/in/ege-ozmeral · github.com/wittyphysicist