main.cpp

#include <vector>
#include <memory>
#include <thread>
#include <array>
#include <cstring>
#include "raylib.h"
#include "C:/raylib/raylib/src/raymath.h"
#include "C:/raylib/raylib/src/rlgl.h"

constexpr int SCREEN_SIZE = 800;
constexpr int MAX_BODIES = 3000;
float G = 1000.0f;
constexpr float THETA = 1.0f;
constexpr float COR = 0.8f;
constexpr float DT = 1.0f / 60.0f;
constexpr float DAMPING = 0.99f;

bool spawn_on_load = true;
bool collisionsEnabled = true;

// Adjust to reduce force spikes
constexpr float SOFTENING = 1000.0f;

// Color palette

Color BG_COLOR = RAYWHITE;
Color GRID_LINE_COLOR = LIGHTGRAY;
Color GRID_OUTLINE_COLOR = DARKGRAY;
Color TEXT_COLOR = BLACK;
Color INPUT_TEXT_COLOR = RED;
Color HINT_TEXT_COLOR = GRAY;
Color VALUE_COLOR = DARKBLUE;

float minX = -5000, minY = -5000, maxX = 5000, maxY = 5000;

Camera2D camera = {0};

// SIMD-friendly SoA data layout
struct alignas(32) BodyData
{
    float x[MAX_BODIES];
    float y[MAX_BODIES];
    float vx[MAX_BODIES];
    float vy[MAX_BODIES];
    float mass[MAX_BODIES];
    float radius[MAX_BODIES];
    float kinEnergy[MAX_BODIES];
    float potEnergy[MAX_BODIES];
    Color color[MAX_BODIES];
    int count = 0;

    void ClampVelocity(int idx)
    {
        const float MAX_SPEED = 500.0f;
        float speed_sq = vx[idx] * vx[idx] + vy[idx] * vy[idx];
        if (speed_sq > MAX_SPEED * MAX_SPEED)
        {
            float scale = MAX_SPEED / sqrtf(speed_sq);
            vx[idx] *= scale;
            vy[idx] *= scale;
        }
    }
};

struct Rect
{
    float x, y, width, height;
};
struct QuadNode
{
    Rect boundary;
    float mass = 0;
    Vector2 com{0};
    bool leaf = true;
    std::array<int, 4> children{-1, -1, -1, -1};
};

std::vector<QuadNode> nodePool;
int nodePointer = 0;

void ResetNodePool()
{
    nodePool.clear();
    nodePool.reserve(MAX_BODIES * 16);
    nodePool.push_back(QuadNode{Rect{0, 0, SCREEN_SIZE, SCREEN_SIZE}});
    nodePointer = 0;
}

void SubdivideNode(int parentIdx)
{
    QuadNode &parent = nodePool[parentIdx];
    const float hw = parent.boundary.width / 2;
    const float hh = parent.boundary.height / 2;
    const float x = parent.boundary.x;
    const float y = parent.boundary.y;

    // Reserve space first to prevent reallocation
    const size_t new_size = nodePool.size() + 4;
    if (new_size > nodePool.capacity())
    {
        nodePool.reserve(nodePool.capacity() * 2);
    }

    parent.children = {
        static_cast<int>(nodePool.size()),
        static_cast<int>(nodePool.size()) + 1,
        static_cast<int>(nodePool.size()) + 2,
        static_cast<int>(nodePool.size()) + 3};

    nodePool.emplace_back(QuadNode{Rect{x + hw, y, hw, hh}});
    nodePool.emplace_back(QuadNode{Rect{x, y, hw, hh}});
    nodePool.emplace_back(QuadNode{Rect{x, y + hh, hw, hh}});
    nodePool.emplace_back(QuadNode{Rect{x + hw, y + hh, hw, hh}});

    parent.leaf = false;
}

