diff --git a/.editorconfig b/.editorconfig
index 0d79c83faacfc50ea18516e732524f5b95b89bda..486152c65ed246f1220b51b87e7b981c56b050fc 100644
--- a/.editorconfig
+++ b/.editorconfig
@@ -16,3 +16,4 @@ indent_size = 4
 
 [Makefile]
 indent_style = tab
+indent_size = 8
diff --git a/.vscode/launch.json b/.vscode/launch.json
index c6284236b9d4f7b67788f25bcb59f4b543e8c3fd..c9f4370a837d6bfa622a9e01da59547c00d56244 100644
--- a/.vscode/launch.json
+++ b/.vscode/launch.json
@@ -5,7 +5,7 @@
           "name": "(gdb) Launch",
           "type": "cppdbg",
           "request": "launch",
-          "program": "${workspaceFolder}/target/program",
+          "program": "${workspaceFolder}/target/debug/bin/planets",
           "args": [],
           "stopAtEntry": false,
           "cwd": "${fileDirname}",
diff --git a/Makefile b/Makefile
index 176b1510240c94d3326e3c406caf356e3e1c84fe..75dfc72981aea039e1532b6adf6f25a2976eba6f 100644
--- a/Makefile
+++ b/Makefile
@@ -1,41 +1,74 @@
-#The compiler
-CC:=gcc
-#The flags passed to the compiler
-CFLAGS:=-g -Ofast -Wall -Wextra -fsanitize=address -fsanitize=leak -std=gnu11
-#The flags passed to the linker
-LDFLAGS:=-lm -lSDL2
-#Path to the lib Vec2
-VPATH:=vec2 gfx planet
-
-# Source: https://yuukidach.github.io/p/makefile-for-projects-with-subdirectories/
-SRC_DIR = src
-TARGET_DIR = target
-SRCS = $(shell find $(SRC_DIR) -type f -name '*.c')
-OBJS = $(patsubst $(SRC_DIR)/%.c, $(TARGET_DIR)/%.o, $(SRCS))
-EXECUTABLE = $(TARGET_DIR)/program
-
-$(EXECUTABLE): $(OBJS)
-	$(CC) $(CFLAGS) $^ -o $@ $(LDFLAGS)
-
-$(TARGET_DIR)/%.o: $(SRC_DIR)/%.c
-	@mkdir -p $(@D)
-	$(CC) $(CFLAGS) -c $< -o $@
+# Programs
+SHELL			:= /bin/sh
+CC			:= gcc
+
+# Flags
+CFLAGS			:= -Wall -Wextra -std=gnu11 -DTRAILS_COUNT=2000 -Wno-unused-variable
+RELEASE_FLAGS		:= $(CFLAGS) -Ofast
+DEBUG_FLAGS		:= $(CFLAGS) -O0 -gdwarf-4 -g3 -fsanitize=address -fsanitize=leak
+LDFLAGS			:= -lm -lSDL2
+
+# Directories & Files
+OBJS			:= main.o gfx/gfx.o vec2/vec2.o celestial_body/celestial_body.o
+PROGRAM 		:= planets
+SRC_DIR			:= src
+TARGET_DIR		:= target
+BIN_DIR			:= bin
+OBJ_DIR			:= obj
+
+RELEASE_DIR		:= $(TARGET_DIR)/release
+RELEASE_OBJ		:= $(RELEASE_DIR)/$(OBJ_DIR)
+RELEASE_BIN		:= $(RELEASE_DIR)/$(BIN_DIR)
 
-default: build
+DEBUG_DIR		:= $(TARGET_DIR)/debug
+DEBUG_OBJ		:= $(DEBUG_DIR)/$(OBJ_DIR)
+DEBUG_BIN		:= $(DEBUG_DIR)/$(BIN_DIR)
 
-build: clean
-build: $(EXECUTABLE)
+default: build-debug
+run: run-debug
+build: build-debug
 
-build-debug: clean
-build-debug: CFLAGS += -gdwarf-4 -g3
-build-debug: $(EXECUTABLE)
+# Release
+run-release: build-release
+	$(RELEASE_BIN)/$(PROGRAM)
+
+build-release: TARGET_DIR := $(RELEASE_DIR)
+build-release: ensure-dir $(RELEASE_BIN)/$(PROGRAM)
+
+$(RELEASE_BIN)/$(PROGRAM): $(addprefix $(RELEASE_OBJ)/,$(OBJS))
+	$(CC) -o $@ $^ $(RELEASE_FLAGS) $(LDFLAGS)
+
+$(RELEASE_OBJ)/%.o: $(SRC_DIR)/%.c
+	@mkdir -p $(@D)
+	$(CC) -o $@ -c $^ -Isrc $(RELEASE_FLAGS)
 
-run: build
-	./$(EXECUTABLE)
+# Debug
+run-debug: build-debug
+	$(DEBUG_BIN)/$(PROGRAM)
+
+build-debug: TARGET_DIR := $(DEBUG_DIR)
+build-debug: ensure-dir $(DEBUG_BIN)/$(PROGRAM)
 
 debug: build-debug
 	gdb $(EXECUTABLE)
 
