N-Body System Simulation

Problem Description

The simulation of N-body systems has historically been one of the most important and fascinating problems in astronomy. The problem consists of simulating the dynamics of a set of N bodies that interact gravitationally with each other.

According to the laws formulated by Newton in 1687, the interaction between two bodies of mass $m_1$ and $m_2$ is governed by an attractive force proportional to their masses and inversely proportional to the square of their separation:

\[F = G\frac{m_1m_2}{r^2}, \quad r = |\vec{r}_1-\vec{r}_2|\]

where $G\approx 6.67\times 10^{-11}\ Nm^2kg^{-2}$ is the universal gravitational constant.

Despite its apparent simplicity, an analytical solution has been found only for the two-body case. For three or more bodies, decades of study ultimately showed that no general solution exists. This underscores the need for numerical simulations to predict planetary orbits, model galaxy formation, and tackle many other problems in astronomy.

Numerical Implementation

For this project we restrict ourselves to two dimensions. Consider $N$ bodies with masses $m_i$. We denote their positions, velocities, and accelerations by $\vec{r}_i$, $\vec{v}_i$, and $\vec{a}_i$, respectively. From Newton’s law of gravitation, the net force on body $i$ at any instant is the sum of the forces exerted by all other bodies:

\[\begin{equation*} \vec F_i = -\sum_{\substack{j=1 \\ j \neq i}}^{N} G \frac{m_i m_j}{|\vec r_i - \vec r_j|^3} (\vec r_i - \vec r_j), \quad i = 1, \dots, N \end{equation*}\]

where we use the vector form of Newton’s law.

Since these equations generally have no analytical solution, they must be discretized and solved numerically, step by step. We use an algorithm called Velocity Verlet, which belongs to a family of methods known as symplectic integrators. These numerical methods preserve key physical quantities—such as energy and angular momentum—approximately over long simulations, mirroring the conservation laws expected from any mechanical system.

The equations you need to implement for a single simulation step are:

\[\begin{align*} \vec r_i(t+\Delta t) &= \vec r_i(t) + \vec v_i(t) \Delta t + \frac{1}{2} \vec a_i(t) \Delta t^2, \\ \vec a_i(t+\Delta t) &= -\sum_{\substack{j=1 \\ j \neq i}}^N G \frac{m_j}{\big(|\vec r_j(t+\Delta t)-\vec r_i(t+\Delta t)|^2 + \epsilon^2\big)^{3/2}} (\vec r_i(t+\Delta t)-\vec r_j(t+\Delta t)), \\ \vec v_i(t+\Delta t) &= \vec v_i(t) + \frac{1}{2} \Big(\vec a_i(t) + \vec a_i(t+\Delta t)\Big) \Delta t \end{align*}\]

Here $\Delta t$ is the time step. Smaller values yield more accurate dynamics but cost more computation; a value on the order of $0.001$–$0.01$ is typically acceptable.

We have also introduced a softening parameter $\epsilon = 10^{-12}$, which prevents the force from diverging when two bodies occupy the same position.

Note 1: The acceleration of body $i$ is computed directly, so its mass $m_i$ does not appear in the force sum.

Note 2: Execute the three update steps in the order shown above. Changing the order will break the dynamics.

Energy Conservation

Energy is one of the most important quantities in any mechanical system. For an isolated system, Newton’s laws guarantee that total energy remains constant.

In an N-body system, the kinetic energy is the sum of the individual contributions from each body:

\[K = \sum_{i=1}^N \frac{1}{2}m_i v_i^2\]

The potential energy is obtained by summing over every pair of bodies:

\[U = - \sum_{i < j} G \frac{m_i m_j}{\left| \vec{r}_i - \vec{r}_j \right|}\]

noting that we sum over $i<j$ to avoid counting pairs twice.

Finally, the total mechanical energy is simply the sum of these two terms: $E=K+U$.

Because Velocity Verlet is a symplectic integrator, energy is conserved approximately: you should see the total energy oscillate around its initial value without any sustained drift upward or downward.

