#version 450
// ============================================================================
// sph_entegre.comp — PBF zaman entegrasyonu (iki fazlı; faz = push const).
// ----------------------------------------------------------------------------
// Macklin & Müller 2013, Algorithm 1:
//   faz=0 TAHMİN (substep başı): x0 = x ; v += g·dt ; x += v·dt
//     (kuvvet = yalnız gravite; basınç PBF konum-düzeltmesiyle uygulanır)
//   faz=1 FİNALİZE (düzeltme iterasyonları sonrası):
//     v = (x − x0)/dt  → sonra XSPH viskozite: v_i += c Σ_j (v_j−v_i) W(x_ij)
//     (Schechter&Bridson 2012 / Macklin§4; sıvı metal düşük visk → küçük c)
// Komşu (XSPH): hash grid 3×3×3 tarama (faz=1'de).
// ============================================================================
#include "sph_ortak.glsl"

layout(local_size_x = 64) in;

layout(std430, binding = 0) buffer Konum          { vec4 konum[]; };
layout(std430, binding = 1) buffer Hiz            { vec4 hiz[]; };
layout(std430, binding = 2) buffer EskiKonum      { vec4 eskiKonum[]; };
layout(std430, binding = 3) readonly buffer Lambda{ float lambda[]; };
layout(std430, binding = 4) readonly buffer HucreSay   { uint hucreSay[]; };
layout(std430, binding = 5) readonly buffer ParcaHucre { uint parcaHucre[]; };

void main() {
  uint i = gl_GlobalInvocationID.x;
  if (i >= itme.n) return;
  if (konum[i].w < 0.5) return;

  if (itme.faz == 0u) {
    // TAHMİN: gravite entegrasyonu + öngörülen konum.
    eskiKonum[i].xyz = konum[i].xyz;
    vec3 v = hiz[i].xyz + vec3(itme.gx, itme.gy, itme.gz) * itme.dt;
    hiz[i].xyz = v;
    konum[i].xyz = konum[i].xyz + v * itme.dt;
    return;
  }

  // FİNALİZE: hızı konum değişiminden türet.
  vec3 xi = konum[i].xyz;
  vec3 v = (xi - eskiKonum[i].xyz) / itme.dt;

  // XSPH viskozite (yumuşatma) — hız gürültüsünü kırar, dolum sakinleşir.
  if (itme.xsphC > 0.0) {
    ivec3 c0 = hashHucre(xi);
    vec3 dv = vec3(0.0);
    // KOHESYON (akışkan-akışkan çekim) — PBF'de yüzey geriliminin karşılığı.
    // Akinci2013 "cohesion/surface tension model": komşu akışkan partikülleri
    // birbirine ÇEKİLİR (yoğunluk gradyanı yönünde). σ = 0 iken sıvı düz
    // duzarlara yayılır (sise benzer) ve jet parçalanmaz. σ > 0 iken serbest
    // yüzey KÜRESELLEŞİR → sise kalkar, dar boğazda jet KOPAR, damla oluşur,
    // duvarda menisk yapar (gerçek erimiş metal).
    // Kuvvet: F = -σ · m · Σ_j r̂_ij · C(r_ij),  spf = -(h−r)³/r³ (Basri 2013)
    // Kapsama (renk alanından türetilen düşük yoğunluk) uygulanmaz — basit ve
    // kararlı; kullanıcı katsayısıyla ölçekler.
    vec3 fc = vec3(0.0);
    for (int dz = -1; dz <= 1; ++dz)
    for (int dy = -1; dy <= 1; ++dy)
    for (int dx = -1; dx <= 1; ++dx) {
      uint cell = hashDuz(c0 + ivec3(dx, dy, dz));
      uint say = min(hucreSay[cell], itme.bucketK);
      for (uint s = 0u; s < say; ++s) {
        uint j = parcaHucre[cell * itme.bucketK + s];
        if (j == i || konum[j].w < 0.5) continue;
        vec3 rij = xi - konum[j].xyz;
        float r2 = dot(rij, rij);
        float r = sqrt(r2);
        if (r >= itme.h || r < 1e-6) continue;
        vec3 vj = (konum[j].xyz - eskiKonum[j].xyz) / itme.dt;
        dv += (vj - v) * poly6(r2, itme.h) * itme.kutle / itme.rho0;
        if (itme.sigmaC > 0.0) {
          // Sere kuvveti: r arttikça zayıflar, r→h'de 0. Kısa menzilde baskın
          // (yüzey kümelenmesi), uzakta ihmal (hacim içi kendini götürür).
          float t = r / itme.h;                    // 0..1
          float spf = (1.0 - t);  spf = spf * spf * spf;
          fc += (rij / r) * spf;
        }
      }
    }
    v += itme.xsphC * dv;
    // Kohesyon ivmesi (çekim yönüne BAĞIMSIZ — simetrik kuvvet).
    if (itme.sigmaC > 0.0) v += itme.sigmaC * fc * itme.dt;
  }
  hiz[i].xyz = v;
}