"""E06 암반 비탈면 — Markland 판정 · Hoek-Bray 평면 · 쐐기 안전율 (요구 345 · 351 · PLAN 536-1).

근거 = 지식DB `original/안정성평가/03`(Markland 부등식 · Hoek-Bray 평면 · 쐐기) · 06(상용 SW) ·
Rocscience RocPlane · SWedge 이론서(힘 벡터 · 미끄럼 방향 8 경우 · 수직력 · 지지력은 수동력).
코드는 직접 짬(Dips · RocPlane · SWedge 코드 안 씀) · NumPy · SciPy 만.

좌표 = 동(+x) · 북(+y) · 위(+z). 면 = (경사방향 α · 경사 ψ) · 선 = (추향 · 선주각).
단위 = 길이 m · 단위중량 kN/m³ · 점착력 kPa · 각 ° · 힘 kN(평면 = 단위 폭당 kN/m).

⚠ 지식DB 03 의 평면 볼트 항(수직 T cos(θ+ψp) · 활동 −T sin(θ+ψp))은 θ = 볼트와 면 법선 사이 각일
때 맞음 — 여기는 θ = 볼트 경사(수평 아래)로 받아 힘 벡터로 풂(수직 T sin(θ+ψp) · 활동 −T cos(θ+ψp)).
"""

from __future__ import annotations

import math
from typing import Any

import numpy as np
from scipy.optimize import minimize_scalar

GAMMA_W = 9.81  # 물 단위중량 kN/m³ (E03 Types 와 같음)
LATERAL_PLANAR_DEG = 20.0  # 평면 측방 한계(Dips 기본 · Wyllie & Mah)
LATERAL_TOPPLING_DEG = 20.0  # 전도 측방 한계(Dips 기본 · 지식DB 03 은 30 — 화면에서 바꿈)


# ---- 방향 ----------------------------------------------------------------------
def ang_diff(a: float, b: float) -> float:
    """두 방위각 사이 최소 각거리 0~180."""
    d = abs((a - b) % 360.0)
    return min(d, 360.0 - d)


def plane_normal(dip_dir: float, dip: float) -> np.ndarray:
    """면의 위쪽 단위 법선."""
    a, p = math.radians(dip_dir), math.radians(dip)
    return np.array([math.sin(a) * math.sin(p), math.cos(a) * math.sin(p), math.cos(p)])


def line_vector(trend: float, plunge: float) -> np.ndarray:
    """선(추향 · 선주각) → 아래로 향한 단위 벡터."""
    t, p = math.radians(trend), math.radians(plunge)
    return np.array([math.sin(t) * math.cos(p), math.cos(t) * math.cos(p), -math.sin(p)])


def vector_line(v: np.ndarray) -> tuple[float, float]:
    """벡터 → (추향 · 선주각) 하반구(위로 향하면 뒤집음)."""
    v = np.asarray(v, float)
    if v[2] > 0:
        v = -v
    v = v / np.linalg.norm(v)
    trend = math.degrees(math.atan2(v[0], v[1])) % 360.0
    return trend, math.degrees(math.asin(max(-1.0, min(1.0, -v[2]))))


def pole(dip_dir: float, dip: float) -> tuple[float, float]:
    """면의 극점(하반구) = (α + 180 · 90 − ψ)."""
    return (dip_dir + 180.0) % 360.0, 90.0 - dip


def intersection(a: dict, b: dict) -> tuple[float, float] | None:
    """두 면 교선(추향 · 선주각) — 나란하면 None."""
    v = np.cross(
        plane_normal(a["dip_direction"], a["dip"]), plane_normal(b["dip_direction"], b["dip"])
    )
    if np.linalg.norm(v) < 1e-9:
        return None
    return vector_line(v)


def apparent_dip(dip: float, dip_dir: float, trend: float) -> float:
    """면 (α · ψ) 의 trend 방향 겉보기 경사 — 반대쪽(90° 넘음)이면 음수."""
    return math.degrees(
        math.atan(math.tan(math.radians(dip)) * math.cos(math.radians(trend - dip_dir)))
    )


def strike_to_dip_direction(strike: float) -> float:
    """오른손 법칙 주향 → 경사방향."""
    return (strike + 90.0) % 360.0