void ComputeForces(BodyData &bodies, int start, int end, int nodeIdx)
{
    std::vector<int> stack{nodeIdx};
    const float SOFTENING = 100.0f; // Add force softening

    for (int i = start; i < end; ++i)
    {
        bodies.potEnergy[i] = 0.0f; // Reset potential energy

        float fx = 0, fy = 0;
        stack = {0};

        while (!stack.empty())
        {
            int currentIdx = stack.back();
            stack.pop_back();
            QuadNode &current = nodePool[currentIdx];

            if (!current.leaf)
            {
                float dx = current.com.x - bodies.x[i];
                float dy = current.com.y - bodies.y[i];
                float dist_sq = dx * dx + dy * dy + SOFTENING;
                float dist = sqrtf(dist_sq);
                float ratio = current.boundary.width / dist;

                if (ratio < THETA)
                {
                    float F = G * bodies.mass[i] * current.mass / dist_sq;
                    fx += F * dx / dist;
                    fy += F * dy / dist;
                }
                else
                {
                    for (int child : current.children)
                        if (child != -1)
                            stack.push_back(child);
                }
            }
            else if (current.mass > 0)
            {
                float dx = current.com.x - bodies.x[i];
                float dy = current.com.y - bodies.y[i];
                float dist_sq = dx * dx + dy * dy + SOFTENING;
                float dist = sqrtf(dist_sq);
                float F = G * bodies.mass[i] * current.mass / dist_sq;

                bodies.potEnergy[i] += F;
                fx += F * dx / dist;
                fy += F * dy / dist;
            }
        }

        bodies.vx[i] += fx * DT / bodies.mass[i];
        bodies.vy[i] += fy * DT / bodies.mass[i];
        bodies.ClampVelocity(i); // Velocity clamping

        // Kinetic energy
        bodies.kinEnergy[i] = 0.5f * bodies.mass[i] * (bodies.vx[i] * bodies.vx[i] + bodies.vy[i] * bodies.vy[i]);
    }
}

void BuildTree(BodyData &bodies)
{
    ResetNodePool();

    for (int i = 0; i < bodies.count; ++i)
    {
        int currentIdx = 0;
        int depth = 0;

        while (true)
        {
            QuadNode &current = nodePool[currentIdx]; // Safe reference after reservation

            if (current.leaf)
            {
                if (current.mass == 0)
                {
                    // Store body in empty leaf
                    current.mass = bodies.mass[i];
                    current.com = {bodies.x[i], bodies.y[i]};
                    break;
                }
                else
                {
                    // Subdivide and reinsert existing mass first
                    const Vector2 oldCOM = current.com;
                    const float oldMass = current.mass; // Now safe

                    SubdivideNode(currentIdx);
                    current.mass = 0;
                    current.com = {0, 0};

                    // Reinsert old mass
                    const float cx = current.boundary.x + current.boundary.width / 2;
                    const float cy = current.boundary.y + current.boundary.height / 2;
                    int quadrant = (oldCOM.x > cx) | ((oldCOM.y > cy) << 1);
                    int childIdx = nodePool[currentIdx].children[quadrant];
                    nodePool[childIdx].mass = oldMass;
                    nodePool[childIdx].com = oldCOM;
                }
            }

            // Update current node's mass and COM
            current.mass += bodies.mass[i];
            current.com.x = (current.com.x * (current.mass - bodies.mass[i]) +
                             bodies.x[i] * bodies.mass[i]) /
                            current.mass;
            current.com.y = (current.com.y * (current.mass - bodies.mass[i]) +
                             bodies.y[i] * bodies.mass[i]) /
                            current.mass;

            // Find child quadrant for current body
            const float cx = current.boundary.x + current.boundary.width / 2;
            const float cy = current.boundary.y + current.boundary.height / 2;
            int quadrant = (bodies.x[i] > cx) | ((bodies.y[i] > cy) << 1);
            currentIdx = current.children[quadrant];

            if (++depth > 30)
            { // Prevent infinite loops
                break;
            }
        }
    }
}

constexpr int GRID_SIZE = 20;
constexpr int MAX_BODIES_PER_CELL = 1024;

const float WORLD_WIDTH = maxX - minX;
const float WORLD_HEIGHT = maxY - minY;
const float CELL_SIZE_X = WORLD_WIDTH / GRID_SIZE;
const float CELL_SIZE_Y = WORLD_HEIGHT / GRID_SIZE;
struct SpatialGrid
{
    int count[GRID_SIZE][GRID_SIZE] = {0};
    int indices[GRID_SIZE][GRID_SIZE][MAX_BODIES_PER_CELL];
};

