Merge remote-tracking branch 'origin/sub_laptop_4' into sub_laptop_5

E03_StabSection/E03_StabSection_Engine_Describe.py

@@ -0,0 +1,114 @@
+"""E03 결과 꾸미기 — 임계면 글상자 · 절편 자유물체도 · F–λ 곡선 · Safety Map 상위 면.
+
+좌표는 원래 방향(x = B06 offset · 좌 +)으로 되돌려 냄. 출력 칸은 덧붙이기만(E03 화면 · E05 보고서가 읽음).
+"""
+
+import numpy as np
+
+from E03_StabSection import E03_StabSection_Engine_GLE as gle
+from E03_StabSection import E03_StabSection_Engine_LEM as lem
+from E03_StabSection.E03_StabSection_Engine_Surface import (
+    circle_slices,
+    evaluate,
+    poly_slices,
+    three_point_center,
+)
+
+#: Safety Map 면 수 · 면마다 꼭짓점 수(원호)
+SAFETY_MAP_TOP = 50
+_MAP_POINTS = 25
+
+
+def _curve(sl, center, method, f=None) -> dict:
+    """F_m · F_f 대 λ (화면 그래프) — Spencer = 상수 f · M-P = 반정현."""
+    if method == "spencer":
+        grid, fm, ff = lem.lambda_curve(
+            sl, center, np.linspace(*lem.LAMBDA_RANGE, gle.CURVE_POINTS)
+        )
+        fm, ff = fm[:, 0], ff[:, 0]
+    else:
+        grid, fm, ff = gle.lambda_curves(sl, center, f)
+        fm, ff = fm[0], ff[0]
+    num = lambda v: [float(x) if np.isfinite(x) else None for x in v]  # noqa: E731
+    return {"lambda": num(grid), "Fm": num(fm), "Ff": num(ff)}
+
+
+def describe(prep, kind, geom, method, n, min_depth, seed_F=None) -> dict:
+    """임계면 하나 → 안전율 · 좌표 · 절편표(자유물체도 힘) · 수렴 · (Spencer · M-P) λ · F–λ 곡선."""
+    s = prep.direction
+    if kind == "circle":
+        xc, zc, R, xE, xX = geom
+        sl = circle_slices(prep, [xc], [zc], [R], [xE], [xX], n, min_depth)
+        center = np.array([[xc, zc]])
+        pts = np.column_stack([sl["xe"][0], sl["ze"][0]])
+    else:
+        vx, vz = geom
+        sl = poly_slices(prep, vx[None, :], vz[None, :], n, min_depth)
+        center = three_point_center(vx[None, :], vz[None, :])
+        pts = np.column_stack([vx, vz])
+    seed = None if seed_F is None else np.array([seed_F], float)
+    res = evaluate(sl, center, (method,), mp_seed=seed)[method]
+    F = float(res["F"][0])
+    lam = float(res["lambda"][0]) if method in ("spencer", "morgenstern_price") else 0.0
+    f = gle.interslice_function(sl["xe"]) if method == "morgenstern_price" else None
+    rows = lem.slice_table(sl, 0, F, lam, method, f) if np.isfinite(F) else []
+    for r in rows:
+        r["material"] = (
+            prep.materials[r["material"]]["id"] if r["material"] < len(prep.materials) else None
+        )
+        # 행 순서 = 머리 → 끝 · E · X 의 left = 머리 쪽 경계 · right = 끝 쪽 경계(해석 틀 그대로)
+        r["x_left"], r["x_right"] = sorted((s * r["x_left"], s * r["x_right"]))
+        r["x_mid"] = s * r["x_mid"]
+    out = {
+        "F": F if np.isfinite(F) else None,
+        "kind": kind,
+        "points": [[float(s * x), float(z)] for x, z in (pts if s > 0 else pts[::-1])],
+        "entry": [float(s * sl["xe"][0, 0]), float(sl["ze"][0, 0])],
+        "exit": [float(s * sl["xe"][0, -1]), float(sl["ze"][0, -1])],
+        "iterations": int(res["iterations"][0]),
+        "converged": bool(res["converged"][0]),
+        "slices": rows,
+    }
+    if kind == "circle":
+        out["center"] = [float(s * center[0, 0]), float(center[0, 1])]
+        out["radius"] = float(geom[2])
+    else:
+        out["moment_center"] = [float(s * center[0, 0]), float(center[0, 1])]
+    if method in ("spencer", "morgenstern_price"):
+        out["lambda"] = lam
+        out["theta_deg"] = float(np.degrees(np.arctan(lam)))
+        if np.isfinite(F):
+            out["lambda_curve"] = _curve(sl, center, method, f)
+    if method == "morgenstern_price":
+        out["interslice_function"] = "half_sine"
+    if method == "janbu":
+        out["F0"] = float(res["F0"][0]) if np.isfinite(res["F0"][0]) else None
+        out["f0"] = float(res["f0"][0])
+    if prep.crack_depth > 0:
+        out["crack_water_force"] = float(sl["Pw"][0])
+    return out
+
+
+def safety_map(prep, circles, F, method, top=SAFETY_MAP_TOP) -> dict:
+    """Safety Map — 거친 탐색 원호 중 안전율 낮은 top 개(F · 중심 · 반지름 · 진입 · 진출 · 꼭짓점)."""
+    s = prep.direction
+    xc, zc, R, xE, xX = circles
+    order = np.argsort(np.where(np.isfinite(F), F, np.inf))[:top]
+    order = order[np.isfinite(F[order])]
+    t = np.linspace(0.0, 1.0, _MAP_POINTS)
+    out = []
+    for k in order:
+        x = xE[k] + (xX[k] - xE[k]) * t
+        z = zc[k] - np.sqrt(np.clip(R[k] ** 2 - (x - xc[k]) ** 2, 0, None))
+        pts = np.column_stack([s * x, z])
+        out.append(
+            {
+                "F": float(F[k]),
+                "center": [float(s * xc[k]), float(zc[k])],
+                "radius": float(R[k]),
+                "entry": [float(s * xE[k]), float(z[0])],
+                "exit": [float(s * xX[k]), float(z[-1])],
+                "points": [[float(a), float(b)] for a, b in (pts if s > 0 else pts[::-1])],
+            }
+        )
+    return {"method": method, "surfaces": out}

E03_StabSection/E03_StabSection_Engine_GLE.py