In addition to energy, isolated mechanical systems conserve two other quantities: linear momentum and angular momentum.

Project Setup

Create a subdirectory inside your existing exercises repository
cd ~/src/exercises
mkdir nbody

Register the new subdirectory with CMake

Append add_subdirectory(nbody) to the top-level CMakeLists.txt:
add_subdirectory(nbody)

Add nbody/CMakeLists.txt for the nbody exercise
add_runnable(nbody main.cpp)

Add a stub main.cpp

int main() {
  return 0;
}

Step-by-step Implementation

Your task is to write a program that simulates the dynamics of an N-body system. The program should accept the number of bodies $N$, their masses, and optionally the number of simulation steps.

Assign initial positions and velocities either randomly (choosing reasonable ranges) or explicitly. If you choose explicit initialization, do it in code. (We did not cover file I/O.)

At each step, compute the total energy of the system to verify that it remains approximately constant.

Finally, output results either to the console or through a graphical visualization (see the SFML hints below).

Define `Vec2d<T>`

Create nbody/vec2d.h

Define template <typename T> struct Vec2d { T x, y; };.

Add non-member friend operators / functions (hidden friends) as needed. Start with adding operator+. See User-defined literals and operators.

Add friend constexpr bool operator==(Vec2d, Vec2d) = default; inside the struct. See Defaulted comparisons.

Mark all functions constexpr. See constexpr specifier.

Write constexpr unit tests

Create nbody/vec2d_test.cpp with the following testing harness:
#include "vec2d.h"

#define VERIFY(expr)                                 \
  do {                                               \
    if (not (expr)) throw "failure at: `" #expr "`"; \
  } while (false)

consteval {
  Vec2d<float> a = {};
  VERIFY(a.x == 0);
  VERIFY(a.y == 0);

  // add more ...
}
Register the test file in nbody/CMakeLists.txt:
add_runnable(nbody main.cpp)

add_library(tests STATIC vec2d_test.cpp)
The mere act of compiling vec2d_test.cpp is enough to exercise the consteval block. If any VERIFY fires, compilation fails with a message naming the failing expression. Extend the test suite as you add each new operation to Vec2d.

Clang and C++26

Clang doesn’t understand the consteval block yet. A simple alternative is to use:
static_assert([] {
  // add your consteval code here
  return true;
}());

Design notes

constexpr testing: Wrapping tests in consteval { ... } runs them at compile time. We could place the tests in the header and thus force a check on every use of the header. However, compile-time costs are real. Have a separate build target for this can be more efficient in the long run.

The VERIFY macro throws a descriptive string on failure, which the compiler prints as part of the diagnostics.

Since invoking undefined behavior inside a constant expression is ill-formed, this approach catches UB early.

add_library(tests STATIC ...): Adding test files as a static library ensures they are compiled and linked into the build without producing a separate executable. Any uncaught constexpr exceptions or static_assert failures surface immediately during compilation.

Hidden friend functions: Keeping operators and utility functions outside the class maintains a minimal interface and follows modern C++ best practices. They also enable argument-dependent lookup (ADL), so once you add a hidden friend called distance, distance(a, b) just works without explicit namespace qualification.

The VERIFY macro

In general, our goal is to avoid macros. However, in tests, macros are much less problematic.

With GCC trunk (to be GCC 17), we should be able to provide a verify_equal function that uses constexpr std::format to create a helpful error string on failure.

Commit to git

Define `Body<T>`

Create nbody/body.h

Include "vec2d.h".

Define template <typename T> struct Body { Vec2d<T> r, v, a; T m; };.

Write void print_state(std::span<const Body<T>> bodies) using structured bindings to iterate and print each body.

In main.cpp, introduce a type alias using T = float;. Initialize three bodies with the figure-eight initial conditions given below, and call print_state<T>(bodies).
r_1 = {-0.97000436, 0.24308753};
r_2 = {0.97000436, -0.24308753};
r_3 = {0, 0};
v_1 = {0.4662036850, 0.4323657300};
v_2 = {0.4662036850, 0.4323657300};
v_3 = {-0.93240737, -0.86473146};