void UpdateSpatialGrid(BodyData &bodies, SpatialGrid &grid)
{
    memset(grid.count, 0, sizeof(grid.count));

    for (int i = 0; i < bodies.count; ++i)
    {
        // Convert world position to grid coordinates
        int x = static_cast<int>((bodies.x[i] - minX) / CELL_SIZE_X);
        int y = static_cast<int>((bodies.y[i] - minY) / CELL_SIZE_Y);

        // Clamp to grid boundaries
        x = Clamp(x, 0, GRID_SIZE - 1);
        y = Clamp(y, 0, GRID_SIZE - 1);

        if (grid.count[x][y] < MAX_BODIES_PER_CELL)
        {
            grid.indices[x][y][grid.count[x][y]++] = i;
        }
    }
}

void ResolveCollisions(BodyData &bodies, SpatialGrid &grid)
{
    UpdateSpatialGrid(bodies, grid);

    for (int cx = 0; cx < GRID_SIZE; ++cx)
    {
        for (int cy = 0; cy < GRID_SIZE; ++cy)
        {

            int cellIndex = cy * GRID_SIZE + cx;
            int *cell = grid.indices[cx][cy];
            int cellCount = std::min(grid.count[cx][cy], MAX_BODIES_PER_CELL); // Safety check

            // Check current cell
            for (int i = 0; i < cellCount; ++i)
            {
                int a = cell[i];
                if (a >= bodies.count)
                    continue; // Boundary check

                for (int j = i + 1; j < cellCount; ++j)
                {
                    int b = cell[j];
                    if (b >= bodies.count)
                        continue; // Boundary check

                    float dx = bodies.x[b] - bodies.x[a];
                    float dy = bodies.y[b] - bodies.y[a];
                    float distSq = dx * dx + dy * dy;
                    float minDist = bodies.radius[a] + bodies.radius[b];

                    if (distSq < minDist * minDist)
                    {
                        // Collision response
                        float dist = sqrtf(distSq);
                        float nx = dx / dist;
                        float ny = dy / dist;

                        // Relative velocity
                        float rvx = bodies.vx[b] - bodies.vx[a];
                        float rvy = bodies.vy[b] - bodies.vy[a];

                        // Impulse calculation
                        float impulse = (1.0f + COR) * (rvx * nx + rvy * ny);
                        impulse /= (1.0f / bodies.mass[a] + 1.0f / bodies.mass[b]);

                        // Apply impulses
                        bodies.vx[a] += impulse * nx / bodies.mass[a];
                        bodies.vy[a] += impulse * ny / bodies.mass[a];
                        bodies.vx[b] -= impulse * nx / bodies.mass[b];
                        bodies.vy[b] -= impulse * ny / bodies.mass[b];

                        // Position correction
                        float penetration = minDist - dist;
                        bodies.x[a] -= nx * penetration * 0.5f;
                        bodies.y[a] -= ny * penetration * 0.5f;
                        bodies.x[b] += nx * penetration * 0.5f;
                        bodies.y[b] += ny * penetration * 0.5f;

                        // Damping
                        bodies.vx[a] *= DAMPING;
                        bodies.vy[a] *= DAMPING;
                        bodies.vx[b] *= DAMPING;
                        bodies.vy[b] *= DAMPING;
                    }
                }
            }

            // Check adjacent cells with boundary validation
            const int neighbors[][2] = {{1, 0}, {0, 1}, {1, 1}, {1, -1}};
            for (auto [dx, dy] : neighbors)
            {
                int nx = cx + dx;
                int ny = cy + dy;

                if (nx >= 0 && nx < GRID_SIZE && ny >= 0 && ny < GRID_SIZE)
                {

                    int neighborCellIndex = ny * GRID_SIZE + nx;
                    int *neighborCell = grid.indices[nx][ny];
                    int neighborCellCount = std::min(grid.count[nx][ny], MAX_BODIES_PER_CELL);

                    for (int i = 0; i < cellCount; ++i)
                    {
                        int a = cell[i];
                        if (a >= bodies.count)
                            continue;

                        for (int j = 0; j < neighborCellCount; ++j)
                        {
                            int b = neighborCell[j];
                            if (b >= bodies.count)
                                continue;

                            float dx = bodies.x[b] - bodies.x[a];
                            float dy = bodies.y[b] - bodies.y[a];
                            float distSq = dx * dx + dy * dy;
                            float minDist = bodies.radius[a] + bodies.radius[b];

                            if (distSq < minDist * minDist)
                            {
                                float dist = sqrtf(distSq);
                                float nx = dx / dist;
                                float ny = dy / dist;

                                // Relative velocity
                                float rvx = bodies.vx[b] - bodies.vx[a];
                                float rvy = bodies.vy[b] - bodies.vy[a];

                                // Impulse calculation
                                float impulse = (1.0f + COR) * (rvx * nx + rvy * ny);
                                impulse /= (1.0f / bodies.mass[a] + 1.0f / bodies.mass[b]);

                                // Apply impulses
                                bodies.vx[a] += impulse * nx / bodies.mass[a];
                                bodies.vy[a] += impulse * ny / bodies.mass[a];
                                bodies.vx[b] -= impulse * nx / bodies.mass[b];
                                bodies.vy[b] -= impulse * ny / bodies.mass[b];

                                // Position correction
                                float penetration = minDist - dist;
                                bodies.x[a] -= nx * penetration * 0.5f;
                                bodies.y[a] -= ny * penetration * 0.5f;
                                bodies.x[b] += nx * penetration * 0.5f;
                                bodies.y[b] += ny * penetration * 0.5f;
                            }
                        }
                    }
                }
            }
        }
    }
}