@@ -0,0 +1,198 @@
+"""E03 엄밀해 — Morgenstern-Price(GLE) · 절편간 힘 X = λ·f(x)·E · f = 반정현(기본) 또는 상수.
+
+식 = LEM 과 같은 절편 평형(`E03_StabSection_Engine_LEM` 머리글) — 절편간 힘 함수만 경계마다 다름:
+  절편 i 연직 평형 → N·(m_α − λ·f_{i+1}·A) = W − c'·sinα/F − X_i + λ·f_{i+1}·(E_i + B0)
+  수평 평형 → E_{i+1} = E_i + B0 + N·A · X_{i+1} = λ·f_{i+1}·E_{i+1} · 머리 E_0 = 균열 수압 · X_0 = 0
+  f ≡ 1 이면 Spencer 와 같은 식(시험으로 확인). F_m · F_f 는 LEM 과 같은 전체 평형식.
+풀이 = λ 훑기(F_m(λ) · F_f(λ) 곡선 — 화면 그래프 자료) → 교점 → (F, λ) 2 차 뉴턴(수치 야코비안).
+계산은 임계면 몇 개만(탐색 고리 밖) — 면을 λ 갯수만큼 복제해 한 번에 셈. 외부 코드 없음(NumPy 만).
+"""
+
+import numpy as np
+
+from E03_StabSection import E03_StabSection_Engine_LEM as lem
+
+FUNCTIONS = ("half_sine", "constant")
+_TINY = 1e-9
+CURVE_POINTS = 21
+#: M-P λ 범위 — 반정현 평균 2/π 라 평균 X/E 가 Spencer 범위(LEM.LAMBDA_RANGE)와 같게 π/2 배
+LAMBDA_RANGE = tuple(v * np.pi / 2 for v in lem.LAMBDA_RANGE)
+
+
+def interslice_function(xe: np.ndarray, kind: str = "half_sine") -> np.ndarray:
+    """절편 경계마다 f(x) — 반정현 sin(π·(x − x_진입)/(x_진출 − x_진입)) · 상수 1."""
+    if kind == "constant":
+        return np.ones_like(xe)
+    if kind != "half_sine":
+        raise ValueError(f"절편간 힘 함수 모름: {kind}")
+    x0, x1 = xe[:, :1], xe[:, -1:]
+    return np.sin(np.pi * (xe - x0) / np.maximum(x1 - x0, _TINY))
+
+
+def forces(sl: dict, F, lam, f):
+    """절편마다 N · 분모 · 경계 E · X (꼴 (S, n) · (S, n) · (S, n+1) · (S, n+1)) — 머리부터 차례로."""
+    p = lem._pre(sl)
+    F = np.asarray(F, float)
+    lam = np.asarray(lam, float)
+    g = 1.0 / F[:, None]
+    m = sl["ca"] + g * p["sat"]
+    A = sl["sa"] - g * p["tca"]
+    B0 = p["khws"] - g * p["cuca"]
+    rhs = sl["W"] - g * p["cusa"]
+    S, n = rhs.shape
+    N = np.empty((S, n))
+    den = np.empty((S, n))
+    E = np.empty((S, n + 1))
+    X = np.zeros((S, n + 1))
+    E[:, 0] = sl["Pw"]
+    for i in range(n):
+        lf = lam * f[:, i + 1]
+        d = m[:, i] - lf * A[:, i]
+        den[:, i] = d
+        N[:, i] = (rhs[:, i] - X[:, i] + lf * (E[:, i] + B0[:, i])) / np.where(
+            np.abs(d) < _TINY, _TINY, d
+        )
+        E[:, i + 1] = E[:, i] + B0[:, i] + N[:, i] * A[:, i]
+        X[:, i + 1] = lf * E[:, i + 1]
+    return N, den, E, X
+
+
+def _fm_ff(sl, center, F, lam, f):
+    N, den, *_ = forces(sl, F, lam, f)
+    at, gn, d0, ma = lem._mom(sl, center)
+    p = lem._pre(sl)
+    fm = lem._ratio(ma + (at * N).sum(1), d0 + (gn * N).sum(1))
+    ff = lem._ratio(p["rc0"] + (p["tca"] * N).sum(1), p["fd0"] + (sl["sa"] * N).sum(1), neg=False)
+    return fm, ff, den
+
+
+def _repeat(sl: dict, k: int) -> dict:
+    """면마다 k 번 복제(λ 훑기 · 야코비안을 한 번에) — 미리 묶은 칸은 버림."""
+    S = len(sl["W"])
+    return {
+        key: (np.repeat(v, k, axis=0) if isinstance(v, np.ndarray) and v.shape[:1] == (S,) else v)
+        for key, v in sl.items()
+        if key not in ("_pre", "_mom")
+    }
+
+
+def _solve_f(fun, F0, iters=60):
+    """F = fun(F) — 수치 미분 뉴턴 · 안 되면 반쯤 섞은 고정점(LEM `_iterate` 와 같은 안전장치)."""
+    F = np.asarray(F0, float).copy()
+    done = np.zeros(len(F), bool)
+    for _ in range(iters):
+        f1 = fun(F)
+        h1 = f1 - F
+        dlt = 1e-6 * np.maximum(np.abs(F), 0.1)
+        dh = ((fun(F + dlt) - (F + dlt)) - h1) / dlt
+        newton = F - h1 / np.where(np.abs(dh) < 1e-12, -1.0, dh)
+        use = np.isfinite(newton) & (dh < -1e-9) & (newton > 0.5 * F) & (newton < 2.0 * F)
+        Fn = np.where(use, newton, 0.5 * (F + f1))
+        Fn = np.where(np.isfinite(Fn) & (Fn > 0), Fn, np.nan)
+        conv = np.abs(h1) < lem.TOL * np.maximum(1.0, np.abs(F))
+        F = np.where(done | conv, F, Fn)
+        done |= conv | ~np.isfinite(F)
+        if done.all():
+            break
+    return np.where(done, F, np.nan)
+
+
+def lambda_curves(sl, center, f, grid=None, F0=None):
+    """λ 마다 F_m(λ) · F_f(λ) — 꼴 (G,) · (S, G) · (S, G)."""
+    grid = np.linspace(*LAMBDA_RANGE, CURVE_POINTS) if grid is None else np.asarray(grid)
+    S, G = len(sl["W"]), len(grid)
+    rep = _repeat(sl, G)
+    rc = np.repeat(center, G, axis=0)
+    rf = np.repeat(f, G, axis=0)
+    lam = np.tile(grid, S)
+    start = np.full(S * G, 1.5) if F0 is None else np.repeat(F0, G)
+    fm = _solve_f(lambda F: _fm_ff(rep, rc, F, lam, rf)[0], start)
+    ff = _solve_f(lambda F: _fm_ff(rep, rc, F, lam, rf)[1], start)
+    return grid, fm.reshape(S, G), ff.reshape(S, G)
+
+
+def _crossing(grid, fm, ff):
+    """F_m − F_f 부호가 바뀌는 곳 중 λ = 0 가까운 쪽 → (F, λ) · 없으면 nan."""
+    S = len(fm)
+    F, L = np.full(S, np.nan), np.full(S, np.nan)
+    best = np.full(S, np.inf)
+    d = fm - ff
+    for k in range(len(grid) - 1):
+        d0, d1 = d[:, k], d[:, k + 1]
+        cross = np.isfinite(d0) & np.isfinite(d1) & (np.sign(d0) != np.sign(d1))
+        t = np.where(cross, d0 / np.where(d0 == d1, 1, d0 - d1), 0)
+        lc = grid[k] + t * (grid[k + 1] - grid[k])
+        fc = fm[:, k] + t * (fm[:, k + 1] - fm[:, k])
+        better = cross & (np.abs(lc) < best)
+        best = np.where(better, np.abs(lc), best)
+        F, L = np.where(better, fc, F), np.where(better, lc, L)
+    return F, L
+
+
+def _newton2(sl, center, f, F, lam, iters=40):
+    """(F_m − F, F_f − F) = 0 — 세 벌 복제(F · F+h · λ+h)로 수치 야코비안 · 걸음 제한 LEM Spencer 와 같음."""
+    S = len(F)
+    rep, rc, rf = _repeat(sl, 3), np.repeat(center, 3, axis=0), np.repeat(f, 3, axis=0)
+    F, lam = F.copy(), lam.copy()
+    done = ~np.isfinite(F)
+    it = np.zeros(S, int)
+    for _ in range(iters):
+        hF = 1e-6 * np.maximum(np.abs(F), 0.1)
+        hl = 1e-6
+        Fs = np.column_stack([F, F + hF, F]).ravel()
+        ls = np.column_stack([lam, lam, lam + hl]).ravel()
+        fm, ff, _ = _fm_ff(rep, rc, np.where(np.isfinite(Fs), Fs, 1.0), ls, rf)
+        rm = (fm - Fs).reshape(S, 3)
+        rf_ = (ff - Fs).reshape(S, 3)
+        j11, j12 = (rm[:, 1] - rm[:, 0]) / hF, (rm[:, 2] - rm[:, 0]) / hl
+        j21, j22 = (rf_[:, 1] - rf_[:, 0]) / hF, (rf_[:, 2] - rf_[:, 0]) / hl
+        det = j11 * j22 - j12 * j21
+        det = np.where(np.abs(det) < 1e-14, np.nan, det)
+        dF = -(j22 * rm[:, 0] - j12 * rf_[:, 0]) / det
+        dl = -(-j21 * rm[:, 0] + j11 * rf_[:, 0]) / det
+        scale = np.minimum(
+            1.0,
+            np.minimum(
+                0.5 * np.abs(F) / np.maximum(np.abs(dF), _TINY),
+                0.3 / np.maximum(np.abs(dl), _TINY),
+            ),
+        )
+        F = np.where(done, F, F + scale * dF)
+        lam = np.where(done, lam, lam + scale * dl)
+        it += ~done
+        small = (np.abs(dF) < lem.TOL * np.maximum(1, np.abs(F))) & (np.abs(dl) < lem.TOL)
+        done |= small | ~np.isfinite(F) | (F < lem.F_MIN)
+        if done.all():
+            break
+    return F, lam, it
+
+
+def morgenstern_price(sl: dict, center: np.ndarray, kind="half_sine", F0=None, curve=True):
+    """M-P (F, λ) · 수렴 · (curve) F_m · F_f–λ 곡선. F0 = 출발 참고값(Spencer · Bishop) — 엉뚱한 근 거름."""
+    f = interslice_function(sl["xe"], kind)
+    grid, fm, ff = lambda_curves(sl, center, f, F0=F0)
+    F, lam = _crossing(grid, fm, ff)
+    F, lam, it = _newton2(sl, center, f, F, lam)
+    ok = np.isfinite(F)
+    safe = np.where(ok & (F > 0), F, 1.0)
+    rm, rff, den = _fm_ff(sl, center, safe, np.where(ok, lam, 0.0), f)
+    p = lem._pre(sl)
+    N = forces(sl, safe, np.where(ok, lam, 0.0), f)[0]
+    conv = ok & (np.abs(rm - F) < 1e-4) & (np.abs(rff - F) < 1e-4)
+    conv &= (lam >= LAMBDA_RANGE[0]) & (lam <= LAMBDA_RANGE[1])
+    if F0 is not None:
+        conv &= (F > 0.5 * F0) & (F < 2.0 * F0)
+    malpha_ok = np.all(den >= lem.MALPHA_MIN, axis=1)
+    R = p["cu"] + N * sl["tanphi"]
+    good = sl["valid"] & conv & (F > lem.F_MIN) & malpha_ok & (R.sum(1) > 0)
+    out = {
+        "F": np.where(good, F, np.nan),
+        "lambda": np.where(good, lam, np.nan),
+        "iterations": it,
+        "converged": conv,
+        "malpha_ok": malpha_ok,
+        "function": kind,
+    }
+    if curve:
+        out["curve"] = (grid, fm, ff)
+    return out

E03_StabSection/E03_StabSection_Engine_LEM.py

@@ -25,6 +25,8 @@ LAMBDA_RANGE = (-1.25, 0.25)
 #: 이보다 작은 F = 1/F 항이 터진 엉뚱한 근(실제 비탈 F 는 이 위) — 버림.
 F_MIN = 0.05
 TOL = 1e-6
+#: Janbu 보정 f0 식의 d/L 상한(지식DB 02 §2.3)
+JANBU_DL_MAX = 0.35
 MAX_ITER = 100
 
 
@@ -116,6 +118,12 @@ class Prepared:
         # 빈 칸 = 아주 낮은 두께 0 구간(nan 처리 없이 바로 셈)
         for name in ("bot_a", "bot_b", "top_a", "top_b"):
             setattr(self, name, np.nan_to_num(getattr(self, name), nan=_PAD_Z))
+        # 절편 계산이 한 번에 모아 읽는 표 (띠, 6, 구간) — 바닥 · 바닥 기울기 · 윗면 · 윗면 기울기 · γt · γsat − γt
+        gt, gs = self.gamma_t[self.mat], self.gamma_sat[self.mat]
+        self.table = np.stack(
+            [self.bot_a, self.bot_b - self.bot_a, self.top_a, self.top_b - self.top_a, gt, gs - gt],
+            axis=1,
+        )
 
     # ---- 조회 -------------------------------------------------------------
     def surface(self, x: np.ndarray) -> np.ndarray:
@@ -147,25 +155,36 @@ def slices(prep: Prepared, xe: np.ndarray, ze: np.ndarray) -> dict:
     xm = 0.5 * (xe[:, :-1] + xe[:, 1:])
     zb = 0.5 * (ze[:, :-1] + ze[:, 1:])
     alpha = np.arctan2(ze[:, :-1] - ze[:, 1:], np.maximum(b, _TINY))
-    bot, top, mat = prep.columns(xm)
+    j = np.clip(np.searchsorted(prep.X, xm, side="right") - 1, 0, len(prep.xa) - 1)
+    t = ((xm - prep.xa[j]) / (prep.xb[j] - prep.xa[j]))[..., None]
+    G = prep.table[j]
+    mat = prep.mat[j]
+    bot = G[..., 0, :] + t * G[..., 1, :]
+    top = G[..., 2, :] + t * G[..., 3, :]
+    gt, dgs = G[..., 4, :], G[..., 5, :]
     lo = np.maximum(bot, zb[..., None])
     thick = np.clip(top - lo, 0.0, None)
     zt = prep.surface(xm)
     if prep.gw_mode == "line":
         zw = np.minimum(prep.water(xm), zt)
         wet = np.clip(np.minimum(top, zw[..., None]) - lo, 0.0, None)
-        dry = thick - wet
         u = GAMMA_W * np.clip(zw - zb, 0.0, None)
-        wet_mid = lo + 0.5 * wet
-        dry_mid = lo + wet + 0.5 * dry
-        g_dry, g_wet = prep.gamma_t[mat] * dry, prep.gamma_sat[mat] * wet
+        weight_per_b = np.sum(gt * thick + dgs * wet, axis=-1)
     else:
-        g_dry, g_wet = prep.gamma_t[mat] * thick, 0.0 * thick
-        dry_mid = wet_mid = lo + 0.5 * thick
+        wet = None
         u = np.zeros_like(zb)
-    weight_per_b = np.sum(g_dry + g_wet, axis=-1)
-    moment_z = np.sum(g_dry * dry_mid + g_wet * wet_mid, axis=-1)
-    zg = np.where(weight_per_b > 0, moment_z / np.maximum(weight_per_b, _TINY), zb)
+        weight_per_b = np.sum(gt * thick, axis=-1)
+    if prep.kh > 0:  # 무게중심 높이는 지진 관성력 팔에만 쓰임
+        if wet is None:
+            moment_z = np.sum(gt * thick * (lo + 0.5 * thick), axis=-1)
+        else:
+            dry = thick - wet
+            moment_z = np.sum(
+                gt * dry * (lo + wet + 0.5 * dry) + (gt + dgs) * wet * (lo + 0.5 * wet), axis=-1
+            )
+        zg = np.where(weight_per_b > 0, moment_z / np.maximum(weight_per_b, _TINY), zb)
+    else:
+        zg = zb
     Ws = weight_per_b * b
     if prep.gw_mode == "ru":
         u = prep.ru * weight_per_b
@@ -225,77 +244,130 @@ def add_crack_water(sl: dict, depth: np.ndarray, fill: float) -> None:
 
 
 # ---- 방법 -------------------------------------------------------------------
-def _parts(sl: dict, F, lam):
-    """주어진 F · λ 의 절편 수직력 N · 저항 R · 분모(m_α − λA).
-
-    유효 수직력 음(인장) 절편도 그대로 둠(Slide2 와 같음 — 0 으로 자르면 ACADS · Arai & Tagyo 대조가 3~6 % 벗어남).
-    대신 저항 합이 0 이하인 면은 `_finish` 가 버림.
-    """
-    F = np.asarray(F, float)[:, None]
-    lam = np.asarray(lam, float)[:, None]
-    W = sl["W"].copy()
-    W[:, 0] += lam[:, 0] * sl["Pw"]  # 머리 경계 수압은 전단 없음 → X_0 = 0
-    t = sl["tanphi"]
-    N, den, A, B0 = _normal(sl, F, lam, W, t)
-    R = sl["c"] * sl["l"] + (N - sl["u"] * sl["l"]) * t
-    return N, R, den, A, B0
+# 한 벌 식을 F 와 무관한 칸(_pre) · 모멘트 중심 칸(_mom) 으로 미리 묶어 두고, F · λ 마다는 N 만 새로 셈.
+# g = 1/F 로 놓으면 N = (P0 − g·P1)/(Q0 + g·Q1) — dN/dg · dN/dλ 가 닫힌 꼴이라 뉴턴 야코비안을 해석적으로 씀.
+def _pre(sl: dict) -> dict:
+    p = sl.get("_pre")
+    if p is None:
+        t, sa, ca = sl["tanphi"], sl["sa"], sl["ca"]
+        cu = sl["c"] * sl["l"] - sl["u"] * sl["l"] * t
+        khws = sl["kh"] * sl["Ws"]
+        lw = khws.copy()
+        lw[:, 0] += sl["Pw"]  # 머리 경계 수압은 전단 없음 → X_0 = 0
+        p = {
+            "cu": cu,
+            "cusa": cu * sa,
+            "cuca": cu * ca,
+            "sat": sa * t,
+            "tca": t * ca,
+            "khws": khws,
+            "lw": lw,
+            "fd0": khws.sum(1) + sl["Pw"],
+            "rc0": (cu * ca).sum(1),
+        }
+        sl["_pre"] = p
+    return p
 
 
-def _normal(sl, F, lam, W, t):
-    sa, ca = sl["sa"], sl["ca"]
-    cu = sl["c"] * sl["l"] - sl["u"] * sl["l"] * t
-    m = ca + sa * t / F
-    A = sa - t * ca / F
-    B0 = sl["kh"] * sl["Ws"] - cu * ca / F
-    den = m - lam * A
-    N = (W + lam * B0 - cu * sa / F) / np.where(np.abs(den) < _TINY, _TINY, den)
-    return N, den, A, B0
-
-
-def _moment_terms(sl: dict, center: np.ndarray, N: np.ndarray):
+def _mom(sl: dict, center: np.ndarray):
+    """모멘트 팔(a) · 분모 N 계수(gN) · 상수 칸 — 같은 중심 배열이면 다시 쓰고 아니면 새로 셈."""
+    cached = sl.get("_mom")
+    if cached is not None and cached[0] is center:
+        return cached[1]
+    p = _pre(sl)
     xo, zo = center[:, 0:1], center[:, 1:2]
-    a = (sl["xm"] - xo) * sl["sa"] + (sl["zb"] - zo) * sl["ca"]
-    D = (
-        (sl["xw"] - xo) * sl["W"]
-        + (sl["zg"] - zo) * sl["kh"] * sl["Ws"]
-        - (sl["xm"] - xo) * N * sl["ca"]
-        + (sl["zb"] - zo) * N * sl["sa"]
-    ).sum(1) + (sl["zp"] - center[:, 1]) * sl["Pw"]
-    return a, D
+    dx, dz = sl["xm"] - xo, sl["zb"] - zo
+    a = dx * sl["sa"] + dz * sl["ca"]
+    gn = dz * sl["sa"] - dx * sl["ca"]
+    d0 = ((sl["xw"] - xo) * sl["W"] + (sl["zg"] - zo) * p["khws"]).sum(1)
+    d0 = d0 + (sl["zp"] - center[:, 1]) * sl["Pw"]
+    out = (a * sl["tanphi"], gn, d0, (a * p["cu"]).sum(1))
+    sl["_mom"] = (center, out)
+    return out
 
 
-def _force_den(sl: dict, N: np.ndarray):
-    return (N * sl["sa"] + sl["kh"] * sl["Ws"]).sum(1) + sl["Pw"]
+def _normal(sl, F, lam, grad=0):
+    """절편 수직력 N · 분모(m_α − λA) · grad ≥ 1 dN/dg · grad 2 dN/dλ — F · λ 꼴 (S,)."""
+    p = _pre(sl)
+    g = 1.0 / np.asarray(F, float)[:, None]
+    L = np.asarray(lam, float)[:, None]
+    P1 = p["cusa"] + L * p["cuca"]
+    Q1 = p["sat"] + L * p["tca"]
+    den = sl["ca"] - L * sl["sa"] + g * Q1
+    safe = np.where(np.abs(den) < _TINY, _TINY, den)
+    N = (sl["W"] + L * p["lw"] - g * P1) / safe
+    if not grad:
+        return N, den
+    dg = -(P1 + N * Q1) / safe
+    if grad == 1:
+        return N, den, dg, None
+    dl = (p["lw"] - g * p["cuca"] + N * (sl["sa"] - g * p["tca"])) / safe
+    return N, den, dg, dl
 
 
-def _fm(sl, center, F, lam):
-    N, R, den, *_ = _parts(sl, F, lam)
-    a, D = _moment_terms(sl, center, N)
-    return (a * R).sum(1) / np.where(np.abs(D) < _TINY, -_TINY, D), den
+def _parts(sl: dict, F, lam):
+    """주어진 F · λ 의 N · 저항 R · 분모 · A · B0(절편표용).
+
+    유효 수직력 음(인장) 절편도 그대로 둠(Slide2 와 같음 — 0 으로 자르면 ACADS · Arai & Tagyo
+    대조가 3~6 % 벗어남). 대신 저항 합이 0 이하인 면은 `_finish` 가 버림.
+    """
+    p = _pre(sl)
+    F = np.asarray(F, float)
+    N, den = _normal(sl, F, lam)
+    g = 1.0 / F[:, None]
+    A = sl["sa"] - p["tca"] * g
+    B0 = p["khws"] - p["cuca"] * g
+    return N, p["cu"] + N * sl["tanphi"], den, A, B0
 
 
-def _ff(sl, F, lam):
-    N, R, den, *_ = _parts(sl, F, lam)
-    return (R * sl["ca"]).sum(1) / np.where(
-        np.abs(_force_den(sl, N)) < _TINY, _TINY, _force_den(sl, N)
-    ), den
+def _ratio(num, den, dnum=None, dden=None, neg=True):
+    """f = num/den (작은 분모는 ∓TINY) · 도함수 있으면 (f, df)."""
+    d = np.where(np.abs(den) < _TINY, -_TINY if neg else _TINY, den)
+    f = num / d
+    return f if dnum is None else (f, (dnum - f * dden) / d)
+
+
+def _balance(cf, cd, c0, d0, nd, F, neg):
+    """f = (c0 + Σcf·N)/(d0 + Σcd·N) 와 df/dF · df/dλ (nd = _normal 결과 · 없는 도함수는 None)."""
+    N, _den, dg, dl = nd
+    num, den = c0 + (cf * N).sum(1), d0 + (cd * N).sum(1)
+    if dg is None:
+        return _ratio(num, den, neg=neg), None, None
+    f, fg = _ratio(num, den, (cf * dg).sum(1), (cd * dg).sum(1), neg)
+    fl = None if dl is None else _ratio(num, den, (cf * dl).sum(1), (cd * dl).sum(1), neg)[1]
+    return f, -fg / np.asarray(F, float) ** 2, fl
+
+
+def _fm(sl, center, F, lam, grad=0, nd=None):
+    """모멘트 평형 F_m = Σ a·R / D → (F_m, dF_m/dF, dF_m/dλ, 분모)."""
+    at, gn, d0, ma = _mom(sl, center)
+    nd = nd or _normal(sl, F, lam, grad)
+    if len(nd) == 2:
+        nd = (*nd, None, None)
+    return (*_balance(at, gn, ma, d0, nd, F, True), nd[1])
+
+
+def _ff(sl, F, lam, grad=0, nd=None):
+    """힘 평형 F_f = Σ R·cosα / (Σ(N·sinα + kh·Ws) + Pw) → (F_f, dF_f/dF, dF_f/dλ, 분모)."""
+    p = _pre(sl)
+    nd = nd or _normal(sl, F, lam, grad)
+    if len(nd) == 2:
+        nd = (*nd, None, None)
+    return (*_balance(p["tca"], sl["sa"], p["rc0"], p["fd0"], nd, F, False), nd[1])
 
 
 def _iterate(fun, F0):
-    """F = fun(F) 풀이(묶음 전체) — 수치 미분 뉴턴(h = fun(F) − F) · 안 되면 반쯤 섞은 고정점 걸음.
+    """F = fun(F) 풀이(묶음 전체) — 해석 도함수 뉴턴 · 안 되면 반쯤 섞은 고정점 걸음.
 
