//main.cpp
#include <glad/glad.h>
#include <GLFW/glfw3.h>
#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/gtc/type_ptr.hpp>

#include <vector>
#include <cmath>
#include <iostream>
#include <random>
#include <algorithm>
#include <omp.h> // OpenMP

// ==========================================
// 配置参数 (Configuration)
// ==========================================
const unsigned int SCR_WIDTH = 1280;
const unsigned int SCR_HEIGHT = 720;

// ==========================================
// 流体参数 (Tait Equation Setup)
// ==========================================

const int PARTICLE_COUNT = 4000;
const float SMOOTHING_RADIUS = 0.2f;
const float MASS = 0.1f;                 // 质量
const float REST_DENSITY = 120.0f;     

// --- Tait方程参数 ---
const float TAIT_EXPONENT = 7.0f;        // 水的典型指数，决定了"硬度"的爆发力
const float STIFFNESS = 800.0f;          // 基础刚度 (B)

// --- 动态参数 ---
const float VISCOSITY = 0.5f;            // 粘度
const float TIME_STEP = 0.0025f;         // 步长
const int SUB_STEPS = 8;                 // 子步数

const glm::vec3 GRAVITY(0.0f, -9.8f, 0.0f);
const float PARTICLE_RADIUS = 0.04f;

// --- 边界反弹 ---
const float BOUND_DAMPING = -0.8f;       // 保留80%动能，制造飞溅
const float WALL_FRICTION = 0.95f;       // 墙壁摩擦

// 容器边界
const glm::vec3 CONTAINER_SIZE(4.0f, 6.0f, 4.0f); // 宽, 高, 深 (以原点为中心)

constexpr float PI_F = 3.14159265358979323846f;
// ==========================================
//  辅助类: 摄像机 (Camera)
// ==========================================
class Camera {
public:
    glm::vec3 Position;
    glm::vec3 Front;
    glm::vec3 Up;
    glm::vec3 Right;
    glm::vec3 WorldUp;
    float Yaw = -90.0f;
    float Pitch = 0.0f;
    float MovementSpeed = 2.5f;
    float MouseSensitivity = 0.1f;

    Camera(glm::vec3 position = glm::vec3(0.0f, 0.0f, 8.0f)) 
        : Position(position), Front(0.0f, 0.0f, -1.0f), WorldUp(0.0f, 1.0f, 0.0f) {
        updateCameraVectors();
    }

    glm::mat4 GetViewMatrix() const {
        return glm::lookAt(Position, Position + Front, Up);
    }

    void ProcessKeyboard(int direction, float deltaTime) {
        float velocity = MovementSpeed * deltaTime;
        if (direction == GLFW_KEY_W) Position += Front * velocity;
        if (direction == GLFW_KEY_S) Position -= Front * velocity;
        if (direction == GLFW_KEY_A) Position -= Right * velocity;
        if (direction == GLFW_KEY_D) Position += Right * velocity;
    }

    void ProcessMouseMovement(float xoffset, float yoffset) {
        xoffset *= MouseSensitivity;
        yoffset *= MouseSensitivity;
        Yaw += xoffset;
        Pitch += yoffset;
        if (Pitch > 89.0f) Pitch = 89.0f;
        if (Pitch < -89.0f) Pitch = -89.0f;
        updateCameraVectors();
    }

private:
    void updateCameraVectors() {
        glm::vec3 front;
        front.x = cos(glm::radians(Yaw)) * cos(glm::radians(Pitch));
        front.y = sin(glm::radians(Pitch));
        front.z = sin(glm::radians(Yaw)) * cos(glm::radians(Pitch));
        Front = glm::normalize(front);
        Right = glm::normalize(glm::cross(Front, WorldUp));
        Up = glm::normalize(glm::cross(Right, Front));
    }
};

// ==========================================
// SPH 求解器 (Physics System)
// ==========================================
class SPHSystem {
public:
    // 使用 SoA (Structure of Arrays) 布局优化缓存
    std::vector<glm::vec3> positions;
    std::vector<glm::vec3> velocities;
    std::vector<glm::vec3> forces;
    std::vector<float> densities;
    std::vector<float> pressures;

    // 空间哈希网格优化
    // key: 网格索引 (x + y*width + z*width*depth), value: 粒子索引列表
    std::vector<std::vector<int>> grid; 
    int gridResX, gridResY, gridResZ;
    glm::vec3 gridMin, gridMax;
    float cellSize;

