#include "fiko/dokum/porozite.hpp"

#include <algorithm>
#include <cmath>
#include <deque>

namespace fiko::dokum {
namespace {

struct Yon {
  i32 di, dj, dk;
};
constexpr Yon kYonlar[6] = {{1, 0, 0}, {-1, 0, 0}, {0, 1, 0},
                            {0, -1, 0}, {0, 0, 1}, {0, 0, -1}};

// Alaşım-aile kritik Niyama [(°C·s)^½/mm] — likidusa göre aile seçimi.
// Çelik ~1.0 (Niyama 1982 orijinal eşiği), dökme demir ~0.6 (grafit genleşmesi
// büzülmeyi kısmen telafi eder), Al ~0.5 (Carlson-Beckermann 2009 Ny* eğrisi
// eşdeğeri). İnce ayar PoroziteAyar::nyOlcek ile.
f64 nyKritik(const Alasim& a) {
  if (a.tLikidus > 1350.0) return 1.0;
  if (a.tLikidus > 1000.0) return 0.6;
  return 0.5;
}

}  // namespace

f64 cekmePorozitesi(f64 niyama, const Alasim& alasim, f64 nyOlcek) {
  // Tavan: toplam katılaşma büzülmesi β/(1+β) hacim % (kütle korunumu).
  const f64 tavan =
      std::max(3.0, 100.0 * alasim.cekmeOrani / (1.0 + alasim.cekmeOrani));
  const f64 nyStar = std::max(0.0, niyama) * std::max(nyOlcek, 1e-9) /
                     nyKritik(alasim);
  // Parçalı-lineer Ny* → porozite % (Carlson-Beckermann 2009 eğri vekili).
  if (nyStar >= 2.0) return 0.0;
  if (nyStar >= 1.0) return 0.3 * (2.0 - nyStar);
  if (nyStar >= 0.7) return 0.3 + (1.0 - nyStar) / 0.3 * 0.7;
  if (nyStar >= 0.1) return 1.0 + (0.7 - nyStar) / 0.6 * 2.0;
  return 3.0 + (0.1 - nyStar) / 0.1 * (tavan - 3.0);
}

PoroziteSonuc poroziteHesapla(const VoxelKafes& kafes, const KatilasmaSonuc& kat,
                              const HavaAnalizSonuc* hava, const Alasim& alasim,
                              const PoroziteAyar& ayar) {
  PoroziteSonuc s;
  const u32 N = kafes.hucreSayi();
  if (!kafes.gecerliMi() || !kat.hazir || kat.niyama.size() != N ||
      kafes.tip.size() != N)
    return s;
  const bool havaVar = hava && hava->hazir && hava->lastAir.size() == N &&
                       hava->kabarcikRisk.size() == N;

  const u32 nx = kafes.nx(), ny = kafes.ny(), nz = kafes.nz();
  auto metal = [&](u32 idx) { return kafes.tip[idx] != HucreTip::Dis; };
  auto koord = [&](u32 idx, u32& i, u32& j, u32& k) {
    i = idx % nx;
    j = (idx / nx) % ny;
    k = idx / (nx * ny);
  };
  auto komsu = [&](u32 idx, const Yon& y) -> u32 {
    u32 i, j, k;
    koord(idx, i, j, k);
    const i64 ni = i64(i) + y.di, nj = i64(j) + y.dj, nk = i64(k) + y.dk;
    if (ni < 0 || nj < 0 || nk < 0 || ni >= i64(nx) || nj >= i64(ny) ||
        nk >= i64(nz))
      return static_cast<u32>(-1);
    return kafes.indeks(u32(ni), u32(nj), u32(nk));
  };

  // ---- Per-hücre bileşenler: çekme (Niyama) + gaz (hava analizi) ----------
  s.porozite.assign(N, 0.0f);
  s.cekmePor.assign(N, 0.0f);
  s.gazPor.assign(N, 0.0f);
  const f64 hucreMm3 = kafes.hucreHacim() * 1e9;  // m³ → mm³
  u32 metalSayi = 0;
  for (u32 c = 0; c < N; ++c) {
    if (!metal(c)) continue;
    ++metalSayi;
    const f64 cekme = cekmePorozitesi(kat.niyama[c], alasim, ayar.nyOlcek);
    f64 gaz = 0;
    if (havaVar) {
      // Dolmayan hücre (lastAir=1.0 tam işareti) = hapsolmuş hava vekili
      // (HavaCebi hücre listesi taşımaz — per-hücre lastAir alanı kullanılır).
      gaz = ayar.gazKatsayi * f64(hava->kabarcikRisk[c]);
      if (hava->lastAir[c] > 0.999f) gaz += ayar.kapaliCepPor;
    }
    const f64 toplam = std::clamp(cekme + gaz, 0.0, 100.0);
    s.cekmePor[c] = f32(cekme);
    s.gazPor[c] = f32(std::min(gaz, 100.0));
    s.porozite[c] = f32(toplam);
    s.toplamPorHacimMm3 += toplam * 0.01 * hucreMm3;
    s.maksPor = std::max(s.maksPor, toplam);
  }
  s.parcaHacimMm3 = f64(metalSayi) * hucreMm3;
  s.poroziteOrani =
      s.parcaHacimMm3 > 0 ? 100.0 * s.toplamPorHacimMm3 / s.parcaHacimMm3 : 0;

  // ---- Baloncuk bölgeleri: porozite > esik bağlı bileşenleri --------------
  const f32 esik = f32(ayar.esik);
  auto aday = [&](u32 c) { return metal(c) && s.porozite[c] > esik; };
  std::vector<u8> ziyaret(N, 0);
  s.bolgeId.assign(N, -1);
  for (u32 c = 0; c < N; ++c) {
    if (ziyaret[c] || !aday(c)) continue;
    const i32 yeniId = static_cast<i32>(s.bolgeler.size());
    std::vector<u32> hucreler;
    std::deque<u32> kuyruk{c};
    ziyaret[c] = 1;
    while (!kuyruk.empty()) {
      const u32 a = kuyruk.front();
      kuyruk.pop_front();
      hucreler.push_back(a);
      s.bolgeId[a] = yeniId;
      for (const Yon& y : kYonlar) {
        const u32 nb = komsu(a, y);
        if (nb == static_cast<u32>(-1) || ziyaret[nb] || !aday(nb)) continue;
        ziyaret[nb] = 1;
        kuyruk.push_back(nb);
      }
    }
    PoroziteBolge b;
    Vec3d toplamM{0, 0, 0};
    f64 porToplam = 0, cekmeT = 0, gazT = 0;
    for (u32 idx : hucreler) {
      u32 i, j, k;
      koord(idx, i, j, k);
      toplamM = toplamM + kafes.merkez(i, j, k);
      porToplam += s.porozite[idx];
      cekmeT += s.cekmePor[idx];
      gazT += s.gazPor[idx];
      b.maksPor = std::max(b.maksPor, f64(s.porozite[idx]));
    }
    const f64 n = f64(hucreler.size());
    b.merkez = Vec3f{f32(toplamM.x / n), f32(toplamM.y / n), f32(toplamM.z / n)};
    b.hucreSayi = u32(hucreler.size());
    b.hacimMm3 = n * hucreMm3;
    b.ortalamaPor = porToplam / n;
    const f64 kucuk = std::min(cekmeT, gazT), buyuk = std::max(cekmeT, gazT);
    b.tip = (buyuk > 0 && kucuk > 0.25 * buyuk) ? PoroziteTip::Karisik
            : (cekmeT >= gazT ? PoroziteTip::Cekme : PoroziteTip::Gaz);
    s.bolgeler.push_back(b);
  }

  s.hazir = true;
  return s;
}

f64 poroziteOtoEsik(const std::vector<f32>& porozite, f64 yuzdelik,
                    f64 asgari, f64 azami) {
  if (azami < asgari) std::swap(asgari, azami);
  std::vector<f32> dolu;
  dolu.reserve(porozite.size());
  for (f32 p : porozite)
    if (p > 0.0f && std::isfinite(p)) dolu.push_back(p);
  if (dolu.empty()) return asgari;
  const f64 oran = std::clamp(yuzdelik, 0.0, 1.0);
  size_t n = static_cast<size_t>(oran * (dolu.size() - 1));
  std::nth_element(dolu.begin(), dolu.begin() + n, dolu.end());
  return std::clamp(static_cast<f64>(dolu[n]), asgari, azami);
}

}  // namespace fiko::dokum