-.PHONY: clean
+$(DEBUG_BIN)/$(PROGRAM): $(addprefix $(DEBUG_OBJ)/,$(OBJS))
+	$(CC) -o $@ $^ $(DEBUG_FLAGS) $(LDFLAGS)
+
+$(DEBUG_OBJ)/%.o: $(SRC_DIR)/%.c
+	@mkdir -p $(@D)
+	$(CC) -o $@ -c $^ -Isrc $(DEBUG_FLAGS)
+
+.PHONY: clean clean-release clean-debug ensure-dir
 clean:
-	rm -rf $(TARGET_DIR)
+	@rm -rf $(TARGET_DIR)
+
+clean-release:
+	@rm -rf $(RELEASE_DIR)
+
+clean-debug:
+	@rm -rf $(DEBUG_DIR)
+
+ensure-dir:
+	@mkdir -p $(TARGET_DIR)/$(BIN_DIR)
+	@mkdir -p $(TARGET_DIR)/$(OBJ_DIR)
diff --git a/README.md b/README.md
index 5de477790ee9657faafe0829f74a6d8975b6f807..1c0ccd3eaa774b3457c226281d0fc18441285447 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,7 @@
-# TP Physique Planètes
-## Developpement
+# Orbital System Simulation
+## Requirements
+- SDL2
+
 ## Project Structure
 
 - `.devcontainer/`: Contains the configuration for the development container.
@@ -10,33 +12,54 @@
   - `vec2/`: Contains the 2D vector manipulation code.
 - `target/`: Contains the compiled binaries.
 
-## Requirements
-- SDL2
+## Compilation Flags
+| Name           | Default | Description                                                                         |
+| -------------- | ------- | ----------------------------------------------------------------------------------- |
+| `TRAILS_COUNT` | `2000`  | Display the last n positions as trails, if not defined, does not display any trails |
+
+## Memory Leaks
+Due to how SDL2 works, it is possible that ASAN detects memory leaks when running `mode_display`. This is not a real memory leak in the sense that the memory is not continually growing over time, but rather, it's a one-time allocation that is cleaned up by the operating system when the process ends. These allocations are made by SDL2 itself and are outside of our control.
 
-### Building
+## Building
 To build the project, run the following command:
 ```sh
-make build
+make build # or make build-debug
 ```
-*Note: this will completely rebuild the project*
 
 ## Running
-To run the project, run the following command:
+To build and run the project, run the following command:
 ```sh
-make run
+make run # or make run-debug
 ```
-*Note: this will also rebuild the project from scratch*
 
 ## Debugging
-By default the built binaries do not come with debug symbols. To build with debug symbols use:
+The `make build` command automatically builds with debug flags. If you wish to be dropped directly into `gdb` you can use
 ```sh
-make build-debug
+make debug
 ```
-*Note: this will also rebuild the project from scratch*
+
+## Release
+By default the built binaries are built in debug mode. To build with release flags use:
+```sh
+make build-release
+```
+
+*Note: you can also use `run-release` to build and run in release mode*
+
+
+*Note: all of these commands will also rebuild the project from scratch*
 
 VS Code debug configuration files are also provided for ease of use.
 
 ## Development Environement
 A devcontainer is provided for convenience. The devcontainer comes preconfigured for VSCode but should also work in other IDEs (provided that you configure them correctly).
 
-To be able to display the program's window on your desktop environement, the devcontainer excepts your X11 session to be at `/tmp/.X11-unix` and accessible from `local` (using `xhost +local:` on your local host). This assumses that you are using Linux.
+### GUI Forwarding
+#### Linux
+To be able to display the program's window on your desktop environement, the devcontainer excepts your X11 session to be at `/tmp/.X11-unix` and accessible from `local` (using `xhost +local:` on your local host).
+
+#### Windows
+To be able to display the program's window on your desktop environement on Windows you must have an X11 session. This is usually achieved by installing and running `vcxsrv`.
+Once installed you need to make some changes to `.devcontainer/devcontainer.json`:
+- Comment out `mounts`
+- Change `containerEnv.DISPLAY` to `${host.docker.internal:0.0}` (this assumes that your container can resolve the DNS, otherwise use your internal host IP)
diff --git a/src/celestial_body/celestial_body.c b/src/celestial_body/celestial_body.c
index 870516aaee3dc36a1378fe4b049189634df64a31..7c739e6c91ed5ea937bc23b2bb4156d868289904 100644
--- a/src/celestial_body/celestial_body.c
+++ b/src/celestial_body/celestial_body.c
@@ -1,9 +1,37 @@
 #include "celestial_body.h"
 #include <stdlib.h>
 