    SPHSystem() {
        positions.resize(PARTICLE_COUNT);
        velocities.resize(PARTICLE_COUNT);
        forces.resize(PARTICLE_COUNT);
        densities.resize(PARTICLE_COUNT);
        pressures.resize(PARTICLE_COUNT);

        // 初始化网格参数
        cellSize = SMOOTHING_RADIUS;
        gridMin = -CONTAINER_SIZE / 2.0f - glm::vec3(1.0f);
        gridMax = CONTAINER_SIZE / 2.0f + glm::vec3(1.0f);
        gridResX = (int)ceil((gridMax.x - gridMin.x) / cellSize);
        gridResY = (int)ceil((gridMax.y - gridMin.y) / cellSize);
        gridResZ = (int)ceil((gridMax.z - gridMin.z) / cellSize);
        grid.resize(gridResX * gridResY * gridResZ);

        reset();
    }

    void reset() {
        std::mt19937 rng(std::random_device{}());
        std::uniform_real_distribution<float> dist(-0.5f, 0.5f);

        // 将粒子生成在空中的一个立方体区域内
        float spawnSize = 2.0f;
        int particlesPerAxis = (int)ceil(cbrt(PARTICLE_COUNT));
        float spacing = spawnSize / particlesPerAxis;

        int idx = 0;
        for(int x=0; x<particlesPerAxis && idx < PARTICLE_COUNT; ++x) {
            for(int y=0; y<particlesPerAxis && idx < PARTICLE_COUNT; ++y) {
                for(int z=0; z<particlesPerAxis && idx < PARTICLE_COUNT; ++z) {
                    // 添加一些随机抖动避免由于完美对齐导致的数值不稳定
                    float jitterX = dist(rng) * 0.01f;
                    float jitterY = dist(rng) * 0.01f;
                    float jitterZ = dist(rng) * 0.01f;

                    // 初始位置在容器上方
                    positions[idx] = glm::vec3(
                        (x * spacing) - spawnSize/2.0f + jitterX,
                        (y * spacing) + 2.0f + jitterY, 
                        (z * spacing) - spawnSize/2.0f + jitterZ
                    );
                    // velocities[idx] = glm::vec3(0.0f);
                    // 在生成位置后，给一个随机微小初速度
                    velocities[idx] = glm::vec3(
                        dist(rng) * 2.0f, // X方向有点初速度
                        dist(rng) * 2.0f, 
                        dist(rng) * 2.0f
                    );
                    idx++;
                }
            }
        }
    }

    // SPH 核函数: Poly6 (用于密度计算)
    float poly6Kernel(float r2) {
        float h2 = SMOOTHING_RADIUS * SMOOTHING_RADIUS;
        if (r2 < 0 || r2 > h2) return 0.0f;
        float diff = h2 - r2;
        return (315.0f / (64.0f * PI_F * pow(SMOOTHING_RADIUS, 9))) * diff * diff * diff;
    }

    // SPH 核函数梯度: Spiky (用于压力计算)
    glm::vec3 spikyKernelGradient(glm::vec3 r, float dist) {
        if (dist <= 0 || dist > SMOOTHING_RADIUS) return glm::vec3(0.0f);
        float t = SMOOTHING_RADIUS - dist;
        float scalar = -45.0f / (PI_F * pow(SMOOTHING_RADIUS, 6)) * t * t;
        return glm::normalize(r) * scalar;
    }

    // SPH 核函数拉普拉斯: Viscosity (用于粘性计算)
    float viscosityKernelLaplacian(float dist) {
        if (dist <= 0 || dist > SMOOTHING_RADIUS) return 0.0f;
        return 45.0f / (PI_F * pow(SMOOTHING_RADIUS, 6)) * (SMOOTHING_RADIUS - dist);
    }

    // 计算网格索引
    int getGridIndex(const glm::vec3& pos) {
        int x = (int)((pos.x - gridMin.x) / cellSize);
        int y = (int)((pos.y - gridMin.y) / cellSize);
        int z = (int)((pos.z - gridMin.z) / cellSize);

        // 边界钳制，防止越界
        x = std::clamp(x, 0, gridResX - 1);
        y = std::clamp(y, 0, gridResY - 1);
        z = std::clamp(z, 0, gridResZ - 1);

        return x + y * gridResX + z * gridResX * gridResY;
    }