# ---- Markland 판정 -------------------------------------------------------------
def kinematic(
    slope: dict,
    joints: list[dict],
    lateral_planar: float = LATERAL_PLANAR_DEG,
    lateral_toppling: float = LATERAL_TOPPLING_DEG,
) -> dict[str, list]:
    """slope = {dip_direction, dip} · joints = [{id, dip_direction, dip, phi}, …].

    평면 = Δ(αf, αp) ≤ 측방 · φ < ψp < 겉보기 사면 경사 / 쐐기 = 교선 Δ < 90 · φ < 선주각 <
    겉보기 사면 경사(φ = 두 면 중 작은 값) /
    전도 = Δ(αf + 180, αp) ≤ 측방 · ψp > 90 − ψf + φ (Goodman).
    """
    af, pf = float(slope["dip_direction"]), float(slope["dip"])
    planar, toppling, wedge = [], [], []
    for j in joints:
        aj, pj, phi = float(j["dip_direction"]), float(j["dip"]), float(j["phi"])
        if ang_diff(af, aj) <= lateral_planar and phi < pj < apparent_dip(pf, af, aj):
            planar.append(j["id"])
        if ang_diff(af + 180.0, aj) <= lateral_toppling and pj > 90.0 - pf + phi:
            toppling.append(j["id"])
    for i, a in enumerate(joints):
        for b in joints[i + 1 :]:
            line = intersection(a, b)
            if line is None:
                continue
            trend, plunge = line
            phi = min(float(a["phi"]), float(b["phi"]))
            if ang_diff(af, trend) < 90.0 and phi < plunge < apparent_dip(pf, af, trend):
                wedge.append(
                    {
                        "ids": [a["id"], b["id"]],
                        "trend": round(trend, 2),
                        "plunge": round(plunge, 2),
                    }
                )
    return {"planar": planar, "wedge": wedge, "toppling": toppling}


# ---- 평면 파괴 (Hoek-Bray · 2D 단위 폭) -----------------------------------------
def _plane_block(H: float, pf: float, pp: float, ps: float, crack_x: float | None):
    """사면 끝(0,0) · x = 산 쪽 수평 · z = 위. 덩이 다각형 · 활동면 끝 x · 균열 깊이.

    윗면 = 마루(H cotψf, H)에서 ψs 로 오름. 균열 없음 = 활동면이 윗면과 만나는 점까지."""
    tf, tp, ts = (math.tan(math.radians(v)) for v in (pf, pp, ps))
    xc = H / tf  # 마루

    def top(x: float) -> float:
        return x * tf if x <= xc else H + (x - xc) * ts

    if crack_x is None:
        if ts >= tp - 1e-12:
            return None  # 활동면이 윗면과 안 만남 = 덩이 안 생김
        xe = (H - xc * ts) / (tp - ts)
        poly = [(0.0, 0.0), (xc, H), (xe, xe * tp)]
        return poly, xe, 0.0
    if crack_x < xc or top(crack_x) - crack_x * tp <= 1e-9:
        return None
    poly = [(0.0, 0.0), (xc, H), (crack_x, top(crack_x)), (crack_x, crack_x * tp)]
    return poly, crack_x, top(crack_x) - crack_x * tp


def _area(poly) -> float:
    xs, zs = np.array(poly).T
    return 0.5 * abs(float(np.dot(xs, np.roll(zs, -1)) - np.dot(zs, np.roll(xs, -1))))


