main.cpp

#include <iostream>
#include <vector>
#include <cmath>
#include <cstdlib>
#include <cstdio>
#include <GLFW/glfw3.h>

#include "imgui.h"
#include "imgui_impl_glfw.h"
#include "imgui_impl_opengl2.h"

struct Particle {
    float x, y;
    float vx, vy;
};

// Global Simulation State
float wallMomentumSum = 0.0f;
float timeAccumulator = 0.0f;
float gravity = 0.0f;          
float restitution = 1.0f;      
float globalRadius = 0.02f;
float particleMass = 1.0f;
float timeScale = 1.0f;
float boxWidth = 1.0f;
float boxHeight = 1.0f;
float tempScaleFactor = 18750000.0f;

// Colors
float particleColor[3] = {1.0f, 0.4f, 0.2f};
float bgColor[3] = {0.1f, 0.15f, 0.2f};

int main() {
    if (!glfwInit()) return -1;

    GLFWwindow* window = glfwCreateWindow(1024, 768, "Thermodynamics Sandbox", NULL, NULL);
    if (!window) {
        glfwTerminate();
        return -1;
    }
    glfwMakeContextCurrent(window);

    IMGUI_CHECKVERSION();
    ImGui::CreateContext();
    ImGuiIO& io = ImGui::GetIO(); (void)io;
    ImGui::StyleColorsDark();
    ImGui_ImplGlfw_InitForOpenGL(window, true);
    ImGui_ImplOpenGL2_Init();

    std::vector<Particle> particles;
    for (int i = 0; i < 50; ++i) {
        float x = ((rand() % 1000) / 500.0f) * 0.8f;
        float y = ((rand() % 1000) / 500.0f) * 0.8f;
        float vx = (((rand() % 1000) / 500.0f) - 1.0f) * 0.004f;
        float vy = (((rand() % 1000) / 500.0f) - 1.0f) * 0.004f;
        particles.push_back({x, y, vx, vy});
    }

    double lastTime = glfwGetTime();
    float currentPressure = 0.0f;
    bool lastMouseState = false;

    int particleCountInput = 50;
    float targetTempInput = 300.0f;
    float targetPressureInput = 10.0f;
    float targetEnergyInput = 1000.0f;

    while (!glfwWindowShouldClose(window)) {
        double currentTime = glfwGetTime();
        float rawDelta = static_cast<float>(currentTime - lastTime);
        lastTime = currentTime;
        
        float deltaTime = rawDelta * timeScale;
        timeAccumulator += rawDelta; 

        glfwPollEvents();

        ImGui_ImplOpenGL2_NewFrame();
        ImGui_ImplGlfw_NewFrame();
        ImGui::NewFrame();

        bool currentMouseState = (glfwGetMouseButton(window, GLFW_MOUSE_BUTTON_LEFT) == GLFW_PRESS);
        if (!io.WantCaptureMouse && currentMouseState && !lastMouseState) {
            double mx, my;
            int ww, wh;
            glfwGetCursorPos(window, &mx, &my);
            glfwGetWindowSize(window, &ww, &wh);
            
            float nx = (mx / (ww / 2.0f)) - 1.0f;
            float ny = 1.0f - (my / (wh / 2.0f));

            if (nx >= -boxWidth && nx <= boxWidth && ny >= -boxHeight && ny <= boxHeight) {
                for (int i = 0; i < 5; ++i) {
                    float vx = (((rand() % 1000) / 500.0f) - 1.0f) * 0.004f;
                    float vy = (((rand() % 1000) / 500.0f) - 1.0f) * 0.004f;
                    particles.push_back({nx, ny, vx, vy});
                }
                particleCountInput = static_cast<int>(particles.size());
            }
        }
        lastMouseState = currentMouseState;

        glClearColor(bgColor[0], bgColor[1], bgColor[2], 1.0f);
        glClear(GL_COLOR_BUFFER_BIT);

        wallMomentumSum = 0.0f;

        // Physics Update
        for (auto& p : particles) {
            p.vy -= gravity * deltaTime;
            p.x += p.vx * timeScale;
            p.y += p.vy * timeScale;

            if (p.x > boxWidth - globalRadius) p.x = boxWidth - globalRadius;
            if (p.x < -boxWidth + globalRadius) p.x = -boxWidth + globalRadius;
            if (p.y > boxHeight - globalRadius) p.y = boxHeight - globalRadius;
            if (p.y < -boxHeight + globalRadius) p.y = -boxHeight + globalRadius;

            if (p.x >= boxWidth - globalRadius || p.x <= -boxWidth + globalRadius) {
                p.vx *= -restitution;
                wallMomentumSum += std::abs(2.0f * particleMass * p.vx);
            }
            if (p.y >= boxHeight - globalRadius || p.y <= -boxHeight + globalRadius) {
                p.vy *= -restitution;
                wallMomentumSum += std::abs(2.0f * particleMass * p.vy);
            }
        }

        // Collisions
        for (size_t i = 0; i < particles.size(); ++i) {
            for (size_t j = i + 1; j < particles.size(); ++j) {
                float dx = particles[j].x - particles[i].x;
                float dy = particles[j].y - particles[i].y;
                float dist = std::sqrt(dx * dx + dy * dy);
                float minDist = globalRadius * 2.0f;

                if (dist < minDist) {
                    if (dist == 0.0f) { particles[j].x += 0.001f; continue; }
                    float overlap = 0.5f * (minDist - dist);
                    particles[i].x -= overlap * (dx / dist);
                    particles[i].y -= overlap * (dy / dist);
                    particles[j].x += overlap * (dx / dist);
                    particles[j].y += overlap * (dy / dist);

                    float nx = dx / dist;
                    float ny = dy / dist;
                    float kx = particles[i].vx - particles[j].vx;
                    float ky = particles[i].vy - particles[j].vy;
                    float p_impulse = (nx * kx + ny * ky);

                    particles[i].vx -= p_impulse * nx;
                    particles[i].vy -= p_impulse * ny;
                    particles[j].vx += p_impulse * nx;
                    particles[j].vy += p_impulse * ny;
                }
            }
        }

        if (timeAccumulator >= 0.5f) {
            currentPressure = (wallMomentumSum / timeAccumulator) * 100.0f;
            timeAccumulator = 0.0f;
        }

        // Calculations
        float totalKE = 0.0f;
        float totalSpeedSq = 0.0f;
        for (const auto& p : particles) {
            float vSq = p.vx * p.vx + p.vy * p.vy;
            totalKE += 0.5f * particleMass * vSq;
            totalSpeedSq += vSq;
        }

        size_t N = particles.size();
        float avgKE = N == 0 ? 0.0f : totalKE / N;
        float temperature = avgKE * tempScaleFactor; 
        
        float internalEnergy = totalKE * tempScaleFactor; 
        float vRMS = N == 0 ? 0.0f : std::sqrt(totalSpeedSq / N) * 1000.0f;
        float containerArea = (2.0f * boxWidth) * (2.0f * boxHeight);
        float meanFreePath = (N == 0 || globalRadius == 0.0f) ? 0.0f : 
                             containerArea / (2.828f * (2.0f * globalRadius) * N);
        float idealStateRatio = (temperature == 0.0f || N == 0) ? 0.0f : 
                                (currentPressure * containerArea) / (N * temperature);

        // Draw Container Box
        glColor3f(0.3f, 0.5f, 0.7f);
        glBegin(GL_LINE_LOOP);
        glVertex2f(-boxWidth, -boxHeight);
        glVertex2f(boxWidth, -boxHeight);
        glVertex2f(boxWidth, boxHeight);
        glVertex2f(-boxWidth, boxHeight);
        glEnd();

        // Draw Particles
        glPointSize(globalRadius * 300.0f);
        glBegin(GL_POINTS);
        glColor3f(particleColor[0], particleColor[1], particleColor[2]);
        for (const auto& p : particles) {
            glVertex2f(p.x, p.y);
        }
        glEnd();

        // Expanded GUI
        ImGui::SetNextWindowSize(ImVec2(360, 720), ImGuiCond_FirstUseEver);
        ImGui::Begin("Thermodynamics Dashboard");
        
        ImGui::TextColored(ImVec4(0.2f, 0.8f, 1.0f, 1.0f), "State Telemetry");
        ImGui::Text("Particle Count (N): %zu", N);
        ImGui::Text("Temperature (T): %.1f K", temperature);
        ImGui::Text("Pressure (P): %.2f P", currentPressure);
        ImGui::Text("Container Area (A): %.2f u²", containerArea);

        ImGui::Separator();
        ImGui::TextColored(ImVec4(0.2f, 0.8f, 1.0f, 1.0f), "Derived Quantities");
        ImGui::Text("Internal Energy (U): %.1f J", internalEnergy);
        ImGui::Text("Avg Kinetic Energy: %.2f e-6 J", avgKE * 1e6f);
        ImGui::Text("RMS Speed (v_rms): %.2f u/s", vRMS);
        ImGui::Text("Mean Free Path (Lambda): %.3f u", meanFreePath);
        ImGui::Text("Ideal Gas Ratio (PA/NT): %.4f", idealStateRatio);

        ImGui::Separator();
        ImGui::TextColored(ImVec4(0.2f, 0.8f, 1.0f, 1.0f), "State Adjustments");
        ImGui::PushItemWidth(100.0f);
        
        // Count Setter
        if (ImGui::InputInt("Count", &particleCountInput)) {
            if (particleCountInput < 0) particleCountInput = 0;
            int currentN = static_cast<int>(particles.size());
            if (particleCountInput > currentN) {
                for (int i = 0; i < particleCountInput - currentN; ++i) {
                    float x = ((rand() % 1000) / 500.0f) * (boxWidth * 0.8f);
                    float y = ((rand() % 1000) / 500.0f) * (boxHeight * 0.8f);
                    float vx = (((rand() % 1000) / 500.0f) - 1.0f) * 0.004f;
                    float vy = (((rand() % 1000) / 500.0f) - 1.0f) * 0.004f;
                    particles.push_back({x, y, vx, vy});
                }
            } else if (particleCountInput < currentN) {
                particles.resize(particleCountInput);
            }
        }

        // Temperature Setter
        ImGui::InputFloat("Temp(K)", &targetTempInput, 10.0f, 50.0f, "%.1f"); ImGui::SameLine();
        if (ImGui::Button("Set T") && temperature > 0.001f && targetTempInput > 0.0f) {
            float scale = std::sqrt(targetTempInput / temperature);
            for (auto& p : particles) { p.vx *= scale; p.vy *= scale; }
        }

        // Internal Energy Setter
        ImGui::InputFloat("Energy(J)", &targetEnergyInput, 100.0f, 500.0f, "%.1f"); ImGui::SameLine();
        if (ImGui::Button("Set U") && internalEnergy > 0.001f && targetEnergyInput > 0.0f) {
            float scale = std::sqrt(targetEnergyInput / internalEnergy);
            for (auto& p : particles) { p.vx *= scale; p.vy *= scale; }
        }

        // Pressure Setter
        ImGui::InputFloat("Pressure", &targetPressureInput, 1.0f, 5.0f, "%.2f"); ImGui::SameLine();
        if (ImGui::Button("Set P") && currentPressure > 0.001f && targetPressureInput > 0.0f) {
            float scale = std::sqrt(targetPressureInput / currentPressure);
            for (auto& p : particles) { p.vx *= scale; p.vy *= scale; }
        }

        ImGui::Separator();
        ImGui::TextColored(ImVec4(0.2f, 0.8f, 1.0f, 1.0f), "Environment & Container");
        ImGui::SliderFloat("Time Scale", &timeScale, 0.0f, 5.0f, "%.2fx");
        ImGui::SliderFloat("Box Width", &boxWidth, 0.1f, 1.0f, "%.2f");
        ImGui::SliderFloat("Box Height", &boxHeight, 0.1f, 1.0f, "%.2f");

        ImGui::Separator();
        ImGui::TextColored(ImVec4(0.2f, 0.8f, 1.0f, 1.0f), "Physics Constants");
        ImGui::InputFloat("Mass", &particleMass, 0.1f, 1.0f, "%.2f");
        ImGui::InputFloat("Radius", &globalRadius, 0.005f, 0.01f, "%.3f");
        ImGui::InputFloat("Gravity", &gravity, 0.01f, 0.1f, "%.3f");
        ImGui::InputFloat("Restitution", &restitution, 0.05f, 0.1f, "%.2f");
        ImGui::InputFloat("Temp Scalar", &tempScaleFactor, 100000.0f, 1000000.0f, "%.0f");

        ImGui::Separator();
        ImGui::TextColored(ImVec4(0.2f, 0.8f, 1.0f, 1.0f), "Visuals");
        ImGui::ColorEdit3("Particle Color", particleColor, ImGuiColorEditFlags_NoInputs);
        ImGui::ColorEdit3("Background", bgColor, ImGuiColorEditFlags_NoInputs);
        ImGui::PopItemWidth();

        if (ImGui::Button("Reset All", ImVec2(100, 30))) {
            particles.clear();
            particleCountInput = 0;
            boxWidth = 1.0f;
            boxHeight = 1.0f;
            timeScale = 1.0f;
            particleMass = 1.0f;
        }
        ImGui::End();

        ImGui::Render();
        ImGui_ImplOpenGL2_RenderDrawData(ImGui::GetDrawData());

        glfwSwapBuffers(window);
    }

    ImGui_ImplOpenGL2_Shutdown();
    ImGui_ImplGlfw_Shutdown();
    ImGui::DestroyContext();
    glfwTerminate();
    return 0;
}