    void update() {
        // 1. 构建空间网格
        // 清空网格
        for(auto& cell : grid) cell.clear();

        // 将粒子填入网格
        for (int i = 0; i < PARTICLE_COUNT; ++i) {
            grid[getGridIndex(positions[i])].push_back(i);
        }

        // 2. 计算密度和压力
        #pragma omp parallel for schedule(static)
        for (int i = 0; i < PARTICLE_COUNT; ++i) {
            densities[i] = 0.0f;

            // 获取当前粒子所在及其周围的网格 (3x3x3)
            int cx = (int)((positions[i].x - gridMin.x) / cellSize);
            int cy = (int)((positions[i].y - gridMin.y) / cellSize);
            int cz = (int)((positions[i].z - gridMin.z) / cellSize);

            for (int dx = -1; dx <= 1; ++dx) {
                for (int dy = -1; dy <= 1; ++dy) {
                    for (int dz = -1; dz <= 1; ++dz) {
                        int nx = cx + dx;
                        int ny = cy + dy;
                        int nz = cz + dz;

                        if (nx >= 0 && nx < gridResX && ny >= 0 && ny < gridResY && nz >= 0 && nz < gridResZ) {
                            int cellIdx = nx + ny * gridResX + nz * gridResX * gridResY;
                            for (int j : grid[cellIdx]) {
                                glm::vec3 diff = positions[i] - positions[j];
                                float r2 = glm::dot(diff, diff);
                                if (r2 < SMOOTHING_RADIUS * SMOOTHING_RADIUS) {
                                    densities[i] += MASS * poly6Kernel(r2);
                                }
                            }
                        }
                    }
                }
            }
            // 状态方程截断负压
            // 如果密度低于静态密度 (REST_DENSITY)，压力设为 0，而不是负数。
            densities[i] = std::max(densities[i], 0.001f); // 防止除零

            // Tait Equation
            // P = B * ((rho / rho0)^gamma - 1)
            // 这种方程在密度接近1时压力很小，但稍微压缩一点压力就暴增
            float densityRatio = densities[i] / REST_DENSITY;

            // 截断负压：只允许 densityRatio > 1 时产生正压
            if (densityRatio > 1.0f) {
                pressures[i] = STIFFNESS * (std::pow(densityRatio, TAIT_EXPONENT) - 1.0f);
            } else {
                pressures[i] = 0.0f;
            }

        }

        // 3. 计算合力 (压力 + 粘性 + 重力)
        #pragma omp parallel for schedule(static)
        for (int i = 0; i < PARTICLE_COUNT; ++i) {
            glm::vec3 pressureForce(0.0f);
            glm::vec3 viscosityForce(0.0f);

            int cx = (int)((positions[i].x - gridMin.x) / cellSize);
            int cy = (int)((positions[i].y - gridMin.y) / cellSize);
            int cz = (int)((positions[i].z - gridMin.z) / cellSize);

            for (int dx = -1; dx <= 1; ++dx) {
                for (int dy = -1; dy <= 1; ++dy) {
                    for (int dz = -1; dz <= 1; ++dz) {
                        int nx = cx + dx; 
                        int ny = cy + dy;
                        int nz = cz + dz;

                        if (nx >= 0 && nx < gridResX && ny >= 0 && ny < gridResY && nz >= 0 && nz < gridResZ) {
                            int cellIdx = nx + ny * gridResX + nz * gridResX * gridResY;
                            for (int j : grid[cellIdx]) {
                                if (i == j) continue;

                                glm::vec3 diff = positions[i] - positions[j];
                                float dist = glm::length(diff);

                                // 最小距离钳制
                                // 防止两个粒子重合时产生 NaN 或者 Inf 的力
                                if (dist < 0.0001f) dist = 0.0001f; 

                                if (dist < SMOOTHING_RADIUS && dist > 0.0001f) {
                                    // 压力 (对称)
                                    glm::vec3 grad = spikyKernelGradient(diff, dist);
                                    pressureForce -= MASS * (pressures[i] + pressures[j]) / (2.0f * densities[j]) * grad;

                                    // 粘滞力
                                    float lap = viscosityKernelLaplacian(dist);
                                    viscosityForce += VISCOSITY * MASS * (velocities[j] - velocities[i]) / densities[j] * lap;
                                }
                            }
                        }
                    }
                }
            }

            glm::vec3 gravityForce = GRAVITY * densities[i];
            forces[i] = pressureForce + viscosityForce + gravityForce;
        }

        // 4. 积分 (Integration) 与 边界碰撞
        #pragma omp parallel for schedule(static)
        for (int i = 0; i < PARTICLE_COUNT; ++i) {
            // F = ma => a = F/rho
            glm::vec3 acceleration = forces[i] / densities[i];

            velocities[i] += acceleration * TIME_STEP;

            // 放宽速度限制,只拦截极端错误的数值爆炸
            float vLength = glm::length(velocities[i]);
            if (vLength > 20.0f) {
                velocities[i] = glm::normalize(velocities[i]) * 20.0f;
            }

            positions[i] += velocities[i] * TIME_STEP;

            // 边界检测
            // 容器半长宽 (Box Half-Size)
            glm::vec3 bounds = CONTAINER_SIZE / 2.0f;
            // 边界检测需考虑粒子半径 (r)，防止穿模
            float r = PARTICLE_RADIUS;
            bool collided = false;

            // X轴碰撞 增加切向摩擦力，让水流撞墙后产生翻滚效果
            if (positions[i].x < -bounds.x + r) {
                positions[i].x = -bounds.x + r;
                velocities[i].x *= BOUND_DAMPING;
                velocities[i].y *= WALL_FRICTION; // 切向摩擦
                velocities[i].z *= WALL_FRICTION;
            }
            else if (positions[i].x > bounds.x - r) {
                positions[i].x = bounds.x - r;
                velocities[i].x *= BOUND_DAMPING;
                velocities[i].y *= WALL_FRICTION;
                velocities[i].z *= WALL_FRICTION;
            }

            // Y轴碰撞 (地板)
            if (positions[i].y < -bounds.y + r) {
                positions[i].y = -bounds.y + r;
                velocities[i].y *= BOUND_DAMPING; // 强力反弹
                velocities[i].x *= WALL_FRICTION; // 地面摩擦
                velocities[i].z *= WALL_FRICTION;
            }

            // Z轴碰撞
            if (positions[i].z < -bounds.z + r) {
                positions[i].z = -bounds.z + r;
                velocities[i].z *= BOUND_DAMPING;
                velocities[i].x *= WALL_FRICTION;
                velocities[i].y *= WALL_FRICTION;
            }
            else if (positions[i].z > bounds.z - r) {
                positions[i].z = bounds.z - r;
                velocities[i].z *= BOUND_DAMPING;
                velocities[i].x *= WALL_FRICTION;
                velocities[i].y *= WALL_FRICTION;
            }
        }
    }
};

