Aislo/common_util/common_util_mass_haul_curve.ts

/* =============================================================================
 * common_util_mass_haul_curve.ts
 * 누가토량 곡선의 기하 — 포물선 보간 · 수평선 교점 · 극값 추출.
 *
 * ── 측점 사이를 무엇으로 잇는가: **직선**(2026-08-02 사용자 확정) ──────────
 * 이론상으로는 포물선이다 — 평균단면법 `V = (A₁+A₂)/2 × L`이 "구간 내 단면적이 선형으로
 * 변한다"를 전제하므로, 그 적분인 누가토량은 2차식이 된다. 참고 도면의 종무대 평균운반거리에
 * `20/√2 = 14.14m`가 12번 찍히는 것도 포물선 + "종거 1/2 수평선" 규칙으로만 설명된다
 * (직선이면 `20/2 = 10m`).
 *
 * 그럼에도 **직선을 쓴다.** 측점을 기준으로 절·성토량을 읽고 판단하는 것이 실무이고, 도면도
 * 직선으로 그린다는 사용자(임도설계 실무자) 판단이다. 대가는 평균운반거리가 포물선 대비
 * 약 40% 짧게 나오는 것이며, 이는 20m·70m 장비 경계 판정에 직접 영향을 준다.
 *
 * 되돌릴 일이 생기면 `PARABOLIC_MASS_CURVE` 하나만 켜면 된다 — 그리기·재기가 함께 바뀐다.
 *
 * 토량 분배 계산(`common_util_mass_haul_balance`)이 700줄을 넘겨 곡선 기하만 여기로 떼어냈다.
 * ========================================================================== */

import type { MassHaulPoint } from "./common_util_mass_haul";

export const EPSILON = 1e-9;

/**
 * 누가토량 곡선을 측점 사이에서 포물선으로 볼지 여부. **false = 직선**(사용자 확정).
 * 이 값 하나가 곡선 렌더링(`curvePath`)과 수평선 교점(`crossFrom`)을 함께 바꾼다 —
 * 둘이 어긋나면 수평선 끝이 곡선에서 떠 보인다(2026-08-02에 이미 겪은 문제).
 */
export const PARABOLIC_MASS_CURVE = false;

/**
 * 구간 하나의 2차식 계수. 측점 사이에서 단면적이 선형으로 변하므로 누가토량은
 *   `V(t) = V₀ + a₀·t + (a₁ − a₀)·t² / (2Δ)`   (t = 구간 시작으로부터의 거리)
 * 이다. `V(Δ) = V₀ + (a₀ + a₁)Δ/2`라 평균단면법 결과와 정확히 같다.
 *
 * 순단면적이 없는 데이터(옛 payload·종단 기준 개략값)는 구간 평균 단면적을 양 끝에 같이
 * 넣어 자동으로 직선으로 떨어진다.
 */
export interface Segment {
  x0: number;
  span: number;
  v0: number;
  a0: number;
  curvature: number;
}

export function segmentOf(points: MassHaulPoint[], index: number): Segment | null {
  const a = points[index - 1];
  const b = points[index];
  const span = b.chainage_m - a.chainage_m;
  if (!(span > 0)) return null;
  const average = (b.cumulative_volume_m3 - a.cumulative_volume_m3) / span;
  const start = Number.isFinite(a.net_area_m2) ? a.net_area_m2 : average;
  const end = Number.isFinite(b.net_area_m2) ? b.net_area_m2 : average;
  // 포물선이 **구간 끝값을 정확히 통과하도록** 두 단면적을 같은 양만큼 옮긴다.
  // `V(Δ) = V₀ + (a₀+a₁)Δ/2`이므로 평균이 어긋나면 끝값이 어긋나고, 그러면 수평선 교점을
  // 못 찾아 블록이 엉뚱하게 길어진다(도면 곡선 검증에서 성토 1455㎥가 한 블록에 묻혔다).
  // 곡률 `(a₁−a₀)`는 그대로라 포물선 모양은 보존된다.
  const shift = average - (start + end) / 2;
  return {
    x0: a.chainage_m,
    span,
    v0: a.cumulative_volume_m3,
    // 직선 모드에서는 구간 평균 단면적 하나로 기울기를 잡는다 → 곡률 0 = 직선.
    a0: PARABOLIC_MASS_CURVE ? start + shift : average,
    curvature: PARABOLIC_MASS_CURVE ? (end - start) / (2 * span) : 0,
  };
}

