<!DOCTYPE html>
<html>
<head>
    <title>PIC Fluid Simulation</title>
    <style>
        canvas { border: 1px solid black; }
        body { margin: 0; display: flex; justify-content: center; background: #f0f0f0; }
    </style>
</head>
<body>
    <canvas id="fluidCanvas"></canvas>

    <script>
        const canvas = document.getElementById('fluidCanvas');
        const ctx = canvas.getContext('2d');

        // Simulation parameters
        const SIM_WIDTH = 800;
        const SIM_HEIGHT = 600;
        const H = 20; // Grid spacing
        const DT = 0.1;
        const GRAVITY = 0.5;
        const NUM_PARTICLES = 5000;

        // Grid dimensions
        const COLS = Math.floor(SIM_WIDTH / H);
        const ROWS = Math.floor(SIM_HEIGHT / H);

        // Initialize grid
        let grid = {
            u: new Float32Array((COLS+1) * ROWS),   // Horizontal velocity
            v: new Float32Array(COLS * (ROWS+1)),   // Vertical velocity
            type: new Array(COLS * ROWS).fill(0),   // 0=air, 1=fluid, 2=solid
            numParts: new Float32Array(COLS * ROWS) // Particle count
        };

        // Initialize particles
        let particles = Array(NUM_PARTICLES).fill().map(() => ({
            pos: [H/2 + Math.random()*(SIM_WIDTH-H), H/2 + Math.random()*(SIM_HEIGHT/2)],
            vel: [0, 0]
        }));

        // Mouse interaction
        let mouse = { x: -100, y: -100, pressed: false };
        canvas.onmousemove = e => {
            const rect = canvas.getBoundingClientRect();
            mouse.x = e.clientX - rect.left;
            mouse.y = e.clientY - rect.top;
        };
        canvas.onmousedown = () => mouse.pressed = true;
        canvas.onmouseup = () => mouse.pressed = false;

        // Bilinear interpolation weights
        function getWeights(x, y) {
            const ix = Math.floor(x / H);
            const iy = Math.floor(y / H);
            const dx = (x - ix * H) / H;
            const dy = (y - iy * H) / H;
            return {
                ix, iy,
                weights: [
                    (1 - dx) * (1 - dy),
                    dx * (1 - dy),
                    (1 - dx) * dy,
                    dx * dy
                ]
            };
        }

        // Transfer particle velocities to grid
        function particlesToGrid() {
            grid.u.fill(0);
            grid.v.fill(0);
            grid.numParts.fill(0);

            particles.forEach(p => {
                // Horizontal velocity (u component)
                const { ix, iy, weights } = getWeights(p.pos[0], p.pos[1] - H/2);
                for (let j = 0; j < 4; j++) {
                    const wx = j % 2;
                    const wy = Math.floor(j / 2);
                    const idx = (ix + wx) + (iy + wy) * (COLS + 1);
                    if (idx >= 0 && idx < grid.u.length) {
                        grid.u[idx] += weights[j] * p.vel[0];
                        grid.numParts[ix + iy * COLS] += weights[j];
                    }
                }

                // Vertical velocity (v component) - FIXED HERE
                const { ix: ixv, iy: iyv, weights: weightsV } = getWeights(p.pos[0] - H/2, p.pos[1]);
                for (let j = 0; j < 4; j++) {
                    const wx = j % 2;
                    const wy = Math.floor(j / 2);
                    const idx = (ixv + wx) + (iyv + wy) * COLS;
                    if (idx >= 0 && idx < grid.v.length) {
                        grid.v[idx] += weightsV[j] * p.vel[1];
                    }
                }
            });

            // Normalize grid velocities
            for (let i = 0; i < grid.u.length; i++) {
                const count = grid.numParts[Math.floor(i/(COLS+1)) * COLS + i%(COLS+1)];
                if (count > 0) grid.u[i] /= count;
            }
            for (let i = 0; i < grid.v.length; i++) {
                const count = grid.numParts[i % COLS + Math.floor(i/COLS) * COLS];
                if (count > 0) grid.v[i] /= count;
            }
        }

        // Projection step to enforce incompressibility
        function project() {
            const div = new Float32Array(COLS * ROWS);
            const p = new Float32Array(COLS * ROWS);

            // Calculate divergence
            for (let i = 0; i < COLS; i++) {
                for (let j = 0; j < ROWS; j++) {
                    if (grid.type[i + j * COLS] === 2) continue;
                    div[i + j * COLS] = (grid.u[i+1 + j*(COLS+1)] - grid.u[i + j*(COLS+1)] +
                                         grid.v[i + (j+1)*COLS] - grid.v[i + j*COLS]) / H;
                }
            }

            // Solve pressure using Gauss-Seidel
            for (let iter = 0; iter < 20; iter++) {
                for (let i = 0; i < COLS; i++) {
                    for (let j = 0; j < ROWS; j++) {
                        if (grid.type[i + j * COLS] !== 1) continue;
                        let sum = 0, count = 0;
                        if (i > 0) { sum += p[i-1 + j*COLS]; count++; }
                        if (i < COLS-1) { sum += p[i+1 + j*COLS]; count++; }
                        if (j > 0) { sum += p[i + (j-1)*COLS]; count++; }
                        if (j < ROWS-1) { sum += p[i + (j+1)*COLS]; count++; }
                        p[i + j*COLS] = (sum - div[i + j*COLS] * H*H) / count;
                    }
                }
            }

            // Apply pressure gradient
            for (let i = 0; i < COLS; i++) {
                for (let j = 0; j < ROWS; j++) {
                    if (grid.type[i + j * COLS] !== 1) continue;
                    grid.u[i + j*(COLS+1)] -= (p[i + j*COLS] - (i>0 ? p[i-1 + j*COLS] : 0)) / H;
                    grid.v[i + j*COLS] -= (p[i + j*COLS] - (j>0 ? p[i + (j-1)*COLS] : 0)) / H;
                }
            }
        }

        // Transfer grid velocities back to particles
        function gridToParticles() {
            particles.forEach(p => {
                // Horizontal velocity
                const { ix, iy, weights } = getWeights(p.pos[0], p.pos[1] - H/2);
                p.vel[0] = 0;
                for (let j = 0; j < 4; j++) {
                    const wx = j % 2;
                    const wy = Math.floor(j / 2);
                    const idx = (ix + wx) + (iy + wy) * (COLS + 1);
                    if (idx >= 0 && idx < grid.u.length) {
                        p.vel[0] += grid.u[idx] * weights[j];
                    }
                }

                // Vertical velocity - FIXED HERE
                const { ix: ixv, iy: iyv, weights: weightsV } = getWeights(p.pos[0] - H/2, p.pos[1]);
                p.vel[1] = 0;
                for (let j = 0; j < 4; j++) {
                    const wx = j % 2;
                    const wy = Math.floor(j / 2);
                    const idx = (ixv + wx) + (iyv + wy) * COLS;
                    if (idx >= 0 && idx < grid.v.length) {
                        p.vel[1] += grid.v[idx] * weightsV[j];
                    }
                }
            });
        }

        // Main simulation loop
        function update() {
            // Add mouse interaction
            if (mouse.pressed) {
                particles.push({
                    pos: [mouse.x, mouse.y],
                    vel: [Math.random()*2-1, Math.random()*2-1]
                });
            }

            // Update particles
            particles.forEach(p => {
                p.vel[1] += GRAVITY * DT;
                p.pos[0] += p.vel[0] * DT;
                p.pos[1] += p.vel[1] * DT;

                // Basic boundary collision
                if (p.pos[0] < 0) { p.pos[0] = 0; p.vel[0] *= -0.5; }
                if (p.pos[0] > SIM_WIDTH) { p.pos[0] = SIM_WIDTH; p.vel[0] *= -0.5; }
                if (p.pos[1] > SIM_HEIGHT) { p.pos[1] = SIM_HEIGHT; p.vel[1] *= -0.5; }
            });

            particlesToGrid();
            project();
            gridToParticles();

            // Draw
            ctx.fillStyle = '#87CEEB';
            ctx.fillRect(0, 0, SIM_WIDTH, SIM_HEIGHT);

            particles.forEach(p => {
                ctx.fillStyle = `rgba(30, 144, 255, ${Math.min(1, 0.1 + p.pos[1]/SIM_HEIGHT)})`;
                ctx.beginPath();
                ctx.arc(p.pos[0], p.pos[1], 2, 0, Math.PI*2);
                ctx.fill();
            });

            requestAnimationFrame(update);
        }

        // Start simulation
        canvas.width = SIM_WIDTH;
        canvas.height = SIM_HEIGHT;
        update();
    </script>
</body>
</html>