void Draw2DGrid()
{
    for (float y = minY; y <= maxY; y += CELL_SIZE_Y)
    {
        if (y == minY || y == maxY)
        {
            DrawLine(minX, y, maxX, y, GRID_OUTLINE_COLOR);
            continue;
        }
        DrawLine(minX, y, maxX, y, GRID_LINE_COLOR);
    }

    for (float x = minX; x <= maxX; x += CELL_SIZE_X)
    {
        if (x == minX or x == maxX)
        {
            DrawLine(x, minY, x, maxY, GRID_OUTLINE_COLOR);
            continue;
        }
        DrawLine(x, minY, x, maxY, GRID_LINE_COLOR);
    }
}

int main()
{
    SetWindowState(FLAG_MSAA_4X_HINT);
    InitWindow(SCREEN_SIZE, SCREEN_SIZE, "Optimized Barnes-Hut");
    SetWindowState(FLAG_WINDOW_RESIZABLE);
    SetTargetFPS(60);

    ChangeDirectory(GetApplicationDirectory());

    Texture2D particleTex = LoadTexture("Assets/particle.png");
    SetTextureWrap(particleTex, TEXTURE_WRAP_CLAMP);

    camera.target = {0, 0};
    camera.offset = {SCREEN_SIZE / 2, SCREEN_SIZE / 2};
    camera.zoom = 0.3f;

    BodyData bodies;
    SpatialGrid spatialGrid;

    // Initialize bodies
    if (spawn_on_load)
    {
        for (int i = 0; i < 1500; ++i)
        {
            bodies.x[i] = GetRandomValue(minX + 500, maxX - 500);
            bodies.y[i] = GetRandomValue(minY + 500, maxY - 500);
            bodies.vx[i] = GetRandomValue(-50, 50);
            bodies.vy[i] = GetRandomValue(-50, 50);
            bodies.mass[i] = 5;
            bodies.radius[i] = sqrtf(bodies.mass[i]) * 10.0f;
            bodies.color[i] = BLACK;
            bodies.kinEnergy[i] = 0.0f;
            bodies.potEnergy[i] = 0.0f;
            bodies.count++;
        }
    }

    float mouse_down_time = 0.0f;
    Vector2 mouse_init_pos = Vector2Zero();

    while (!WindowShouldClose())
    {
        if (GetMouseWheelMove() != 0)
        {

            Vector2 mouseWorldPos = GetScreenToWorld2D(GetMousePosition(), camera);

            camera.offset = GetMousePosition();

            camera.target = mouseWorldPos;

            // Update zoom level
            camera.zoom += GetMouseWheelMove() / 20.0f;

            float scaleFactor = 1.0f + (0.25f * fabsf(GetMouseWheelMove()));

            if (GetMouseWheelMove() < 0)
                scaleFactor = 1.0f / scaleFactor;

            camera.zoom = Clamp(camera.zoom * scaleFactor, 0.0825f, 64.0f);
        }

        if (IsKeyPressed(KEY_LEFT_BRACKET))
        {
            G = G / 2.0f;
        }

        if (IsKeyPressed(KEY_RIGHT_BRACKET))
        {
            G *= 2.0f;
        }

        if (IsKeyPressed(KEY_P))
        {
            (spawn_on_load) ? spawn_on_load = false : spawn_on_load = true;
        }

        if (IsKeyPressed(KEY_R))
        {

            for (int i = 0; i < MAX_BODIES; ++i)
            {
                bodies.x[i] = 0;
                bodies.y[i] = 0;
                bodies.vx[i] = 0;
                bodies.vy[i] = 0;
                bodies.mass[i] = 0;
                bodies.radius[i] = 0;
                bodies.color[i] = BLACK;
                bodies.kinEnergy[i] = 0.0f;
                bodies.potEnergy[i] = 0.0f;
                bodies.count = 0;
            }

            if (spawn_on_load == true)
            {
                for (int i = 0; i < 1500; ++i)
                {
                    bodies.x[i] = GetRandomValue(minX + 500, maxX - 500);
                    bodies.y[i] = GetRandomValue(minY + 500, maxY - 500);
                    bodies.vx[i] = GetRandomValue(-50, 50);
                    bodies.vy[i] = GetRandomValue(-50, 50);
                    bodies.mass[i] = 5;
                    bodies.radius[i] = sqrtf(bodies.mass[i]) * 10.0f;
                    bodies.color[i] = BLACK;
                    bodies.kinEnergy[i] = 0.0f;
                    bodies.potEnergy[i] = 0.0f;
                    bodies.count++;
                }
            }
        }

        if (IsMouseButtonPressed(MOUSE_LEFT_BUTTON))
        {
            mouse_down_time = 0.0f;
            mouse_init_pos = GetScreenToWorld2D(GetMousePosition(), camera);
        }

        if (IsMouseButtonDown(MOUSE_LEFT_BUTTON))
        {
            mouse_down_time += GetFrameTime();
            mouse_down_time = Clamp(mouse_down_time, 0.0f, 5.0f);
        }

        if (IsMouseButtonReleased(MOUSE_LEFT_BUTTON))
        {
            Vector2 mousePos = GetScreenToWorld2D(GetMousePosition(), camera);

            if (bodies.count < MAX_BODIES)
            {
                bodies.x[bodies.count] = mouse_init_pos.x;
                bodies.y[bodies.count] = mouse_init_pos.y;
                bodies.vx[bodies.count] = Vector2Subtract(mousePos, mouse_init_pos).x * 5;
                bodies.vy[bodies.count] = Vector2Subtract(mousePos, mouse_init_pos).y * 5;
                bodies.mass[bodies.count] = 1 + pow(2.4, mouse_down_time);
                bodies.radius[bodies.count] = sqrtf(bodies.mass[bodies.count]) * 10.0f;
                bodies.color[bodies.count] = BLACK;
                bodies.kinEnergy[bodies.count] = 0.0f;
                bodies.potEnergy[bodies.count] = 0.0f;
                bodies.count++;
            }
        }

        if (IsKeyPressed(KEY_C))
        {
            (collisionsEnabled) ? collisionsEnabled = false : collisionsEnabled = true;
        }

        if (IsMouseButtonDown(MOUSE_BUTTON_RIGHT))
        {
            Vector2 delta = GetMouseDelta();
            delta = Vector2Scale(delta, -1.0f / camera.zoom);
            camera.target = Vector2Add(camera.target, delta);
        }

        G = Clamp(G, 0.0f, 409600.0f);

        camera.target.x += (IsKeyDown(KEY_RIGHT) - IsKeyDown(KEY_LEFT)) * 30;
        camera.target.y += (IsKeyDown(KEY_DOWN) - IsKeyDown(KEY_UP)) * 30;

        // Physics update
        BuildTree(bodies);

        // Parallel force computation
        const int numThreads = std::thread::hardware_concurrency();
        std::vector<std::thread> threads;
        int chunk = bodies.count / numThreads;
        for (int t = 0; t < numThreads; t++)
        {
            int start = t * chunk;
            int end = (t == numThreads - 1) ? bodies.count : start + chunk;
            threads.emplace_back(ComputeForces, std::ref(bodies), start, end, 0);
        }
        for (auto &t : threads)
            t.join();

        // Update particles
        for (int i = 0; i < bodies.count; ++i)
        {
            bodies.x[i] += bodies.vx[i] * DT;
            bodies.y[i] += bodies.vy[i] * DT;

            // Boundary checks
            if (bodies.x[i] < minX + bodies.radius[i])
                bodies.vx[i] *= -COR;
            if (bodies.x[i] > maxX - bodies.radius[i])
                bodies.vx[i] *= -COR;
            if (bodies.y[i] < minY + bodies.radius[i] || bodies.y[i] > maxY - bodies.radius[i])
                bodies.vy[i] *= -COR;

            bodies.x[i] = Clamp(bodies.x[i], minX + bodies.radius[i], maxX - bodies.radius[i]);
            bodies.y[i] = Clamp(bodies.y[i], minY + bodies.radius[i], maxY - bodies.radius[i]);

            // Update colors

            float totalEnergy = bodies.kinEnergy[i] * 0.01 - bodies.potEnergy[i] * 1000000000; // Note the minus sign!

            float normalizedEnergy = (totalEnergy * 0.18f) + // 1800/G where G=10000 -> 0.18
                                     (bodies.kinEnergy[i] / 150000000000.0f);

            normalizedEnergy = Clamp(normalizedEnergy * 0.01f, 0.0f, 1.0f); // Adjusted scaling

            Color gradient[] = {
                {0, 12, 255, 255},    // Deep blue (cold, low energy)
                {0, 120, 255, 255},   // Medium blue
                {0, 255, 255, 255},   // Cyan
                {50, 255, 150, 255},  // Blue-green
                {100, 255, 100, 255}, // Green
                {255, 255, 0, 255},   // Yellow
                {255, 120, 0, 255},   // Orange
                {255, 12, 0, 255}     // Red (hot, high energy)
            };
            const int gradientStops = sizeof(gradient) / sizeof(gradient[0]);

            float position = normalizedEnergy * (gradientStops - 1);
            int idx1 = Clamp((int)position, 0, gradientStops - 2);
            int idx2 = idx1 + 1;
            float t = position - idx1;

            t = t * t * (3.0f - 2.0f * t); // Smoothstep function

            Color finalColor;
            finalColor.r = gradient[idx1].r + t * (gradient[idx2].r - gradient[idx1].r);
            finalColor.g = gradient[idx1].g + t * (gradient[idx2].g - gradient[idx1].g);
            finalColor.b = gradient[idx1].b + t * (gradient[idx2].b - gradient[idx1].b);
            finalColor.a = 255;

            const float gamma = 2.2f;
            finalColor.r = powf(finalColor.r / 255.0f, gamma) * 255;
            finalColor.g = powf(finalColor.g / 255.0f, gamma) * 255;
            finalColor.b = powf(finalColor.b / 255.0f, gamma) * 255;

            bodies.color[i] = finalColor;
        }

        // Handle collisions
        if (collisionsEnabled)
        {
            ResolveCollisions(bodies, spatialGrid);
        }

        // Rendering
        BeginDrawing();
        ClearBackground(BG_COLOR);

        BeginMode2D(camera);

        Draw2DGrid();

        for (int i = 0; i < bodies.count; i++)
        {
            DrawTextureEx(particleTex, {bodies.x[i] - (particleTex.width * (2.0f * bodies.radius[i]) / particleTex.width) / 2.0f, bodies.y[i] - (particleTex.height * (2.0f * bodies.radius[i]) / particleTex.width) / 2.0f}, 0, (2.0f * bodies.radius[i]) / particleTex.width, bodies.color[i]);
        }

        if (IsMouseButtonDown(MOUSE_LEFT_BUTTON))
        {
            DrawTextureEx(particleTex, {mouse_init_pos.x - (particleTex.width * (2.0f * (sqrtf(1 + pow(2, mouse_down_time)) * 10.0f) / particleTex.width) / 2.0f), mouse_init_pos.y - (particleTex.height * (2.0f * (sqrtf(1 + pow(2.4, mouse_down_time)) * 10.0f) / particleTex.width) / 2.0f)}, 0, (2.0f * sqrtf(1 + pow(2.4, mouse_down_time)) * 10.0f) / particleTex.width, BLACK);
        }

        EndMode2D();

        DrawText("N-Body Simulation |", 20, 10, 20, TEXT_COLOR);
        DrawFPS(220, 10);

        DrawText("Camera Position:", 20, 50, 20, TEXT_COLOR);
        DrawText(TextFormat("(%.0f, %.0f)", camera.target.x, camera.target.y), 250, 50, 20, VALUE_COLOR);
        DrawText("Hold", 20, 70, 20, HINT_TEXT_COLOR);
        DrawText("Right Click", 70, 70, 20, INPUT_TEXT_COLOR);
        DrawText("or use the", 185, 70, 20, HINT_TEXT_COLOR);
        DrawText("Arrow Keys", 300, 70, 20, INPUT_TEXT_COLOR);
        DrawText("to pan the camera", 425, 70, 20, HINT_TEXT_COLOR);

        DrawText("Camera Zoom:", 20, 110, 20, TEXT_COLOR);
        DrawText(TextFormat("%.2f", camera.zoom), 250, 110, 20, VALUE_COLOR);
        DrawText("Use the", 20, 130, 20, HINT_TEXT_COLOR);
        DrawText("Mouse Wheel", 108, 130, 20, INPUT_TEXT_COLOR);
        DrawText("to zoom in/out", 245, 130, 20, HINT_TEXT_COLOR);

        DrawText("Particles:", 20, 170, 20, TEXT_COLOR);
        DrawText(TextFormat("%i", bodies.count), 250, 170, 20, VALUE_COLOR);
        DrawText("Left Click", 20, 190, 20, INPUT_TEXT_COLOR);
        DrawText("to spawn particles (Hold to set velocity, long hold to set mass)", 130, 190, 20, HINT_TEXT_COLOR);

        DrawText("Strength of Gravity:", 20, 230, 20, TEXT_COLOR);
        DrawText(TextFormat("%.1f", G), 250, 230, 20, VALUE_COLOR);
        DrawText("Press", 20, 250, 20, HINT_TEXT_COLOR);
        DrawText("]", 90, 250, 20, INPUT_TEXT_COLOR);
        DrawText("to double Gravity", 105, 250, 20, HINT_TEXT_COLOR);
        DrawText("Press", 20, 270, 20, HINT_TEXT_COLOR);
        DrawText("[", 90, 270, 20, INPUT_TEXT_COLOR);
        DrawText("to halve Gravity", 105, 270, 20, HINT_TEXT_COLOR);

        DrawText("Press", 20, 310, 20, HINT_TEXT_COLOR);
        DrawText("R", 90, 310, 20, INPUT_TEXT_COLOR);
        DrawText("to reset the simulation.", 110, 310, 20, HINT_TEXT_COLOR);

        DrawText("Spawn bundle of particles", 20, 350, 20, TEXT_COLOR);
        DrawText("on start/restart:", 20, 370, 20, TEXT_COLOR);
        if (spawn_on_load == true)
        {
            DrawText("ENABLED", 250, 370, 20, GREEN);
        }
        else
        {
            DrawText("DISABLED", 250, 370, 20, RED);
        }
        DrawText("Press", 20, 390, 20, HINT_TEXT_COLOR);
        DrawText("P", 90, 390, 20, INPUT_TEXT_COLOR);
        DrawText("to toggle spawn_on_load", 110, 390, 20, HINT_TEXT_COLOR);

        DrawText("Collisions enabled:", 20, 430, 20, TEXT_COLOR);
        if (collisionsEnabled == true)
        {
            DrawText("ENABLED", 250, 430, 20, GREEN);
        }
        else
        {
            DrawText("DISABLED", 250, 430, 20, RED);
        }
        DrawText("Press", 20, 450, 20, HINT_TEXT_COLOR);
        DrawText("C", 90, 450, 20, INPUT_TEXT_COLOR);
        DrawText("to toggle collisions", 110, 450, 20, HINT_TEXT_COLOR);

        EndDrawing();
    }

    UnloadTexture(particleTex);
    CloseWindow();

    return 0;
}