Type aliases and std::span

The using T = float; alias means every function call that requires an explicit template argument uses T rather than repeating float everywhere. When you later experiment with double or SIMD types, you change exactly one line.

Using std::span instead of const std::vector& decouples your algorithms from the container.

Commit to git

Compute gravitational forces

Implement the force computation

Refer back to the Numerical Implementation section for the formula.

Add hidden friends to Vec2d as needed.

Implement void compute_accelerations(std::span<Body<T>> bodies) that computes the pairwise gravitational acceleration for every body. Store the result in each body’s a field.

Choose a value for $G$. Since we are working in simulation units, pick $G = 1$ for simplicity.

Call compute_accelerations<T>(bodies) from main() and print the resulting accelerations to verify.

Commit to git

Implement Velocity Verlet integration

Implement a single simulation step

Refer back to the Numerical Implementation section for the full set of equations. For a single time step with step size $\Delta t$, update in this exact order:

Update positions using current velocities and accelerations.

Recompute accelerations at the new positions.

Update velocities using the average of old and new accelerations.

Split the work into three focused functions:

void update_positions(T dt, std::span<Body<T>> bodies) (only touches r fields).

Refactor compute_accelerations into void update_a_and_v(T dt, std::span<Body<T>> bodies). This function recomputes accelerations at the new positions and updates velocities using the average of old and new accelerations. Can you split the function into an update for a and another function for v? Consider the trade-offs and document your choice. The goal here is readability and reduction of cognitive load.

void do_timestep(T dt, std::span<Body<T>> bodies) (composes the two calls above).

Run a loop of several hundred steps and print the final positions. Verify that bodies move and interact.

Why split the integrator?

Separating position updates from force computation makes each piece independently testable. It also implicitly documents the code, via the names of the functions. This is an important aspect of writing self-documenting code.

Would you tend to avoid the function calls for performance? The compiler agrees with you and will inline function calls, so we don’t have to. We focus on creating readable code first, before we restructure our code for optimization. And we’ll only do the latter when we have data that supports such restructuring.

For an intuitive introduction to Velocity Verlet, see this video.

Commit to git

Readability checkpoint: comparing implementations

Readability checkpoint

Consider two approaches for the pairwise force computation loop:
// Approach A: nested indexed loops
for (std::size_t i = 0; i < bodies.size(); ++i)
  for (std::size_t j = 0; j < bodies.size(); ++j)
    if (i != j)
    {
      // ... compute and accumulate force on bodies[i].a
    }

// Approach B: range-based for
for (auto& body_i : bodies)
  for (const auto& body_j : bodies)
    if (&body_i != &body_j)
    {
      // ... compute and accumulate force on body_i.a
    }
Discuss in pairs (5 min):

Which version makes it clearer that we skip self-interaction?

Which would you rather debug at 2 AM when energy drifts unexpectedly?

Which one gives the compiler better semantic information about what you’re trying to do?

What is the mental load of each? Count: how many variables do you need to track mentally while reading the inner loop body?

There isn’t always a single correct answer — both produce working code. The skill is learning to articulate the trade-offs.

A third way: structured bindings

Both approaches above carry friction. Approach A forces you to track indices, bounds, and field names simultaneously. Approach B hides the index, but still leaves .r, .m, .a scattered through the loop body—and the self-interaction guard (&body_i != &body_j) reads like pointer arithmetic rather than physics.

C++17 introduced structured bindings, which let you unpack a struct into named variables. Combined with the _ discard pattern, they give you the control of indices with the clarity of semantics:
for (std::size_t i = 0; i < bodies.size(); ++i)
{
  const auto& [ri, _, ai, mi] = bodies[i];
  Vec2d<T> a = {};
  for (std::size_t j = 0; j < bodies.size(); ++j)
  {
    if (i == j) continue;
    const auto& [rj, _, _, mj] = bodies[j];
    
    auto dr = ri - rj;
    auto dist2 = dr.x * dr.x + dr.y * dr.y;
    // ... a -= mj / pow(dist2 + epsilon2, 1.5f) * dr;
  }
  bodies[i].v += (ai + a) * (dt / 2);
  bodies[i].a = a;
}
Notice what changes:

You keep the index because you need it for the i == j guard.

You unpack bodies[i] once, giving yourself ri, ai, mi. The _ explicitly discards velocity where it isn’t needed yet.

The force computation reads like the physics equation, not array subscripting.

Field access repetition vanishes. Cognitive load drops significantly.

This is the kind of abstraction that can reduce mental load without costing performance. The compiler sees exactly the same memory accesses; you just gave them meaningful names.

Commit to git (if needed)

Energy conservation check

Compute and monitor total energy

Refer back to the Energy Conservation section for the formulas.

You will need a dot product in Vec2d for computing $v_i^2$.

Implement constexpr T kinetic_energy(std::span<const Body<T>> bodies).

Implement constexpr T potential_energy(std::span<const Body<T>> bodies). Remember to sum over $i < j$ to avoid counting pairs twice.

Print the total energy $E = K + U$ at regular intervals during the simulation (e.g., every 100 steps).

Verify that $E$ oscillates around its initial value without sustained drift.

Why this matters

Velocity Verlet is a symplectic integrator, meaning it preserves the geometric structure of Hamiltonian mechanics. You should see energy oscillate around its initial value without any sustained upward or downward drift. Significant drift indicates a bug in the integration order or a $\Delta t$ that is too large.

Commit to git

Compile-time simulation tests

Refactor VERIFY into a shared header

Move the VERIFY macro from vec2d_test.cpp into its own nbody/verify.h. Include it from both test files instead of duplicating the definition.

Write constexpr simulation tests

Create nbody/nbody_simulation_test.cpp and add it to the tests static library in CMakeLists.txt:
add_library(tests STATIC vec2d_test.cpp nbody_simulation_test.cpp)
Again, use consteval {} blocks. Each block is a self-contained test:
#include "nbody_simulation.h"
#include "verify.h"

consteval
{
  // one test body here
}
Design each test around a specific physical property or edge case. Useful categories include:

Empty system. Zero bodies should yield zero kinetic and potential energy, and stepping should be a no-op.

Single body. A lone body with no forces should translate linearly: position shifts by velocity times $\Delta t$, kinetic energy stays constant, potential energy stays zero.

Two-body attraction. Two bodies at rest should develop opposing velocities after one step. Kinetic energy should increase while potential energy becomes more negative.

Multi-body energy. Three or more bodies should still show positive kinetic energy after stepping from rest.

Symmetry preservation. A symmetric configuration (mirror-image positions and velocities) should maintain its symmetry after a timestep. Verify that corresponding coordinates remain negations of each other.

Multi-step energy conservation. Run 100 or more timesteps and check that total energy remains within a small tolerance band of its initial value. Bodies should also move toward each other gravitationally.

Inside each consteval block, construct a std::vector<Body<float>> with the scenario you want to test, call do_timestep (or a loop of them), and VERIFY the expected invariants.

Commit to git

Verification: the figure-eight orbit

Validate against a known solution

A three-body system with equal masses, initialized with specific conditions, produces a famous periodic figure-eight orbit discovered by Chenciner and Montgomery in 2000. Use these initial conditions:
\[\begin{align*} \vec r_1&=(-0.97000436, 0.24308753) \\ \vec r_2&=(0.97000436, -0.24308753) \\ \vec r_3&=(0,0) \\ \vec v_1&=(0.4662036850, 0.4323657300)\\ \vec v_2&=(0.4662036850, 0.4323657300)\\ \vec v_3&=(-0.93240737, -0.86473146) \end{align*}\]

Ensure your main() initializes exactly these three bodies with mass $m=1$.

Run the simulation for many steps (try $\Delta t = 0.001$ and several thousand steps).