// ==========================================
// 渲染器 (Shader & Setup)
// ==========================================

// 顶点着色器: 支持实例化渲染 (Instancing)
const char* vertexShaderSource = R"(
#version 330 core
layout (location = 0) in vec3 aPos;      // 基础球体模型顶点
layout (location = 1) in vec3 aNormal;   // 法线
layout (location = 2) in vec3 aOffset;   // 实例偏移量 (每个粒子的位置)

out vec3 Normal;
out vec3 FragPos;

uniform mat4 view;
uniform mat4 projection;
uniform float radius;

void main()
{
    // 计算实际世界坐标：模型偏移 + 顶点缩放
    vec3 worldPos = aOffset + aPos * radius; 
    FragPos = worldPos;
    Normal = aNormal; // 简化，假设球体不旋转

    gl_Position = projection * view * vec4(worldPos, 1.0);
}
)";

// 片段着色器: 简单的蓝色水体效果
const char* fragmentShaderSource = R"(
#version 330 core
out vec4 FragColor;

in vec3 Normal;
in vec3 FragPos;

uniform vec3 viewPos;

void main()
{
    // 环境光
    float ambientStrength = 0.3;
    vec3 ambient = ambientStrength * vec3(0.2, 0.5, 1.0);

    // 漫反射
    vec3 norm = normalize(Normal);
    vec3 lightDir = normalize(vec3(0.5, 1.0, 0.5)); // 假定一个平行光
    float diff = max(dot(norm, lightDir), 0.0);
    vec3 diffuse = diff * vec3(0.0, 0.6, 1.0);

    // 镜面光
    float specularStrength = 0.8;
    vec3 viewDir = normalize(viewPos - FragPos);
    vec3 reflectDir = reflect(-lightDir, norm);  
    float spec = pow(max(dot(viewDir, reflectDir), 0.0), 32);
    vec3 specular = specularStrength * vec3(1.0, 1.0, 1.0);  

    vec3 result = ambient + diffuse + specular;
    FragColor = vec4(result, 0.85); // 略微透明
}
)";