def plane_fs(
    *,
    H: float,
    slope_dip: float,
    plane_dip: float,
    c: float,
    phi: float,
    gamma: float,
    upper_dip: float = 0.0,
    crack: dict | None = None,
    water_fill: float = 0.0,
    kh: float = 0.0,
    bolt_t: float = 0.0,
    bolt_plunge: float = 0.0,
    gamma_w: float = GAMMA_W,
) -> dict[str, Any] | None:
    """평면 파괴 안전율 — 덩이가 안 생기면 None.

    crack = None | {"mode": "depth", "z": m} | {"mode": "critical"}(안전율 최소 자리 찾음) ·
    물 = 균열 있음: zw = fill·z · U = ½γw zw A · V = ½γw zw² / 균열 없음: 가운데 최대 수두
    fill·h/2(h = 활동면 오름) · U = ½γw (fill·h/2) A (Hoek · RocPlane) ·
    kh = 바깥 수평 · 볼트 = 수동력(산 쪽 · 수평 아래 bolt_plunge°)."""
    if not 0 < plane_dip < slope_dip < 90 or H <= 0 or upper_dip >= slope_dip:
        return None
    mode = (crack or {}).get("mode", "none")
    if mode == "critical":
        return _critical_crack(locals())
    tp = math.tan(math.radians(plane_dip))
    xc = H / math.tan(math.radians(slope_dip))
    crack_x = None
    if mode == "depth":
        z = float(crack["z"])
        crack_x = _crack_x_for_depth(H, slope_dip, plane_dip, upper_dip, z)
        if crack_x is None:
            return None
    elif mode == "x":
        crack_x = xc + float(crack["b"])
    block = _plane_block(H, slope_dip, plane_dip, upper_dip, crack_x)
    if block is None:
        return None
    poly, x_end, z = block
    psi = math.radians(plane_dip)
    A = x_end / math.cos(psi)
    W = gamma * _area(poly)
    if crack_x is not None:
        zw = water_fill * z
        U, V = 0.5 * gamma_w * zw * A, 0.5 * gamma_w * zw * zw
    else:
        zw = water_fill * (x_end * tp) / 2.0
        U, V = 0.5 * gamma_w * zw * A, 0.0
    d = np.array([-math.cos(psi), -math.sin(psi)])  # 아래로 미끄러짐(바깥 = −x)
    n = np.array([-math.sin(psi), math.cos(psi)])  # 면에서 덩이 쪽
    active = np.array([-kh * W - V, -W]) + U * n
    drive = float(active @ d)
    if drive <= 0:
        return None
    bolt = bolt_t * np.array(
        [math.cos(math.radians(bolt_plunge)), -math.sin(math.radians(bolt_plunge))]
    )
    N = max(0.0, -float((active + bolt) @ n))
    resist = c * A + N * math.tan(math.radians(phi)) - float(bolt @ d)
    return {
        "F": resist / drive,
        "W": W,
        "A": A,
        "U": U,
        "V": V,
        "z": z,
        "zw": zw,
        "crack_b": None if crack_x is None else crack_x - xc,
        "N": N,
        "drive": drive,
        "resist": resist,
        "polygon": [[round(x, 3), round(y, 3)] for x, y in poly],
    }


def _crack_x_for_depth(H, pf, pp, ps, z) -> float | None:
    """윗면 아래 균열 깊이 z 가 되는 자리(마루 뒤) — 깊이가 마루 틈보다 크면 None."""
    tf, tp, ts = (math.tan(math.radians(v)) for v in (pf, pp, ps))
    xc = H / tf
    gap = H - xc * tp  # 마루에서 윗면 − 활동면
    if z <= 0 or z >= gap:
        return None
    if ts >= tp:  # 틈이 커지기만 함 = 뒤로 갈수록 깊어짐
        return xc + (z - gap) / (ts - tp) if ts > tp else None
    return xc + (gap - z) / (tp - ts)


def plane_critical_depth_hoek(H: float, slope_dip: float, plane_dip: float) -> float:
    """Hoek-Bray 임계 인장균열 깊이(윗면 수평 · 건조) z = H(1 − √(cotψf tanψp))."""
    return H * (
        1.0 - math.sqrt(math.tan(math.radians(plane_dip)) / math.tan(math.radians(slope_dip)))
    )


def _critical_crack(kw: dict) -> dict | None:
    """균열 자리를 마루 뒤로 옮겨 안전율 최소 — 균열 없음보다 크면 균열 없음."""
    args = {k: kw[k] for k in kw if k not in ("crack", "mode")}
    tf, tp, ts = (math.tan(math.radians(kw[k])) for k in ("slope_dip", "plane_dip", "upper_dip"))
    xc = kw["H"] / tf
    x_far = (kw["H"] - xc * ts) / (tp - ts) if ts < tp else xc + 10 * kw["H"]
    span = x_far - xc

    def fs_at(b: float) -> float:
        r = plane_fs(**args, crack={"mode": "x", "b": b})
        return r["F"] if r else math.inf

    best = minimize_scalar(fs_at, bounds=(span * 1e-3, span * (1 - 1e-3)), method="bounded")
    found = plane_fs(**args, crack={"mode": "x", "b": float(best.x)})
    plain = plane_fs(**args, crack=None)
    if found is None or (plain is not None and plain["F"] <= found["F"]):
        return plain
    return found


# ---- 쐐기 파괴 (Hoek-Bray · SWedge 꼴 사면체) ------------------------------------
def _meet(planes: list[tuple[np.ndarray, float]]) -> np.ndarray | None:
    """세 면 n·x = d 의 교점."""
    m = np.array([p[0] for p in planes])
    if abs(np.linalg.det(m)) < 1e-10:
        return None
    return np.linalg.solve(m, np.array([p[1] for p in planes]))