-    단순 고정점은 급한 머리 절편에서 진동해 느림(도함수 ≈ −1) — 뉴턴은 몇 번에 끝남.
+    fun(F) → (f, df/dF, …). 단순 고정점은 급한 머리 절편에서 진동해 느림(도함수 ≈ −1) — 뉴턴은 몇 번에 끝남.
     """
     F = np.asarray(F0, float).copy()
     done = np.zeros(len(F), bool)
     iters = np.zeros(len(F), int)
-    den = None
     for _ in range(MAX_ITER):
-        f1, den = fun(F)
+        f1, df = fun(F)[:2]
         h1 = f1 - F
-        delta = 1e-6 * np.maximum(np.abs(F), 0.1)
-        f2, _ = fun(F + delta)
-        dh = ((f2 - (F + delta)) - h1) / delta
+        dh = df - 1.0
         newton = F - h1 / np.where(np.abs(dh) < 1e-12, -1.0, dh)
         relaxed = 0.5 * (F + f1)
         use = np.isfinite(newton) & (dh < -1e-9) & (newton > 0.5 * F) & (newton < 2.0 * F)
@@ -307,17 +379,20 @@ def _iterate(fun, F0):
         done |= conv | ~np.isfinite(F)
         if done.all():
             break
-    f1, den = fun(np.where(np.isfinite(F), F, 1.0))
-    return F, iters, done & np.isfinite(F), den
+    return F, iters, done & np.isfinite(F)
 
 
 def ordinary(sl: dict, center: np.ndarray) -> dict:
     """Fellenius — 절편간 힘 무시 · 바닥 수직 분해 N = W·cosα − kh·Ws·sinα · 모멘트 평형."""
-    N = sl["W"] * sl["ca"] - sl["kh"] * sl["Ws"] * sl["sa"]
-    R = sl["c"] * sl["l"] + (N - sl["u"] * sl["l"]) * sl["tanphi"]
-    a, D = _moment_terms(sl, center, N)
-    F = (a * R).sum(1) / np.where(np.abs(D) < _TINY, -_TINY, D)
-    ok = sl["valid"] & np.isfinite(F) & (F > F_MIN) & (R.sum(1) > 0)
+    N = sl["W"] * sl["ca"] - _pre(sl)["khws"] * sl["sa"]
+    at, gn, d0, ma = _mom(sl, center)
+    F = _ratio(ma + (at * N).sum(1), d0 + (gn * N).sum(1))
+    ok = (
+        sl["valid"]
+        & np.isfinite(F)
+        & (F > F_MIN)
+        & ((_pre(sl)["cu"] + N * sl["tanphi"]).sum(1) > 0)
+    )
     # 같은 꼴 거르기 — 진출부 급경사(m_α < 0.2) 면은 Bishop 처럼 버림
     m = sl["ca"] + sl["sa"] * sl["tanphi"] / np.where(ok, F, 1.0)[:, None]
     ok &= np.all(m >= MALPHA_MIN, axis=1)
@@ -327,7 +402,7 @@ def ordinary(sl: dict, center: np.ndarray) -> dict:
 def bishop(sl: dict, center: np.ndarray, F0=None) -> dict:
     F0 = np.full(len(sl["W"]), 1.5) if F0 is None else F0
     lam = np.zeros(len(F0))
-    F, it, conv, _ = _iterate(lambda f: _fm(sl, center, f, lam), F0)
+    F, it, conv = _iterate(lambda f: _fm(sl, center, f, lam, 1), F0)
     return _finish(sl, F, it, conv, lam)
 
 
@@ -335,7 +410,7 @@ def janbu(sl: dict, F0=None) -> dict:
     """Janbu 간편(F_0) + 보정 f_0 — 지식DB 02 §2.3 · b1 = 0.31(c=0) · 0.69(φ=0) · 0.50."""
     F0 = np.full(len(sl["W"]), 1.5) if F0 is None else F0
     lam = np.zeros(len(F0))
-    F, it, conv, _ = _iterate(lambda f: _ff(sl, f, lam), F0)
+    F, it, conv = _iterate(lambda f: _ff(sl, f, lam, 1), F0)
     out = _finish(sl, F, it, conv, lam)
     out["F0"] = out["F"]
     out["f0"] = janbu_f0(sl)
@@ -350,7 +425,8 @@ def janbu_f0(sl: dict) -> np.ndarray:
     dist = np.abs(dz * (xe - xe[:, :1]) - dx * (ze - ze[:, :1]))[:, :] / np.maximum(
         L[:, None], _TINY
     )
-    ratio = dist.max(1) / np.maximum(L, _TINY)
+    # d/L ≤ 0.35 에서만 식이 맞음(지식DB 02 §2.3) — 넘으면 0.35 값(f0 곡선 꼭대기 근처)
+    ratio = np.minimum(dist.max(1) / np.maximum(L, _TINY), JANBU_DL_MAX)
     c_all0 = np.all(sl["c"] <= 0, axis=1)
     phi_all0 = np.all(sl["tanphi"] <= 0, axis=1)
     b1 = np.where(c_all0, 0.31, np.where(phi_all0, 0.69, 0.50))
@@ -358,7 +434,7 @@ def janbu_f0(sl: dict) -> np.ndarray:
 
 
 def spencer(sl: dict, center: np.ndarray, F0=None, scan: bool = True) -> dict:
-    """Spencer — 모멘트 · 힘 두 평형을 함께 맞추는 (F, λ) · 뉴턴(수치 야코비안) · 안 맞으면 λ 훑기.
+    """Spencer — 모멘트 · 힘 두 평형을 함께 맞추는 (F, λ) · 뉴턴(해석 야코비안) · 안 맞으면 λ 훑기.
 
     scan=False — 탐색 고리 안(시험 면 대부분이 버려짐)에서는 훑기를 건너뜀(시간 대부분이 여기).
     """
@@ -381,27 +457,27 @@ def spencer(sl: dict, center: np.ndarray, F0=None, scan: bool = True) -> dict:
     return out
 
 
-def _residual(sl, center, F, lam):
-    N, R, den, *_ = _parts(sl, F, lam)
-    a, D = _moment_terms(sl, center, N)
-    fd = _force_den(sl, N)
-    rm = ((a * R).sum(1) - F * D) / np.maximum(np.abs(D), _TINY)
-    rf = ((R * sl["ca"]).sum(1) - F * fd) / np.maximum(np.abs(fd), _TINY)
-    return rm, rf
-
-
 def _spencer_newton(sl, center, F, lam, iters=40):
+    """(F_m − F, F_f − F) = 0 뉴턴 — 걸음 크기 제한(F 절반 · λ 0.3) · 수렴 = 두 잔차 < 1e-4 · λ 범위 안.
+
+    끝난 면은 빼고(남은 면이 절반 아래면 묶음을 줄임) · F < F_MIN 으로 빠진 면은 버림(어차피 `_finish` 가 버림).
+    """
     F, lam = F.copy(), lam.copy()
     it = np.zeros(len(F), int)
     done = np.zeros(len(F), bool)
+    idx, w, wc = np.arange(len(F)), sl, center
     for _ in range(iters):
-        rm, rf = _residual(sl, center, F, lam)
-        hF = 1e-6 * np.maximum(F, 0.1)
-        hl = 1e-6
-        rmF, rfF = _residual(sl, center, F + hF, lam)
-        rml, rfl = _residual(sl, center, F, lam + hl)
-        j11, j12 = (rmF - rm) / hF, (rml - rm) / hl
-        j21, j22 = (rfF - rf) / hF, (rfl - rf) / hl
+        act = ~done[idx]
+        if not act.any():
+            break
+        if act.sum() < 0.5 * len(idx):
+            idx, w, wc = idx[act], _subset(w, act), wc[act]
+        f, lm, dn = F[idx], lam[idx], done[idx]
+        nd = _normal(w, f, lm, 2)
+        fm, fmF, fml, _ = _fm(w, wc, f, lm, nd=nd)
+        ff, ffF, ffl, _ = _ff(w, f, lm, nd=nd)
+        rm, rf = fm - f, ff - f
+        j11, j12, j21, j22 = fmF - 1.0, fml, ffF - 1.0, ffl
         det = j11 * j22 - j12 * j21
         det = np.where(np.abs(det) < 1e-14, np.nan, det)
         dF = -(j22 * rm - j12 * rf) / det
@@ -409,34 +485,41 @@ def _spencer_newton(sl, center, F, lam, iters=40):
         scale = np.minimum(
             1.0,
             np.minimum(
-                0.5 * F / np.maximum(np.abs(dF), _TINY), 0.3 / np.maximum(np.abs(dl), _TINY)
+                0.5 * f / np.maximum(np.abs(dF), _TINY), 0.3 / np.maximum(np.abs(dl), _TINY)
             ),
         )
-        F = np.where(done, F, F + scale * dF)
-        lam = np.where(done, lam, lam + scale * dl)
-        it += ~done
-        small = (np.abs(dF) < TOL * np.maximum(1, np.abs(F))) & (np.abs(dl) < TOL)
-        done |= small | ~np.isfinite(F)
-        if done.all():
-            break
-    rm, rf = _residual(sl, center, F, lam)
+        f = np.where(dn, f, f + scale * dF)
+        F[idx], lam[idx] = f, np.where(dn, lm, lm + scale * dl)
+        it[idx] += ~dn
+        small = (np.abs(dF) < TOL * np.maximum(1, np.abs(f))) & (np.abs(dl) < TOL)
+        done[idx] = dn | small | ~np.isfinite(f) | (f < F_MIN)
+    nd = _normal(sl, F, lam)
+    rm = _fm(sl, center, F, lam, nd=nd)[0] - F
+    rf = _ff(sl, F, lam, nd=nd)[0] - F
     conv = done & np.isfinite(F) & (F > 0) & (np.abs(rm) < 1e-4) & (np.abs(rf) < 1e-4)
     conv &= (lam >= LAMBDA_RANGE[0]) & (lam <= LAMBDA_RANGE[1])
     return F, lam, it, conv
 
 
-def _spencer_scan(sl, center):
-    """λ 를 훑어 F_m(λ) − F_f(λ) 부호가 바뀌는 곳(λ=0 가까운 쪽)을 출발점으로."""
-    grid = np.linspace(LAMBDA_RANGE[0], LAMBDA_RANGE[1], 16)
+def lambda_curve(sl, center, grid=None):
+    """λ 마다 모멘트 · 힘 평형 F_m(λ) · F_f(λ) — 꼴 (len(grid), S) 둘(화면 F–λ 그래프 · 훑기)."""
+    grid = np.linspace(LAMBDA_RANGE[0], LAMBDA_RANGE[1], 16) if grid is None else grid
     S = len(sl["W"])