export function evaluate(segment: Segment, t: number): number {
  return segment.v0 + segment.a0 * t + segment.curvature * t * t;
}

/** 구간 2차식이 `[tLo, tHi]` 안에서 `level`을 처음 지나는 t. 없으면 null. */
function solveSegment(segment: Segment, level: number, tLo: number, tHi: number): number | null {
  const c = segment.v0 - level;
  const roots: number[] = [];
  if (Math.abs(segment.curvature) < 1e-12) {
    if (Math.abs(segment.a0) > EPSILON) roots.push(-c / segment.a0);
  } else {
    const discriminant = segment.a0 * segment.a0 - 4 * segment.curvature * c;
    if (discriminant < 0) return null;
    const root = Math.sqrt(discriminant);
    roots.push((-segment.a0 - root) / (2 * segment.curvature));
    roots.push((-segment.a0 + root) / (2 * segment.curvature));
  }
  const inside = roots
    .filter((t) => t >= tLo - 1e-9 && t <= tHi + 1e-9)
    .sort((left, right) => left - right);
  return inside.length ? Math.min(Math.max(inside[0], tLo), tHi) : null;
}

/**
 * `[fromX, toX]`를 왼쪽부터 훑어 곡선이 `level`을 **처음** 지나는 x. 없으면 null.
 * 평형선·경계현·반종거 수평선의 교점이 전부 이 함수 하나로 나온다.
 */
export function crossFrom(
  points: MassHaulPoint[],
  fromX: number,
  toX: number,
  level: number,
): number | null {
  if (!(toX > fromX)) return null;
  for (let index = 1; index < points.length; index += 1) {
    const segment = segmentOf(points, index);
    if (!segment) continue;
    const x0 = Math.max(segment.x0, fromX);
    const x1 = Math.min(segment.x0 + segment.span, toX);
    if (!(x1 > x0)) continue;
    const t = solveSegment(segment, level, x0 - segment.x0, x1 - segment.x0);
    if (t !== null) return segment.x0 + t;
    // 구간 끝에서 부호가 바뀌었는데 2차 근을 못 찾으면(수치 한계) 직선으로 떨어뜨린다.
    // 교점을 놓치면 블록 경계가 통째로 밀리므로 여기서 반드시 하나는 돌려줘야 한다.
    const y0 = evaluate(segment, x0 - segment.x0);
    const y1 = evaluate(segment, x1 - segment.x0);
    if ((y0 - level) * (y1 - level) <= 0) {
      if (Math.abs(y1 - y0) < EPSILON) return x0;
      return x0 + ((level - y0) / (y1 - y0)) * (x1 - x0);
    }
  }
  return null;
}

/** 곡선의 극값 인덱스 수열(양 끝점 포함). 기울기 부호가 바뀌는 자리가 극값이다. */
export function extremaIndices(points: MassHaulPoint[]): number[] {
  const list = [0];
  let sign = 0;
  for (let index = 1; index < points.length; index += 1) {
    const delta = points[index].cumulative_volume_m3 - points[index - 1].cumulative_volume_m3;
    if (Math.abs(delta) < EPSILON) continue;
    const next = delta > 0 ? 1 : -1;
    if (sign !== 0 && next !== sign && list[list.length - 1] !== index - 1) list.push(index - 1);
    sign = next;
  }
  if (list[list.length - 1] !== points.length - 1) list.push(points.length - 1);
  return list;
}

/**
 * 잔진동 제거. 진폭이 `minSwing` 미만인 극값은 **짝이 되는 이웃과 함께** 지운다 —
 * 하나만 지우면 극대·극소 교대가 깨져 뒤 계산이 방향을 잃는다. 양 끝점은 남긴다.
 */
export function pruneExtrema(
  points: MassHaulPoint[],
  indices: number[],
  minSwing: number,
): number[] {
  const valueOf = (index: number): number => points[index].cumulative_volume_m3;
  let list = indices.slice();
  let changed = true;
  while (changed && list.length > 2) {
    changed = false;
    for (let i = 1; i < list.length - 1; i += 1) {
      const left = Math.abs(valueOf(list[i]) - valueOf(list[i - 1]));
      const right = Math.abs(valueOf(list[i]) - valueOf(list[i + 1]));
      if (Math.min(left, right) >= minSwing) continue;
      const partner = left <= right ? i - 1 : i + 1;
      const drop = new Set<number>([i]);
      if (partner > 0 && partner < list.length - 1) drop.add(partner);
      list = list.filter((_, position) => !drop.has(position));
      changed = true;
      break;
    }
  }
  return list;
}