// 构建简单球体网格
void createSphere(std::vector<float>& vertices, std::vector<unsigned int>& indices, int sectors, int stacks) {
    float x, y, z, xy;  
    float nx, ny, nz, lengthInv = 1.0f;         

    float sectorStep = 2 * PI_F / sectors;
    float stackStep = PI_F / stacks;
    float sectorAngle, stackAngle;

    for(int i = 0; i <= stacks; ++i) {
        stackAngle = PI_F / 2 - i * stackStep;
        xy = cos(stackAngle);
        z = sin(stackAngle);

        for(int j = 0; j <= sectors; ++j) {
            sectorAngle = j * sectorStep;
            x = xy * cos(sectorAngle);
            y = xy * sin(sectorAngle);

            // 顶点
            vertices.push_back(x);
            vertices.push_back(z); // 注意: 这里的生成逻辑y朝上，我们交换y/z适应习惯或Shader
            vertices.push_back(y);

            // 法线 (球体中心在原点，法线就是位置)
            vertices.push_back(x);
            vertices.push_back(z);
            vertices.push_back(y);
        }
    }

    int k1, k2;
    for(int i = 0; i < stacks; ++i) {
        k1 = i * (sectors + 1);
        k2 = k1 + sectors + 1;
        for(int j = 0; j < sectors; ++j, ++k1, ++k2) {
            if(i != 0) {
                indices.push_back(k1);
                indices.push_back(k2);
                indices.push_back(k1 + 1);
            }
            if(i != (stacks - 1)) {
                indices.push_back(k1 + 1);
                indices.push_back(k2);
                indices.push_back(k2 + 1);
            }
        }
    }
}

// 全局变量用于回调
Camera camera(glm::vec3(0.0f, 2.0f, 10.0f));
float lastX = SCR_WIDTH / 2.0f;
float lastY = SCR_HEIGHT / 2.0f;
bool firstMouse = true;
bool resetSim = false;

void mouse_callback(GLFWwindow* window, double xposIn, double yposIn) {
    float xpos = static_cast<float>(xposIn);
    float ypos = static_cast<float>(yposIn);
    if (firstMouse) { lastX = xpos; lastY = ypos; firstMouse = false; }
    float xoffset = xpos - lastX;
    float yoffset = lastY - ypos; 
    lastX = xpos; lastY = ypos;
    camera.ProcessMouseMovement(xoffset, yoffset);
}

void key_callback(GLFWwindow* window, int key, int scancode, int action, int mods) {
    if (key == GLFW_KEY_ESCAPE && action == GLFW_PRESS)
        glfwSetWindowShouldClose(window, true);
    if (key == GLFW_KEY_R && action == GLFW_PRESS)
        resetSim = true;
}

// 编译 Shader 辅助函数
unsigned int compileShader(const char* vSource, const char* fSource) {
    unsigned int vertex, fragment;
    int success;
    char infoLog[512];

    vertex = glCreateShader(GL_VERTEX_SHADER);
    glShaderSource(vertex, 1, &vSource, NULL);
    glCompileShader(vertex);
    glGetShaderiv(vertex, GL_COMPILE_STATUS, &success);
    if(!success) { glGetShaderInfoLog(vertex, 512, NULL, infoLog); std::cout << "Vertex Shader Error: " << infoLog << std::endl; }

    fragment = glCreateShader(GL_FRAGMENT_SHADER);
    glShaderSource(fragment, 1, &fSource, NULL);
    glCompileShader(fragment);
    glGetShaderiv(fragment, GL_COMPILE_STATUS, &success);
    if(!success) { glGetShaderInfoLog(fragment, 512, NULL, infoLog); std::cout << "Fragment Shader Error: " << infoLog << std::endl; }

    unsigned int ID = glCreateProgram();
    glAttachShader(ID, vertex);
    glAttachShader(ID, fragment);
    glLinkProgram(ID);
    glGetProgramiv(ID, GL_LINK_STATUS, &success);
    if(!success) { glGetProgramInfoLog(ID, 512, NULL, infoLog); std::cout << "Link Error: " << infoLog << std::endl; }

    glDeleteShader(vertex);
    glDeleteShader(fragment);
    return ID;
}


