using OpenTK.Graphics.OpenGL4;
using OpenTK.Windowing.Common;
using OpenTK.Windowing.Desktop;
using OpenTK.Mathematics;
using OpenTK.Windowing.GraphicsLibraryFramework;

namespace fluidsC_;

public class FluidWindow3D_GPU : GameWindow
{
    private int _fluidSizeX = 64;
    private int _fluidSizeY = 64;
    private int _fluidSizeZ = 64;

    // GPU buffers - double buffered
    private int[] _densityTextures = new int[2];
    private int[] _velocityTextures = new int[2];
    private int[] _pressureTextures = new int[2];
    private int _divergenceTexture;
    private int _obstacleTexture;

    // Shader programs
    private int _advectionShader;
    private int _forceShader;
    private int _divergenceShader;
    private int _pressureShader;
    private int _gradientShader;
    private int _renderShader;
    private int _clearShader;
    private int _obstacleShader;

    private int _vertexArrayObject;
    private int _pointCount;

    private Camera _camera;
    private bool _firstMove = true;
    private Vector2 _lastPos;
    private float _time = 0;

    private int _currentDensity = 0;
    private int _currentVelocity = 0;
    private int _currentPressure = 0;

    // Simulation parameters
    private float _diffusion = 0.0001f;
    private float _viscosity = 0.0001f;

    public FluidWindow3D_GPU(GameWindowSettings gameWindowSettings, NativeWindowSettings nativeWindowSettings)
        : base(gameWindowSettings, nativeWindowSettings)
    {
        _camera = new Camera(Vector3.UnitZ * 3, Size.X / (float)Size.Y);
    }

    protected override void OnLoad()
    {
        base.OnLoad();

        Console.WriteLine($"OpenGL Version: {GL.GetString(StringName.Version)}");
        Console.WriteLine($"GLSL Version: {GL.GetString(StringName.ShadingLanguageVersion)}");

        GL.ClearColor(0.02f, 0.02f, 0.05f, 1.0f);
        GL.Enable(EnableCap.DepthTest);
        GL.Enable(EnableCap.Blend);
        GL.BlendFunc(BlendingFactor.SrcAlpha, BlendingFactor.OneMinusSrcAlpha);
        GL.Enable(EnableCap.ProgramPointSize);

        CreateTextures();
        SetupShaders();
        SetupPointCloud();
        InitializeObstacles();

        // Add initial fluid
        AddInitialFluid();

        Console.WriteLine("🎯 WORKING GPU Fluid Simulation - CLICK AND DRAG!");
        Console.WriteLine("Left Click: Add Fluid | Right Click: Add Force | R: Reset");
    }

    private void CreateTextures()
    {
        for (int i = 0; i < 2; i++)
        {
            _densityTextures[i] = Create3DTexture(_fluidSizeX, _fluidSizeY, _fluidSizeZ, false);
            _velocityTextures[i] = Create3DTexture(_fluidSizeX, _fluidSizeY, _fluidSizeZ, true);
            _pressureTextures[i] = Create3DTexture(_fluidSizeX, _fluidSizeY, _fluidSizeZ, false);
        }
        _divergenceTexture = Create3DTexture(_fluidSizeX, _fluidSizeY, _fluidSizeZ, false);
        _obstacleTexture = Create3DTexture(_fluidSizeX, _fluidSizeY, _fluidSizeZ, false);

        ClearAllTextures();
    }

    private int Create3DTexture(int width, int height, int depth, bool rgb)
    {
        int texture = GL.GenTexture();
        GL.BindTexture(TextureTarget.Texture3D, texture);

        GL.TexParameter(TextureTarget.Texture3D, TextureParameterName.TextureMinFilter, (int)TextureMinFilter.Linear);
        GL.TexParameter(TextureTarget.Texture3D, TextureParameterName.TextureMagFilter, (int)TextureMagFilter.Linear);
        GL.TexParameter(TextureTarget.Texture3D, TextureParameterName.TextureWrapS, (int)TextureWrapMode.ClampToEdge);
        GL.TexParameter(TextureTarget.Texture3D, TextureParameterName.TextureWrapT, (int)TextureWrapMode.ClampToEdge);
        GL.TexParameter(TextureTarget.Texture3D, TextureParameterName.TextureWrapR, (int)TextureWrapMode.ClampToEdge);

        var format = rgb ? PixelInternalFormat.Rgb32f : PixelInternalFormat.R32f;
        var pixelFormat = rgb ? PixelFormat.Rgb : PixelFormat.Red;

        GL.TexImage3D(TextureTarget.Texture3D, 0, format, width, height, depth, 0, pixelFormat, PixelType.Float, IntPtr.Zero);

        return texture;
    }