-celestial_body_t create_celestial_body(double mass, double semi_major_axis, double eccentricity, double size, uint32_t color) {
-    vec2 position = get_initial_position(semi_major_axis);
+double meters_to_au(double meters)
+{
+    return meters / AU;
+}
+
+celestial_body_t *create_celestial_body(system_t *system, char *name, celestial_body_t *central_body, double mass, double semi_major_axis, double eccentricity, uint32_t size, uint32_t color)
+{
+    vec2 position = get_initial_position(semi_major_axis, eccentricity, central_body);
+
+    #ifdef TRAILS_COUNT
+    vec2 *history = malloc(TRAILS_COUNT * sizeof(vec2));
+    if (history == NULL)
+    {
+        fprintf(stderr, "Failed to allocate memory for the history of the celestial body %s\n", name);
+        exit(EXIT_FAILURE);
+    }
+
+    for (int i = 0; i < TRAILS_COUNT; i++)
+    {
+        history[i] = vec2_create_zero();
+    }
+    #endif
+
+    if (central_body != NULL)
+    {
+        central_body->is_central_body = true;
+    }
+
     celestial_body_t celestial_body = {
+        .central_body = central_body,
+        .is_central_body = false,
         .mass = mass,
         .semi_major_axis = semi_major_axis,
         .eccentricity = eccentricity,
@@ -12,88 +40,216 @@ celestial_body_t create_celestial_body(double mass, double semi_major_axis, doub
 
         .size = size,
         .color = color,
+        .name = name,
+
+        #ifdef TRAILS_COUNT
+        .history = history,
+        #endif
     };
 
-    return celestial_body;
+    system->celestial_bodies[system->celestial_bodies_count] = celestial_body;
+    system->celestial_bodies_count++;
+    return &system->celestial_bodies[system->celestial_bodies_count - 1];
 }
 
-system_t create_system(double delta_t) {
+system_t create_system(double delta_t)
+{
     system_t system = {
-        .central_body = NULL,
-        .celestial_bodies_count = 3,
-        .celestial_bodies = malloc(3 * sizeof(celestial_body_t)),
+        .display_central_body = NULL,
+        .celestial_bodies_count = 0,
+        .celestial_bodies = malloc(30 * sizeof(celestial_body_t)),
     };
 
-    celestial_body_t sun = create_celestial_body(M_SUN, 0.0, 0.0, 35.0, 0xFFFF00);
-    celestial_body_t earth = create_celestial_body(M_EARTH, A_EARTH, E_EARTH, 20, 0x0000FF);
-    celestial_body_t mars = create_celestial_body(M_MARS, A_MARS, A_MARS, 15, 0xFF0000);
-
-    system.celestial_bodies[0] = sun;
-    system.celestial_bodies[1] = earth;
-    system.celestial_bodies[2] = mars;
+    celestial_body_t *sun = create_celestial_body(&system, "Sun", NULL, M_SUN, 0.0, 0.0, 30, 0xFFE484);
+    celestial_body_t *mercury = create_celestial_body(&system, "Mercury", sun, M_MERCURY, A_MERCURY, E_MERCURY, 10, 0x8C8887);
+    celestial_body_t *venus = create_celestial_body(&system, "Venus", sun, M_VENUS, A_VENUS, E_VENUS, 15, 0xECE9E0);
+    celestial_body_t *earth = create_celestial_body(&system, "Earth", sun, M_EARTH, A_EARTH, E_EARTH, 20, 0x4E6E94);
+    celestial_body_t *moon = create_celestial_body(&system, "Moon", earth, M_MOON, A_MOON, E_MOON, 10, 0xAAAAAA);
+    celestial_body_t *mars = create_celestial_body(&system, "Mars", sun, M_MARS, A_MARS, E_MARS, 15, 0xD39564);
+    celestial_body_t *tschai = create_celestial_body(&system, "Tschaï", sun, M_TSCHAI, A_TSCHAI, E_TSCHAI, 15, 0xFC9CFF);
+    celestial_body_t *albuchef = create_celestial_body(&system, "Albuchef", sun, M_ALBUCHEF, A_ALBUCHEF, E_ALBUCHEF, 25, 0x4DD679);
+    celestial_body_t *malaspinas = create_celestial_body(&system, "Malaspinas", albuchef, M_MALASPINAS, A_MALASPINAS, E_MALASPINAS, 15, 0xE8C4B7);
 
+    system.display_central_body = sun;
+    update_system_with_function(&system, delta_t, &get_first_position);
     return system;
 }
 
-void show_system(struct gfx_context_t *ctxt, system_t *system) {
+void show_system(struct gfx_context_t *ctxt, system_t *system, celestial_body_t *central_body)
+{
+    double radius = 0.0;
+    for (uint32_t i = 0; i < system->celestial_bodies_count; i++)
+    {
+        celestial_body_t *body = &system->celestial_bodies[i];
+        if (body->central_body != central_body)
+        {
+            continue;
+        }
+
+        double r = body->semi_major_axis * (1 + body->eccentricity);
+        if (r > radius)
+        {
+            radius = r;
+        }
+    }
+
+    radius *= DISPLAY_SCALE;
     for (uint32_t i = 0; i < system->celestial_bodies_count; i++)
     {
         celestial_body_t *body = &system->celestial_bodies[i];
-        vec2 position = vec2_mul(DISPLAY_SCALE, body->position);
+        if (body->central_body != central_body && body != central_body)
+        {
+            continue;
+        }
+
+        #ifdef TRAILS_COUNT
+        coordinates prev = {
+            .row = 0,
+            .column = 0,
+        };
+
+        // When zooming in, draw less trails
+        int start = central_body->central_body == NULL ? 0 : (TRAILS_COUNT * 0.95);
+        for (int j = start; j < TRAILS_COUNT; j++)
+        {
+            vec2 position = body->history[j];
+            if (body->central_body != central_body && body->central_body != NULL)
+            {
+                position = vec2_add(position, body->central_body->position);
+                position = vec2_sub(position, central_body->position);
+            }
+
+            position = vec2_mul(1.0 / radius, position);
+            coordinates coords = vec2_to_coordinates(position, ctxt->width, ctxt->height);
+
+            double alpha = ((j - start) / (double)(TRAILS_COUNT - start)) * 0.75;
+            uint32_t color = gfx_rgba_to_rgb(body->color, 0x000000, alpha);
+            if (
+                prev.row != 0 && prev.column != 0)
+            {
+                draw_line(ctxt, prev.column, prev.row, coords.column, coords.row, color);
+            }
+
+            prev = coords;
+        }
+        #endif
+
+        vec2 position = body->position;
+        position = vec2_sub(position, central_body->position);
+        position = vec2_mul(1.0 / radius, position);
         coordinates coords = vec2_to_coordinates(position, ctxt->width, ctxt->height);
         draw_full_circle(ctxt, coords.column, coords.row, body->size, body->color);
     }
 }
 
-void update_system(system_t *system, double delta_t) {
-    // TODO: Compute all the bodies and only then update their positions
-    return;
+void update_system_with_function(system_t *system, double delta_t, vec2 (*position_fn)(celestial_body_t *body, system_t *system, double delta_t))
+{
+    uint32_t bodies_count = system->celestial_bodies_count;
+    celestial_body_t *bodies = system->celestial_bodies;
+
+    // All functions are instantenous so we cannot
+    // partially update the system states otherwise
+    // the computation would be (ever so slightly) wrong
+    vec2 next_positions[bodies_count];
+    for (uint32_t i = 0; i < bodies_count; i++)
+    {
+        celestial_body_t *body = &bodies[i];
+        vec2 next_position = position_fn(body, system, delta_t);
+        next_positions[i] = next_position;
+    }
+
+    for (uint32_t i = 0; i < bodies_count; i++)
+    {
+        celestial_body_t *body = &bodies[i];
+        body->previous_position = body->position;
+        body->position = next_positions[i];
+
+        #ifdef TRAILS_COUNT
+        for (int j = 0; j < (TRAILS_COUNT - 1); j++)
+        {
+            body->history[j] = body->history[j + 1];
+        }
+
+        vec2 offset = vec2_create_zero();
+        if (body->central_body != NULL)
+        {
+            offset = body->central_body->position;
+        }
+
+        body->history[TRAILS_COUNT - 1] = vec2_sub(body->position, offset);
+        #endif
+    }
+}
+
+void update_system(system_t *system, double delta_t)
+{
+    update_system_with_function(system, delta_t, &get_next_position);
 }
 
-void free_system(system_t *system) {
+void free_system(system_t *system)
+{
+#ifdef TRAILS_COUNT
+    for (uint32_t i = 0; i < system->celestial_bodies_count; i++)
+    {
+        free(system->celestial_bodies[i].history);
+    }
+#endif
+
     free(system->celestial_bodies);
 }
 
-vec2 get_gravitational_force(celestial_body_t *a, celestial_body_t *b) {
+vec2 get_gravitational_force(celestial_body_t *a, celestial_body_t *b)
+{
     vec2 r = vec2_sub(b->position, a->position);
     double n = vec2_norm(r);
 
     double ma = a->mass;
     double mb = b->mass;
 
-    double factor = (ma * mb) / powf(n, 3.0);
+    double factor = (ma * mb) / pow(n, 3.0);
     return vec2_mul(G * factor, r);
 }
 
-vec2 get_initial_position(double semi_major_axis) {
-    return vec2_create(semi_major_axis, 0.0);
+vec2 get_initial_position(double semi_major_axis, double eccentricity, celestial_body_t *central_body)
+{
+    vec2 position = vec2_create(semi_major_axis * (1 - eccentricity), 0.0);
+    if (central_body != NULL)
+    {
+        position = vec2_add(position, central_body->position);
+    }
+
+    return position;
 }
 
-vec2 get_initial_velocity(celestial_body_t *body, celestial_body_t *central_body) {
-    if (body == central_body)
+vec2 get_initial_velocity(celestial_body_t *body)
+{
+    celestial_body_t *central_body = body->central_body;
+    if (central_body == NULL)
     {
         return vec2_create_zero();
     }
 
     double e = body->eccentricity;
     double a = body->semi_major_axis;
+    double perihelion = a * (1 - e);
+
+    vec2 rp = vec2_sub(central_body->position, body->position);
+    rp = vec2_create(-rp.y, rp.x);
+
+    double gm = G * central_body->mass;
+    double f = sqrt((gm * (1 + e)) / perihelion);
 
-    // Position should always be set to the perpendicular vector
-    // upon creation. Given that no update has been performed yet,
-    // it is safe to assume that the position is still the initial
-    // one.
-    vec2 rp = vec2_create(-body->position.y, body->position.x);
+    vec2 velocity = vec2_mul(1.0 / vec2_norm(rp), rp);
+    velocity = vec2_mul(f, velocity);
 
-    // It's a bit long, but coming up with names for all the
-    // intermediate variables would make it even longer.
-    vec2 velocity = vec2_mul(
-        sqrt((G * central_body->mass * (1 + e)) / (a * (1 - e))) / vec2_norm(rp),
-        rp);
+    vec2 offset = get_initial_velocity(central_body);
+    velocity = vec2_add(velocity, offset);
 
     return velocity;
 }
 
-vec2 get_instantaneous_forces(celestial_body_t *body, system_t *system) {
+vec2 get_instantaneous_forces(celestial_body_t *body, system_t *system)
+{
     vec2 total_gravitational_force = vec2_create_zero();
     for (uint32_t i = 0; i < system->celestial_bodies_count; i++)
     {
@@ -110,17 +266,31 @@ vec2 get_instantaneous_forces(celestial_body_t *body, system_t *system) {
     return total_gravitational_force;
 }
 
-vec2 get_instantaneous_acceleration(celestial_body_t *body, system_t *system) {
+vec2 get_instantaneous_acceleration(celestial_body_t *body, system_t *system)
+{
     vec2 instantaneous_forces = get_instantaneous_forces(body, system);
     return vec2_mul(1.0 / body->mass, instantaneous_forces);
 }
 
-vec2 get_first_position(celestial_body_t *body, system_t *system, double delta_t) {
-    // TODO
-    return vec2_create_zero();
+vec2 get_first_position(celestial_body_t *body, system_t *system, double delta_t)
+{
+    vec2 initial_position = body->position;
+    vec2 velocity = get_initial_velocity(body);
+    vec2 acceleration = get_instantaneous_acceleration(body, system);
+
+    vec2 displacement = vec2_mul(pow(delta_t, 2.0) / 2.0, acceleration);
+    vec2 next_position = vec2_add(vec2_add(initial_position, vec2_mul(delta_t, velocity)), displacement);
+    return next_position;
 }
 
-vec2 get_next_position(celestial_body_t *body, system_t *system, double delta_t) {
-    // TODO
-    return vec2_create_zero();
+vec2 get_next_position(celestial_body_t *body, system_t *system, double delta_t)
+{
+    vec2 acceleration = get_instantaneous_acceleration(body, system);
+
+    vec2 current_position = vec2_mul(2.0, body->position);
+    vec2 previous_position = body->previous_position;
+    vec2 displacement = vec2_mul(pow(delta_t, 2.0), acceleration);
+
+    vec2 next_position = vec2_add(vec2_sub(current_position, previous_position), displacement);
+    return next_position;
 }
diff --git a/src/celestial_body/celestial_body.h b/src/celestial_body/celestial_body.h
index 44df39aaca87217daad3c0b5d0014feb20996729..e8e601edcb8899b75f217687efe0cb51cc79f4bb 100644
--- a/src/celestial_body/celestial_body.h
+++ b/src/celestial_body/celestial_body.h
@@ -4,29 +4,72 @@
 #include "../vec2/vec2.h"
 #include "../gfx/gfx.h"
 
-// Display scale to fit our system on the window
-#define DISPLAY_SCALE (2e-12)
+#define DISPLAY_SCALE 1.1 // 110% of the maximum semi-major axis
 
 // Gravitational constant
 #define G 6.67408e-11 // [m^3 kg^-1 s^-2]
 
+// Astronomic unit
+#define AU (1.495978707e+11)
+
 // Sun
-#define M_SUN 1.989e30 // Mass [kg]
+// source: https://nssdc.gsfc.nasa.gov/planetary/factsheet/sunfact.html
+#define M_SUN 1988500e24 // Mass [kg]
+
+// Mercury
+// source: https://nssdc.gsfc.nasa.gov/planetary/factsheet/mercuryfact.html
+#define M_MERCURY 0.33010e24 // Mass [kg]
+#define E_MERCURY 0.2056     // Eccentricity
+#define A_MERCURY 57.909e9   // Semi-major axis [m]
+
+// Venus
+// source: https://nssdc.gsfc.nasa.gov/planetary/factsheet/venusfact.html
+#define M_VENUS 4.8673e24
+#define E_VENUS 0.0068
+#define A_VENUS 108.210e9
 
 // Earth
-#define M_EARTH 5.972e24 // Mass [kg]
-#define E_EARTH 0.01671022 // Eccentricity
-#define A_EARTH 149597870700 // Semi-major axis [m]
+// source: https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
+#define M_EARTH 5.9722e24
+#define E_EARTH 0.0167
+#define A_EARTH 149.598e9
+
+// Moon
+// source: https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html
+#define M_MOON 0.07346e24
+#define E_MOON 0.0549
+#define A_MOON 0.384e9
 
 // Mars
-#define M_MARS 6.39e23 // Mass [kg]
-#define E_MARS 0.0934 // Eccentricity
-#define A_MARS 227939200000 // Semi-major axis [m]
+// source: https://nssdc.gsfc.nasa.gov/planetary/factsheet/marsfact.html
+#define M_MARS 0.64169e24
+#define E_MARS 0.0935
+#define A_MARS 227.956e9
+
+// Fictional planets
+#define M_TSCHAI 2.9862e24
+#define E_TSCHAI 0.42
+#define A_TSCHAI 420.302e9
+
+#define M_ALBUCHEF 7.3461e24
+#define E_ALBUCHEF 0.0549
+#define A_ALBUCHEF 384.400e9
+
+#define M_MALASPINAS 4.7998e22
+#define E_MALASPINAS 0.0094
+#define A_MALASPINAS 671.100e6
 
 /// @brief Structure representing a celestial body
 typedef struct _celestial_body
 {
     /// @section Physical properties
+    /// @brief Central body of the body
+    /// @example The Sun for the Earth, the Earth for the Moon, ...
+    struct _celestial_body *central_body;
+
+    /// @brief Whether the body is the central body another body
+    bool is_central_body;
+
     /// @brief Mass of the body [kg]
     double mass;
 
@@ -44,33 +87,52 @@ typedef struct _celestial_body
 
     /// @section Display properties
     /// @brief Screen size of the body [px]
-    double size;
+    uint32_t size;
 
     /// @brief Color of the body [RGB]
     uint32_t color;
+
+    /// @brief Name of the body
+    char *name;
+
+    #ifdef TRAILS_COUNT
+    vec2 *history;
+    #endif
 } celestial_body_t;
 
 /// @brief Structure representing a planetary system
 typedef struct _system
 {
-    /// @brief Central body of the system
-    /// @example The sun
-    celestial_body_t *central_body;
+    /// @brief Body at the center of the screen (and not neccessarely the central body of the system)
+    celestial_body_t *display_central_body;
 
     /// @brief Number of bodies currently in the system (including the central body)
     uint32_t celestial_bodies_count;
 
+    /// @brief Number of bodies that can be stored in the system
+    uint32_t celestial_bodies_capacity;
+
     /// @brief Array of bodies
     /// @note For convenience, the central body is included in this array
     celestial_body_t *celestial_bodies;
 } system_t;
 
+/// @brief Convert meters to AU (Astronomical Unit)
+/// @param meters Distance in meters
+/// @return Distance in AU
+double meters_to_au(double meters);
+
 /// @brief Create a celestial body
+/// @param system System in which to create the body
+/// @param central_body Central body of the body
+/// @param name Name of the body
 /// @param mass Mass of the body [kg]
 /// @param semi_major_axis Semi-major axis of the body [m]
 /// @param eccentricity Eccentricity of the body's orbit (between 0 and 1)
+/// @param size Screen size of the body [px]
+/// @param color Color of the body [RGB]
 /// @return The created celestial body
-celestial_body_t create_celestial_body(double mass, double semi_major_axis, double eccentricity, double size, uint32_t color);
+celestial_body_t *create_celestial_body(system_t *system, char *name, celestial_body_t *central_body, double mass, double semi_major_axis, double eccentricity, uint32_t size, uint32_t color);
 
 /// @brief Create a planetary system
 /// @param delta_t Initial time step [s]
@@ -80,7 +142,14 @@ system_t create_system(double delta_t);
 /// @brief Show the system on the screen
 /// @param ctxt Graphics context
 /// @param system System to show
-void show_system(struct gfx_context_t *ctxt, system_t *system);
+/// @param central_body Body at the center of the screen
+void show_system(struct gfx_context_t *ctxt, system_t *system, celestial_body_t *central_body);
+
+/// @brief Update the system by one step using a function
+/// @param system System to update
+/// @param delta_t Time step [s]
+/// @param position_fn Function to retrieve the next position of a body
+void update_system_with_function(system_t *system, double delta_t, vec2 (*position_fn)(celestial_body_t *body, system_t *system, double delta_t));
 
 /// @brief Update the system by one step
 /// @param system System to update
@@ -100,16 +169,17 @@ vec2 get_gravitational_force(celestial_body_t *a, celestial_body_t *b);
 
 /// @brief Get the initial position of a body
 /// @param semi_major_axis Semi-major axis of the body's orbit
+/// @param eccentricity Eccentricity of the body's orbit (between 0 and 1)
+/// @param central_body Central body of the body
 /// @return Initial position of the body
 /// @note This is dervied from Figure 2 from the assignment
-vec2 get_initial_position(double semi_major_axis);
+vec2 get_initial_position(double semi_major_axis, double eccentricity, celestial_body_t *central_body);
 
 /// @brief Get the initial velocity of a body
 /// @param body Body for which to get the initial velocity
-/// @param central_body Body around which the body orbits
 /// @return Initial velocity of the body
 /// @note This is formula (6) from the assignment
-vec2 get_initial_velocity(celestial_body_t *body, celestial_body_t *central_body);
+vec2 get_initial_velocity(celestial_body_t *body);
 
 /// @brief Get the instantaneous forces applied on a body
 /// @param body Body for which to get the forces
diff --git a/src/gfx/gfx.c b/src/gfx/gfx.c
index 3821077bbbf7964c49e73b9b94d6dc123bcb9d9f..8c65325b108d531e082933840b9820320c9255fc 100644
--- a/src/gfx/gfx.c
+++ b/src/gfx/gfx.c
@@ -10,6 +10,29 @@
 #include "gfx.h"
 #include <assert.h>
 
+/// @brief Convert an RGBA color to RGB
+/// @param foreground Color of the foreground
+/// @param background Color of the background
+/// @param alpha Transparency of the foreground (between 0 and 1)
+/// @return
+uint32_t gfx_rgba_to_rgb(uint32_t foreground, uint32_t background, double alpha)
+{
+    double fr = ((foreground >> 16) & 0xFF) / 255.0;
+    double fg = ((foreground >> 8) & 0xFF) / 255.0;
+    double fb = (foreground & 0xFF) / 255.0;
+
+    double br = ((background >> 16) & 0xFF) / 255.0;
+    double bg = ((background >> 8) & 0xFF) / 255.0;
+    double bb = (background & 0xFF) / 255.0;
+
+    double r = (1.0 - alpha) * br + alpha * fr;
+    double g = (1.0 - alpha) * bg + alpha * fg;
+    double b = (1.0 - alpha) * bb + alpha * fb;
+
+    uint32_t color = ((uint32_t)(r * 255.0) << 16) | ((uint32_t)(g * 255.0) << 8) | (uint32_t)(b * 255.0);
+    return color;
+}
+
 /// Create a fullscreen graphic window.
 /// @param title Title of the window.
 /// @param width Width of the window in pixels.
@@ -107,6 +130,45 @@ SDL_Keycode gfx_keypressed()
     return 0;
 }
 
+void draw_line(struct gfx_context_t *ctxt, int32_t x1, int32_t y1, int32_t x2, int32_t y2, uint32_t color)
+{
+    // Calculate differences and absolute values
+    int dx = abs(x2 - x1);
+    int dy = abs(y2 - y1);
+    int sx = (x1 < x2) ? 1 : -1;
+    int sy = (y1 < y2) ? 1 : -1;
+
+    // Initial decision parameter
+    int decision = (dx > dy ? dx : -dy) / 2;
+    int error;
+
+    // Start drawing the line
+    while (1)
+    {
+        // Plot the current pixel
+        gfx_putpixel(ctxt, x1, y1, color);
+
+        // Check if we reached the end point
+        if (x1 == x2 && y1 == y2)
+        {
+            break;
+        }
+
+        // Update decision parameter and position
+        error = decision;
+        if (error > -dx)
+        {
+            decision -= dy;
+            x1 += sx;
+        }
+        if (error < dy)
+        {
+            decision += dx;
+            y1 += sy;
+        }
+    }
+}
+
 /// Draw a full circle using Andres's discrete circle algorithm.
 /// @param ctxt Graphic context to clear.
 /// @param c_column X coordinate of the circle center.
diff --git a/src/gfx/gfx.h b/src/gfx/gfx.h
index 20b7815e4c5ef7a8daa871153c3f4affeeb448b6..9eab186526829b30a21d55c42ad02d6d6ae04c48 100644
--- a/src/gfx/gfx.h
+++ b/src/gfx/gfx.h
@@ -29,6 +29,7 @@ struct gfx_context_t
     uint32_t height;
 };
 
+extern uint32_t gfx_rgba_to_rgb(uint32_t foreground, uint32_t background, double alpha);
 extern void gfx_putpixel(
     struct gfx_context_t *ctxt, uint32_t column, uint32_t row, uint32_t color);
 extern void gfx_clear(struct gfx_context_t *ctxt, uint32_t color);
@@ -37,5 +38,6 @@ extern void gfx_destroy(struct gfx_context_t *ctxt);
 extern void gfx_present(struct gfx_context_t *ctxt);
 extern SDL_Keycode gfx_keypressed();
 extern void draw_full_circle(struct gfx_context_t *ctxt, uint32_t c_column, uint32_t c_row, uint32_t r, uint32_t color);
+void draw_line(struct gfx_context_t *ctxt, int32_t x1, int32_t y1, int32_t x2, int32_t y2, uint32_t color);
 
 #endif
diff --git a/src/main.c b/src/main.c
index a969730b7a3525a460a524de64ebdffc4767c24a..3c9150a93c16d5c6eef0ab036b69bafaf3745504 100644
--- a/src/main.c
+++ b/src/main.c
@@ -8,7 +8,9 @@
 #define SCREEN_WIDTH 1000
 #define SCREEN_HEIGHT 1000
 
-#define DELTA_T (60 * 60 * 6) // 6 hours
+// Update the system by `UPDATE_STEP` seconds every `UPDATE_RATE`
+#define UPDATE_RATE 24            // [Hz]
+#define UPDATE_STEP (60 * 60 * 4) // 4 hours
 
 /// @brief Memory area containing the currently processed event
 /// @note This variable is allocated once and kept until the end of the program
@@ -16,21 +18,69 @@
 ///       This is possible because events are handled one by one by our code.
 SDL_Event *event;
 
+/// @brief Flag indicating if the user has exited the program
+bool has_exited = false;
+
 /// @brief Check if the user has exited the program
 /// @return true if the user has exited the program, false otherwise
-bool has_exited()
+void handle_events(system_t *system)
 {
     while (SDL_PollEvent(event) == 1)
     {
         switch (event->type)
         {
         case SDL_QUIT:
-            return true;
+            has_exited = true;
+            return;
 
         case SDL_KEYDOWN:
-            if (event->key.keysym.sym == SDLK_ESCAPE)
+            switch (event->key.keysym.sym)
+            {
+            case SDLK_ESCAPE:
+                has_exited = true;
+                return;
+            case SDLK_RIGHT:
             {
-                return true;
+                // Find where our current central body is
+                int i;
+                for (i = 0; &system->celestial_bodies[i] != system->display_central_body; i++);
+
+                for (uint32_t j = i + 1; j < system->celestial_bodies_count; j++)
+                {
+                    // Find the next body which is a central body
+                    celestial_body_t *body = &system->celestial_bodies[j];
+                    if (body->is_central_body)
+                    {
+                        system->display_central_body = body;
+                        break;
+                    }
+                }
+                break;
+            }
+
+            case SDLK_LEFT:
+            {
+                // Find where our current central body is
+                int i;
+                for (i = 0; &system->celestial_bodies[i] != system->display_central_body; i++);
+
+                if (i == 0)
+                {
+                    break;
+                }
+
+                for (uint32_t j = i - 1; j > 0; j--)
+                {
+                    // Find the previous body which is a central body
+                    celestial_body_t *body = &system->celestial_bodies[j];
+                    if (body->is_central_body)
+                    {
+                        system->display_central_body = body;
+                        break;
+                    }
+                }
+                break;
+            }
             }
             break;
 
@@ -38,11 +88,9 @@ bool has_exited()
             break;
         }
     }
-
-    return false;
 }
 
-int main()
+bool mode_display(system_t *system)
 {
     srand(time(NULL));
     struct gfx_context_t *ctxt =
@@ -50,25 +98,88 @@ int main()
     if (!ctxt)
     {
         fprintf(stderr, "Graphics initialization failed!\n");
-        return EXIT_FAILURE;
+        return false;
     }
 
     // Events are handled one by one so we can allocate once
     // and keep this alloc until the end of the program
     event = malloc(sizeof(SDL_Event));
 
-    system_t system = create_system(DELTA_T);
-
-    while (!has_exited())
+    while (!has_exited)
     {
-        update_system(&system, DELTA_T);
+        update_system(system, UPDATE_STEP);
+
+        handle_events(system);
         gfx_present(ctxt);
         gfx_clear(ctxt, COLOR_BLACK);
-        show_system(ctxt, &system);
+        show_system(ctxt, system, system->display_central_body);
+
+        // Avoid using 100% CPU for nothing
+        SDL_Delay((1.0 / UPDATE_RATE) * 10);
     }
 
-    free_system(&system);
     gfx_destroy(ctxt);
     free(event);
-    return EXIT_SUCCESS;
+    return true;
+}
+
+void mode_avg_dist(system_t *system, celestial_body_t *body, double sample_rate, double samples)
+{
+    double avg_distances[system->celestial_bodies_count];
+    for (int i = 0; i < samples; i++)
+    {
+        update_system(system, sample_rate);
+        for (uint32_t j = 0; j < system->celestial_bodies_count; j++)
+        {
+            celestial_body_t *other_body = &system->celestial_bodies[j];
+            // Skip the body we're looking at
+            if (other_body == body)
+            {
+                continue;
+            }
+
+            if (other_body == body->central_body)
+            {
+                continue;
+            }
+
+            // Initialize the average distance to 0
+            if (i == 0)
+            {
+                avg_distances[j] = 0;
+            }
+
+            vec2 vec = vec2_sub(other_body->position, body->position);
+            double distance = vec2_norm(vec);
+            avg_distances[j] += distance;
+        }
+    }
+
+    printf("Average distance from %s to:\n", body->name);
+    for (uint32_t i = 0; i < system->celestial_bodies_count; i++)
+    {
+        celestial_body_t *other_body = &system->celestial_bodies[i];
+        if (other_body == body || other_body->central_body != body->central_body)
+        {
+            continue;
+        }
+
+        avg_distances[i] /= samples;
+
+        // Convert to AU (Astronomical Unit) to make the distances
+        // more readable
+        double au = meters_to_au(avg_distances[i]);
+        double rounded = round(au * 100) / 100;
+        printf("\t%s: %f [AU]\n", other_body->name, rounded);
+    }
+}
+
+int main()
+{
+    system_t system = create_system(UPDATE_STEP);
+    mode_avg_dist(&system, &system.celestial_bodies[3], UPDATE_STEP, 10000);
+
+    bool is_display_success = mode_display(&system);
+    free_system(&system);
+    return is_display_success ? EXIT_SUCCESS : EXIT_FAILURE;
 }