numerics
Loading...
Searching...
No Matches
nbody.hpp
Go to the documentation of this file.
1/// @file nbody.hpp
2/// @brief Gravitational N-body simulation using num::ode_verlet / num::ode_rk4.
3///
4/// NBodySim wraps the ODE module to provide a per-frame step() interface.
5/// Each call to step(dt) runs exactly one time step of the chosen integrator:
6///
7/// Verlet (symplectic, default) -- ode_verlet: O(h²), bounded energy error
8/// RK4 (non-symplectic) -- ode_rk4: O(h⁴), secular energy drift
9///
10/// This lets the visualisation directly compare the two on the same scenario.
11#pragma once
12
13#include "ode/ode.hpp"
14#include <vector>
15#include <utility>
16#include <cstdint>
17#include <cmath>
18#include <string>
19
20namespace nbody {
21
22// ── Body descriptor ───────────────────────────────────────────────────────────
23
24struct Body {
25 double mass;
26 float display_radius;
27 uint32_t color; ///< RGBA packed (raylib Color layout)
28 float phys_radius = 0.0f; ///< physics-unit radius for merge detection (0 = disabled)
29};
30
31// ── Scenarios ─────────────────────────────────────────────────────────────────
32
33enum class Scenario {
34 Figure8, ///< Chenciner-Montgomery figure-8 choreography (3 equal masses)
35 SolarSystem, ///< Sun + 4 planets on circular Keplerian orbits
36 BinaryPlus, ///< Equal-mass binary + one test particle on a wide orbit
37 Galaxy, ///< Random N bodies — gravity only, bodies merge on contact
38};
39
40inline const char* scenario_name(Scenario s) {
41 switch (s) {
42 case Scenario::Figure8: return "Figure-8 Choreography";
43 case Scenario::SolarSystem: return "Keplerian Orbits";
44 case Scenario::BinaryPlus: return "Binary + Test Particle";
45 case Scenario::Galaxy: return "Galaxy Collapse";
46 }
47 return "?";
48}
49
50// ── Simulation ────────────────────────────────────────────────────────────────
51
52struct NBodySim {
53 std::vector<Body> bodies;
54 num::Vector q; ///< Positions: [x0,y0, x1,y1, ...]
55 num::Vector v; ///< Velocities: [vx0,vy0, vx1,vy1, ...]
56 double t = 0.0;
57 double G = 1.0;
58 double E0 = 0.0; ///< Initial total energy (for drift tracking)
59 double soft = 1e-3; ///< Softening length (avoids singularities)
60
61 bool use_verlet = true; ///< true = Verlet (symplectic), false = RK4
62 bool enable_merges = false; ///< true = check_merges() is active (Galaxy scenario)
63
64 /// Reset to a preset scenario.
65 void reset(Scenario s);
66
67 /// Advance by exactly one time step dt.
68 void step(double dt);
69
70 /// Check all pairs for contact; merge overlapping bodies.
71 /// Returns the list of (removed_idx, swapped_from_idx) operations in order
72 /// so the caller can mirror them on any parallel index array (e.g. trails).
73 std::vector<std::pair<int,int>> check_merges();
74
75 int n() const { return static_cast<int>(bodies.size()); }
76
77 double kinetic_energy() const;
78 double potential_energy() const;
79 double total_energy() const { return kinetic_energy() + potential_energy(); }
80
81 /// Relative energy drift since reset: (E - E0) / |E0|
82 double energy_drift() const {
83 return E0 != 0.0 ? (total_energy() - E0) / std::abs(E0) : 0.0;
84 }
85
86 /// Fill acc from current positions (gravitational acceleration on each body).
87 void make_accel(const num::Vector& pos, num::Vector& acc) const;
88};
89
90} // namespace nbody
nbody
Definition main.cpp:27
nbody::scenario_name
const char * scenario_name(Scenario s)
Definition nbody.hpp:40
nbody::Scenario
Scenario
Definition nbody.hpp:33
nbody::Scenario::BinaryPlus
@ BinaryPlus
Equal-mass binary + one test particle on a wide orbit.
nbody::Scenario::Galaxy
@ Galaxy
Random N bodies — gravity only, bodies merge on contact.
nbody::Scenario::Figure8
@ Figure8
Chenciner-Montgomery figure-8 choreography (3 equal masses)
nbody::Scenario::SolarSystem
@ SolarSystem
Sun + 4 planets on circular Keplerian orbits.
ode.hpp
ODE time integrators: Euler, RK4, adaptive RK45 (Dormand-Prince), and symplectic Velocity Verlet / Yo...
nbody::Body::phys_radius
float phys_radius
physics-unit radius for merge detection (0 = disabled)
Definition nbody.hpp:28
nbody::Body::mass
double mass
Definition nbody.hpp:25
nbody::Body::display_radius
float display_radius
Definition nbody.hpp:26
nbody::Body::color
uint32_t color
RGBA packed (raylib Color layout)
Definition nbody.hpp:27
nbody::NBodySim::kinetic_energy
double kinetic_energy() const
Definition main.cpp:70
nbody::NBodySim::potential_energy
double potential_energy() const
Definition main.cpp:77
nbody::NBodySim::step
void step(double dt)
Advance by exactly one time step dt.
Definition main.cpp:89
nbody::NBodySim::q
num::Vector q
Positions: [x0,y0, x1,y1, ...].
Definition nbody.hpp:54
nbody::NBodySim::n
int n() const
Definition nbody.hpp:75
nbody::NBodySim::use_verlet
bool use_verlet
true = Verlet (symplectic), false = RK4
Definition nbody.hpp:61
nbody::NBodySim::make_accel
void make_accel(const num::Vector &pos, num::Vector &acc) const
Fill acc from current positions (gravitational acceleration on each body).
Definition main.cpp:53
nbody::NBodySim::reset
void reset(Scenario s)
Reset to a preset scenario.
Definition main.cpp:165
nbody::NBodySim::total_energy
double total_energy() const
Definition nbody.hpp:79
nbody::NBodySim::v
num::Vector v
Velocities: [vx0,vy0, vx1,vy1, ...].
Definition nbody.hpp:55
nbody::NBodySim::check_merges
std::vector< std::pair< int, int > > check_merges()
Definition main.cpp:118
nbody::NBodySim::energy_drift
double energy_drift() const
Relative energy drift since reset: (E - E0) / |E0|.
Definition nbody.hpp:82
nbody::NBodySim::bodies
std::vector< Body > bodies
Definition nbody.hpp:53
nbody::NBodySim::E0
double E0
Initial total energy (for drift tracking)
Definition nbody.hpp:58
nbody::NBodySim::soft
double soft
Softening length (avoids singularities)
Definition nbody.hpp:59
nbody::NBodySim::enable_merges
bool enable_merges
true = check_merges() is active (Galaxy scenario)
Definition nbody.hpp:62