    private void SetupShaders()
    {
        // PROPER FLUID SIMULATION SHADERS
        _advectionShader = CreateComputeShader(@"
            #version 430
            layout(local_size_x = 4, local_size_y = 4, local_size_z = 4) in;
            
            layout(r32f, binding = 0) uniform image3D densityIn;
            layout(rgba32f, binding = 1) uniform image3D velocityIn;
            layout(r32f, binding = 2) uniform image3D densityOut;
            layout(rgba32f, binding = 3) uniform image3D velocityOut;
            layout(r32f, binding = 4) uniform image3D obstacles;
            
            uniform float dt;
            uniform float dissipation;
            
            vec3 sampleVelocity(vec3 pos) {
                ivec3 coord = ivec3(clamp(pos, vec3(0), vec3(imageSize(velocityIn)) - vec3(1)));
                return imageLoad(velocityIn, coord).rgb;
            }
            
            float sampleDensity(vec3 pos) {
                ivec3 coord = ivec3(clamp(pos, vec3(0), vec3(imageSize(densityIn)) - vec3(1)));
                return imageLoad(densityIn, coord).r;
            }
            
            void main() {
                ivec3 coord = ivec3(gl_GlobalInvocationID.xyz);
                ivec3 size = imageSize(densityIn);
                
                if(coord.x >= size.x || coord.y >= size.y || coord.z >= size.z) 
                    return;
                
                if(imageLoad(obstacles, coord).r > 0.5) {
                    imageStore(densityOut, coord, vec4(0.0));
                    imageStore(velocityOut, coord, vec4(0.0));
                    return;
                }
                
                // Semi-Lagrangian advection
                vec3 vel = sampleVelocity(vec3(coord));
                vec3 prevPos = vec3(coord) - vel * dt * 10.0;
                
                // Advect density
                float density = sampleDensity(prevPos);
                imageStore(densityOut, coord, vec4(density * dissipation, 0.0, 0.0, 1.0));
                
                // Advect velocity (self-advection)
                vec3 prevVel = sampleVelocity(prevPos);
                imageStore(velocityOut, coord, vec4(prevVel, 1.0));
            }");

        _forceShader = CreateComputeShader(@"
            #version 430
            layout(local_size_x = 4, local_size_y = 4, local_size_z = 4) in;
            
            layout(rgba32f, binding = 0) uniform image3D velocity;
            layout(r32f, binding = 1) uniform image3D density;
            layout(r32f, binding = 2) uniform image3D obstacles;
            
            uniform vec3 force;
            uniform vec3 position;
            uniform float radius;
            uniform float densityAmount;
            
            void main() {
                ivec3 coord = ivec3(gl_GlobalInvocationID.xyz);
                
                if(imageLoad(obstacles, coord).r > 0.5) 
                    return;
                
                vec3 cellPos = vec3(coord);
                vec3 diff = cellPos - position;
                float dist = length(diff);
                
                if(dist < radius) {
                    float influence = (1.0 - dist/radius) * (1.0 - dist/radius);
                    
                    // Add force to velocity
                    vec3 currentVel = imageLoad(velocity, coord).rgb;
                    vec3 newVel = currentVel + force * influence * 0.5;
                    imageStore(velocity, coord, vec4(newVel, 1.0));
                    
                    // Add density
                    if(densityAmount > 0.0) {
                        float currentDensity = imageLoad(density, coord).r;
                        float newDensity = currentDensity + densityAmount * influence;
                        imageStore(density, coord, vec4(newDensity, 0.0, 0.0, 1.0));
                    }
                }
            }");

        _divergenceShader = CreateComputeShader(@"
            #version 430
            layout(local_size_x = 4, local_size_y = 4, local_size_z = 4) in;
            
            layout(rgba32f, binding = 0) uniform image3D velocity;
            layout(r32f, binding = 1) uniform image3D divergence;
            layout(r32f, binding = 2) uniform image3D obstacles;
            
            void main() {
                ivec3 coord = ivec3(gl_GlobalInvocationID.xyz);
                ivec3 size = imageSize(velocity);
                
                if(coord.x < 1 || coord.x >= size.x-1 || 
                   coord.y < 1 || coord.y >= size.y-1 ||
                   coord.z < 1 || coord.z >= size.z-1) 
                    return;
                
                if(imageLoad(obstacles, coord).r > 0.5) {
                    imageStore(divergence, coord, vec4(0.0));
                    return;
                }
                
                // Calculate divergence
                float left = imageLoad(velocity, coord - ivec3(1,0,0)).r;
                float right = imageLoad(velocity, coord + ivec3(1,0,0)).r;
                float down = imageLoad(velocity, coord - ivec3(0,1,0)).g;
                float up = imageLoad(velocity, coord + ivec3(0,1,0)).g;
                float back = imageLoad(velocity, coord - ivec3(0,0,1)).b;
                float front = imageLoad(velocity, coord + ivec3(0,0,1)).b;
                
                float div = (right - left + up - down + front - back) * 0.5;
                imageStore(divergence, coord, vec4(div, 0.0, 0.0, 1.0));
            }");

        _pressureShader = CreateComputeShader(@"
            #version 430
            layout(local_size_x = 4, local_size_y = 4, local_size_z = 4) in;
            
            layout(r32f, binding = 0) uniform image3D pressureIn;
            layout(r32f, binding = 1) uniform image3D pressureOut;
            layout(r32f, binding = 2) uniform image3D divergence;
            layout(r32f, binding = 3) uniform image3D obstacles;
            
            void main() {
                ivec3 coord = ivec3(gl_GlobalInvocationID.xyz);
                ivec3 size = imageSize(pressureIn);
                
                if(coord.x < 1 || coord.x >= size.x-1 || 
                   coord.y < 1 || coord.y >= size.y-1 ||
                   coord.z < 1 || coord.z >= size.z-1) 
                    return;
                
                if(imageLoad(obstacles, coord).r > 0.5) {
                    imageStore(pressureOut, coord, vec4(0.0));
                    return;
                }
                
                // Jacobi iteration for pressure
                float left = imageLoad(pressureIn, coord - ivec3(1,0,0)).r;
                float right = imageLoad(pressureIn, coord + ivec3(1,0,0)).r;
                float down = imageLoad(pressureIn, coord - ivec3(0,1,0)).r;
                float up = imageLoad(pressureIn, coord + ivec3(0,1,0)).r;
                float back = imageLoad(pressureIn, coord - ivec3(0,0,1)).r;
                float front = imageLoad(pressureIn, coord + ivec3(0,0,1)).r;
                
                float div = imageLoad(divergence, coord).r;
                float pressure = (left + right + down + up + back + front - div) / 6.0;
                imageStore(pressureOut, coord, vec4(pressure, 0.0, 0.0, 1.0));
            }");

        _gradientShader = CreateComputeShader(@"
            #version 430
            layout(local_size_x = 4, local_size_y = 4, local_size_z = 4) in;
            
            layout(rgba32f, binding = 0) uniform image3D velocity;
            layout(r32f, binding = 1) uniform image3D pressure;
            layout(r32f, binding = 2) uniform image3D obstacles;
            
            void main() {
                ivec3 coord = ivec3(gl_GlobalInvocationID.xyz);
                ivec3 size = imageSize(velocity);
                
                if(coord.x < 1 || coord.x >= size.x-1 || 
                   coord.y < 1 || coord.y >= size.y-1 ||
                   coord.z < 1 || coord.z >= size.z-1) 
                    return;
                
                if(imageLoad(obstacles, coord).r > 0.5) {
                    imageStore(velocity, coord, vec4(0.0));
                    return;
                }
                
                // Subtract pressure gradient
                float pressureLeft = imageLoad(pressure, coord - ivec3(1,0,0)).r;
                float pressureRight = imageLoad(pressure, coord + ivec3(1,0,0)).r;
                float pressureDown = imageLoad(pressure, coord - ivec3(0,1,0)).r;
                float pressureUp = imageLoad(pressure, coord + ivec3(0,1,0)).r;
                float pressureBack = imageLoad(pressure, coord - ivec3(0,0,1)).r;
                float pressureFront = imageLoad(pressure, coord + ivec3(0,0,1)).r;
                
                vec3 grad = vec3(
                    pressureRight - pressureLeft,
                    pressureUp - pressureDown, 
                    pressureFront - pressureBack
                ) * 0.5;
                
                vec3 currentVel = imageLoad(velocity, coord).rgb;
                vec3 newVel = currentVel - grad;
                imageStore(velocity, coord, vec4(newVel, 1.0));
            }");

        _clearShader = CreateComputeShader(@"
            #version 430
            layout(local_size_x = 8, local_size_y = 8, local_size_z = 8) in;
            layout(r32f, binding = 0) uniform image3D textureToClear;
            
            uniform float value;
            
            void main() {
                ivec3 coord = ivec3(gl_GlobalInvocationID.xyz);
                ivec3 size = imageSize(textureToClear);
                if(coord.x >= size.x || coord.y >= size.y || coord.z >= size.z) 
                    return;
                
                imageStore(textureToClear, coord, vec4(value, 0.0, 0.0, 1.0));
            }");

        _obstacleShader = CreateComputeShader(@"
            #version 430
            layout(local_size_x = 8, local_size_y = 8, local_size_z = 8) in;
            layout(r32f, binding = 0) uniform image3D obstacles;
            
            uniform vec3 position;
            uniform float radius;
            uniform float setValue;
            
            void main() {
                ivec3 coord = ivec3(gl_GlobalInvocationID.xyz);
                ivec3 size = imageSize(obstacles);
                if(coord.x >= size.x || coord.y >= size.y || coord.z >= size.z) 
                    return;
                
                vec3 cellPos = vec3(coord);
                float dist = length(cellPos - position);
                
                if(dist <= radius) {
                    imageStore(obstacles, coord, vec4(setValue, 0.0, 0.0, 1.0));
                }
            }");

        // RENDER SHADER
        _renderShader = CreateRenderShader();
    }

    private int CreateComputeShader(string source)
    {
        var shader = GL.CreateShader(ShaderType.ComputeShader);
        GL.ShaderSource(shader, source);
        GL.CompileShader(shader);

        GL.GetShader(shader, ShaderParameter.CompileStatus, out int success);
        if (success == 0)
        {
            var infoLog = GL.GetShaderInfoLog(shader);
            Console.WriteLine($"Compute Shader Compilation Failed:\n{infoLog}");
        }

        var program = GL.CreateProgram();
        GL.AttachShader(program, shader);
        GL.LinkProgram(program);

        GL.DeleteShader(shader);
        return program;
    }

    private int CreateRenderShader()
    {
        const string vertexShaderSource = @"
            #version 330 core
            layout (location = 0) in vec3 aPosition;
            out vec3 fragTexCoord;
            uniform mat4 model;
            uniform mat4 view;
            uniform mat4 projection;
            void main()
            {
                gl_Position = projection * view * model * vec4(aPosition, 1.0);
                fragTexCoord = (aPosition + 1.0) * 0.5;
                gl_PointSize = 4.0;
            }";

        const string fragmentShaderSource = @"
            #version 330 core
            in vec3 fragTexCoord;
            out vec4 fragColor;
            uniform sampler3D densityTexture;
            uniform sampler3D obstacleTexture;
            uniform float time;
            void main()
            {
                vec3 texCoord = fragTexCoord;
                if(texCoord.x < 0.0 || texCoord.x > 1.0 || 
                   texCoord.y < 0.0 || texCoord.y > 1.0 || 
                   texCoord.z < 0.0 || texCoord.z > 1.0) 
                    discard;
                
                float density = texture(densityTexture, texCoord).r;
                float obstacle = texture(obstacleTexture, texCoord).r;
                
                if(obstacle > 0.5) {
                    fragColor = vec4(0.3, 0.3, 0.3, 1.0);
                } else if(density > 0.01) {
                    vec3 color = mix(vec3(0.1, 0.3, 0.8), vec3(0.0, 0.8, 1.0), density * 2.0);
                    color += 0.1 * sin(time + fragTexCoord.y * 10.0);
                    fragColor = vec4(color, density * 0.9);
                } else {
                    discard;
                }
            }";

        var vertexShader = GL.CreateShader(ShaderType.VertexShader);
        GL.ShaderSource(vertexShader, vertexShaderSource);
        GL.CompileShader(vertexShader);

        var fragmentShader = GL.CreateShader(ShaderType.FragmentShader);
        GL.ShaderSource(fragmentShader, fragmentShaderSource);
        GL.CompileShader(fragmentShader);

        var program = GL.CreateProgram();
        GL.AttachShader(program, vertexShader);
        GL.AttachShader(program, fragmentShader);
        GL.LinkProgram(program);

        GL.DeleteShader(vertexShader);
        GL.DeleteShader(fragmentShader);

        return program;
    }

    private void SetupPointCloud()
    {
        var vertices = new List<Vector3>();

        // Reduced point count for performance but still good coverage
        for (int x = 0; x < _fluidSizeX; x += 1)
        {
            for (int y = 0; y < _fluidSizeY; y += 1)
            {
                for (int z = 0; z < _fluidSizeZ; z += 1)
                {
                    float px = (x / (float)_fluidSizeX) * 2.0f - 1.0f;
                    float py = (y / (float)_fluidSizeY) * 2.0f - 1.0f;
                    float pz = (z / (float)_fluidSizeZ) * 2.0f - 1.0f;
                    vertices.Add(new Vector3(px, py, pz));
                }
            }
        }

        _pointCount = vertices.Count;
        Console.WriteLine($"Rendering {_pointCount} points");

        _vertexArrayObject = GL.GenVertexArray();
        int vbo = GL.GenBuffer();

        GL.BindVertexArray(_vertexArrayObject);
        GL.BindBuffer(BufferTarget.ArrayBuffer, vbo);
        GL.BufferData(BufferTarget.ArrayBuffer, vertices.Count * 3 * sizeof(float), vertices.ToArray(), BufferUsageHint.StaticDraw);

        GL.VertexAttribPointer(0, 3, VertexAttribPointerType.Float, false, 3 * sizeof(float), 0);
        GL.EnableVertexAttribArray(0);
    }

    private void InitializeObstacles()
    {
        // Simple box boundaries
        AddBoxObstacle(new Vector3(0, 0, 0), new Vector3(_fluidSizeX, 2, _fluidSizeZ)); // Floor
        AddBoxObstacle(new Vector3(0, 0, 0), new Vector3(2, _fluidSizeY, _fluidSizeZ)); // Left
        AddBoxObstacle(new Vector3(_fluidSizeX - 2, 0, 0), new Vector3(2, _fluidSizeY, _fluidSizeZ)); // Right
        AddBoxObstacle(new Vector3(0, 0, 0), new Vector3(_fluidSizeX, _fluidSizeY, 2)); // Back
        AddBoxObstacle(new Vector3(0, 0, _fluidSizeZ - 2), new Vector3(_fluidSizeX, _fluidSizeY, 2)); // Front

        // Central obstacle
        AddSphereObstacle(new Vector3(_fluidSizeX / 2, _fluidSizeY / 4, _fluidSizeZ / 2), 6.0f);
    }

    private void AddInitialFluid()
    {
        // Add initial falling fluid
        Vector3 sourcePos = new Vector3(_fluidSizeX / 2, _fluidSizeY - 10, _fluidSizeZ / 2);
        AddForce(sourcePos, new Vector3(0.0f, -2.0f, 0.0f), 8.0f, 10.0f);

        // Add some random fluid blobs
        Random rand = new Random();
        for (int i = 0; i < 5; i++)
        {
            Vector3 randomPos = new Vector3(
                rand.Next(10, _fluidSizeX - 10),
                rand.Next(10, _fluidSizeY - 10),
                rand.Next(10, _fluidSizeZ - 10)
            );
            AddForce(randomPos, Vector3.Zero, 4.0f, 3.0f);
        }
    }

    private void AddSphereObstacle(Vector3 center, float radius)
    {
        GL.UseProgram(_obstacleShader);
        BindImageTexture(0, _obstacleTexture, TextureAccess.WriteOnly);

        GL.Uniform3(GL.GetUniformLocation(_obstacleShader, "position"), center);
        GL.Uniform1(GL.GetUniformLocation(_obstacleShader, "radius"), radius);
        GL.Uniform1(GL.GetUniformLocation(_obstacleShader, "setValue"), 1.0f);

        DispatchCompute(_fluidSizeX, _fluidSizeY, _fluidSizeZ, 8, 8, 8);
        GL.MemoryBarrier(MemoryBarrierFlags.ShaderImageAccessBarrierBit);
    }

    private void AddBoxObstacle(Vector3 min, Vector3 max)
    {
        Vector3 center = (min + max) * 0.5f;
        Vector3 size = max - min;
        AddSphereObstacle(center, Math.Max(size.X, Math.Max(size.Y, size.Z)) * 0.5f);
    }

    protected override void OnRenderFrame(FrameEventArgs e)
    {
        base.OnRenderFrame(e);

        GL.Clear(ClearBufferMask.ColorBufferBit | ClearBufferMask.DepthBufferBit);

        StepSimulation((float)e.Time);
        RenderFluid();

        SwapBuffers();
    }

    private void StepSimulation(float dt)
    {
        dt = Math.Min(dt, 0.016f); // Cap at ~60fps

        // Full fluid simulation steps
        Advect(dt, 0.99f);
        AddExternalForces(dt);
        Project(dt);
        AddDemoSources();

        // Swap buffers
        _currentDensity = 1 - _currentDensity;
        _currentVelocity = 1 - _currentVelocity;
        _currentPressure = 1 - _currentPressure;
    }

    private void Advect(float dt, float dissipation)
    {
        GL.UseProgram(_advectionShader);

        BindImageTexture(0, _densityTextures[_currentDensity], TextureAccess.ReadOnly);
        BindImageTexture(1, _velocityTextures[_currentVelocity], TextureAccess.ReadOnly);
        BindImageTexture(2, _densityTextures[1 - _currentDensity], TextureAccess.WriteOnly);
        BindImageTexture(3, _velocityTextures[1 - _currentVelocity], TextureAccess.WriteOnly);
        BindImageTexture(4, _obstacleTexture, TextureAccess.ReadOnly);

        GL.Uniform1(GL.GetUniformLocation(_advectionShader, "dt"), dt);
        GL.Uniform1(GL.GetUniformLocation(_advectionShader, "dissipation"), dissipation);

        DispatchCompute(_fluidSizeX, _fluidSizeY, _fluidSizeZ, 4, 4, 4);
        GL.MemoryBarrier(MemoryBarrierFlags.ShaderImageAccessBarrierBit);
    }

    private void AddExternalForces(float dt)
    {
        // Gravity
        Vector3 gravity = new Vector3(0.0f, -9.8f, 0.0f) * dt;

        GL.UseProgram(_forceShader);
        BindImageTexture(0, _velocityTextures[_currentVelocity], TextureAccess.ReadWrite);
        BindImageTexture(1, _densityTextures[_currentDensity], TextureAccess.ReadWrite);
        BindImageTexture(2, _obstacleTexture, TextureAccess.ReadOnly);

        GL.Uniform3(GL.GetUniformLocation(_forceShader, "force"), gravity);
        GL.Uniform3(GL.GetUniformLocation(_forceShader, "position"), new Vector3(_fluidSizeX / 2, _fluidSizeY / 2, _fluidSizeZ / 2));
        GL.Uniform1(GL.GetUniformLocation(_forceShader, "radius"), _fluidSizeX * 2.0f); // Affect entire field
        GL.Uniform1(GL.GetUniformLocation(_forceShader, "densityAmount"), 0.0f);

        DispatchCompute(_fluidSizeX, _fluidSizeY, _fluidSizeZ, 4, 4, 4);
        GL.MemoryBarrier(MemoryBarrierFlags.ShaderImageAccessBarrierBit);
    }

    private void Project(float dt)
    {
        // Calculate divergence
        GL.UseProgram(_divergenceShader);
        BindImageTexture(0, _velocityTextures[_currentVelocity], TextureAccess.ReadOnly);
        BindImageTexture(1, _divergenceTexture, TextureAccess.WriteOnly);
        BindImageTexture(2, _obstacleTexture, TextureAccess.ReadOnly);
        DispatchCompute(_fluidSizeX, _fluidSizeY, _fluidSizeZ, 4, 4, 4);
        GL.MemoryBarrier(MemoryBarrierFlags.ShaderImageAccessBarrierBit);

        // Solve pressure (Jacobi iterations)
        for (int i = 0; i < 10; i++)
        {
            GL.UseProgram(_pressureShader);
            BindImageTexture(0, _pressureTextures[_currentPressure], TextureAccess.ReadOnly);
            BindImageTexture(1, _pressureTextures[1 - _currentPressure], TextureAccess.WriteOnly);
            BindImageTexture(2, _divergenceTexture, TextureAccess.ReadOnly);
            BindImageTexture(3, _obstacleTexture, TextureAccess.ReadOnly);
            DispatchCompute(_fluidSizeX, _fluidSizeY, _fluidSizeZ, 4, 4, 4);
            GL.MemoryBarrier(MemoryBarrierFlags.ShaderImageAccessBarrierBit);
            _currentPressure = 1 - _currentPressure;
        }

        // Subtract pressure gradient
        GL.UseProgram(_gradientShader);
        BindImageTexture(0, _velocityTextures[_currentVelocity], TextureAccess.ReadWrite);
        BindImageTexture(1, _pressureTextures[_currentPressure], TextureAccess.ReadOnly);
        BindImageTexture(2, _obstacleTexture, TextureAccess.ReadOnly);
        DispatchCompute(_fluidSizeX, _fluidSizeY, _fluidSizeZ, 4, 4, 4);
        GL.MemoryBarrier(MemoryBarrierFlags.ShaderImageAccessBarrierBit);
    }

    private void AddDemoSources()
    {
        // Continuous falling source
        Vector3 sourcePos = new Vector3(_fluidSizeX / 2, _fluidSizeY - 8, _fluidSizeZ / 2);
        AddForce(sourcePos, new Vector3(0.0f, -3.0f, 0.0f), 6.0f, 2.0f);

        // Random turbulence
        if ((int)(_time * 2) % 4 == 0)
        {
            Vector3 randomPos = new Vector3(
                _fluidSizeX * 0.3f + MathF.Sin(_time) * 5.0f,
                _fluidSizeY - 5,
                _fluidSizeZ * 0.7f + MathF.Cos(_time) * 5.0f
            );
            AddForce(randomPos, new Vector3(0.0f, -2.0f, 0.0f), 4.0f, 1.5f);
        }
    }

    private void AddForce(Vector3 position, Vector3 force, float radius, float densityAmount)
    {
        GL.UseProgram(_forceShader);

        BindImageTexture(0, _velocityTextures[_currentVelocity], TextureAccess.ReadWrite);
        BindImageTexture(1, _densityTextures[_currentDensity], TextureAccess.ReadWrite);
        BindImageTexture(2, _obstacleTexture, TextureAccess.ReadOnly);

        GL.Uniform3(GL.GetUniformLocation(_forceShader, "force"), force);
        GL.Uniform3(GL.GetUniformLocation(_forceShader, "position"), position);
        GL.Uniform1(GL.GetUniformLocation(_forceShader, "radius"), radius);
        GL.Uniform1(GL.GetUniformLocation(_forceShader, "densityAmount"), densityAmount);

        DispatchCompute(_fluidSizeX, _fluidSizeY, _fluidSizeZ, 4, 4, 4);
        GL.MemoryBarrier(MemoryBarrierFlags.ShaderImageAccessBarrierBit);
    }

    private void RenderFluid()
    {
        GL.UseProgram(_renderShader);

        var model = Matrix4.Identity;
        var view = _camera.GetViewMatrix();
        var projection = _camera.GetProjectionMatrix();

        GL.UniformMatrix4(GL.GetUniformLocation(_renderShader, "model"), false, ref model);
        GL.UniformMatrix4(GL.GetUniformLocation(_renderShader, "view"), false, ref view);
        GL.UniformMatrix4(GL.GetUniformLocation(_renderShader, "projection"), false, ref projection);
        GL.Uniform1(GL.GetUniformLocation(_renderShader, "time"), _time);

        GL.ActiveTexture(TextureUnit.Texture0);
        GL.BindTexture(TextureTarget.Texture3D, _densityTextures[_currentDensity]);
        GL.Uniform1(GL.GetUniformLocation(_renderShader, "densityTexture"), 0);

        GL.ActiveTexture(TextureUnit.Texture1);
        GL.BindTexture(TextureTarget.Texture3D, _obstacleTexture);
        GL.Uniform1(GL.GetUniformLocation(_renderShader, "obstacleTexture"), 1);

        GL.BindVertexArray(_vertexArrayObject);
        GL.DrawArrays(PrimitiveType.Points, 0, _pointCount);
    }

    private void BindImageTexture(int unit, int texture, TextureAccess access)
    {
        GL.BindImageTexture(unit, texture, 0, false, 0, access, SizedInternalFormat.R32f);
    }

    private void DispatchCompute(int sizeX, int sizeY, int sizeZ, int localSizeX, int localSizeY, int localSizeZ)
    {
        int groupsX = (sizeX + localSizeX - 1) / localSizeX;
        int groupsY = (sizeY + localSizeY - 1) / localSizeY;
        int groupsZ = (sizeZ + localSizeZ - 1) / localSizeZ;

        GL.DispatchCompute(groupsX, groupsY, groupsZ);
    }

    private void ClearAllTextures()
    {
        GL.UseProgram(_clearShader);

        foreach (var texture in _densityTextures)
        {
            BindImageTexture(0, texture, TextureAccess.WriteOnly);
            GL.Uniform1(GL.GetUniformLocation(_clearShader, "value"), 0.0f);
            DispatchCompute(_fluidSizeX, _fluidSizeY, _fluidSizeZ, 8, 8, 8);
        }

        foreach (var texture in _velocityTextures)
        {
            BindImageTexture(0, texture, TextureAccess.WriteOnly);
            GL.Uniform1(GL.GetUniformLocation(_clearShader, "value"), 0.0f);
            DispatchCompute(_fluidSizeX, _fluidSizeY, _fluidSizeZ, 8, 8, 8);
        }

        foreach (var texture in _pressureTextures)
        {
            BindImageTexture(0, texture, TextureAccess.WriteOnly);
            GL.Uniform1(GL.GetUniformLocation(_clearShader, "value"), 0.0f);
            DispatchCompute(_fluidSizeX, _fluidSizeY, _fluidSizeZ, 8, 8, 8);
        }

        GL.MemoryBarrier(MemoryBarrierFlags.AllBarrierBits);

        // Re-add initial fluid
        AddInitialFluid();
    }

    protected override void OnUpdateFrame(FrameEventArgs e)
    {
        base.OnUpdateFrame(e);

        _time += (float)e.Time;

        _camera.UpdateVectors();

        var input = KeyboardState;
        float cameraSpeed = _camera.Speed * (float)e.Time;

        if (input.IsKeyDown(Keys.W))
            _camera.Position += _camera.Front * cameraSpeed;
        if (input.IsKeyDown(Keys.S))
            _camera.Position -= _camera.Front * cameraSpeed;
        if (input.IsKeyDown(Keys.A))
            _camera.Position -= _camera.Right * cameraSpeed;
        if (input.IsKeyDown(Keys.D))
            _camera.Position += _camera.Right * cameraSpeed;
        if (input.IsKeyDown(Keys.Space))
            _camera.Position += _camera.Up * cameraSpeed;
        if (input.IsKeyDown(Keys.LeftShift))
            _camera.Position -= _camera.Up * cameraSpeed;

        Title = $"🎯 WORKING Fluid Simulation - FPS: {1.0 / e.Time:0} - CLICK TO INTERACT!";
    }

    protected override void OnMouseMove(MouseMoveEventArgs e)
    {
        base.OnMouseMove(e);

        if (IsFocused)
        {
            if (_firstMove)
            {
                _lastPos = new Vector2(e.X, e.Y);
                _firstMove = false;
            }
            else
            {
                var deltaX = e.X - _lastPos.X;
                var deltaY = e.Y - _lastPos.Y;
                _lastPos = new Vector2(e.X, e.Y);

                _camera.Yaw += deltaX * _camera.Sensitivity;
                _camera.Pitch -= deltaY * _camera.Sensitivity;

                _camera.Pitch = Math.Clamp(_camera.Pitch, -MathHelper.PiOver2 + 0.1f, MathHelper.PiOver2 - 0.1f);
            }
        }
    }

    protected override void OnMouseDown(MouseButtonEventArgs e)
    {
        base.OnMouseDown(e);

        var mousePos = new Vector2(MouseState.X, MouseState.Y);

        // Simple direct mapping
        Vector3 fluidPos = new Vector3(
            (mousePos.X / Size.X) * _fluidSizeX,
            (1.0f - mousePos.Y / Size.Y) * _fluidSizeY,
            _fluidSizeZ / 2
        );

        // Clamp to reasonable bounds
        fluidPos.X = Math.Clamp(fluidPos.X, 5, _fluidSizeX - 5);
        fluidPos.Y = Math.Clamp(fluidPos.Y, 5, _fluidSizeY - 5);
        fluidPos.Z = Math.Clamp(fluidPos.Z, 5, _fluidSizeZ - 5);

        if (e.Button == MouseButton.Left)
        {
            AddForce(fluidPos, Vector3.Zero, 8.0f, 15.0f);
            Console.WriteLine($"💧 Added DENSITY at: {fluidPos}");
        }
        else if (e.Button == MouseButton.Right)
        {
            // Add velocity away from center for explosion effect
            Vector3 center = new Vector3(_fluidSizeX / 2, _fluidSizeY / 2, _fluidSizeZ / 2);
            Vector3 dir = Vector3.Normalize(fluidPos - center);
            AddForce(fluidPos, dir * 20.0f, 6.0f, 0.0f);
            Console.WriteLine($"💨 Added VELOCITY at: {fluidPos}");
        }
    }

    protected override void OnKeyDown(KeyboardKeyEventArgs e)
    {
        base.OnKeyDown(e);

        if (e.Key == Keys.Escape)
            Close();
        else if (e.Key == Keys.R)
        {
            ClearAllTextures();
            Console.WriteLine("🔄 Simulation RESET!");
        }
        else if (e.Key == Keys.O)
        {
            var mousePos = new Vector2(MouseState.X, MouseState.Y);
            Vector3 fluidPos = new Vector3(
                (mousePos.X / Size.X) * _fluidSizeX,
                (1.0f - mousePos.Y / Size.Y) * _fluidSizeY,
                _fluidSizeZ / 2
            );
            AddSphereObstacle(fluidPos, 5.0f);
            Console.WriteLine($"🧱 Added OBSTACLE at: {fluidPos}");
        }
    }

    protected override void OnResize(ResizeEventArgs e)
    {
        base.OnResize(e);
        GL.Viewport(0, 0, Size.X, Size.Y);
        _camera.AspectRatio = Size.X / (float)Size.Y;
    }
}