"""E03 활동면 → 절편 · 방법 묶음 평가 — 원호(중심 · 반지름) · 꺾은선(꼭짓점) · 인장균열 · 모멘트 중심.

활동 방향 +x 로 맞춘 `Prepared` 위의 좌표. 탐색(`Engine_Search`) · 결과 꾸미기(`Engine_Describe`)가 같이 씀.
"""

import numpy as np

from E03_StabSection import E03_StabSection_Engine_GLE as gle
from E03_StabSection import E03_StabSection_Engine_LEM as lem

METHODS = ("ordinary", "bishop", "janbu", "spencer", "morgenstern_price")
_SURFACE_TOL = 1e-3  # 활동면이 지표면 위로 나온 높이 허용(m)


# ---- 면 → 절편 ---------------------------------------------------------------
def _crack_start(prep, xE, xX, base_fn):
    """인장균열 — 머리 쪽에서 활동면 깊이가 균열 깊이에 닿는 x (없으면 원래 진입점)."""
    if prep.crack_depth <= 0:
        return xE, np.zeros(len(xE))
    t = np.linspace(0.0, 1.0, 201)
    xs = xE[:, None] + (xX - xE)[:, None] * t
    depth = prep.surface(xs) - base_fn(xs)
    reach = depth >= prep.crack_depth
    k = np.where(reach.any(1), np.argmax(reach, 1), 0)
    k0 = np.maximum(k - 1, 0)
    rows = np.arange(len(xE))
    d0, d1 = depth[rows, k0], depth[rows, k]
    w = np.where(d1 > d0, (prep.crack_depth - d0) / np.maximum(d1 - d0, 1e-12), 0.0)
    x_new = np.where(k > 0, xs[rows, k0] + np.clip(w, 0, 1) * (xs[rows, k] - xs[rows, k0]), xE)
    return x_new, np.where(k > 0, prep.crack_depth, 0.0)


def _finish_slices(prep, xe, ze, crack, min_depth):
    sl = lem.slices(prep, xe, ze)
    if prep.crack_depth > 0:
        lem.add_crack_water(sl, crack, prep.crack_water)
    height = sl["zt"] - sl["zb"]
    sl["valid"] &= np.all(height > -_SURFACE_TOL, axis=1)
    sl["valid"] &= np.nanmax(height, axis=1) >= max(min_depth, 1e-3)
    return sl


def circle_slices(prep, xc, zc, R, xE, xX, n, min_depth=0.0):
    xc, zc, R = (np.asarray(v, float) for v in (xc, zc, R))

    def base(x):
        return zc[:, None] - np.sqrt(np.clip(R[:, None] ** 2 - (x - xc[:, None]) ** 2, 0, None))

    xE, crack = _crack_start(prep, np.asarray(xE, float), np.asarray(xX, float), base)
    t = np.linspace(0.0, 1.0, n + 1)
    xe = xE[:, None] + (xX - xE)[:, None] * t
    sl = _finish_slices(prep, xe, base(xe), crack, min_depth)
    sl["valid"] &= (xE >= xc - R) & (xX <= xc + R)
    return sl


def poly_base(vx, vz, x):
    """꼭짓점(S, V) 꺾은선을 x(S, m) 에서 — 행마다 선형 보간(양 끝 밖은 끝 마디 연장)."""
    z = None
    for k in range(vx.shape[1] - 1):
        x0, x1 = vx[:, k : k + 1], vx[:, k + 1 : k + 2]
        t = (x - x0) / np.where(x1 - x0 > 1e-12, x1 - x0, 1e-12)
        seg = vz[:, k : k + 1] + t * (vz[:, k + 1 : k + 2] - vz[:, k : k + 1])
        z = seg if z is None else np.where(x >= x0, seg, z)
    return z


def poly_slices(prep, vx, vz, n, min_depth=0.0):
    vx, vz = np.asarray(vx, float), np.asarray(vz, float)

    def base(x):
        return poly_base(vx, vz, x)

    xE, crack = _crack_start(prep, vx[:, 0], vx[:, -1], base)
    # 절편 경계 = 고른 간격 + 안쪽 꼭짓점(꺾임이 절편 안에 묻히지 않게 — SLOPE/W 와 같은 생각)
    t = np.linspace(0.0, 1.0, n + 1 - (vx.shape[1] - 2))
    xe = xE[:, None] + (vx[:, -1] - xE)[:, None] * t
    inner = np.clip(vx[:, 1:-1], xE[:, None], vx[:, -1:])
    xe = np.sort(np.concatenate([xe, inner], axis=1), axis=1)
    sl = _finish_slices(prep, xe, base(xe), crack, min_depth)
    sl["valid"] &= np.all(np.diff(vx, axis=1) > 1e-6, axis=1)
    return sl