-    diff = np.full((len(grid), S), np.nan)
     fm_all = np.full((len(grid), S), np.nan)
+    ff_all = np.full((len(grid), S), np.nan)
     F0 = np.full(S, 1.5)
     for k, g in enumerate(grid):
         lam = np.full(S, g)
-        fm, *_ = _iterate(lambda f: _fm(sl, center, f, lam), F0)
-        ff, *_ = _iterate(lambda f: _ff(sl, f, lam), F0)
-        diff[k], fm_all[k] = fm - ff, fm
+        fm_all[k] = _iterate(lambda f: _fm(sl, center, f, lam, 1), F0)[0]
+        ff_all[k] = _iterate(lambda f: _ff(sl, f, lam, 1), F0)[0]
+    return grid, fm_all, ff_all
+
+
+def _spencer_scan(sl, center):
+    """λ 를 훑어 F_m(λ) − F_f(λ) 부호가 바뀌는 곳(λ=0 가까운 쪽)을 출발점으로."""
+    grid, fm_all, ff_all = lambda_curve(sl, center)
+    S = len(sl["W"])
+    diff = fm_all - ff_all
     best_F, best_l = np.full(S, 1.5), np.zeros(S)
     best_gap = np.full(S, np.inf)
     for k in range(len(grid) - 1):
@@ -453,15 +536,18 @@ def _spencer_scan(sl, center):
 
 
 def _subset(sl: dict, mask: np.ndarray) -> dict:
-    out = {
+    """면 묶음 일부 — 미리 묶은 칸(_pre · _mom)은 버림(부분 꼴이 달라 다시 셈)."""
+    return {
         k: (v[mask] if isinstance(v, np.ndarray) and v.shape[:1] == mask.shape else v)
         for k, v in sl.items()
+        if k not in ("_pre", "_mom")
     }
-    return out
 
 
 def _finish(sl, F, it, conv, lam):
-    _, R, den, *_ = _parts(sl, np.where(np.isfinite(F) & (F > 0), F, 1.0), lam)
+    Fs = np.where(np.isfinite(F) & (F > 0), F, 1.0)
+    N, den = _normal(sl, Fs, lam)
+    R = _pre(sl)["cu"] + N * sl["tanphi"]
     malpha_ok = np.all(den >= MALPHA_MIN, axis=1)
     ok = sl["valid"] & conv & np.isfinite(F) & (F > F_MIN) & malpha_ok & (R.sum(1) > 0)
     return {
@@ -472,19 +558,33 @@ def _finish(sl, F, it, conv, lam):
     }
 
 
-def slice_table(sl: dict, k: int, F: float, lam: float = 0.0, method: str = "bishop") -> list[dict]:
-    """한 면(k 번째)의 절편표 — 폭 · 높이 · 경사 · 무게 · 수압 · 강도 · N' · 동원 전단 · (Spencer) E · X."""
+def slice_table(
+    sl: dict, k: int, F: float, lam: float = 0.0, method: str = "bishop", f=None
+) -> list[dict]:
+    """한 면(k 번째)의 절편표 + 자유물체도 힘 — 행 = 머리 → 끝 · left = 머리 쪽 경계.
+
+    폭 · 높이 · 경사 · 무게 · 수압 · 강도 · N(전) · N' · 동원 전단 · 저항 · 지진력 · 경계 E · X.
+    M-P 는 f(경계마다 절편간 힘 함수)를 받아 차례 계산(`Engine_GLE.forces`).
+    """
     one = _subset(sl, np.arange(len(sl["W"])) == k)
     if method == "ordinary":
         N = one["W"] * one["ca"] - one["kh"] * one["Ws"] * one["sa"]
         E = np.zeros((1, N.shape[1] + 1))
+        X = E
+    elif method == "morgenstern_price":
+        from E03_StabSection.E03_StabSection_Engine_GLE import forces
+
+        N, _den, E, X = forces(one, np.array([F]), np.array([lam]), f[k : k + 1])
     else:
         N, R, _den, A, B0 = _parts(one, np.array([F]), np.array([lam]))
         dE = B0 + N * A
         E = np.concatenate([one["Pw"][:, None], one["Pw"][:, None] + np.cumsum(dE, axis=1)], axis=1)
+        X = lam * E
+        X[:, 0] = 0.0
     R = one["c"] * one["l"] + (N - one["u"] * one["l"]) * one["tanphi"]
     rows = []
     for i in range(N.shape[1]):
+        ul = float(one["u"][0, i] * one["l"][0, i])
         rows.append(
             {
                 "x_left": float(one["xe"][0, i]),
@@ -499,10 +599,20 @@ def slice_table(sl: dict, k: int, F: float, lam: float = 0.0, method: str = "bis
                 "material": int(one["mat"][0, i]),
                 "c": float(one["c"][0, i]),
                 "phi_deg": float(np.degrees(np.arctan(one["tanphi"][0, i]))),
-                "normal_eff": float(N[0, i] - one["u"][0, i] * one["l"][0, i]),
+                "normal_eff": float(N[0, i]) - ul,
                 "shear_mobilized": float(R[0, i] / F),
                 "E_left": float(E[0, i]),
-                "X_left": float(lam * E[0, i]),
+                "X_left": float(X[0, i]),
+                # 덧(2차) — 자유물체도
+                "x_mid": float(one["xm"][0, i]),
+                "z_base": float(one["zb"][0, i]),
+                "z_top": float(one["zt"][0, i]),
+                "normal_total": float(N[0, i]),
+                "pore_force": ul,
+                "shear_strength": float(R[0, i]),
+                "seismic_force": float(one["kh"] * one["Ws"][0, i]),
+                "E_right": float(E[0, i + 1]),
+                "X_right": float(X[0, i + 1]),
             }
         )
     return rows

E03_StabSection/E03_StabSection_Engine_Search.py

@@ -2,8 +2,8 @@
 
 활동 방향 +x 로 맞춘 `Prepared` 위에서 찾고 좌표를 원래 방향으로 되돌려 냄.
 다듬기 · 비원호 = 차분진화(best/1/bin · 무리 전체를 한 번에 평가 · seed 고정 → 같은 입력 같은 답)
-+ 비원호 마지막 국부 다듬기 = SciPy Nelder-Mead. SciPy differential_evolution 은 개체마다 파이썬 돌연변이 고리라
-한 비탈 시간의 30 % 를 먹어 직접 묶음 꼴로 씀.
++ 비원호 마지막 국부 다듬기 = 묶음 Nelder-Mead(후보 넷을 한 번에). SciPy 최적화는 점마다 파이썬 고리라 느려
+직접 묶음 꼴로 씀. 방법마다 찾은 면을 모든 방법으로 다시 셈(교차) → 한 방법 탐색이 놓친 골짜기를 메움.
 상용 SW 대응: Grid & Radius(SLOPE/W · Slide) · Entry & Exit(SLOPE/W) · Auto Refine(Slide) ·
 Block Search + Optimize(Slide · SLOPE/W 꼭짓점 최적화).
 """
@@ -11,181 +11,38 @@ Block Search + Optimize(Slide · SLOPE/W 꼭짓점 최적화).
 import time
 
 import numpy as np
-from scipy.optimize import minimize
 
 from E03_StabSection import E03_StabSection_Engine_LEM as lem
+from E03_StabSection.E03_StabSection_Engine_Describe import describe, safety_map
+from E03_StabSection.E03_StabSection_Engine_Surface import (  # noqa: F401 — 밖에서 Search 로 부르던 이름
+    METHODS,
+    circle_cross_surface,
+    circle_from_chord,
+    circle_slices,
+    evaluate,
+    objective,
+    poly_base,
+    poly_slices,
+    three_point_center,
+)
 
-METHODS = ("ordinary", "bishop", "janbu", "spencer")
 _BIG = 1e3
-_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) -> dict:
-    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" in methods or "spencer" in 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)
-    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
+# 탐색 예산(실험 손잡이)
+LS_ITER = 30
+LS_ITER_SPENCER = 60  # Spencer 비원호는 풀리는 면이 조각나 더 걸음
+POLISH = "nm"
+NM_ITER = 120
+PATTERN_ITER = 30
+DE_MAXITER = 60
+DE_PATIENCE = 12
+SPENCER_SEED_POOL = 300
+#: 교차 재평가에 넣는 방법별 마지막 무리 상위 원 수
+CROSS_TOP = 10
+#: Spencer 원호 씨앗 — Bishop 임계 둘레 구름 크기
+CLOUD_SIZE = 200
 
 
 # ---- 원호 ---------------------------------------------------------------------
-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
-
-
 def entry_exit_circles(prep, limits, n_points=25, thetas_deg=(10, 20, 35, 50, 65, 80, 100, 120)):
     xl, xr = limits
     pts = np.linspace(xl, xr, n_points)
@@ -237,9 +94,9 @@ def _circle_eval(prep, xc, zc, R, xE, xX, n, min_depth, methods, chunk=4000):
     return res
 
 
-def refine_circle(prep, limits, method, seeds, n, min_depth, seed=0, maxiter=60):
+def refine_circle(prep, limits, method, seeds, n, min_depth, seed=0, maxiter=None):
     """원호 다듬기 — (진입 x, 진출 x, 중심각) 차분진화 · 첫 무리 = 거친 탐색의 좋은 면."""
-    fun = _objective(method)
+    fun = objective(method)
     lo = np.array([limits[0], limits[0], np.radians(3)])
     hi = np.array([limits[1], limits[1], np.radians(170)])
 
@@ -252,9 +109,35 @@ def refine_circle(prep, limits, method, seeds, n, min_depth, seed=0, maxiter=60)
             F[ok] = np.where(np.isfinite(f), f, _BIG)
         return F
 
-    x, fmin, _pop, _f = _de(obj, lo, hi, _population(seeds, lo, hi, 30, seed), maxiter, seed)
+    init = _population(seeds, lo, hi, 30, seed)
+    x, fmin, pop, f = _de(obj, lo, hi, init, maxiter or DE_MAXITER, seed, patience=DE_PATIENCE)
+    top = pop[np.argsort(f)[:CROSS_TOP]]
+    top = top[obj(top) < _BIG]
+    xc, zc, R, _ = circle_from_chord(prep, top[:, 0], top[:, 1], top[:, 2])
+    near = [
+        (float(a), float(b), float(c), float(d), float(e))
+        for a, b, c, d, e in zip(xc, zc, R, top[:, 0], top[:, 1])
+    ]
     xc, zc, R, _ = circle_from_chord(prep, x[0:1], x[1:2], x[2:3])
-    return float(fmin), (float(xc[0]), float(zc[0]), float(R[0]), float(x[0]), float(x[1]))
+    best = (float(xc[0]), float(zc[0]), float(R[0]), float(x[0]), float(x[1]))
+    return float(fmin), best, near, top
+
+
+def _cloud(prep, limits, method, p0, n, min_depth, seed, k=CLOUD_SIZE, keep=10):
+    """원 (진입 x, 진출 x, 중심각) 둘레 구름 k 개를 method 로 한 번에 셈 → 풀리는 좋은 것 keep 개."""
+    rng = np.random.default_rng(seed)
+    span = limits[1] - limits[0]
+    lo = np.array([limits[0], limits[0], np.radians(3)])
+    hi = np.array([limits[1], limits[1], np.radians(170)])
+    p = np.clip(p0 + rng.normal(size=(k, 3)) * [0.03 * span, 0.03 * span, np.radians(4)], lo, hi)
+    xc, zc, R, ok = circle_from_chord(prep, p[:, 0], p[:, 1], p[:, 2])
+    F = np.full(k, np.inf)
+    if ok.any():
+        sl = circle_slices(prep, xc[ok], zc[ok], R[ok], p[ok, 0], p[ok, 1], n, min_depth)
+        f = evaluate(sl, np.column_stack([xc[ok], zc[ok]]), (method,), scan=False)[method]["F"]
+        F[ok] = np.where(np.isfinite(f), f, np.inf)
+    order = np.argsort(F)[:keep]
+    return p[order[np.isfinite(F[order])]].reshape(-1, 3)
 
 
 def _de(obj, lo, hi, init, maxiter, seed, cr=0.7, tol=1e-6, patience=25):
@@ -310,6 +193,66 @@ def _local_search(obj, starts, lo, hi, step, maxiter, seed, lam=24, starts_keep=
     return x[i].copy(), float(f[i]), x, f
 
 
+def _nelder_mead(obj, x0, lo, hi, maxiter, xatol=1e-4, fatol=1e-6):
+    """묶음 Nelder-Mead — 반사 · 확장 · 바깥 · 안쪽 수축 넷을 한 번에 셈(판단은 표준 NM 그대로 · 줄어듦도 한 번에).
+
+    SciPy NM 은 점마다 따로 불러 한 면 계산의 파이썬 부담이 큼 — 부르는 수를 줄임. 첫 꼭짓점 = SciPy 와 같은 5 % 흔들기.
+    """
+    D = len(x0)
+    sim = np.tile(np.asarray(x0, float), (D + 1, 1))
+    for k in range(D):
+        sim[k + 1, k] = sim[k + 1, k] * 1.05 if sim[k + 1, k] != 0 else 0.00025
+    sim = np.clip(sim, lo, hi)
+    fs = obj(sim)
+    for _ in range(maxiter):
+        order = np.argsort(fs, kind="stable")
+        sim, fs = sim[order], fs[order]
+        if np.max(np.abs(sim[1:] - sim[0])) <= xatol and np.max(np.abs(fs[1:] - fs[0])) <= fatol:
+            break
+        xbar = sim[:-1].mean(0)
+        xw = sim[-1]
+        cand = np.clip(
+            np.vstack(
+                [2 * xbar - xw, 3 * xbar - 2 * xw, 1.5 * xbar - 0.5 * xw, 0.5 * xbar + 0.5 * xw]
+            ),
+            lo,
+            hi,
+        )
+        fr, fe, fc, fcc = obj(cand)
+        if fr < fs[0]:
+            sim[-1], fs[-1] = (cand[1], fe) if fe < fr else (cand[0], fr)
+        elif fr < fs[-2]:
+            sim[-1], fs[-1] = cand[0], fr
+        elif fr < fs[-1] and fc <= fr:
+            sim[-1], fs[-1] = cand[2], fc
+        elif fr >= fs[-1] and fcc < fs[-1]:
+            sim[-1], fs[-1] = cand[3], fcc
+        else:
+            sim[1:] = np.clip(sim[0] + 0.5 * (sim[1:] - sim[0]), lo, hi)
+            fs[1:] = obj(sim[1:])
+    k = int(np.argmin(fs))
+    return sim[k].copy(), float(fs[k])
+
+
+def _pattern(obj, x, fx, lo, hi, step, iters):
+    """묶음 좌표 패턴 탐색 — 축마다 ±걸음 2D 점을 한 번에 셈 · 나아지면 옮기고 그 쪽으로 한 번 더 · 아니면 걸음 ½."""
+    h = step.copy()
+    for _ in range(iters):
+        cand = np.clip(np.vstack([x + np.diag(h), x - np.diag(h)]), lo, hi)
+        f = obj(cand)
+        k = int(np.argmin(f))
+        if f[k] < fx - 1e-9:
+            move = cand[k] - x
+            x2 = np.clip(cand[k] + move, lo, hi)
+            f2 = float(obj(x2[None, :])[0])
+            x, fx = (x2, f2) if f2 < f[k] else (cand[k], float(f[k]))
+        else:
+            h *= 0.5
+            if np.all(h < step * 1e-3):
+                break
+    return x, fx
+
+
 def _population(seeds, lo, hi, size, seed, spread=0.0):
     """첫 무리 = 씨앗 + (spread > 0 이면 씨앗 둘레 흔든 것 80 % · 고른 무작위 20 %)."""
     rng = np.random.default_rng(seed)
@@ -352,10 +295,10 @@ def layer_bottom_depth(prep, x, material):
 
 
 def optimize_polyline(
-    prep, limits, method, seeds, n, min_depth, m=4, seed=0, maxiter=60, polish=True
+    prep, limits, method, seeds, n, min_depth, m=4, seed=0, maxiter=None, polish=True
 ):
     """비원호 꼭짓점 최적화 — 차분진화 뒤 가장 좋은 면을 Nelder-Mead 로 국부 다듬기(상용 SAO · Optimize 대응)."""
-    fun = _objective(method)
+    fun = objective(method)
     xl, xr = limits
     depth_max = max(prep.surface(np.linspace(xl, xr, 50)).max() - prep.bottom_z, 1.0)
     lo = np.array([xl, xl] + [0.02] * m + [0.0] * m)
@@ -373,16 +316,15 @@ def optimize_polyline(
         return F
 
     step = np.array([0.5, 0.5] + [0.03] * m + [0.3] * m)
-    x, fmin, pop, f = _local_search(obj, np.asarray(seeds, float), lo, hi, step, maxiter, seed)
-    if polish and fmin < _BIG:
-        r = minimize(
-            lambda q: float(obj(np.clip(q, lo, hi))[0]),
-            x,
-            method="Nelder-Mead",
-            options={"maxfev": 200, "xatol": 1e-4, "fatol": 1e-6},
-        )
-        if r.fun < fmin:
-            x, fmin = np.clip(r.x, lo, hi), float(r.fun)
+    x, fmin, pop, f = _local_search(
+        obj, np.asarray(seeds, float), lo, hi, step, maxiter or LS_ITER, seed
+    )
+    if polish and fmin < _BIG and POLISH == "pattern":
+        x, fmin = _pattern(obj, x, fmin, lo, hi, step * 0.1, PATTERN_ITER)
+    elif polish and fmin < _BIG and POLISH == "nm":
+        xn, fn = _nelder_mead(obj, x, lo, hi, NM_ITER)
+        if fn < fmin:
+            x, fmin = xn, fn
     order = np.argsort(f)
     vx, vz = _decode(prep, x[None, :], m)
     return fmin, x, (vx[0], vz[0]), [x, *pop[order]]
@@ -452,53 +394,6 @@ def block_seeds(prep, limits, m, n_points=12, angles=((55, 35), (65, 25))):
     return out
 
 
-# ---- 결과 꾸미기 ----------------------------------------------------------------
-def describe(prep, kind, geom, method, n, min_depth) -> dict:
-    s = prep.direction
-    if kind == "circle":
-        xc, zc, R, xE, xX = geom
-        sl = circle_slices(prep, [xc], [zc], [R], [xE], [xX], n, min_depth)
-        center = np.array([[xc, zc]])
-        xs = sl["xe"][0]
-        pts = np.column_stack([xs, sl["ze"][0]])
-    else:
-        vx, vz = geom
-        sl = poly_slices(prep, vx[None, :], vz[None, :], n, min_depth)
-        center = three_point_center(vx[None, :], vz[None, :])
-        pts = np.column_stack([vx, vz])
-    res = evaluate(sl, center, (method,))[method]
-    F = float(res["F"][0])
-    lam = float(res.get("lambda", [0.0])[0]) if method == "spencer" else 0.0
-    rows = lem.slice_table(sl, 0, F, lam, method) if np.isfinite(F) else []
-    for r in rows:
-        r["material"] = (
-            prep.materials[r["material"]]["id"] if r["material"] < len(prep.materials) else None
-        )
-        r["x_left"], r["x_right"] = sorted((s * r["x_left"], s * r["x_right"]))
-    out = {
-        "F": F if np.isfinite(F) else None,
-        "kind": kind,
-        "points": [[float(s * x), float(z)] for x, z in (pts if s > 0 else pts[::-1])],
-        "entry": [float(s * sl["xe"][0, 0]), float(sl["ze"][0, 0])],
-        "exit": [float(s * sl["xe"][0, -1]), float(sl["ze"][0, -1])],
-        "iterations": int(res["iterations"][0]),
-        "converged": bool(res["converged"][0]),
-        "slices": rows,
-    }
-    if kind == "circle":
-        out["center"] = [float(s * center[0, 0]), float(center[0, 1])]
-        out["radius"] = float(geom[2])
-    if method == "spencer":
-        out["lambda"] = lam
-        out["theta_deg"] = float(np.degrees(np.arctan(lam)))
-    if method == "janbu":
-        out["F0"] = float(res["F0"][0]) if np.isfinite(res["F0"][0]) else None
-        out["f0"] = float(res["f0"][0])
-    if prep.crack_depth > 0:
-        out["crack_water_force"] = float(sl["Pw"][0])
-    return out
-
-
 # ---- 한 비탈 해석 ----------------------------------------------------------------
 def analyze(
     model: dict,
@@ -526,14 +421,15 @@ def analyze(
         lim = tuple(sorted((s * float(limits[0]), s * float(limits[1]))))
     results: dict = {m: {} for m in methods}
     counts = {}
-    contour = None
+    contour = smap = None
     best_circles: dict = {}
     if circular:
         a = entry_exit_circles(prep, lim, entry_exit_points)
         gx, gz = grid_centers(prep, lim, grid[0], grid[1])
         b = grid_radius_circles(prep, lim, gx, gz, grid[2])
         coarse = tuple(np.concatenate([a[k], b[k]]) for k in range(5))
-        cheap = tuple(m for m in methods if m != "spencer") or ("bishop",)
+        cheap = tuple(m for m in methods if m not in ("spencer", "morgenstern_price"))
+        cheap = cheap or ("bishop",)
         F = _circle_eval(prep, *coarse, coarse_slices, min_depth, cheap)
         counts["circles"] = int(len(coarse[0]))
         counts["circles_valid"] = int(np.sum(np.isfinite(F[cheap[0]])))
@@ -553,17 +449,44 @@ def analyze(
                 for row in cmin.reshape(len(gz), len(gx))
             ],
         }
-        for m in methods:
+        smap = safety_map(prep, coarse, F[key], key)
+        if "spencer" in methods:
+            # Spencer 씨앗 — 거친 탐색 상위 원에서 Spencer 를 따로 셈(풀리는 원이 조각나 있어 Bishop 순위만으론 놓침)
+            top = np.argsort(np.where(np.isfinite(F[key]), F[key], np.inf))[:SPENCER_SEED_POOL]
+            top = top[np.isfinite(F[key][top])]
+            F["spencer"] = np.full(len(coarse[0]), np.nan)
+            if len(top):
+                c5 = [v[top] for v in coarse[:5]]
+                sl = circle_slices(prep, *c5, coarse_slices, min_depth)
+                center = np.column_stack([c5[0], c5[1]])
+                F["spencer"][top] = evaluate(sl, center, ("spencer",), scan=False)["spencer"]["F"]
+        searched = [m for m in methods if m != "morgenstern_price"]
+        pool, tops = [], {}
+        for m in searched:
             score = F[m] if m in F else F[key]
             order = np.argsort(np.where(np.isfinite(score), score, np.inf))[:30]
             order = order[np.isfinite(score[order])]
             if not len(order):
                 continue
             seeds = np.column_stack([coarse[3][order], coarse[4][order], _theta(coarse, order)])
-            fmin, geom = refine_circle(prep, lim, m, seeds, n_slices, min_depth, seed)
+            if m == "spencer" and "bishop" in tops:
+                # Spencer 가 풀리는 원은 Bishop 임계 둘레에 좁게 · 조각나 있음 — 그 둘레 구름에서 풀리는 원을 씨앗에 보탬
+                cloud = _cloud(prep, lim, m, tops["bishop"][0], n_slices, min_depth, seed)
+                seeds = np.vstack([seeds[:20], cloud, tops["bishop"]])
+            fmin, geom, near, top = refine_circle(prep, lim, m, seeds, n_slices, min_depth, seed)
             if fmin < _BIG:
-                results[m]["circular"] = describe(prep, "circle", geom, m, n_slices, min_depth)
                 best_circles[m] = geom
+                pool += [geom, *near]
+                tops[m] = top
+        # 방법마다 찾은 원(+ 마지막 무리 상위)을 모든 방법으로 다시 셈 → 방법별 가장 낮은 면
+        # (한 방법 탐색이 놓친 골짜기 메움 · Spencer 는 풀리는 원이 조각나 있어 특히 필요)
+        for m in methods:
+            # M-P(차례 계산 · 느림)는 다른 방법이 고른 원에서만
+            geoms = list(best_circles.values()) if m == "morgenstern_price" else pool
+            g = _cross(prep, "circle", geoms, m, n_slices, min_depth)
+            if g is not None:
+                best_circles[m] = g
+                results[m]["circular"] = _describe(prep, "circle", g, m, n_slices, min_depth)
     if noncircular:
         mv = polyline_vertices
         seeds = []
@@ -576,35 +499,73 @@ def analyze(
             seeds.append(circle_to_params(prep, circle, mv))
         seeds += block_seeds(prep, lim, mv)
         seeds = _best_seeds(prep, seeds, mv, n_slices, min_depth, 40)
-        # 힘 평형 방법(Janbu) 먼저 · Spencer 는 그 마지막 무리에서 이어 찾음 · 모멘트 방법(Fellenius · Bishop)은
-        # 비원호 중심이 임의라 따로 최적화하지 않고 찾은 면에서 값만 셈(상용 SW 도 비원호 최적화 = Janbu · Spencer · GLE)
+        # 힘 평형 방법(Janbu) 먼저 · Spencer 는 그 마지막 무리에서 이어 찾음 · 모멘트 방법(Fellenius ·
+        # Bishop)은 비원호 중심이 임의라 따로 최적화하지 않고 찾은 면에서 값만 셈(참고 칸 · 상용 SW 도
+        # 비원호 최적화 = Janbu · Spencer · GLE) · M-P 는 찾은 면에서 값만
         found = []
         for m in [k for k in ("janbu", "spencer") if k in methods]:
+            iters = LS_ITER_SPENCER if m == "spencer" else LS_ITER
             fmin, _x, geom, top = optimize_polyline(
-                prep, lim, m, seeds, n_slices, min_depth, mv, seed
+                prep, lim, m, seeds, n_slices, min_depth, mv, seed, iters
             )
             seeds = [*top, *seeds]
             if fmin < _BIG:
-                results[m]["noncircular"] = describe(prep, "polyline", geom, m, n_slices, min_depth)
                 found.append(geom)
-        for m in [k for k in ("ordinary", "bishop") if k in methods]:
-            cands = [describe(prep, "polyline", g, m, n_slices, min_depth) for g in found]
-            cands = [c for c in cands if c["F"] is not None]
-            if cands:
-                results[m]["noncircular"] = min(cands, key=lambda c: c["F"])
+        for m in methods:
+            g = _cross(prep, "polyline", found, m, n_slices, min_depth)
+            if g is not None:
+                results[m]["noncircular"] = _describe(prep, "polyline", g, m, n_slices, min_depth)
     for m in methods:
-        cands = [v for v in results[m].values() if v.get("F") is not None]
+        # 비원호 면 위 Fellenius · Bishop = 모멘트 중심 임의 → 참고 칸(최소 F 후보에서 뺌 · Slide2 와 같음)
+        if m in ("ordinary", "bishop") and results[m].get("noncircular"):
+            results[m]["noncircular"]["reference_only"] = True
+        cands = [
+            v
+            for k, v in results[m].items()
+            if v.get("F") is not None and not v.get("reference_only")
+        ]
         results[m]["min"] = min(cands, key=lambda v: v["F"]) if cands else None
     return {
         "direction": s,
         "limits": sorted([s * lim[0], s * lim[1]]),
         "methods": results,
         "contour": contour,
+        "safety_map": smap,
         "counts": counts,
         "elapsed_s": round(time.perf_counter() - t0, 3),
     }
 
 
+def _stack(prep, kind, geoms, n, min_depth):
+    if kind == "circle":
+        xc, zc, R, xE, xX = (np.array([g[k] for g in geoms], float) for k in range(5))
+        return circle_slices(prep, xc, zc, R, xE, xX, n, min_depth), np.column_stack([xc, zc])
+    vx = np.array([g[0] for g in geoms], float)
+    vz = np.array([g[1] for g in geoms], float)
+    return poly_slices(prep, vx, vz, n, min_depth), three_point_center(vx, vz)
+
+
+def _cross(prep, kind, geoms, method, n, min_depth):
+    """찾은 면들 위 method 안전율을 한 번에 → 가장 낮은 면(없으면 None)."""
+    if not geoms:
+        return None
+    sl, center = _stack(prep, kind, geoms, n, min_depth)
+    F = evaluate(sl, center, (method,))[method]["F"]
+    if not np.isfinite(F).any():
+        return None
+    return geoms[int(np.nanargmin(F))]
+
+
+def _describe(prep, kind, geom, method, n, min_depth):
+    """M-P 는 같은 면 Spencer 값을 출발 참고값으로(엉뚱한 근 거름)."""
+    seed = None
+    if method == "morgenstern_price":
+        sl, center = _stack(prep, kind, [geom], n, min_depth)
+        sp = evaluate(sl, center, ("spencer",))["spencer"]["F"][0]
+        seed = float(sp) if np.isfinite(sp) else None
+    return describe(prep, kind, geom, method, n, min_depth, seed)
+
+
 def _theta(coarse, order):
     xc, zc, R, xE, xX = (v[order] for v in coarse)
     zE = zc - np.sqrt(np.clip(R**2 - (xE - xc) ** 2, 0, None))
@@ -630,5 +591,5 @@ def evaluate_specified(
         pts = sorted((s * float(x), float(z)) for x, z in polyline)
         geom, kind = (np.array([p[0] for p in pts]), np.array([p[1] for p in pts])), "polyline"
     for m in methods:
-        out[m] = describe(prep, kind, geom, m, n_slices, 0.0)
+        out[m] = _describe(prep, kind, geom, m, n_slices, 0.0)
     return out

E03_StabSection/E03_StabSection_Engine_Station.py

@@ -25,10 +25,16 @@ def analyze_station(
     methods=METHODS,
     noncircular: bool = True,
     min_depth: float = MIN_DEPTH_M,
+    case_workers: int = 1,
     **model_kw,
 ) -> dict:
+    """측점 하나 — case_workers > 1 이면 비탈 × 조건 경우를 프로세스로 나눠 돎(측점 하나만 다시 셀 때).
+
+    서버 프로세스 안에서는 1(차례로) 그대로 둠 — 자식이 main.py 를 다시 읽음(`Worker` 머리글).
+    """
     t0 = time.perf_counter()
     out = {"chainage_m": cross.get("chainage_m"), "label": cross.get("label"), "cases": []}
+    jobs = []
     for cond in conditions:
         model = section_model_from_b06(
             cross,
@@ -41,17 +47,24 @@ def analyze_station(
             **model_kw,
         )
         out.setdefault("areas", model["areas"])
-        for side in model["sides"]:
-            res = analyze(
-                model,
-                direction=side["direction"],
-                limits=side["limits"],
-                methods=methods,
-                noncircular=noncircular,
-                min_depth=min_depth,
-            )
-            mins = {m: (res["methods"][m]["min"] or {}).get("F") for m in methods}
-            out["cases"].append({"condition": cond, **side, "min_F": mins, "result": res})
+        jobs += [(cond, model, side) for side in model["sides"]]
+    kw = {"methods": methods, "noncircular": noncircular, "min_depth": min_depth}
+    if case_workers > 1 and len(jobs) > 1:
+        from concurrent.futures import ProcessPoolExecutor
+
+        with ProcessPoolExecutor(max_workers=min(case_workers, len(jobs))) as pool:
+            futures = [
+                pool.submit(analyze, m, direction=sd["direction"], limits=sd["limits"], **kw)
+                for _, m, sd in jobs
+            ]
+            results = [f.result() for f in futures]
+    else:
+        results = [
+            analyze(m, direction=sd["direction"], limits=sd["limits"], **kw) for _, m, sd in jobs
+        ]
+    for (cond, _m, side), res in zip(jobs, results):
+        mins = {m: (res["methods"][m]["min"] or {}).get("F") for m in methods}
+        out["cases"].append({"condition": cond, **side, "min_F": mins, "result": res})
     out["elapsed_s"] = round(time.perf_counter() - t0, 3)
     return out
 

E03_StabSection/E03_StabSection_Engine_Surface.py

@@ -0,0 +1,185 @@
+"""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

resources/tester/test_532_2_stab_lem_benchmarks.py

@@ -180,9 +180,8 @@ def test_acads_3a_weak_seam_noncircular_search():
     near(nonc(r, "spencer"), 1.258)
 
 
-def test_acads_4_inclined_seam_water_surcharge():
+def _acads4_model():
     # S2VM #9 (ACADS 4) — 기운 약층 · 지하수위선 · 등분포 20 kPa + 선형 20→40 kPa
-    # 최적화 비원호 Spencer 0.707 · Janbu 보정 0.699 (약층 좌표 = 저해상 그림 판독 — 출처 ±5 % 권함)
     ext = Polygon([[20, 15], [84, 15], [84, 40], [67.5, 40], [43, 27.75], [20, 27.75]])
     seam = Polygon([[20, 18.26], [84, 36.2], [84, 36.8], [20, 18.88]])
     regions = [{"material_id": "seam", "polygon": list(ext.intersection(seam).exterior.coords)}] + [
@@ -208,7 +207,12 @@ def test_acads_4_inclined_seam_water_surcharge():
             ],
         },
     }