int main() {
    // 初始化 GLFW
    glfwInit();
    glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
    glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3);
    glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);

    GLFWwindow* window = glfwCreateWindow(SCR_WIDTH, SCR_HEIGHT, "Fluid Simulation (SPH)", NULL, NULL);
    if (window == NULL) {
        std::cout << "Failed to create GLFW window" << std::endl;
        glfwTerminate();
        return -1;
    }
    glfwMakeContextCurrent(window);
    glfwSetCursorPosCallback(window, mouse_callback);
    glfwSetKeyCallback(window, key_callback);
    glfwSetInputMode(window, GLFW_CURSOR, GLFW_CURSOR_DISABLED);

    // 初始化 GLAD
    if (!gladLoadGLLoader((GLADloadproc)glfwGetProcAddress)) {
        std::cout << "Failed to initialize GLAD" << std::endl;
        return -1;
    }

    glEnable(GL_DEPTH_TEST);
    glEnable(GL_BLEND);
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

    // 编译 Shader
    unsigned int shaderProgram = compileShader(vertexShaderSource, fragmentShaderSource);

    // 创建球体模型数据
    std::vector<float> sphereVertices;
    std::vector<unsigned int> sphereIndices;
    createSphere(sphereVertices, sphereIndices, 12, 12); // 低多边形球体以提升性能

    // 配置 OpenGL Buffer
    unsigned int VAO, VBO, EBO, instanceVBO;
    glGenVertexArrays(1, &VAO);
    glGenBuffers(1, &VBO);
    glGenBuffers(1, &EBO);
    glGenBuffers(1, &instanceVBO);

    glBindVertexArray(VAO);

    // 1. 球体模型 VBO
    glBindBuffer(GL_ARRAY_BUFFER, VBO);
    glBufferData(GL_ARRAY_BUFFER, sphereVertices.size() * sizeof(float), sphereVertices.data(), GL_STATIC_DRAW);
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, EBO);
    glBufferData(GL_ELEMENT_ARRAY_BUFFER, sphereIndices.size() * sizeof(unsigned int), sphereIndices.data(), GL_STATIC_DRAW);

    // 属性 0: 顶点位置
    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 6 * sizeof(float), (void*)0);
    glEnableVertexAttribArray(0);
    // 属性 1: 法线
    glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 6 * sizeof(float), (void*)(3 * sizeof(float)));
    glEnableVertexAttribArray(1);

    // 2. 实例位置 VBO 
    glBindBuffer(GL_ARRAY_BUFFER, instanceVBO);
    // 初始分配大小，但数据稍后更新 (DYNAMIC_DRAW)
    glBufferData(GL_ARRAY_BUFFER, PARTICLE_COUNT * sizeof(glm::vec3), NULL, GL_DYNAMIC_DRAW);

    // 属性 2: 实例偏移 (layout 2)
    glVertexAttribPointer(2, 3, GL_FLOAT, GL_FALSE, sizeof(glm::vec3), (void*)0);
    glEnableVertexAttribArray(2);
    // 告诉 OpenGL 这个属性是每实例更新一次，而不是每顶点
    glVertexAttribDivisor(2, 1); 

    glBindVertexArray(0);

    // 初始化仿真系统
    SPHSystem sph;

    // 容器线框数据 (用于可视化边界)
    float containerVertices[] = {
        // 底部矩形
        -CONTAINER_SIZE.x/2, -CONTAINER_SIZE.y/2, -CONTAINER_SIZE.z/2,
         CONTAINER_SIZE.x/2, -CONTAINER_SIZE.y/2, -CONTAINER_SIZE.z/2,
         CONTAINER_SIZE.x/2, -CONTAINER_SIZE.y/2,  CONTAINER_SIZE.z/2,
        -CONTAINER_SIZE.x/2, -CONTAINER_SIZE.y/2,  CONTAINER_SIZE.z/2,
        // 顶部矩形
        -CONTAINER_SIZE.x/2,  CONTAINER_SIZE.y/2, -CONTAINER_SIZE.z/2,
         CONTAINER_SIZE.x/2,  CONTAINER_SIZE.y/2, -CONTAINER_SIZE.z/2,
         CONTAINER_SIZE.x/2,  CONTAINER_SIZE.y/2,  CONTAINER_SIZE.z/2,
        -CONTAINER_SIZE.x/2,  CONTAINER_SIZE.y/2,  CONTAINER_SIZE.z/2,
    };
    unsigned int containerIndices[] = {
        0, 1, 1, 2, 2, 3, 3, 0, // 底
        4, 5, 5, 6, 6, 7, 7, 4, // 顶
        0, 4, 1, 5, 2, 6, 3, 7  // 柱
    };
    unsigned int cVAO, cVBO, cEBO;
    glGenVertexArrays(1, &cVAO);
    glGenBuffers(1, &cVBO);
    glGenBuffers(1, &cEBO);
    glBindVertexArray(cVAO);
    glBindBuffer(GL_ARRAY_BUFFER, cVBO);
    glBufferData(GL_ARRAY_BUFFER, sizeof(containerVertices), containerVertices, GL_STATIC_DRAW);
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, cEBO);
    glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(containerIndices), containerIndices, GL_STATIC_DRAW);
    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(float), (void*)0);
    glEnableVertexAttribArray(0);

    // 时间统计
    float lastFrame = 0.0f;

    // 定义物理子步数
    const int SUB_STEPS = 4; 

    while (!glfwWindowShouldClose(window)) {
        // 计算时间差
        float currentFrame = static_cast<float>(glfwGetTime());
        float deltaTime = currentFrame - lastFrame;
        lastFrame = currentFrame;

        // 输入处理
        if (glfwGetKey(window, GLFW_KEY_W) == GLFW_PRESS) camera.ProcessKeyboard(GLFW_KEY_W, deltaTime);
        if (glfwGetKey(window, GLFW_KEY_S) == GLFW_PRESS) camera.ProcessKeyboard(GLFW_KEY_S, deltaTime);
        if (glfwGetKey(window, GLFW_KEY_A) == GLFW_PRESS) camera.ProcessKeyboard(GLFW_KEY_A, deltaTime);
        if (glfwGetKey(window, GLFW_KEY_D) == GLFW_PRESS) camera.ProcessKeyboard(GLFW_KEY_D, deltaTime);

        // 逻辑重置
        if (resetSim) {
            sph.reset();
            resetSim = false;
        }


        // 无论当前渲染帧率是多少，物理引擎都必须以精细的步伐前进
        // 这能彻底解决粒子穿墙和爆炸的问题
        for (int s = 0; s < SUB_STEPS; ++s) {
            sph.update();
        }

        // 渲染准备
        glClearColor(0.1f, 0.1f, 0.1f, 1.0f);
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

        glm::mat4 projection = glm::perspective(glm::radians(45.0f), (float)SCR_WIDTH / (float)SCR_HEIGHT, 0.1f, 100.0f);
        glm::mat4 view = camera.GetViewMatrix();

        glUseProgram(shaderProgram);
        glUniformMatrix4fv(glGetUniformLocation(shaderProgram, "projection"), 1, GL_FALSE, glm::value_ptr(projection));
        glUniformMatrix4fv(glGetUniformLocation(shaderProgram, "view"), 1, GL_FALSE, glm::value_ptr(view));
        glUniform3fv(glGetUniformLocation(shaderProgram, "viewPos"), 1, glm::value_ptr(camera.Position));
        glUniform1f(glGetUniformLocation(shaderProgram, "radius"), PARTICLE_RADIUS);

        // 更新实例数据 (将最新的粒子位置传给 GPU)
        glBindBuffer(GL_ARRAY_BUFFER, instanceVBO);
        glBufferSubData(GL_ARRAY_BUFFER, 0, PARTICLE_COUNT * sizeof(glm::vec3), sph.positions.data());

        // 绘制流体粒子 (Instanced)
        glBindVertexArray(VAO);
        glDrawElementsInstanced(GL_TRIANGLES, (unsigned int)sphereIndices.size(), GL_UNSIGNED_INT, 0, PARTICLE_COUNT);

        // 绘制容器边界 (简单的白色线框)
        glBindVertexArray(cVAO);
        glUniform1f(glGetUniformLocation(shaderProgram, "radius"), 1.0f); // 容器不需要radius缩放
        glDrawElements(GL_LINES, 24, GL_UNSIGNED_INT, 0);

        glfwSwapBuffers(window);
        glfwPollEvents();
    }

    glfwTerminate();
    return 0;
}