----------------------[ Making a simple physics engine: ]----------------------- Publication date: 1685923200 The simplest way to do a stable physics simulation is with Verlet integration: xn + 1 = 2xn - xn - 1 + a Δt2 To make it easier to see what is going on here, lets rearrange the equation a bit: xn + 1 = xn + (xn - xn - 1) + a Δt2 Instead of storing a velocity, it is computed from the current and last positions. This has the advantage that it is collisions can just update the velocity, and it works. These are some basic functions for working with 2d vectors, because C has no built in vector type: /////////////////////////////// // low level math functions. // /////////////////////////////// typedef struct Vector2 { float x; float y; } Vector2; Vector2 vector_add(Vector2 v1, Vector2 v2) { Vector2 output = {.x = v1.x + v2.x, .y = v1.y + v2.y}; return output; } Vector2 vector_sub(Vector2 v1, Vector2 v2) { Vector2 output = {.x = v1.x - v2.x, .y = v1.y - v2.y}; return output; } Vector2 vector_mul_scaler(Vector2 v1, float s) { Vector2 output = {.x = v1.x * s, .y = v1.y * s}; return output; } float vector_length(Vector2 v1) { return sqrt(v1.x * v1.x + v1.y * v1.y); } For each physics object, track the current position, last position. Acceleration here will be used to sum up the effect of forces acting on the object every simulation step: typedef struct Body { float radius; Vector2 position_old; Vector2 position; Vector2 acceleration; } Body; Time for the actual Verlet integration. The timestep is how much time passes between every simulation step. If the simulation is running at 60 updates per second, it should be 1/60 (~0.017). // Use verlet integration to apply velocity and acceleration to the body void update_position(Body* body, float dt) { // Compute velocity in terms of timestep based on position Vector2 velocity = vector_sub(body->position, body->position_old); // Save the position before updating body->position_old = body->position; // Move the object by the velocity, and the acceleration. body->position = vector_add(body->position,vector_add(velocity, vector_mul_scaler(body->acceleration, dt * dt ))); // Reset acceleration body->acceleration.x = 0; body->acceleration.y = 0; } One physics object is all right, but some more would be nice. This is a struct that holds a collection of physics object (and some rather uninteresting functions to manage memory). Depending on your programming language, this could just become a list of objects: typedef struct World { Body* objects; int size; // How many objects are in the world at the given moment int capacity; // The total amount that can be stored at a given time } World; World world_with_capacity(int capacity) { World w = { .objects = malloc(capacity * sizeof(Body)), .size = 0, .capacity = capacity }; return w; } // Returns 1 if sucessfull, zero if not int insert_object(World* world, Body object) { if (world->capacity > world->size) { world->objects[world->size] = object; world->size++; return 1; } else { return 0; } } // Frees internal datastructures of a world, does not free the passed pointer. void cleanup_world(World* w) { free(w->objects); // Set these to make sure no one tries to access any of the freed data; w->size = 0; w->capacity = 0; w->objects = 0; } void update_positions(World* w, float dt) { for (int i = 0; i < w->size; i++) { update_position(>objects[i], dt); } } void apply_gravity(World* w) { for (int i = 0; i < w->size; i++) { w->objects[i].acceleration.y -= 9.8; } } Now, simply calling `update_positions` and `apply_gravity`, can simulate a whole bunch of objects at once. Objects falling into an infinite void is fine and all, but isn't very intresting. It needs collisions. This is where the magic of Verlet integration comes in, all that has to be done is to check if the object is outside of the bounds, and if it isn't, just move it inside. Verlet integration takes care of all the dynamics, and the object ends up nicely sliding around inside of the sphere: void collide(World* w) { // Bound the center of an object to 9 units away from (0, 0) for (int i = 0; i < w->size; i++) { if (vector_length(w->objects[i].position) > 9) { float position_scale = 9 / vector_length(w->objects[i].position); // position = position * scale w->objects[i].position = vector_mul_scaler(w->objects[i].position,position_scale); } } } It would also be nice of the objects would collide with one another. Continuing the collide() function, There are much faster methods then this \\(O(n^2)\\) one, but it is straightforward: void collide(World* w) { ... // Check object colision, this is super inefficent. for (int i = 0; i < w->size; i++) { // This iterates up to i and not w->size to avoid rendundent checks and checking an object against itself for (int e = 0; e < i; e++) { // The closest the objects can be without intersecting float mindistance = w->objects[i].radius + w->objects[e].radius; Vector2* object1 = >objects[i].position; Vector2* object2 = >objects[e].position; Vector2 difference = vector_sub(*object1, *object2); // The current distance float distance = vector_length(difference); // Check for intersections if (mindistance > distance) { // How much each object has to be moved in order to not be intersecting float delta = (mindistance - distance)/2; // adjustment = normalize(difference) * deta Vector2 adjustment = vector_mul_scaler(vector_mul_scaler(difference, 1.0/distance), delta); // Move the objects by that distance *object1 = vector_add(*object1, adjustment); *object2 = vector_sub(*object2, adjustment); } } } } To use it, just call apply_gravity(), update_positions(), and collide() in a loop. That's it, a working physics engine. Currently, collisions are inelastic, meaning that all kinetic energy is absorbed (no bounciness), but there is no friction. Although at high timesteps objects sliding on a curved surface will slow down, but this is an artifact of limited integration steps. The collision solver can cause problems with large amounts of objects, this can be fixed by calling `collide` multiple times during a single simulation step. Code here Video here