Monitor energy drift. It should remain very small.

Useful reference

A 3D interactive N-body simulator with preset configurations, including the figure-eight orbit, is available at trisolarchaos.com.

Plot the tracks

Add void print_positions_csv(std::span<const Body<T>> bodies, T energy) that emits one row per timestep: x1,y1,x2,y2,...,energy. Redirect output to a file and plot externally to visualize the trajectories.

Plot with gnuplot

plot 'nbody.csv' using 1:2 with lines title "P1", \
              '' using 3:4 with lines title "P2", \
              '' using 5:6 with lines title "P3"

plot 'nbody.csv' using 7 with lines title "Energy"

Commit to git

Code Review & Refactoring

Review

Before we measure performance, take your time to review the code you’ve written.

Your habit should become to

code and design simply,

then refactor mercilessly.

Optimizing unreadable code only accelerates bugs. Use this step to polish your implementation:

Refactoring Goals

1. Develop a Common Vocabulary

Read your function and variable names aloud. Do they match the physics equations at the top of this page? Do the functions make the code read like prose (self-documenting)?

If one partner writes acc and the other writes a, pick one and standardize it across the whole project. Consistent naming reduces cognitive load when you return to this code later.

2. Code and Design Simply

Hunt for magic numbers (e.g., hardcoded $G$, $\epsilon$, or time steps). Extract them into named constexpr variables.

Look for functions doing more than one thing. Can any logic be extracted into a smaller, focused helper?

Remove unused includes, stale comments, or dead code branches.

3. Refactor Mercilessly

Apply your agreed-upon coding standards uniformly. Fix indentation, spacing, and brace placement so the code reads smoothly. (use clang-format)

Ensure every function has a clear single responsibility. If you had to deviate from your goals document “the why”.

4. Integrate Continually & Verify

Run your tests again. Refactoring should never break correctness. If it does, fix the regression immediately.

Commit to git

with a message like refactor: clean up naming, extract constants, simplify integrator

Performance (Day 2)

Measure baseline performance

Increase $N$ to something substantial (e.g., 500 or 1000 bodies) and run enough steps to get measurable wall-clock time.

Build in Release mode: use ccmake in the build directory to set CMAKE_BUILD_TYPE to Release.

Benchmark with hyperfine ./nbody.

Profile with perf stat (prepend taskset 1 to inhibit CPU migrations).

For later

schedtool -F -p 10 -a 1 -e \
perf stat -e fp_arith_inst_retired.scalar_single,\
fp_arith_inst_retired.128b_packed_single,\
fp_arith_inst_retired.256b_packed_single \
sh -c './nbody >/dev/null'

Use *_double instead of *_single depending on the type of T.

CMake defaults

I disagree with some CMake defaults:

Release or RelWithDebInfo should be default

CMAKE_CXX_FLAGS_RELEASE defaults to -O3 -DNDEBUG; consider -O2 -DNDEBUG

CMAKE_CXX_FLAGS_RELWITHDEBINFO defaults to -O2 -g -DNDEBUG; consider -Og -g -DNDEBUG

A useful flag to add nowadays is -fhardened.

Think about data layout

Our struct Body stores all properties of a single body contiguously (Array of Structures). For the force computation, we repeatedly access only position and acceleration components across all bodies. Consider whether a Structure of Arrays layout might improve cache utilization. We will revisit this later.

Commit to git

Switch precision

Change T

Change using T = float; to using T = double; in main.cpp. Rebuild and compare energy drift over the same number of steps. Does the drift amplitude change? What about wall-clock time?

Payoff of templates

Because every algorithm is templated on T, switching precision requires changing exactly one line. Zero algorithm modifications are needed.

Commit to git if needed

Visualization

Visual simulation with SFML

For real-time visualization, integrate SFML 3. Add to nbody/CMakeLists.txt:

find_package(SFML 3 QUIET COMPONENTS system window graphics)
if(NOT SFML_FOUND)
  include(FetchContent)
  FetchContent_Declare(
    SFML
    GIT_REPOSITORY https://github.com/SFML/SFML.git
    GIT_TAG 3.1.0
    GIT_SHALLOW TRUE
    OVERRIDE_FIND_PACKAGE
  )
  option(SFML_BUILD_AUDIO "" OFF)
  option(SFML_BUILD_NETWORK "" OFF)
  FetchContent_MakeAvailable(SFML)
endif()
target_link_libraries(nbody PRIVATE SFML::System SFML::Window SFML::Graphics)

Create a simulate_to_sfml function with a structure like this:

sf::RenderWindow window(sf::VideoMode({1000, 1000}), "Window");
constexpr T dt = 0.0001;
while (window.isOpen())
  {
    while (std::optional event = window.pollEvent())
      { // Close window: exit
        if (event->is<sf::Event::Closed>())
          window.close();
      }
    window.clear();
    for (const auto& [r, ..._, m] : bodies)
      {
        sf::CircleShape point(m);
        // determine x and y
        point.setPosition({x, y});
        window.draw(point);
      }
    window.display();
    do_timestep<T>(dt, bodies);

    //using namespace std::literals;
    //std::this_thread::sleep_for(0.5ms);
  }

Render each body as an sf::CircleShape whose radius encodes its mass. Scale x and y coordinates into the window’s pixel space. Use structured bindings with the discard pattern [r, ..._, m] to unpack only the fields you need.

Call window.clear(), draw all bodies, then window.display() each frame. Note that SFML 3’s pollEvent() returns std::optional.

If the simulation runs too fast, use std::this_thread::sleep_for() between timesteps. Also consider different dt, G, and epsilon.

Commit to git

Further Topics (Optional)

Simulating Infinite Space

Change position vectors to use infinite space

In many astrophysical simulations, the region of interest is a small patch of an otherwise vast, homogeneous universe. Rather than simulating an impractically large volume, you can model space as periodic: any body that crosses one edge of the simulation domain re-enters from the opposite side. The gravitational force on body $i$ due to body $j$ is computed using only the nearest image of $j$, known as the minimum-image convention.

Hints

The key insight is that unsigned integer arithmetic wraps naturally modulo $2^{32}$, implementing periodic boundary conditions for free. By encoding positions as scaled unsigned integers, subtraction already produces the correct wrapped displacement. All that remains is choosing the shorter of the two possible directions.

The InfSpace<T> type. Create nbody/infspace.h defining a template struct parameterized by a compile-time floating-point domain length:
template <std::floating_point auto Length>
struct InfSpace
{
  using FpType = decltype(Length);
  static constexpr FpType length = Length;
  static constexpr FpType scale
    = Length / std::bit_floor(std::numeric_limits<unsigned>::max()) / 2;