-    r = analyze(m, direction=-1, methods=("janbu", "spencer"))
+    return m
+
+
+def test_acads_4_inclined_seam_water_surcharge():
+    # S2VM #9 — 최적화 비원호 Spencer 0.707 · Janbu 보정 0.699 (약층 좌표 = 저해상 그림 판독 — 출처 ±5 % 권함)
+    r = analyze(_acads4_model(), direction=-1, methods=("janbu", "spencer"))
     near(nonc(r, "spencer"), 0.707)
     near(nonc(r, "janbu"), 0.699)
 

resources/tester/test_532_2_stab_mp_outputs.py

@@ -0,0 +1,141 @@
+# -*- coding: utf-8 -*-
+"""532-2 2차 E03 2D — Morgenstern-Price(반정현) 대조 · 결과 자료(자유물체도 · F–λ 곡선 · Safety Map).
+
+기대값 = Rocscience Slide2 Slope Stability Verification Manual(지식DB 참고문헌) 원문 GLE · M-P 값.
+- #1 ACADS 1a GLE 0.986 · #3 ACADS 1c GLE 1.374 · #8 ACADS 3b 지정 비원호 GLE 1.262 ·
+  #9 ACADS 4 블록 + 최적화 GLE 0.683 · #21 F&K 지정 원 M-P 건기 2.075 · r_u 1.760(Slide 열).
+- 뺀 것: #7 · #15 — 원문 GLE(1.275 · 0.437)가 같은 표 Spencer(1.246 · 0.409)보다 높음(Slide 탐색이 GLE
+  임계면을 놓친 값 · 같은 문제 #8 지정면 GLE 는 Spencer 보다 낮음) → 대조 기준으로 못 씀.
+
+    ./venv/Scripts/python.exe -m pytest resources/tester/test_532_2_stab_mp_outputs.py -q
+"""
+
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+ROOT = Path(__file__).resolve().parents[2]
+sys.path.insert(0, str(ROOT))
+sys.path.insert(0, str(Path(__file__).resolve().parent))
+
+from E03_StabSection import E03_StabSection_Engine_GLE as gle  # noqa: E402
+from E03_StabSection import E03_StabSection_Engine_LEM as lem  # noqa: E402
+from E03_StabSection.E03_StabSection_Engine_Search import (  # noqa: E402
+    analyze,
+    evaluate_specified,
+)
+from E03_StabSection.E03_StabSection_Engine_Surface import (  # noqa: E402
+    poly_slices,
+    three_point_center,
+)
+from test_532_2_stab_lem_benchmarks import (  # noqa: E402
+    ACADS1C_MATS,
+    ACADS1C_REGIONS,
+    ACADS_EXT,
+    FK_EXT,
+    _acads3_model,
+    _acads4_model,
+    circ,
+    mats,
+    near,
+    nonc,
+)
+
+MP = "morgenstern_price"
+SEAM_3B = [[41.85, 27.75], [44.0, 26.5], [63.5, 27.0], [73.31, 40.0]]
+
+
+def test_mp_constant_function_equals_spencer():
+    # f(x) ≡ 1 이면 M-P 식 = Spencer 식 — 같은 F · λ
+    prep = lem.Prepared(_acads3_model(), -1)
+    pts = sorted((-x, z) for x, z in SEAM_3B)
+    vx, vz = np.array([[p[0] for p in pts]]), np.array([[p[1] for p in pts]])
+    sl = poly_slices(prep, vx, vz, 50)
+    center = three_point_center(vx, vz)
+    sp = lem.spencer(sl, center, np.array([1.3]))
+    mp = gle.morgenstern_price(sl, center, "constant", F0=np.array([1.3]))
+    assert abs(mp["F"][0] - sp["F"][0]) < 1e-6
+    assert abs(mp["lambda"][0] - sp["lambda"][0]) < 1e-5
+
+
+def test_mp_specified_surfaces_slide2():
+    # #8 지정 비원호 GLE 1.262 · #21 지정 원 M-P 건기 2.075 · r_u 0.25 1.760
+    r = evaluate_specified(_acads3_model(), polyline=SEAM_3B, direction=-1, methods=(MP,))
+    near(r[MP]["F"], 1.262, 0.01)
+    m = {"materials": mats(("s", 120, 600, 20)), "regions": [{"material_id": "s", "polygon": FK_EXT}]}
+    near(evaluate_specified(m, circle=[120, 90, 80], direction=1, methods=(MP,))[MP]["F"], 2.075, 0.01)
+    m["groundwater"] = {"mode": "ru", "ru": 0.25}
+    near(evaluate_specified(m, circle=[120, 90, 80], direction=1, methods=(MP,))[MP]["F"], 1.760, 0.01)
+
+
+def test_mp_search_slide2():
+    # #1 원호 GLE 0.986 · #3 원호 GLE 1.374
+    m = {"materials": mats(("s", 20, 3, 19.6)), "regions": [{"material_id": "s", "polygon": ACADS_EXT}]}
+    near(circ(analyze(m, direction=-1, noncircular=False), MP), 0.986, 0.01)
+    r = analyze({"materials": ACADS1C_MATS, "regions": ACADS1C_REGIONS}, direction=-1, noncircular=False)
+    near(circ(r, MP), 1.374, 0.01)
+
+
+def test_mp_acads4_noncircular():
+    # #9 블록 + 최적화 GLE 0.683 (약층 좌표 = 그림 판독 → ±2 %)
+    r = analyze(_acads4_model(), direction=-1, methods=("janbu", "spencer", MP))
+    near(nonc(r, MP), 0.683, 0.02)
+
+
+def _closure(rows, lam, f_last=1.0):
+    weight = sum(r["weight"] for r in rows)
+    for a, b in zip(rows, rows[1:]):
+        assert abs(a["E_right"] - b["E_left"]) < 1e-9 and abs(a["X_right"] - b["X_left"]) < 1e-9
+    assert abs(rows[-1]["E_right"]) < 1e-3 * weight  # 끝 절편에서 힘 평형이 닫힘
+    assert abs(rows[0]["X_left"]) < 1e-12
+    for r in rows:
+        # 연직 평형 — N·cosα + S·sinα = W + (X_right − X_left) (지진 · 균열 없음)
+        a = np.radians(r["alpha_deg"])
+        lhs = r["normal_total"] * np.cos(a) + r["shear_mobilized"] * np.sin(a)
+        assert abs(lhs - (r["weight"] + r["X_right"] - r["X_left"])) < 1e-6 * weight
+
+
+def test_free_body_and_lambda_curve_outputs():
+    m = {"materials": mats(("s", 20, 3, 19.6)), "regions": [{"material_id": "s", "polygon": ACADS_EXT}]}
+    r = evaluate_specified(m, circle=[29.104, 55.527, 30.527], direction=-1, methods=("spencer", MP))
+    for key in ("spencer", MP):
+        out = r[key]
+        rows = out["slices"]
+        assert {"normal_total", "pore_force", "shear_strength", "E_right", "X_right", "x_mid"} <= rows[0].keys()
+        _closure(rows, out["lambda"])
+        curve = out["lambda_curve"]
+        lam, fm, ff = (np.array(curve[k], float) for k in ("lambda", "Fm", "Ff"))
+        assert len(lam) == gle.CURVE_POINTS and np.isfinite(fm).sum() > 10
+        # 곡선 교점 = 해 (F · λ)
+        d = fm - ff
+        k = np.where(np.isfinite(d[:-1]) & np.isfinite(d[1:]) & (np.sign(d[:-1]) != np.sign(d[1:])))[0]
+        assert len(k)
+        j = k[np.argmin(np.abs(lam[k] - out["lambda"]))]
+        t = d[j] / (d[j] - d[j + 1])
+        assert abs(lam[j] + t * (lam[j + 1] - lam[j]) - out["lambda"]) < 0.02
+        assert abs(fm[j] + t * (fm[j + 1] - fm[j]) - out["F"]) < 0.01 * out["F"]
+    assert abs(r[MP]["slices"][-1]["X_right"]) < 1e-9  # 반정현 f = 0 (끝)
+
+
+def test_safety_map_and_reference_only():
+    m = {"materials": mats(("s", 17.64, 9.8, 10)), "regions": [{"material_id": "s", "polygon": [[0, 0], [25, 0], [25, 10], [15, 10], [5, 5], [0, 5]]}]}
+    r = analyze(m, direction=-1)
+    smap = r["safety_map"]["surfaces"]
+    assert 0 < len(smap) <= 50
+    F = [s["F"] for s in smap]
+    assert F == sorted(F) and F[0] >= 0.97 * circ(r, "bishop")
+    assert all(len(s["points"]) == 25 and s["points"][0][0] < s["points"][-1][0] for s in smap)
+    # 비원호 Fellenius · Bishop = 참고 칸 · 최소 F 는 원호에서만
+    for k in ("ordinary", "bishop"):
+        res = r["methods"][k]
+        assert res["noncircular"]["reference_only"] is True
+        assert res["min"]["kind"] == "circle"
+    assert r["methods"]["spencer"]["min"]["F"] <= circ(r, "spencer") + 1e-12
+
+
+if __name__ == "__main__":
+    sys.exit(pytest.main([__file__, "-q"]))