def _tri_area(a, b, c) -> float:
    return 0.5 * float(np.linalg.norm(np.cross(b - a, c - a)))


def wedge_geometry(j1: dict, j2: dict, slope: dict, upper: dict, H: float) -> dict | None:
    """사면체 쐐기 — 꼭짓점 O(두 절리 · 사면 = 사면 끝) · A(1 · 사면 · 윗면) · B(2 · 사면 · 윗면) ·
    C(1 · 2 · 윗면). 사면 높이 H = 사면 위 꼭짓점 A · B 중 높은 쪽의 O 위 높이(SWedge 와 같음 —
    Priest 예제 부피 45.201 m³ 재현). 교선이 사면 · 윗면으로 안 나오면 None(쐐기 안 생김)."""
    n1 = plane_normal(j1["dip_direction"], j1["dip"])
    n2 = plane_normal(j2["dip_direction"], j2["dip"])
    nf = plane_normal(slope["dip_direction"], slope["dip"])
    nu = plane_normal(upper["dip_direction"], upper["dip"])
    # 윗면 = 사면 위로 H 오른 점을 지남(사면 경사방향 반대로 H cotψf 들어간 점)
    back = -line_vector(slope["dip_direction"], 0.0) * (H / math.tan(math.radians(slope["dip"])))
    crest = back + np.array([0.0, 0.0, H])
    du = float(nu @ crest)
    O = np.zeros(3)  # noqa: E741 — 꼭짓점 이름(SWedge 그림과 같음)
    A = _meet([(n1, 0.0), (nf, 0.0), (nu, du)])
    B = _meet([(n2, 0.0), (nf, 0.0), (nu, du)])
    C = _meet([(n1, 0.0), (n2, 0.0), (nu, du)])
    if A is None or B is None or C is None:
        return None
    # C 는 사면 뒤(산 쪽) · 위 — 교선이 사면으로 나와 윗면에 닿음
    if C[2] <= 1e-9 or float(nf @ C) >= -1e-9 or A[2] <= 1e-9 or B[2] <= 1e-9:
        return None
    scale = H / max(A[2], B[2])  # 면이 모두 O 를 지나 닮은꼴 — 높이만 맞춤
    A, B, C = A * scale, B * scale, C * scale
    volume = abs(float(np.linalg.det(np.array([A - O, B - O, C - O])))) / 6.0
    if volume <= 1e-9:
        return None
    centroid = (O + A + B + C) / 4.0
    inward = [n if float(n @ centroid) > 0 else -n for n in (n1, n2)]
    return {
        "vertices": {"O": O, "A": A, "B": B, "C": C},
        "volume": volume,
        "areas": [_tri_area(O, A, C), _tri_area(O, B, C)],
        "face_area": _tri_area(O, A, B),
        "upper_area": _tri_area(A, B, C),
        "normals": inward,
        "Hw": float(C[2]),
        "line_length": float(np.linalg.norm(C)),
    }


def _unit(v: np.ndarray) -> np.ndarray:
    n = float(np.linalg.norm(v))
    return v / n if n > 1e-12 else v * 0.0


def _sliding(A: np.ndarray, n1: np.ndarray, n2: np.ndarray, W: np.ndarray):
    """SWedge 미끄럼 방향 — (갈래, 방향). 갈래 = lift · 1 · 2 · 12 · stable."""
    s1 = _unit(np.cross(np.cross(n1, A), n1))
    s2 = _unit(np.cross(np.cross(n2, A), n2))
    c12 = np.cross(n1, n2)
    s12 = _unit(c12) * np.sign(float(c12 @ A))
    if A @ n1 > 0 and A @ n2 > 0:
        return "lift", _unit(A)
    if A @ n1 <= 0 and s1 @ n2 > 0:
        return "1", s1
    if A @ n2 <= 0 and s2 @ n1 > 0:
        return "2", s2
    if s1 @ n2 <= 0 and s2 @ n1 <= 0:
        return "12", s12
    return "stable", None


def _normals(mode: str, F: np.ndarray, n1: np.ndarray, n2: np.ndarray) -> tuple[float, float]:
    if mode == "1":
        return -float(F @ n1), 0.0
    if mode == "2":
        return 0.0, -float(F @ n2)
    if mode == "12":
        c12, c21 = np.cross(n1, n2), np.cross(n2, n1)
        return (
            -float(np.cross(F, n2) @ c12) / float(c12 @ c12),
            -float(np.cross(F, n1) @ c21) / float(c21 @ c21),
        )
    return 0.0, 0.0