def three_point_center(vx, vz):
    """비원호 모멘트 중심 = 진입 · 가장 깊은 꼭짓점 · 진출을 지나는 원의 중심(곧은 면은 위로 멀리)."""
    p1 = np.column_stack([vx[:, 0], vz[:, 0]])
    p3 = np.column_stack([vx[:, -1], vz[:, -1]])
    chord_z = vz[:, :1] + (vx - vx[:, :1]) * (
        (vz[:, -1:] - vz[:, :1]) / np.maximum(vx[:, -1:] - vx[:, :1], 1e-9)
    )
    k = np.argmax(chord_z - vz, axis=1)
    p2 = np.column_stack([vx[np.arange(len(vx)), k], vz[np.arange(len(vx)), k]])
    ax, ay = p1.T
    bx, by = p2.T
    cx, cy = p3.T
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    ux = (
        (ax**2 + ay**2) * (by - cy) + (bx**2 + by**2) * (cy - ay) + (cx**2 + cy**2) * (ay - by)
    ) / np.where(np.abs(d) < 1e-9, 1, d)
    uy = (
        (ax**2 + ay**2) * (cx - bx) + (bx**2 + by**2) * (ax - cx) + (cx**2 + cy**2) * (bx - ax)
    ) / np.where(np.abs(d) < 1e-9, 1, d)
    far = np.abs(d) < 1e-9
    span = np.hypot(cx - ax, cy - ay)
    ux = np.where(far, 0.5 * (ax + cx), ux)
    uy = np.where(far, np.maximum(ay, cy) + 10 * span, uy)
    return np.column_stack([ux, uy])


# ---- 방법 묶음 평가 ------------------------------------------------------------
def evaluate(sl, center, methods=METHODS, scan=True, mp_seed=None) -> dict:
    """방법별 F 묶음 — 출발값 차례 Fellenius → Bishop → (Janbu · Spencer · M-P).

    M-P = 임계면 몇 개에만(차례 계산이라 느림 · 탐색 고리에서 쓰지 않음) · 출발 참고값 = mp_seed 또는 Bishop.
    """
    out = {}
    base = lem.ordinary(sl, center)
    seed = np.where(np.isfinite(base["F"]), base["F"], 1.5)
    if "ordinary" in methods:
        out["ordinary"] = base
    if {"bishop", "spencer", "morgenstern_price"} & set(methods):
        out["bishop"] = lem.bishop(sl, center, seed)
        seed = np.where(np.isfinite(out["bishop"]["F"]), out["bishop"]["F"], seed)
    if "janbu" in methods:
        out["janbu"] = lem.janbu(sl, seed)
    if "spencer" in methods:
        out["spencer"] = lem.spencer(sl, center, seed, scan)
    if "morgenstern_price" in methods:
        ref = seed if mp_seed is None else np.where(np.isfinite(mp_seed), mp_seed, seed)
        out["morgenstern_price"] = gle.morgenstern_price(sl, center, F0=ref, curve=False)
    return {m: out[m] for m in methods}


def objective(method):
    def fun(sl, center):
        return evaluate(sl, center, (method,), scan=False)[method]["F"]

    return fun


# ---- 원호 기하 ---------------------------------------------------------------
def circle_from_chord(prep, xE, xX, theta):
    """진입 · 진출 x 와 중심각 θ → 중심 · 반지름(중심은 현 위쪽 · 두 끝이 중심보다 아래여야)."""
    zE, zX = prep.surface(xE), prep.surface(xX)
    dx, dz = xX - xE, zX - zE
    c = np.hypot(dx, dz)
    R = c / (2 * np.sin(theta / 2))
    d = R * np.cos(theta / 2)
    nx, nz = -dz / np.maximum(c, 1e-12), dx / np.maximum(c, 1e-12)
    flip = nz < 0
    nx, nz = np.where(flip, -nx, nx), np.where(flip, -nz, nz)
    xc = 0.5 * (xE + xX) + d * nx
    zc = 0.5 * (zE + zX) + d * nz
    ok = (zE <= zc) & (zX <= zc) & (dx > 0) & np.isfinite(zE) & np.isfinite(zX)
    return xc, zc, R, ok


def circle_cross_surface(prep, xc, zc, R):
    """원 아래 반쪽과 지표면의 교점 중 지표면 아래로 들어간 가장 넓은 구간 → (진입, 진출)."""
    sx, sz = prep.surface_x, prep.surface_z
    ax, az, bx, bz = sx[:-1], sz[:-1], sx[1:], sz[1:]
    dx, dz = bx - ax, bz - az
    fx, fz = ax[None, :] - xc[:, None], az[None, :] - zc[:, None]
    a = dx**2 + dz**2
    b = 2 * (fx * dx + fz * dz)
    c = fx**2 + fz**2 - R[:, None] ** 2
    disc = b**2 - 4 * a * c
    root = np.sqrt(np.clip(disc, 0, None))
    roots = []
    for sign in (-1, 1):
        t = (-b + sign * root) / (2 * np.maximum(a, 1e-12))
        x = ax + t * dx
        z = az + t * dz
        keep = (disc >= 0) & (t >= 0) & (t <= 1) & (z <= zc[:, None] + 1e-9)
        roots.append(np.where(keep, x, np.nan))
    xs = np.sort(np.concatenate(roots, axis=1), axis=1)
    mids = 0.5 * (xs[:, :-1] + xs[:, 1:])
    arc = zc[:, None] - np.sqrt(np.clip(R[:, None] ** 2 - (mids - xc[:, None]) ** 2, 0, None))
    below = np.isfinite(mids) & (prep.surface(np.nan_to_num(mids)) > arc)
    width = np.where(below, xs[:, 1:] - xs[:, :-1], -1.0)
    k = np.argmax(width, axis=1)
    rows = np.arange(len(xc))
    ok = width[rows, k] > 0
    return xs[rows, k], xs[rows, k + 1], ok