  unsigned value;
The constructor accepts a floating-point coordinate and folds it into [0, Length) before dividing by scale to produce an unsigned value. Provide an explicit conversion operator back to FpType. See explicit conversions.

The delta abstraction. The critical operation is computing the signed shortest displacement between two periodic coordinates. For a scalar InfSpace:
friend constexpr FpType
delta(InfSpace a, InfSpace b)
{
  unsigned d = a.value - b.value;
  if (d < -d)
    return scale * d;
  else
    return -scale * -d;
}
Unsigned subtraction a.value - b.value wraps automatically. The expression -d also wraps, producing the complementary distance in the other direction. Comparing d < -d selects whichever is shorter. The trick of returning -scale * -d (rather than -(scale * d)) avoids negating the unsigned value directly and instead lets the compiler emit clean code.

Polymorphic Vec2d. Generalize Vec2d so it works with InfSpace coordinates. Add a free-function overload of delta for plain floating-point types constrained by the std::floating_point concept:
template <std::floating_point T>
  constexpr T
  delta(T a, T b) { return a - b; }
Then rewrite distance, distance2, and add a new delta(Vec2d, Vec2d) hidden friend inside Vec2d that delegates to the scalar delta:
friend constexpr auto
delta(Vec2d a, Vec2d b)
{
  auto x = delta(a.x, b.x);
  auto y = delta(a.y, b.y);
  return Vec2d<decltype(x)>(x, y);
}
Update distance and distance2 to call delta(a.x, b.x) instead of a.x - b.x. Use auto return types so they adapt to whatever delta produces.

Using it in Body. Change the position field in Body<T> to use InfSpace:
// Vec2d<T> r;  // old
Vec2d<InfSpace<T(3.5)>> r;  // new — domain size 3.5
Velocity and acceleration remain Vec2d<T>. The simulation force computation replaces (ri - rj) with delta(ri, rj), which returns a Vec2d<float> displacement suitable for force calculation regardless of whether positions are stored as floats or InfSpace.

Verify energy conservation. With periodic boundaries, energy should still oscillate around a constant value. Extend your consteval tests to confirm behavior under wrapping: place two bodies near opposite edges and verify they attract through the boundary.

Collision Handling

Can we replace $\epsilon$?

Instead of relying on the softening parameter $\epsilon$, you can implement explicit collision detection. Assign each body a radius $k_i$. After each simulation step, check every pair $(i, j)$ to see whether their separation satisfies:
\[|\vec{r}_i - \vec{r}_j| \leq k_i+k_j\]
When they do, merge them into a single body with mass $m_i + m_j$. The merged body’s position and velocity are those of the combined center of mass:
\[\begin{align*} \vec{r}_{new} &= \frac{m_i\vec{r}_i+m_j\vec{r}_j}{m_i+m_j} \\ \vec{v}_{new} &= \frac{m_i\vec{v}_i+m_j\vec{v}_j}{m_i+m_j} \end{align*}\]
Note: This models a perfectly inelastic collision. Momentum is conserved, but the total energy of the system can change.

Note: Keep careful track of the number of bodies remaining in the simulation as merges occur.

Acknowledgement

The introduction text is taken from an exercise written by Francesco Giacomini.

References

Argument-dependent lookup (ADL) — cppreference
constexpr specifier — cppreference
consteval specifier — cppreference
Defaulted comparisons — cppreference
explicit conversions — cppreference
std::floating_point concept — cppreference
std::format — cppreference
User-defined literals and operators — cppreference
std::optional — cppreference
std::this_thread::sleep_for — cppreference
std::span — cppreference
Structured bindings — cppreference
Velocity Verlet tutorial video — YouTube
trisolarchaos.com — interactive N-body simulator
GCC -fhardened option

N-Body System Simulation

N-Body System Simulation

Problem Description

Numerical Implementation

Energy Conservation

Project Setup

Step-by-step Implementation

Define `Vec2d<T>`

Define `Body<T>`

Compute gravitational forces

Implement Velocity Verlet integration

Readability checkpoint: comparing implementations

Energy conservation check

Compile-time simulation tests

Verification: the figure-eight orbit

Code Review & Refactoring

1. Develop a Common Vocabulary

2. Code and Design Simply

3. Refactor Mercilessly

4. Integrate Continually & Verify

Performance (Day 2)

Switch precision

Visualization

Further Topics (Optional)

Simulating Infinite Space

Collision Handling

Acknowledgement

References

results matching ""

No results matching ""

N-Body System Simulation

Problem Description

Numerical Implementation

Energy Conservation

Project Setup

Step-by-step Implementation

Define Vec2d<T>

Define Body<T>

Compute gravitational forces

Implement Velocity Verlet integration

Readability checkpoint: comparing implementations

Energy conservation check

Compile-time simulation tests

Verification: the figure-eight orbit

Code Review & Refactoring

1. Develop a Common Vocabulary

2. Code and Design Simply

3. Refactor Mercilessly

4. Integrate Continually & Verify

Performance (Day 2)

Switch precision

Visualization

Further Topics (Optional)

Simulating Infinite Space

Collision Handling

Acknowledgement

References

results matching ""

No results matching ""

Define `Vec2d<T>`

Define `Body<T>`