def wedge_fs(
    *,
    j1: dict,
    j2: dict,
    slope: dict,
    upper: dict,
    H: float,
    gamma: float,
    water: dict | None = None,
    kh: float = 0.0,
    seismic_trend: float | None = None,
    external: list[dict] | None = None,
    bolt: dict | None = None,
    gamma_w: float = GAMMA_W,
) -> dict[str, Any] | None:
    """쐐기 안전율(SWedge 이론서 꼴) — 쐐기가 안 생기면 None.

    j = {dip_direction, dip, c, phi} · water = {"fill": p}(중력식 U = p³γw Hw a/6) |
    {"pressure": [u1, u2]}(평균 수압 kPa) · kh = 수평(기본 교선 추향) · external = [{force, trend,
    plunge}] · bolt = {force, trend, plunge}(수동력). FS = max(지지 전 · 지지 뒤 · 들림)."""
    geo = wedge_geometry(j1, j2, slope, upper, H)
    if geo is None:
        return None
    n1, n2 = geo["normals"]
    a1, a2 = geo["areas"]
    weight = gamma * geo["volume"]
    Wv = np.array([0.0, 0.0, -weight])
    active = Wv.copy()
    water = water or {}
    if "pressure" in water:
        u1, u2 = (float(u) for u in water["pressure"])
    else:
        p = float(water.get("fill", 0.0))
        u1 = u2 = p**3 * gamma_w * geo["Hw"] / 6.0
    Uf = [u1 * a1 * n1, u2 * a2 * n2]
    active += Uf[0] + Uf[1]
    line = vector_line(np.cross(n1, n2))
    if kh:
        trend = line[0] if seismic_trend is None else seismic_trend
        active += kh * weight * line_vector(trend, 0.0)
    for item in external or []:
        active += float(item["force"]) * line_vector(item["trend"], item["plunge"])
    passive = np.zeros(3)
    if bolt and bolt.get("force"):
        passive = float(bolt["force"]) * line_vector(bolt["trend"], bolt["plunge"])
    mode, s = _sliding(active, n1, n2, Wv)
    out = {
        "volume": geo["volume"],
        "W": weight,
        "areas": [a1, a2],
        "face_area": geo["face_area"],
        "upper_area": geo["upper_area"],
        "Hw": geo["Hw"],
        "line_length": geo["line_length"],
        "line": [round(line[0], 2), round(line[1], 2)],
        "water_force": [float(np.linalg.norm(Uf[0])), float(np.linalg.norm(Uf[1]))],
        "mode": mode,
        "vertices": {k: [round(float(x), 3) for x in v] for k, v in geo["vertices"].items()},
    }
    if mode == "stable":
        return {**out, "F": math.inf, "N": [0.0, 0.0], "drive": 0.0, "resist": 0.0}
    drive = float(active @ s)
    if mode == "lift":
        fs = -float(passive @ s) / drive if drive > 0 else math.inf
        return {**out, "F": fs, "N": [0.0, 0.0], "drive": drive, "resist": -float(passive @ s)}
    joints = (j1, j2)

    def resisting(F: np.ndarray) -> tuple[float, list[float]]:
        N = _normals(mode, F, n1, n2)
        total = 0.0
        for k, (n_k, a_k) in enumerate(zip(N, (a1, a2))):
            if str(k + 1) not in mode or n_k <= 0:
                continue  # 떨어진 면 = 저항 없음(SWedge 와 같음)
            tau = float(joints[k]["c"]) + n_k / a_k * math.tan(
                math.radians(float(joints[k]["phi"]))
            )
            total += tau * a_k
        return total, [max(0.0, n) for n in N]

    j_u, n_u = resisting(active)
    fs_u = j_u / drive if drive > 0 else math.inf
    fs, n_used, resist = fs_u, n_u, j_u
    if np.any(passive):
        j_s, n_s = resisting(active + passive)
        fs_s = (-float(passive @ s) + j_s) / drive if drive > 0 else math.inf
        fs_lift = -float(passive @ _unit(active)) / float(np.linalg.norm(active))
        if fs_s > fs:
            fs, n_used, resist = fs_s, n_s, -float(passive @ s) + j_s
        fs = max(fs, fs_lift)
    return {**out, "F": fs, "N": n_used, "drive": drive, "resist": resist}