from __future__ import annotations

from dataclasses import dataclass
import math
import re
import unicodedata
from typing import Iterable


class HydrologyError(ValueError):
    """Raised when hydrology input data is invalid or incomplete."""


def _text_width(value: object) -> int:
    text = str(value)
    width = 0
    for char in text:
        width += 2 if unicodedata.east_asian_width(char) in {"F", "W"} else 1
    return width


def _center_cell(value: object, width: int) -> str:
    text = str(value)
    pad = max(0, width - _text_width(text))
    left = pad // 2
    right = pad - left
    return " " * left + text + " " * right


def _fixed_table(rows: Iterable[Iterable[object]], min_width: int = 8) -> str:
    data = [[str(cell) for cell in row] for row in rows]
    if not data:
        return ""
    col_count = max(len(row) for row in data)
    normalized = [row + [""] * (col_count - len(row)) for row in data]
    widths = [
        max(min_width, max(_text_width(row[col]) for row in normalized) + 2)
        for col in range(col_count)
    ]
    column_gap = "    "
    lines = []
    for row in normalized:
        line = column_gap.join(_center_cell(row[col], widths[col]) for col in range(col_count)).rstrip()
        lines.append(line)
    return "\n".join(lines)


def format_fixed_table(rows: Iterable[Iterable[object]], min_width: int = 8) -> str:
    return _fixed_table(rows, min_width=min_width)


TIE_FREQUENCIES = ("100", "50", "30", "25", "20", "10", "5", "2")
TIE_FREQUENCY_PROBABILITIES = {
    "100": 1.0,
    "50": 2.0,
    "30": 100 / 30,
    "25": 4.0,
    "20": 5.0,
    "10": 10.0,
    "5": 20.0,
    "2": 50.0,
}
TIE_SOIL_CLASS_PARAMS = {
    "II": {
        "wet": {"r": 0.83, "r1": 0.56},
        "medium": {"r": 0.93, "r1": 0.63},
        "dry": {"r": 1.00, "r1": 0.68},
    },
    "III": {
        "wet": {"r": 0.93, "r1": 0.63},
        "medium": {"r": 1.02, "r1": 0.69},
        "dry": {"r": 1.08, "r1": 0.75},
    },
    "IV": {
        "wet": {"r": 0.98, "r1": 0.66},
        "medium": {"r": 1.10, "r1": 0.76},
        "dry": {"r": 1.16, "r1": 0.81},
    },
    "V": {
        "wet": {"r": 1.10, "r1": 0.76},
        "medium": {"r": 1.18, "r1": 0.83},
        "dry": {"r": 1.22, "r1": 0.87},
    },
    "VI": {
        "wet": {"r": 1.22, "r1": 0.87},
        "medium": {"r": 1.25, "r1": 0.90},
        "dry": {"r": 1.27, "r1": 0.92},
    },
}
TIE_Y_TABLE = [
    (0.06, 0.479), (0.07, 0.475), (0.08, 0.472), (0.09, 0.469), (0.10, 0.466),
    (0.12, 0.461), (0.14, 0.456), (0.16, 0.452), (0.18, 0.447), (0.20, 0.443),
    (0.22, 0.439), (0.24, 0.436), (0.26, 0.432), (0.28, 0.429), (0.30, 0.426),
    (0.32, 0.423), (0.34, 0.420), (0.36, 0.417), (0.38, 0.414), (0.40, 0.412),
    (0.42, 0.410), (0.44, 0.407), (0.46, 0.405), (0.48, 0.403), (0.50, 0.401),
    (0.52, 0.399), (0.54, 0.397), (0.56, 0.396), (0.58, 0.394), (0.60, 0.392),
    (0.62, 0.391), (0.64, 0.389), (0.66, 0.388), (0.68, 0.386), (0.70, 0.385),
    (0.72, 0.383), (0.74, 0.382), (0.76, 0.381), (0.78, 0.380), (0.80, 0.378),
    (0.82, 0.377), (0.84, 0.376), (0.86, 0.375), (0.88, 0.374), (0.90, 0.373),
    (0.92, 0.372), (0.94, 0.371), (0.96, 0.370), (0.98, 0.369), (1.00, 0.368),
    (1.02, 0.368), (1.04, 0.367), (1.06, 0.366), (1.08, 0.365), (1.10, 0.364),
    (1.12, 0.363), (1.14, 0.363), (1.16, 0.362), (1.18, 0.361), (1.20, 0.360),
    (1.25, 0.359), (1.30, 0.357), (1.35, 0.356), (1.40, 0.354), (1.45, 0.353),
    (1.50, 0.352), (1.55, 0.351), (1.60, 0.350), (1.65, 0.348), (1.70, 0.347),
]
TIE_GAMMA_TABLE = [
    (0.10, 2.41), (0.15, 2.36), (0.20, 2.31), (0.30, 2.26), (0.40, 2.22),
    (0.50, 2.19), (0.70, 2.15), (1.00, 2.10), (1.50, 2.05), (2.00, 2.01),
    (3.00, 1.97), (4.00, 1.93), (5.00, 1.91), (7.00, 1.87), (10.00, 1.83),
    (20.00, 1.76), (30.00, 1.71), (40.00, 1.68), (50.00, 1.66), (70.00, 1.63),
    (100.00, 1.59), (150.00, 1.56),
]
TIE_CN_R1_VALUES = [0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90]
TIE_CN_TABLE = {
    0.10: [1.04, 1.04, 1.03, 1.03, 1.02, 1.02, 1.01],
    0.12: [1.05, 1.05, 1.04, 1.03, 1.03, 1.02, 1.01],
    0.14: [1.07, 1.06, 1.06, 1.04, 1.03, 1.03, 1.02],
    0.16: [1.08, 1.07, 1.06, 1.05, 1.04, 1.03, 1.02],
    0.18: [1.09, 1.08, 1.07, 1.06, 1.04, 1.03, 1.02],
    0.20: [1.10, 1.09, 1.08, 1.06, 1.05, 1.04, 1.03],
    0.22: [1.11, 1.10, 1.09, 1.07, 1.06, 1.04, 1.03],
    0.24: [1.13, 1.11, 1.10, 1.08, 1.06, 1.05, 1.03],
    0.26: [1.14, 1.12, 1.11, 1.09, 1.07, 1.05, 1.04],
    0.28: [1.16, 1.14, 1.12, 1.10, 1.08, 1.06, 1.04],
    0.30: [1.17, 1.15, 1.13, 1.11, 1.09, 1.06, 1.04],
    0.32: [1.19, 1.17, 1.14, 1.12, 1.09, 1.07, 1.05],
    0.34: [1.21, 1.18, 1.16, 1.13, 1.10, 1.08, 1.05],
    0.36: [1.23, 1.20, 1.17, 1.14, 1.11, 1.08, 1.06],
    0.38: [1.25, 1.21, 1.18, 1.15, 1.12, 1.09, 1.06],
    0.40: [1.27, 1.23, 1.20, 1.17, 1.13, 1.10, 1.07],
    0.42: [1.29, 1.25, 1.22, 1.18, 1.14, 1.11, 1.07],
    0.44: [1.31, 1.27, 1.24, 1.20, 1.16, 1.12, 1.08],
    0.46: [1.34, 1.30, 1.26, 1.21, 1.17, 1.13, 1.08],
    0.48: [1.37, 1.32, 1.28, 1.23, 1.18, 1.14, 1.09],
    0.50: [1.40, 1.36, 1.30, 1.26, 1.20, 1.15, 1.10],
    0.52: [1.43, 1.38, 1.33, 1.27, 1.22, 1.16, 1.11],
    0.54: [1.47, 1.41, 1.35, 1.29, 1.24, 1.18, 1.12],
    0.56: [1.51, 1.45, 1.38, 1.32, 1.26, 1.19, 1.13],
    0.58: [1.55, 1.48, 1.41, 1.34, 1.28, 1.21, 1.14],
    0.60: [1.60, 1.52, 1.45, 1.37, 1.30, 1.22, 1.15],
    0.62: [1.65, 1.57, 1.49, 1.41, 1.32, 1.24, 1.16],
    0.64: [1.71, 1.62, 1.53, 1.44, 1.36, 1.27, 1.18],
    0.66: [1.78, 1.68, 1.58, 1.48, 1.39, 1.29, 1.20],
    0.68: [1.85, 1.74, 1.63, 1.53, 1.43, 1.32, 1.21],
    0.70: [1.93, 1.82, 1.70, 1.58, 1.47, 1.35, 1.23],
}
TIE_P1_GAMMA_VALUES = [1.50, 1.55, 1.60, 1.65, 1.70, 1.75, 1.80, 1.85, 1.90, 1.95, 2.00, 2.05, 2.10, 2.15, 2.20, 2.25, 2.30, 2.35, 2.40, 2.45, 2.50]
TIE_P1_N_VALUES = [0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95]
TIE_P1_TABLE = [
    [0.86, 0.83, 0.79, 0.74, 0.68, 0.62, 0.53, 0.44, 0.33, 0.18],
    [0.84, 0.80, 0.76, 0.71, 0.65, 0.59, 0.50, 0.41, 0.30, 0.16],
    [0.82, 0.78, 0.73, 0.68, 0.63, 0.56, 0.47, 0.38, 0.27, 0.15],
    [0.80, 0.75, 0.71, 0.66, 0.60, 0.53, 0.45, 0.36, 0.25, 0.14],
    [0.78, 0.73, 0.69, 0.63, 0.58, 0.51, 0.43, 0.34, 0.25, 0.13],
    [0.76, 0.71, 0.67, 0.61, 0.55, 0.49, 0.41, 0.33, 0.24, 0.13],
    [0.74, 0.69, 0.64, 0.59, 0.53, 0.46, 0.39, 0.31, 0.22, 0.12],
    [0.72, 0.67, 0.62, 0.57, 0.51, 0.44, 0.37, 0.29, 0.21, 0.11],
    [0.70, 0.65, 0.60, 0.55, 0.49, 0.43, 0.36, 0.28, 0.20, 0.11],
    [0.68, 0.64, 0.59, 0.54, 0.48, 0.41, 0.34, 0.27, 0.19, 0.10],
    [0.67, 0.62, 0.57, 0.52, 0.46, 0.40, 0.33, 0.26, 0.18, 0.10],
    [0.65, 0.61, 0.56, 0.50, 0.45, 0.39, 0.32, 0.25, 0.17, 0.09],
    [0.63, 0.59, 0.54, 0.48, 0.43, 0.37, 0.31, 0.24, 0.16, 0.08],
    [0.62, 0.58, 0.53, 0.47, 0.42, 0.36, 0.30, 0.23, 0.16, 0.08],
    [0.61, 0.56, 0.51, 0.46, 0.41, 0.35, 0.29, 0.22, 0.15, 0.08],
    [0.60, 0.55, 0.50, 0.45, 0.40, 0.34, 0.28, 0.21, 0.14, 0.07],
    [0.58, 0.53, 0.49, 0.44, 0.39, 0.33, 0.27, 0.20, 0.14, 0.07],
    [0.57, 0.52, 0.48, 0.43, 0.38, 0.32, 0.27, 0.20, 0.13, 0.07],
    [0.56, 0.51, 0.47, 0.42, 0.37, 0.31, 0.26, 0.19, 0.13, 0.07],
    [0.55, 0.50, 0.46, 0.41, 0.36, 0.30, 0.25, 0.19, 0.13, 0.07],
    [0.53, 0.49, 0.44, 0.39, 0.34, 0.29, 0.24, 0.18, 0.12, 0.06],
]
TIE_C_S_VALUES = [20, 40, 70, 100, 200]
TIE_C_MOISTURE_CORRECTION = {
    "II": {"wet": 1.12, "medium": 1.00, "dry": 0.87},
    "III": {"wet": 1.26, "medium": 1.00, "dry": 0.80},
    "IV": {"wet": 1.25, "medium": 1.00, "dry": 0.80},
    "V": {"wet": 1.40, "medium": 1.00, "dry": 0.70},
    "VI": {"wet": 1.60, "medium": 1.00, "dry": 0.60},
}
TIE_C_AREA_RANGES = [
    {"高低山": (0.01, 1.0)},
    {"高低山": (1.01, 6.0), "丘陵": (0.01, 1.0)},
    {"高低山": (6.01, 20.0), "丘陵": (1.01, 6.0), "平坦": (0.01, 1.0)},
    {"高低山": (20.01, 50.0), "丘陵": (5.01, 20.0), "平坦": (1.01, 5.0)},
    {"高低山": (50.01, 100.0), "丘陵": (20.01, 50.0), "平坦": (5.01, 20.0)},
    {"丘陵": (50.01, 100.0), "平坦": (20.01, 50.0)},
    {"平坦": (50.01, 1000.0)},
]
TIE_C_TABLE = {
    "II": [
        {"t": 0.1, "n04": [0.78, 0.83, 0.86, 0.88, 0.91], "n07": [0.87, 0.87, 0.90, 0.91, 0.93]},
        {"t": 0.2, "n04": [0.76, 0.81, 0.85, 0.87, 0.90], "n07": [0.80, 0.84, 0.87, 0.89, 0.91]},
        {"t": 0.4, "n04": [0.73, 0.79, 0.83, 0.85, 0.89], "n07": [0.76, 0.81, 0.85, 0.87, 0.90]},
        {"t": 0.6, "n04": [0.72, 0.78, 0.82, 0.84, 0.88], "n07": [0.73, 0.79, 0.83, 0.85, 0.89]},
        {"t": 0.8, "n04": [0.70, 0.77, 0.81, 0.83, 0.87], "n07": [0.71, 0.78, 0.82, 0.84, 0.88]},
        {"t": 1.0, "n04": [0.69, 0.76, 0.81, 0.83, 0.87], "n07": [0.69, 0.76, 0.81, 0.83, 0.87]},
        {"t": 2.5, "n04": [0.65, 0.73, 0.78, 0.81, 0.85], "n07": [0.61, 0.70, 0.76, 0.79, 0.83]},
    ],
    "III": [
        {"t": 0.1, "n04": [0.70, 0.76, 0.80, 0.82, 0.85], "n07": [0.76, 0.80, 0.83, 0.85, 0.88]},
        {"t": 0.2, "n04": [0.67, 0.73, 0.78, 0.80, 0.84], "n07": [0.72, 0.77, 0.81, 0.83, 0.86]},
        {"t": 0.4, "n04": [0.64, 0.71, 0.76, 0.78, 0.82], "n07": [0.67, 0.73, 0.78, 0.80, 0.84]},
        {"t": 0.8, "n04": [0.61, 0.68, 0.73, 0.76, 0.81], "n07": [0.62, 0.69, 0.74, 0.77, 0.81]},
        {"t": 1.0, "n04": [0.60, 0.67, 0.73, 0.76, 0.80], "n07": [0.60, 0.68, 0.73, 0.76, 0.80]},
        {"t": 1.5, "n04": [0.58, 0.66, 0.71, 0.74, 0.79], "n07": [0.56, 0.66, 0.70, 0.73, 0.78]},
        {"t": 3.0, "n04": [0.54, 0.63, 0.69, 0.72, 0.77], "n07": [0.49, 0.59, 0.65, 0.69, 0.75]},
    ],
    "IV": [
        {"t": 0.1, "n04": [0.57, 0.64, 0.68, 0.71, 0.75], "n07": [0.64, 0.69, 0.73, 0.75, 0.79]},
        {"t": 0.2, "n04": [0.54, 0.61, 0.66, 0.69, 0.74], "n07": [0.59, 0.65, 0.70, 0.72, 0.77]},
        {"t": 0.4, "n04": [0.51, 0.58, 0.64, 0.67, 0.72], "n07": [0.54, 0.61, 0.66, 0.69, 0.74]},
        {"t": 0.8, "n04": [0.48, 0.56, 0.61, 0.64, 0.70], "n07": [0.48, 0.56, 0.62, 0.66, 0.70]},
        {"t": 1.0, "n04": [0.46, 0.56, 0.60, 0.64, 0.69], "n07": [0.46, 0.55, 0.60, 0.64, 0.69]},
        {"t": 1.5, "n04": [0.44, 0.53, 0.59, 0.62, 0.68], "n07": [0.43, 0.51, 0.58, 0.61, 0.67]},
        {"t": 3.0, "n04": [0.40, 0.50, 0.56, 0.60, 0.66], "n07": [0.38, 0.46, 0.52, 0.56, 0.63]},
    ],
    "V": [
        {"t": 0.1, "n04": [0.39, 0.46, 0.51, 0.54, 0.59], "n07": [0.46, 0.52, 0.56, 0.59, 0.64]},
        {"t": 0.2, "n04": [0.36, 0.44, 0.49, 0.52, 0.57], "n07": [0.41, 0.48, 0.53, 0.56, 0.60]},
        {"t": 0.4, "n04": [0.33, 0.41, 0.46, 0.49, 0.55], "n07": [0.36, 0.44, 0.49, 0.52, 0.57]},
        {"t": 0.8, "n04": [0.30, 0.38, 0.44, 0.47, 0.53], "n07": [0.31, 0.39, 0.44, 0.48, 0.53]},
        {"t": 1.0, "n04": [0.29, 0.37, 0.43, 0.46, 0.52], "n07": [0.29, 0.37, 0.43, 0.46, 0.52]},
        {"t": 2.0, "n04": [0.26, 0.34, 0.40, 0.44, 0.50], "n07": [0.23, 0.32, 0.38, 0.41, 0.48]},
        {"t": 3.5, "n04": [0.23, 0.31, 0.38, 0.41, 0.48], "n07": [0.18, 0.27, 0.34, 0.37, 0.44]},
    ],
    "VI": [
        {"t": 0.1, "n04": [0.16, 0.21, 0.26, 0.28, 0.33], "n07": [0.21, 0.26, 0.30, 0.33, 0.37]},
        {"t": 0.2, "n04": [0.13, 0.19, 0.23, 0.26, 0.31], "n07": [0.17, 0.23, 0.27, 0.30, 0.34]},
        {"t": 0.4, "n04": [0.11, 0.17, 0.21, 0.24, 0.29], "n07": [0.13, 0.19, 0.23, 0.26, 0.31]},
        {"t": 0.8, "n04": [0.08, 0.14, 0.19, 0.22, 0.27], "n07": [0.09, 0.15, 0.20, 0.22, 0.28]},
        {"t": 1.5, "n04": [0.06, 0.12, 0.17, 0.20, 0.25], "n07": [0.05, 0.11, 0.16, 0.19, 0.24]},
        {"t": 3.0, "n04": [0.03, 0.10, 0.15, 0.18, 0.23], "n07": [0.00, 0.06, 0.12, 0.16, 0.21]},
        {"t": 4.0, "n04": [0.02, 0.09, 0.14, 0.17, 0.22], "n07": [0.00, 0.05, 0.10, 0.13, 0.19]},
    ],
}
SMALL_RIVER_FREQUENCIES = ("10", "5", "2", "1")
SMALL_RIVER_DURATION_RULE_TEXT = (
    "H采用的暴雨历时由面积F控制：F<20km²取10分钟，20~50km²取30分钟，"
    "50~300km²取1小时，300~1000km²取3小时，1000~5000km²取6小时，"
    ">5000km²取24小时。面积接近分界值时，可用相邻两个历时试算后合理选取。"
)
SMALL_RIVER_RAINFALL_SCOPE_TEXT = (
    "公式中的H必须是面积对应历时、同一频率的流域平均雨量(mm)。推荐直接输入各频率流域平均设计H；"
    "若只有流域平均雨量均值H0，可按Cv、Cs=3.5Cv计算频率模比系数Kp，再用H=H0×Kp。"
    "单个雨量站或单点等值线值不能直接代表流域平均雨量。"
)
REGION_FREQUENCIES = ("0.1%", "0.2%", "0.5%", "1%", "2%", "5%", "10%")
INSTANT_UNIT_LINE_FREQUENCIES = TIE_FREQUENCIES
DESIGN_FREQUENCIES_7 = INSTANT_UNIT_LINE_FREQUENCIES
DESIGN_FREQUENCY_PROBABILITIES = {
    "100": 1.0,
    "50": 2.0,
    "30": 100 / 30,
    "25": 4.0,
    "20": 5.0,
    "10": 10.0,
    "5": 20.0,
    "2": 50.0,
    "0.1%": 0.1,
    "0.2%": 0.2,
    "0.5%": 0.5,
    "1%": 1.0,
    "2%": 2.0,
    "5%": 5.0,
    "10%": 10.0,
}


def _require_positive(name: str, value: float) -> float:
    try:
        number = float(value)
    except (TypeError, ValueError) as exc:
        raise HydrologyError(f"{name} 必须是数字") from exc
    if not math.isfinite(number) or number <= 0:
        raise HydrologyError(f"{name} 必须大于 0")
    return number


def _require_non_negative(name: str, value: float) -> float:
    try:
        number = float(value)
    except (TypeError, ValueError) as exc:
        raise HydrologyError(f"{name} 必须是数字") from exc
    if not math.isfinite(number) or number < 0:
        raise HydrologyError(f"{name} 必须大于等于 0")
    return number


def _safe_round(value: float, digits: int = 2) -> float:
    if not math.isfinite(value):
        raise HydrologyError("计算结果不是有效数字，请检查输入参数")
    return round(value, digits)


@dataclass(frozen=True)
class TieyuanInput:
    basin_name: str
    area: float
    river_length_l1: float
    river_slope_j1: float
    rainfall_24h: dict[str, float]
    slope_length_l2: float | None = None
    slope_slope_i2: float | None = None
    a1: float = 0.341
    a2: float = 0.025
    n: float = 0.65
    soil_class: str = "III"
    moisture_condition: str = "medium"
    runoff_coefficient: float | None = 0.86
    runoff_coefficient_auto: bool = False
    c_landform: str = "高低山"
    r_param: float = 1.02
    r1_param: float = 0.69
    beta1: float = 0.65
    rainfall_24h_is_areal: bool = False


@dataclass(frozen=True)
class TieyuanFrequencyParams:
    rainfall_24h: float
    storm_parameter_s: float
    runoff_coefficient: float
    c_table_t: float | None
    loss_factor_k2: float
    cn: float
    n_prime: float
    p1: float
    p: float
    q: float | None


@dataclass(frozen=True)
class TieyuanResult:
    basin_name: str
    a1: float
    a2: float
    n: float
    soil_class: str
    moisture_condition: str
    runoff_coefficient: float | None
    runoff_coefficient_auto: bool
    c_landform: str
    r_param: float
    r1_param: float
    beta1: float
    l2: float
    j2: float
    k1: float
    k2: float
    x: float
    gamma: float
    y: float
    flows: dict[str, float | None]
    mu_values: dict[str, float]
    frequency_params: dict[str, TieyuanFrequencyParams]
    warnings: list[str]


def _linear_interpolate(points: list[tuple[float, float]], x_value: float) -> float:
    points = sorted(points)
    if x_value <= points[0][0]:
        return points[0][1]
    if x_value >= points[-1][0]:
        return points[-1][1]
    for (x0, y0), (x1, y1) in zip(points, points[1:]):
        if x0 <= x_value <= x1:
            if x1 == x0:
                return y0
            ratio = (x_value - x0) / (x1 - x0)
            return y0 + ratio * (y1 - y0)
    return points[-1][1]


def _bilinear_interpolate(x_values: list[float], y_values: list[float], table: list[list[float]], x_value: float, y_value: float) -> float:
    x_value = min(max(x_value, x_values[0]), x_values[-1])
    y_value = min(max(y_value, y_values[0]), y_values[-1])

    def bounds(values: list[float], value: float) -> tuple[int, int]:
        if value <= values[0]:
            return 0, 0
        if value >= values[-1]:
            last = len(values) - 1
            return last, last
        for idx in range(len(values) - 1):
            if values[idx] <= value <= values[idx + 1]:
                return idx, idx + 1
        last = len(values) - 1
        return last, last

    x0_idx, x1_idx = bounds(x_values, x_value)
    y0_idx, y1_idx = bounds(y_values, y_value)
    x0, x1 = x_values[x0_idx], x_values[x1_idx]
    y0, y1 = y_values[y0_idx], y_values[y1_idx]
    q00 = table[x0_idx][y0_idx]
    q01 = table[x0_idx][y1_idx]
    q10 = table[x1_idx][y0_idx]
    q11 = table[x1_idx][y1_idx]
    if x0 == x1 and y0 == y1:
        return q00
    if x0 == x1:
        return q00 + (q01 - q00) * ((y_value - y0) / (y1 - y0))
    if y0 == y1:
        return q00 + (q10 - q00) * ((x_value - x0) / (x1 - x0))
    tx = (x_value - x0) / (x1 - x0)
    ty = (y_value - y0) / (y1 - y0)
    return (q00 * (1 - tx) * (1 - ty)) + (q10 * tx * (1 - ty)) + (q01 * (1 - tx) * ty) + (q11 * tx * ty)


def tieyuan_y_from_ratio(k1_over_k2: float) -> float:
    return _safe_round(_linear_interpolate(TIE_Y_TABLE, k1_over_k2), 3)


def tieyuan_gamma_from_k_sum(k_sum: float) -> float:
    return _safe_round(_linear_interpolate(TIE_GAMMA_TABLE, k_sum), 3)


def tieyuan_cn_from_loss_factor(loss_factor_k2: float, r1: float) -> float:
    x_values = sorted(TIE_CN_TABLE)
    table = [TIE_CN_TABLE[x] for x in x_values]
    return _safe_round(_bilinear_interpolate(x_values, TIE_CN_R1_VALUES, table, loss_factor_k2, r1), 3)


def tieyuan_p1_from_gamma_n(gamma: float, n_prime: float) -> float:
    return _safe_round(_bilinear_interpolate(TIE_P1_GAMMA_VALUES, TIE_P1_N_VALUES, TIE_P1_TABLE, gamma, n_prime), 3)


def tieyuan_runoff_coefficient_from_table(
    soil_class: str,
    moisture_condition: str,
    landform: str,
    area: float,
    n: float,
    storm_parameter_s: float,
) -> tuple[float, float]:
    soil_class = soil_class.strip().upper()
    if soil_class not in TIE_C_TABLE:
        raise HydrologyError("表3-3仅支持土类 II、III、IV、V、VI")
    moisture_condition = moisture_condition.strip().lower()
    if moisture_condition not in TIE_C_MOISTURE_CORRECTION[soil_class]:
        raise HydrologyError("前期土壤水分必须为 wet、medium、dry")
    landform = landform.strip() or "高低山"
    if landform not in ("高低山", "丘陵", "平坦"):
        raise HydrologyError("表3-3地貌类型必须为 高低山、丘陵、平坦")
    area = _require_positive("面积 F", area)
    storm_parameter_s = _require_positive("暴雨参数 S", storm_parameter_s)

    row_index: int | None = None
    for idx, ranges in enumerate(TIE_C_AREA_RANGES):
        if landform not in ranges:
            continue
        low, high = ranges[landform]
        if low <= area <= high:
            row_index = idx
            break
    if row_index is None:
        raise HydrologyError(f"表3-3中未找到 {landform}、F={area:g} km² 对应的面积范围")

    row = TIE_C_TABLE[soil_class][row_index]
    column_key = "n04" if n <= 0.55 else "n07"
    base_c = _linear_interpolate(list(zip(TIE_C_S_VALUES, row[column_key])), storm_parameter_s)
    correction = TIE_C_MOISTURE_CORRECTION[soil_class][moisture_condition]
    c_value = min(1.0, max(0.0, base_c * correction))
    return _safe_round(c_value, 3), row["t"]


def tieyuan_peak_flow_formula(area: float, n: float, runoff_coefficient: float, storm_parameter_s: float, x_factor: float, p1: float, p: float, y: float) -> float:
    k1 = 0.278 * storm_parameter_s * area
    denominator = (x_factor * p1) ** n
    exponent = 1 / (1 - n * y)
    return (k1 * runoff_coefficient * p / denominator) ** exponent


def calculate_tieyuan(data: TieyuanInput) -> TieyuanResult:
    area = _require_positive("面积 F", data.area)
    l1 = _require_positive("河长 L1", data.river_length_l1)
    j1 = _require_positive("河坡 J1", data.river_slope_j1)
    a1 = _require_positive("A1", data.a1)
    a2 = _require_positive("A2", data.a2)
    n = _require_positive("n", data.n)
    soil_class = data.soil_class.strip().upper() or "III"
    if soil_class not in TIE_SOIL_CLASS_PARAMS:
        raise HydrologyError("土类必须为 II、III、IV、V、VI")
    moisture_condition = data.moisture_condition.strip().lower() or "medium"
    if moisture_condition not in TIE_SOIL_CLASS_PARAMS[soil_class]:
        raise HydrologyError("前期土壤水分必须为 wet、medium、dry")
    r_param = _require_positive("R", data.r_param)
    r1_param = _require_positive("r1", data.r1_param)
    runoff_coefficient = None if data.runoff_coefficient_auto else _require_positive("径流系数 C", data.runoff_coefficient)
    c_landform = data.c_landform.strip() or "高低山"
    if c_landform not in ("高低山", "丘陵", "平坦"):
        raise HydrologyError("表3-3地貌类型必须为 高低山、丘陵、平坦")
    beta1 = _require_positive("beta1", data.beta1)

    l2 = _require_positive("坡面长度 L2", data.slope_length_l2) if data.slope_length_l2 is not None else area / (2.7 * l1)
    j2 = _require_positive("山坡平均坡度 I2", data.slope_slope_i2) if data.slope_slope_i2 is not None else j1 * 0.6
    k1 = 0.278 * l1 / (a1 * (j1**0.35))
    k2 = 0.278 * (l2**0.5) * (area**0.5) / (a2 * (j2**0.33))
    x_factor = k1 + 0.95 * k2
    gamma = tieyuan_gamma_from_k_sum(k1 + k2)
    y_value = tieyuan_y_from_ratio(k1 / k2)

    flows: dict[str, float | None] = {}
    mu_values: dict[str, float] = {}
    frequency_params: dict[str, TieyuanFrequencyParams] = {}
    warnings: list[str] = []
    for freq in TIE_FREQUENCIES:
        h24 = float(data.rainfall_24h.get(freq, 0) or 0)
        if h24 <= 0:
            flows[freq] = None
            warnings.append(f"缺少 H₂₄_{freq}，未计算 Q{freq}")
            continue
        storm_parameter_s = h24 * (24 ** (n - 1))
        if data.runoff_coefficient_auto:
            c_value, c_table_t = tieyuan_runoff_coefficient_from_table(
                soil_class,
                moisture_condition,
                c_landform,
                area,
                n,
                storm_parameter_s,
            )
        else:
            c_value = runoff_coefficient
            c_table_t = None
        loss_factor_k2 = r_param * (storm_parameter_s ** (r1_param - 1))
        cn = tieyuan_cn_from_loss_factor(loss_factor_k2, r1_param)
        n_prime = cn * n
        p1 = tieyuan_p1_from_gamma_n(gamma, n_prime)
        p_value = 1 - ((1 - p1) ** gamma)
        flow = tieyuan_peak_flow_formula(area, n, c_value, storm_parameter_s, x_factor, p1, p_value, y_value)
        flows[freq] = _safe_round(flow)
        mu_values[freq] = _safe_round(loss_factor_k2, 3)
        frequency_params[freq] = TieyuanFrequencyParams(
            rainfall_24h=_safe_round(h24, 3),
            storm_parameter_s=_safe_round(storm_parameter_s, 3),
            runoff_coefficient=_safe_round(c_value, 3),
            c_table_t=None if c_table_t is None else _safe_round(c_table_t, 3),
            loss_factor_k2=_safe_round(loss_factor_k2, 3),
            cn=cn,
            n_prime=_safe_round(n_prime, 3),
            p1=p1,
            p=_safe_round(p_value, 3),
            q=flows[freq],
        )

    if "2" not in data.rainfall_24h or float(data.rainfall_24h.get("2", 0) or 0) <= 0:
        warnings.append("Python 版不再用 H₂₄_5 代替 H₂₄_2；需要 H₂₄_2 才计算 Q2")

    return TieyuanResult(
        basin_name=data.basin_name,
        a1=a1,
        a2=a2,
        n=n,
        soil_class=soil_class,
        moisture_condition=moisture_condition,
        runoff_coefficient=runoff_coefficient,
        runoff_coefficient_auto=data.runoff_coefficient_auto,
        c_landform=c_landform,
        r_param=r_param,
        r1_param=r1_param,
        beta1=beta1,
        l2=_safe_round(l2, 3),
        j2=_safe_round(j2, 3),
        k1=_safe_round(k1, 3),
        k2=_safe_round(k2, 3),
        x=_safe_round(x_factor, 3),
        gamma=_safe_round(gamma, 3),
        y=_safe_round(y_value, 3),
        flows=flows,
        mu_values=mu_values,
        frequency_params=frequency_params,
        warnings=warnings,
    )


SMALL_RIVER_COEFFICIENTS: dict[int, dict[str, tuple[float, float, float, float, float, float]]] = {
    1: {"10": (1.272, 1.221, 0.47, 0.9, 0.57, 0.47), "5": (1.272, 1.221, 0.47, 0.9, 0.57, 0.47), "2": (1.272, 1.321, 0.47, 0.9, 0.57, 0.47), "1": (1.272, 1.371, 0.47, 0.9, 0.57, 0.47)},
    2: {"10": (1.074, 0.082, 0.65, 0.9, 0.69, -0.67), "5": (1.074, 0.152, 0.65, 0.9, 0.69, -0.67), "2": (1.074, 0.252, 0.65, 0.9, 0.69, -0.67), "1": (1.074, 0.252, 0.65, 0.9, 0.69, -0.67)},
    3: {"10": (0.707, 0.6, 0.65, 0.9, 0.67, -0.64), "5": (0.707, 0.65, 0.65, 0.9, 0.67, -0.64), "2": (0.707, 0.79, 0.65, 0.9, 0.67, -0.64), "1": (0.707, 0.89, 0.65, 0.9, 0.67, -0.64)},
    4: {"10": (1.144, -0.55, 0.65, 0.9, 0.67, -0.38), "5": (1.144, -0.519, 0.65, 0.9, 0.67, -0.38), "2": (1.144, -0.319, 0.65, 0.9, 0.67, -0.38), "1": (1.144, -0.419, 0.65, 0.9, 0.67, -0.38)},
    5: {"10": (1.508, -0.392, 0.5, 0.9, 0.57, -0.19), "5": (1.508, -0.342, 0.5, 0.9, 0.57, -0.19), "2": (1.508, -0.292, 0.5, 0.9, 0.57, -0.19), "1": (1.508, -0.292, 0.5, 0.9, 0.57, -0.19)},
    6: {"10": (1.53, -0.38, 0.65, 0.9, 0.67, 0.47), "5": (1.53, -0.28, 0.65, 0.9, 0.67, 0.47), "2": (1.53, -0.38, 0.65, 0.9, 0.67, 0.47), "1": (1.53, -0.38, 0.65, 0.9, 0.67, 0.47)},
    7: {"10": (1.662, -1.04, 0.35, 0.9, 0.47, -0.51), "5": (1.662, -1.14, 0.35, 0.9, 0.47, -0.51), "2": (1.662, -0.84, 0.35, 0.9, 0.47, -0.51), "1": (1.662, -0.84, 0.35, 0.9, 0.47, -0.51)},
    8: {"10": (1.207, -1.255, 0.32, 0.9, -0.182, 0.308), "5": (1.207, -1.155, 0.32, 0.9, -0.182, 0.308), "2": (1.207, -1.1, 0.32, 0.9, -0.182, 0.308), "1": (1.207, -1.055, 0.32, 0.9, -0.182, 0.308)},
}

SMALL_RIVER_REGION_NAMES = {
    1: "黄土沟壑区",
    2: "石山林区",
    3: "土山林区",
    4: "草原区",
    5: "土石山沟壑混合区",
    6: "祁连山区",
    7: "荒漠区",
    8: "河西走廊区",
}


def small_river_region_label(region: int) -> str:
    if region not in SMALL_RIVER_REGION_NAMES:
        raise HydrologyError("中小河流分区必须为 1-8")
    return f"{region}. {SMALL_RIVER_REGION_NAMES[region]}"


def small_river_parameter_rows(region: int) -> list[tuple[str, float, float, float, float, float, float]]:
    if region not in SMALL_RIVER_COEFFICIENTS:
        raise HydrologyError("中小河流分区必须为 1-8")
    return [
        (freq, *SMALL_RIVER_COEFFICIENTS[region][freq])
        for freq in SMALL_RIVER_FREQUENCIES
    ]


def rainfall_duration(area: float) -> str:
    area = _require_positive("面积 F", area)
    if area < 20:
        return "10分钟"
    if area <= 50:
        return "30分钟"
    if area <= 300:
        return "1小时"
    if area <= 1000:
        return "3小时"
    if area <= 5000:
        return "6小时"
    return "24小时"


def kp_pearson_iii(cv: float, cs: float, probability_percent: int) -> float:
    cv = _require_positive("Cv", cv)
    if probability_percent not in (1, 2, 5, 10):
        raise HydrologyError("频率必须是 1、2、5、10")
    cs_ratio = cs / cv
    return pearson_iii_kp(cv, cs_ratio, probability_percent)


def small_river_peak_flow(region: int, frequency: str, area: float, rainfall: float, slope: float, shape_factor: float) -> float:
    if region not in SMALL_RIVER_COEFFICIENTS:
        raise HydrologyError("中小河流分区必须为 1-8")
    if frequency not in SMALL_RIVER_COEFFICIENTS[region]:
        raise HydrologyError("中小河流频率必须为 10、5、2、1")
    area = _require_positive("面积 F", area)
    rainfall = _require_positive(f"H{frequency}", rainfall)
    slope = _require_positive("比降 J", slope)
    shape_factor = _require_positive("形状系数 f", shape_factor)

    k, b, n, m, alpha, beta = SMALL_RIVER_COEFFICIENTS[region][frequency]
    base = (area**n) * (rainfall**m) * (slope**alpha) * (shape_factor**beta)
    return 10 ** (k * (math.log(base) / math.log(10)) + b)


@dataclass(frozen=True)
class SmallRiverInput:
    basin_name: str
    region: int
    area: float
    river_length: float
    slope: float
    cv: float
    rainfall: dict[str, float]
    rainfall_is_mean: bool = True


@dataclass(frozen=True)
class SmallRiverResult:
    basin_name: str
    shape_factor: float
    duration: str
    kp_values: dict[str, float | None]
    design_rainfall: dict[str, float]
    flows: dict[str, float]


def calculate_small_river(data: SmallRiverInput) -> SmallRiverResult:
    area = _require_positive("面积 F", data.area)
    length = _require_positive("河长 L", data.river_length)
    slope = _require_positive("比降 J", data.slope)
    shape = area / (length * length)
    cv = _require_positive("Cv", data.cv) if data.rainfall_is_mean else None
    cs = 3.5 * cv if cv is not None else None

    kp_values: dict[str, float | None] = {}
    design_rainfall: dict[str, float] = {}
    flows: dict[str, float] = {}
    for freq in SMALL_RIVER_FREQUENCIES:
        rainfall = _require_positive(f"H{freq}", data.rainfall.get(freq, 0))
        kp = kp_pearson_iii(cv, cs, int(freq)) if cv is not None and cs is not None else None
        kp_values[freq] = _safe_round(kp, 3) if kp is not None else None
        h_design = rainfall * kp if kp is not None else rainfall
        design_rainfall[freq] = _safe_round(h_design, 3)
        flows[freq] = _safe_round(small_river_peak_flow(data.region, freq, area, h_design, slope, shape))

    return SmallRiverResult(
        basin_name=data.basin_name,
        shape_factor=_safe_round(shape, 4),
        duration=rainfall_duration(area),
        kp_values=kp_values,
        design_rainfall=design_rainfall,
        flows=flows,
    )


@dataclass(frozen=True)
class StationFlow:
    name: str
    area: float
    flows: dict[str, float]


def double_log_interpolate_flows(
    station_a: StationFlow,
    station_b: StationFlow,
    design_area: float,
    frequencies: Iterable[str] = TIE_FREQUENCIES,
) -> dict[str, dict[str, float | str | None]]:
    """Interpolate design flows on log(F)-log(Q) coordinates between two stations."""

    area_a = _require_positive(f"{station_a.name} 面积", station_a.area)
    area_b = _require_positive(f"{station_b.name} 面积", station_b.area)
    design_area = _require_positive("设计流域面积 F设", design_area)
    if abs(math.log10(area_b) - math.log10(area_a)) < 1e-12:
        raise HydrologyError("双对数内插法要求两个水文站面积不能相同")

    if area_b < area_a:
        station_a, station_b = station_b, station_a
        area_a, area_b = area_b, area_a

    log_area_a = math.log10(area_a)
    log_area_b = math.log10(area_b)
    log_design_area = math.log10(design_area)
    area_position = "内插" if area_a <= design_area <= area_b else "外推"
    area_weight = (log_design_area - log_area_a) / (log_area_b - log_area_a)

    result: dict[str, dict[str, float | str | None]] = {}
    for freq in frequencies:
        q_a = float(station_a.flows.get(freq, 0) or 0)
        q_b = float(station_b.flows.get(freq, 0) or 0)
        if q_a <= 0 or q_b <= 0:
            result[freq] = {
                "station_a": station_a.name,
                "area_a": area_a,
                "q_a": None if q_a <= 0 else _safe_round(q_a, 3),
                "station_b": station_b.name,
                "area_b": area_b,
                "q_b": None if q_b <= 0 else _safe_round(q_b, 3),
                "n": None,
                "weight": _safe_round(area_weight, 3),
                "position": "缺少流量",
                "q_design": None,
            }
            continue
        log_q_a = math.log10(q_a)
        log_q_b = math.log10(q_b)
        n_value = (log_q_b - log_q_a) / (log_area_b - log_area_a)
        log_q_design = log_q_a + (log_q_b - log_q_a) * area_weight
        result[freq] = {
            "station_a": station_a.name,
            "area_a": area_a,
            "q_a": _safe_round(q_a, 3),
            "station_b": station_b.name,
            "area_b": area_b,
            "q_b": _safe_round(q_b, 3),
            "n": _safe_round(n_value, 3),
            "weight": _safe_round(area_weight, 3),
            "position": area_position,
            "q_design": _safe_round(10 ** log_q_design, 3),
        }
    return result


@dataclass(frozen=True)
class FloodSynthesisComponent:
    name: str
    time_offset_h: float
    flows: list[float]


@dataclass(frozen=True)
class FloodSynthesisPoint:
    time_h: float
    component_flows: dict[str, float]
    total_flow: float


@dataclass(frozen=True)
class FloodSynthesisResult:
    delta_t_h: float
    components: tuple[FloodSynthesisComponent, ...]
    points: tuple[FloodSynthesisPoint, ...]
    peak_time_h: float
    peak_flow_m3s: float
    flood_volume_10000m3: float


@dataclass(frozen=True)
class ManningSectionInput:
    name: str
    roughness_n: float
    slope: float
    design_flow_m3s: float
    points: tuple[tuple[float, float], ...]


@dataclass(frozen=True)
class ManningSectionResult:
    name: str
    water_level: float
    thalweg_x_m: float
    thalweg_elevation_m: float
    area_m2: float
    wetted_perimeter_m: float
    hydraulic_radius_m: float
    top_width_m: float
    max_depth_m: float
    velocity_ms: float
    froude_number: float
    calculated_flow_m3s: float


@dataclass(frozen=True)
class HydraulicSectionInput:
    name: str
    points: tuple[tuple[float, float], ...]
    roughness_n: float
    alpha: float = 1.0


@dataclass(frozen=True)
class HydraulicSectionState:
    name: str
    water_level_m: float
    thalweg_x_m: float
    thalweg_elevation_m: float
    area_m2: float
    wetted_perimeter_m: float
    hydraulic_radius_m: float
    top_width_m: float
    max_depth_m: float
    velocity_ms: float
    velocity_head_m: float
    froude_number: float
    energy_grade_m: float
    flow_modulus_k: float
    friction_slope: float


@dataclass(frozen=True)
class EnergyEquationInput:
    flow_m3s: float
    known_section: HydraulicSectionInput
    target_section: HydraulicSectionInput
    known_water_level_m: float
    reach_length_m: float
    known_is_upstream: bool
    contraction_coeff: float = 0.1
    expansion_coeff: float = 0.3
    flow_regime: str = "subcritical"
    friction_method: str = "average_conveyance"


@dataclass(frozen=True)
class EnergyEquationResult:
    flow_m3s: float
    known_is_upstream: bool
    known_state: HydraulicSectionState
    target_state: HydraulicSectionState
    reach_length_m: float
    friction_loss_m: float
    local_loss_m: float
    total_loss_m: float
    energy_residual_m: float


@dataclass(frozen=True)
class HecReachInput:
    flow_m3s: float
    sections_downstream_to_upstream: tuple[HydraulicSectionInput, ...]
    reach_lengths_m: tuple[float, ...]
    downstream_water_level_m: float
    upstream_water_level_m: float | None = None
    flow_regime: str = "subcritical"
    contraction_coeff: float = 0.1
    expansion_coeff: float = 0.3
    friction_method: str = "average_conveyance"


@dataclass(frozen=True)
class HecProfileStep:
    index: int
    downstream_section: str
    upstream_section: str
    reach_length_m: float
    downstream_water_level_m: float
    upstream_water_level_m: float
    friction_loss_m: float
    local_loss_m: float
    total_loss_m: float
    energy_residual_m: float | None = None


@dataclass(frozen=True)
class HecProfileResult:
    flow_m3s: float
    states_downstream_to_upstream: tuple[HydraulicSectionState, ...]
    steps: tuple[HecProfileStep, ...]
    flow_regime: str
    calculation_direction: str
    boundary_note: str


@dataclass(frozen=True)
class WeirFlowInput:
    weir_type: str
    design_flow_m3s: float
    crest_elevation_m: float
    downstream_water_level_m: float = 0.0
    opening_count: int = 1
    single_width_m: float = 1.0
    sill_height_m: float = 0.0
    abutment_type: int = 1
    pier_type: int = 1
    sill_type: int = 1
    practical_m: float = 0.48
    upstream_weir_height_m: float = 0.5
    orifice_height_m: float = 1.0


@dataclass(frozen=True)
class WeirFlowResult:
    weir_type: str
    design_flow_m3s: float
    upstream_water_level_m: float
    crest_elevation_m: float
    downstream_water_level_m: float
    head_m: float
    downstream_head_m: float
    submergence_ratio: float
    submergence_coeff: float
    discharge_coeff: float
    abutment_coeff: float
    pier_coeff: float
    contraction_coeff: float
    calculated_flow_m3s: float
    classification_note: str
    formula_note: str


@dataclass(frozen=True)
class WeirSectionRatingInput:
    section_name: str
    points: tuple[tuple[float, float], ...]
    thalweg_elevation_m: float
    max_depth_m: float
    weir_type: str = "宽顶堰"
    crest_elevation_m: float = 0.0
    downstream_water_level_m: float = 0.0
    depth_step_m: float = 0.2
    submergence_coeff: float = 1.0
    contraction_coeff: float = 1.0
    discharge_coeff: float = 0.36
    opening_count: int = 1
    gravity: float = 9.81
    orifice_height_m: float = 1.0


@dataclass(frozen=True)
class WeirSectionRatingRow:
    flow_m3s: float
    submergence_coeff: float
    contraction_coeff: float
    discharge_coeff: float
    gravity: float
    depth_m: float
    head_m: float
    downstream_head_m: float
    submergence_ratio: float
    overflow_state: str
    thalweg_elevation_m: float
    water_level_m: float
    top_width_m: float


@dataclass(frozen=True)
class WeirSectionRatingResult:
    section_name: str
    rows: tuple[WeirSectionRatingRow, ...]
    formula_note: str


_BROAD_WEIR_SUBMERGENCE_X = (
    0.80, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89,
    0.90, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98,
)
_BROAD_WEIR_SUBMERGENCE_Y = (
    1.000, 0.995, 0.990, 0.980, 0.970, 0.960, 0.950, 0.930, 0.900, 0.870,
    0.840, 0.820, 0.780, 0.740, 0.700, 0.650, 0.590, 0.500, 0.400,
)
_PRACTICAL_WEIR_SUBMERGENCE_X = (
    0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50,
    0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.00,
)
_PRACTICAL_WEIR_SUBMERGENCE_Y = (
    1.000, 0.996, 0.991, 0.986, 0.981, 0.976, 0.970, 0.963, 0.956, 0.948,
    0.937, 0.923, 0.907, 0.886, 0.860, 0.821, 0.778, 0.709, 0.621, 0.438, 0.000,
)
_PIER_COEFF_X = (0.75, 0.80, 0.85, 0.90, 0.95)
_PIER_COEFF_ROWS = {
    1: (0.80, 0.86, 0.92, 0.98, 1.00),
    2: (0.45, 0.51, 0.57, 0.63, 0.69),
    3: (0.45, 0.51, 0.57, 0.63, 0.69),
    4: (0.25, 0.32, 0.39, 0.46, 0.53),
    5: (0.00, 0.00, 0.00, 0.00, 0.00),
}


def _is_open_practical_weir(weir_type: str) -> bool:
    return weir_type in {"实用堰", "实用堰-敞式"}


def _is_orifice_practical_weir(weir_type: str) -> bool:
    return weir_type in {"实用堰-胸墙孔口式", "实用堰-胸墙孔口式平底闸堰"}


def _linear_interpolate_clamped(x: float, xs: tuple[float, ...], ys: tuple[float, ...]) -> float:
    if x <= xs[0]:
        return ys[0]
    if x >= xs[-1]:
        return ys[-1]
    for idx in range(len(xs) - 1):
        if xs[idx] <= x <= xs[idx + 1]:
            fraction = (x - xs[idx]) / (xs[idx + 1] - xs[idx])
            return ys[idx] + (ys[idx + 1] - ys[idx]) * fraction
    return ys[-1]


def _weir_abutment_coeff(abutment_type: int) -> float:
    if abutment_type == 1:
        return 1.0
    if abutment_type in {2, 3}:
        return 0.7
    return 0.4


def _weir_pier_coeff(pier_type: int, submergence_ratio: float) -> float:
    values = _PIER_COEFF_ROWS.get(int(pier_type), _PIER_COEFF_ROWS[1])
    return _linear_interpolate_clamped(submergence_ratio, _PIER_COEFF_X, values)


def _weir_common_contraction(
    head_m: float,
    opening_count: int,
    single_width_m: float,
    abutment_coeff: float,
    pier_coeff: float,
) -> float:
    if opening_count <= 0:
        raise HydrologyError("孔数必须大于0")
    width = _require_positive("单孔宽B", single_width_m)
    value = 1 - 0.2 * (abutment_coeff + (opening_count - 1) * pier_coeff) * head_m / (opening_count * width)
    return max(0.0, value)


def _broad_weir_discharge_coeff(head_m: float, sill_height_m: float, sill_type: int) -> float:
    ratio = sill_height_m / head_m if head_m > 0 else 0.0
    if int(sill_type) == 1:
        return 0.33 + 0.01 * (3 - ratio) / (0.46 + 0.75 * ratio)
    return 0.36 + 0.01 * (3 - ratio) / (1.2 + 1.5 * ratio)


def _thin_wall_discharge_coeff(head_m: float, upstream_weir_height_m: float) -> float:
    height = _require_positive("上游堰高a", upstream_weir_height_m)
    # Follows the provided Excel formula: (0.405+0.0027/H)*(1+0.55*(H/(H+a))^4).
    return (0.405 + 0.0027 / head_m) * (1 + 0.55 * (head_m / (head_m + height)) ** 4)


def _weir_discharge_for_head(data: WeirFlowInput, head_m: float) -> tuple[float, float, float, float, float, float, float, float]:
    if head_m <= 0:
        return 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0
    weir_type = data.weir_type.strip()
    width = _require_positive("单孔宽/溢流宽", data.single_width_m)
    opening_count = int(data.opening_count)
    if opening_count <= 0:
        raise HydrologyError("孔数必须大于0")
    downstream_head = max(0.0, data.downstream_water_level_m - data.crest_elevation_m)
    submergence_ratio = downstream_head / head_m if head_m > 0 else 0.0
    abutment = 1.0
    pier = 0.0
    contraction = 1.0
    if weir_type == "宽顶堰":
        submergence = 1.0 if submergence_ratio <= 0.80 else _linear_interpolate_clamped(
            submergence_ratio, _BROAD_WEIR_SUBMERGENCE_X, _BROAD_WEIR_SUBMERGENCE_Y
        )
        discharge_coeff = _broad_weir_discharge_coeff(head_m, data.sill_height_m, data.sill_type)
        abutment = _weir_abutment_coeff(data.abutment_type)
        pier = _weir_pier_coeff(data.pier_type, submergence_ratio)
        contraction = _weir_common_contraction(head_m, opening_count, width, abutment, pier)
    elif _is_open_practical_weir(weir_type):
        submergence = 1.0
        discharge_coeff = _require_positive("实用堰流量系数m", data.practical_m)
        abutment = _weir_abutment_coeff(data.abutment_type)
        pier = _weir_pier_coeff(data.pier_type, submergence_ratio)
        contraction = _weir_common_contraction(head_m, opening_count, width, abutment, pier)
    elif _is_orifice_practical_weir(weir_type):
        discharge_coeff = _require_positive("孔口流量系数μ", data.practical_m)
        orifice_height = _require_positive("孔口高度D", data.orifice_height_m)
        effective_head = max(0.0, head_m - orifice_height)
        area = opening_count * width * orifice_height
        flow = discharge_coeff * area * math.sqrt(2 * 9.81 * effective_head) if effective_head > 0 else 0.0
        return flow, downstream_head, submergence_ratio, 1.0, discharge_coeff, abutment, pier, contraction
    elif weir_type == "薄壁堰":
        submergence = 1.0
        discharge_coeff = _thin_wall_discharge_coeff(head_m, data.upstream_weir_height_m)
    else:
        raise HydrologyError(f"不支持的堰型：{data.weir_type}")
    flow = submergence * discharge_coeff * contraction * opening_count * width * math.sqrt(2 * 9.81) * head_m ** 1.5
    return flow, downstream_head, submergence_ratio, submergence, discharge_coeff, abutment, pier, contraction


def _weir_submergence_for_ratio(weir_type: str, ratio: float) -> tuple[float, str]:
    if weir_type == "宽顶堰":
        if ratio <= 0.80:
            return 1.0, "自由溢流"
        return _linear_interpolate_clamped(ratio, _BROAD_WEIR_SUBMERGENCE_X, _BROAD_WEIR_SUBMERGENCE_Y), "淹没溢流"
    if _is_open_practical_weir(weir_type):
        return 1.0, "实用堰出流"
    if _is_orifice_practical_weir(weir_type):
        return 1.0, "孔口出流"
    if weir_type == "薄壁堰":
        return 1.0, "薄壁堰出流"
    return 1.0, "堰流计算"


def calculate_weir_water_level(data: WeirFlowInput) -> WeirFlowResult:
    target_flow = _require_positive("设计流量Q", data.design_flow_m3s)
    _require_positive("单孔宽/溢流宽", data.single_width_m)
    low = 1e-6
    high = max(0.5, (target_flow / max(data.opening_count * data.single_width_m, 1e-6)) ** (2 / 3))
    high_flow = _weir_discharge_for_head(data, high)[0]
    expand_count = 0
    while high_flow < target_flow and expand_count < 80:
        high *= 1.8
        high_flow = _weir_discharge_for_head(data, high)[0]
        expand_count += 1
    if high_flow < target_flow:
        raise HydrologyError("堰流公式无法在合理水头范围内反算水位，请检查流量、堰宽、孔数和系数")
    for _ in range(90):
        mid = (low + high) / 2
        mid_flow = _weir_discharge_for_head(data, mid)[0]
        if mid_flow >= target_flow:
            high = mid
        else:
            low = mid
    head = (low + high) / 2
    flow, downstream_head, sub_ratio, sub_coeff, discharge_coeff, abutment, pier, contraction = _weir_discharge_for_head(data, head)
    if data.weir_type == "宽顶堰":
        classification = "宽顶堰适用条件参考：2.5<δ/H<10。"
        formula = "Q=σ·m·ε·n·B·sqrt(2g)·H^(3/2)"
    elif _is_open_practical_weir(data.weir_type):
        classification = "敞式实用堰：按实用堰公式计算，不区分自由溢流和淹没溢流；流量系数m需按项目取值。"
        formula = "Q=m·ε·n·B·sqrt(2g)·H^(3/2)"
    elif data.weir_type == "实用堰-胸墙孔口式":
        classification = "胸墙孔口式实用堰：按孔口出流计算，μ一般约0.6~0.8，圆滑孔口可取约0.9。"
        formula = "Q=μ·A·sqrt(2g(H-D)), A=n·B·D"
    elif data.weir_type == "实用堰-胸墙孔口式平底闸堰":
        classification = "胸墙孔口式平底闸堰：按孔口出流计算，μ一般约0.6~0.8。"
        formula = "Q=μ·A·sqrt(2g(H-D)), A=n·B·D"
    else:
        classification = "薄壁堰适用条件参考：δ/H<0.67。"
        formula = "Q=m₀·b·sqrt(2g)·H^(3/2)"
    return WeirFlowResult(
        weir_type=data.weir_type,
        design_flow_m3s=_safe_round(target_flow, 3),
        upstream_water_level_m=_safe_round(data.crest_elevation_m + head, 3),
        crest_elevation_m=_safe_round(data.crest_elevation_m, 3),
        downstream_water_level_m=_safe_round(data.downstream_water_level_m, 3),
        head_m=_safe_round(head, 3),
        downstream_head_m=_safe_round(downstream_head, 3),
        submergence_ratio=_safe_round(sub_ratio, 3),
        submergence_coeff=_safe_round(sub_coeff, 3),
        discharge_coeff=_safe_round(discharge_coeff, 6),
        abutment_coeff=_safe_round(abutment, 3),
        pier_coeff=_safe_round(pier, 3),
        contraction_coeff=_safe_round(contraction, 3),
        calculated_flow_m3s=_safe_round(flow, 3),
        classification_note=classification,
        formula_note=formula,
    )


def _section_geometry(points: tuple[tuple[float, float], ...], water_level: float) -> tuple[float, float, float, float]:
    area = 0.0
    wetted = 0.0
    top_width = 0.0
    for (x1, z1), (x2, z2) in zip(points, points[1:]):
        dx = x2 - x1
        if abs(dx) < 1e-12:
            continue
        d1 = water_level - z1
        d2 = water_level - z2
        length = math.hypot(x2 - x1, z2 - z1)
        if d1 <= 0 and d2 <= 0:
            continue
        if d1 >= 0 and d2 >= 0:
            area += (d1 + d2) * abs(dx) / 2
            wetted += length
            top_width += abs(dx)
            continue
        ratio = d1 / (d1 - d2)
        xi = x1 + ratio * (x2 - x1)
        zi = water_level
        if d1 > 0:
            area += d1 * abs(xi - x1) / 2
            wetted += math.hypot(xi - x1, zi - z1)
            top_width += abs(xi - x1)
        else:
            area += d2 * abs(x2 - xi) / 2
            wetted += math.hypot(x2 - xi, z2 - zi)
            top_width += abs(x2 - xi)
    max_depth = max(0.0, max((water_level - z for _x, z in points), default=0.0))
    return area, wetted, top_width, max_depth


def calculate_weir_section_rating(data: WeirSectionRatingInput) -> WeirSectionRatingResult:
    section = tuple(sorted((float(x), float(z)) for x, z in data.points))
    if len(section) < 2:
        raise HydrologyError("堰流断面点至少需要2个")
    if len({x for x, _z in section}) != len(section):
        raise HydrologyError("堰流断面点里程不能重复")
    thalweg = float(data.thalweg_elevation_m)
    max_depth = _require_positive("最大计算水深", data.max_depth_m)
    step = _require_positive("水深步长", data.depth_step_m)
    epsilon = _require_positive("侧收缩系数ε", data.contraction_coeff)
    discharge_coeff = _require_positive("流量系数m", data.discharge_coeff)
    gravity = _require_positive("重力加速度g", data.gravity)
    orifice_height = _require_positive("孔口高度D", data.orifice_height_m)
    crest = float(data.crest_elevation_m) if data.crest_elevation_m else thalweg
    downstream_head = max(0.0, float(data.downstream_water_level_m) - crest)
    opening_count = int(data.opening_count)
    if opening_count <= 0:
        raise HydrologyError("孔数n必须大于0")

    depths: list[float] = []
    depth = 0.0
    guard = 0
    while depth < max_depth - 1e-9 and guard < 10000:
        depths.append(round(depth, 10))
        depth += step
        guard += 1
    depths.append(max_depth)

    rows: list[WeirSectionRatingRow] = []
    for depth in depths:
        water_level = thalweg + depth
        _area, _wetted, top_width, _section_max_depth = _section_geometry(section, water_level)
        head = max(0.0, water_level - crest)
        if head <= 0:
            ratio = 0.0
            sigma = 1.0
            state = "未过堰"
        elif _is_orifice_practical_weir(data.weir_type.strip()):
            ratio = downstream_head / head
            sigma = 1.0
            state = "孔口出流" if head > orifice_height else "未形成孔口出流"
        else:
            ratio = downstream_head / head
            sigma, state = _weir_submergence_for_ratio(data.weir_type.strip(), ratio)
        if head <= 0 or top_width <= 0:
            flow = 0.0
        elif _is_orifice_practical_weir(data.weir_type.strip()):
            effective_head = max(0.0, head - orifice_height)
            area = opening_count * top_width * orifice_height
            flow = discharge_coeff * area * math.sqrt(2 * gravity * effective_head) if effective_head > 0 else 0.0
        else:
            flow = sigma * epsilon * discharge_coeff * opening_count * top_width * math.sqrt(2 * gravity) * head ** 1.5
        rows.append(
            WeirSectionRatingRow(
                flow_m3s=_safe_round(flow, 2),
                submergence_coeff=_safe_round(sigma, 3),
                contraction_coeff=_safe_round(epsilon, 3),
                discharge_coeff=_safe_round(discharge_coeff, 3),
                gravity=_safe_round(gravity, 2),
                depth_m=_safe_round(depth, 2),
                head_m=_safe_round(head, 3),
                downstream_head_m=_safe_round(downstream_head, 3),
                submergence_ratio=_safe_round(ratio, 3),
                overflow_state=state,
                thalweg_elevation_m=_safe_round(thalweg, 3),
                water_level_m=_safe_round(water_level, 3),
                top_width_m=_safe_round(top_width, 2),
            )
        )

    return WeirSectionRatingResult(
        section_name=data.section_name.strip() or "堰流计算断面",
        rows=tuple(rows),
        formula_note="宽顶堰：Q=δ·ε·m·n·B·sqrt(2g)·H^(3/2)；实用堰、薄壁堰按对应公式计算，不按自由/淹没溢流分型。",
    )


def manning_flow_for_level(
    points: Iterable[tuple[float, float]],
    water_level: float,
    roughness_n: float,
    slope: float,
) -> tuple[float, float, float, float, float, float]:
    n_value = _require_positive("糙率 n", roughness_n)
    slope_value = _require_positive("水力坡降 I", slope)
    section = tuple(sorted((float(x), float(z)) for x, z in points))
    if len(section) < 2:
        raise HydrologyError("断面点至少需要 2 个")
    area, wetted, top_width, max_depth = _section_geometry(section, water_level)
    if area <= 0 or wetted <= 0:
        return 0.0, area, wetted, 0.0, top_width, max_depth
    radius = area / wetted
    flow = (1 / n_value) * area * (radius ** (2 / 3)) * math.sqrt(slope_value)
    return flow, area, wetted, radius, top_width, max_depth


def hydraulic_state_for_level(section_input: HydraulicSectionInput, flow_m3s: float, water_level_m: float) -> HydraulicSectionState:
    name = section_input.name.strip() or "未命名断面"
    flow = _require_positive("流量Q", flow_m3s)
    n_value = _require_positive(f"{name} 糙率n", section_input.roughness_n)
    alpha = _require_positive(f"{name} 动能修正系数α", section_input.alpha)
    points = tuple(sorted((float(x), float(z)) for x, z in section_input.points))
    if len(points) < 2:
        raise HydrologyError(f"{name} 至少需要2个断面点")
    if len({x for x, _z in points}) != len(points):
        raise HydrologyError(f"{name} 断面点起点距不能重复")
    area, wetted, top_width, max_depth = _section_geometry(points, water_level_m)
    if area <= 0 or wetted <= 0:
        raise HydrologyError(f"{name} 水位未淹没有效断面，无法计算水力要素")
    radius = area / wetted
    velocity = flow / area
    velocity_head = alpha * velocity * velocity / (2 * 9.80665)
    hydraulic_depth = area / top_width if top_width > 0 else 0.0
    froude = velocity / math.sqrt(9.80665 * hydraulic_depth) if hydraulic_depth > 0 else 0.0
    flow_modulus = area * (radius ** (2 / 3)) / n_value
    friction_slope = (flow / flow_modulus) ** 2
    thalweg_x, thalweg_z = min(points, key=lambda item: item[1])
    return HydraulicSectionState(
        name=name,
        water_level_m=water_level_m,
        thalweg_x_m=thalweg_x,
        thalweg_elevation_m=thalweg_z,
        area_m2=area,
        wetted_perimeter_m=wetted,
        hydraulic_radius_m=radius,
        top_width_m=top_width,
        max_depth_m=max_depth,
        velocity_ms=velocity,
        velocity_head_m=velocity_head,
        froude_number=froude,
        energy_grade_m=water_level_m + velocity_head,
        flow_modulus_k=flow_modulus,
        friction_slope=friction_slope,
    )


def _critical_water_level_raw(points: Iterable[tuple[float, float]], flow_m3s: float) -> tuple[float, float, float, float]:
    flow = _require_positive("流量Q", flow_m3s)
    section = tuple(sorted((float(x), float(z)) for x, z in points))
    if len(section) < 2:
        raise HydrologyError("临界水深计算至少需要2个断面点")
    min_z = min(z for _x, z in section)
    max_z = max(z for _x, z in section)

    def critical_residual(level: float) -> float:
        area, _wetted, top_width, _max_depth = _section_geometry(section, level)
        if area <= 0 or top_width <= 0:
            return float("inf")
        return flow * flow * top_width / (9.80665 * area ** 3) - 1.0

    low = min_z + 1e-6
    high = max_z + max(1.0, (max(x for x, _z in section) - min(x for x, _z in section)) * 0.05)
    low_residual = critical_residual(low)
    high_residual = critical_residual(high)
    expand_count = 0
    while low_residual * high_residual > 0 and expand_count < 80:
        high += max(1.0, high - low)
        high_residual = critical_residual(high)
        expand_count += 1
    if low_residual * high_residual > 0:
        raise HydrologyError("无法在合理范围内求得临界水深")

    for _ in range(90):
        mid = (low + high) / 2
        mid_residual = critical_residual(mid)
        if low_residual * mid_residual <= 0:
            high = mid
            high_residual = mid_residual
        else:
            low = mid
            low_residual = mid_residual

    level = (low + high) / 2
    area, _wetted, top_width, max_depth = _section_geometry(section, level)
    return level, max_depth, area, top_width


def critical_water_level_for_section(points: Iterable[tuple[float, float]], flow_m3s: float) -> tuple[float, float, float, float]:
    level, max_depth, area, top_width = _critical_water_level_raw(points, flow_m3s)
    return _safe_round(level, 3), _safe_round(max_depth, 3), _safe_round(area, 3), _safe_round(top_width, 3)


def normal_water_level_for_section(
    section_input: HydraulicSectionInput,
    flow_m3s: float,
    bed_slope: float,
) -> tuple[float, float, float]:
    slope = _require_positive("河床坡降", bed_slope)
    result = calculate_manning_water_level(
        ManningSectionInput(
            name=section_input.name,
            roughness_n=section_input.roughness_n,
            slope=slope,
            design_flow_m3s=flow_m3s,
            points=section_input.points,
        )
    )
    return result.water_level, result.max_depth_m, result.froude_number


def _energy_loss_between(
    known_state: HydraulicSectionState,
    target_state: HydraulicSectionState,
    reach_length_m: float,
    known_is_upstream: bool,
    contraction_coeff: float,
    expansion_coeff: float,
    friction_method: str = "average_conveyance",
) -> tuple[float, float, float]:
    length = _require_positive("河段长度L", reach_length_m)
    contraction = _require_non_negative("收缩损失系数", contraction_coeff)
    expansion = _require_non_negative("扩张损失系数", expansion_coeff)
    method = friction_method.strip().lower()
    if method in {"average_conveyance", "hec", "平均流量模数"}:
        known_q = known_state.flow_modulus_k * math.sqrt(known_state.friction_slope)
        target_q = target_state.flow_modulus_k * math.sqrt(target_state.friction_slope)
        representative_slope = ((known_q + target_q) / (known_state.flow_modulus_k + target_state.flow_modulus_k)) ** 2
    elif method in {"average_friction_slope", "arithmetic", "平均摩阻坡降"}:
        representative_slope = (known_state.friction_slope + target_state.friction_slope) / 2
    else:
        raise HydrologyError(f"不支持的摩阻坡降计算方法：{friction_method}")
    friction = representative_slope * length
    if known_is_upstream:
        upstream_vh = known_state.velocity_head_m
        downstream_vh = target_state.velocity_head_m
    else:
        upstream_vh = target_state.velocity_head_m
        downstream_vh = known_state.velocity_head_m
    coeff = contraction if downstream_vh > upstream_vh else expansion
    local = coeff * abs(downstream_vh - upstream_vh)
    return friction, local, friction + local


def _normalize_flow_regime(value: str) -> str:
    regime = value.strip().lower()
    if regime in {"亚临界流", "subcritical", "sub"}:
        return "subcritical"
    if regime in {"超临界流", "supercritical", "super"}:
        return "supercritical"
    if regime in {"混合流", "mixed", "mix"}:
        raise HydrologyError("当前标准步法不支持混合流自动计算。请先划分亚临界、超临界河段；水跃和过临界控制断面应采用专门模型复核。")
    raise HydrologyError(f"不支持的计算流态：{value}")


def _find_energy_root(residual, low: float, high: float, expand_high: bool) -> tuple[float, float]:
    for _ in range(30):
        previous_level = low
        previous_value = residual(previous_level)
        for index in range(1, 321):
            level = low + (high - low) * index / 320
            value = residual(level)
            if abs(value) <= 1e-10:
                return level, level
            if previous_value * value < 0:
                return previous_level, level
            previous_level = level
            previous_value = value
        if not expand_high:
            break
        high += max(1.0, high - low)
    raise HydrologyError("所选流态分支内没有可收敛的水位解")


def calculate_energy_equation_water_level(data: EnergyEquationInput) -> EnergyEquationResult:
    flow = _require_positive("流量Q", data.flow_m3s)
    regime = _normalize_flow_regime(data.flow_regime)
    known_state = hydraulic_state_for_level(data.known_section, flow, data.known_water_level_m)
    if regime == "subcritical" and known_state.froude_number >= 1.0:
        raise HydrologyError(
            f"已知断面Fr={known_state.froude_number:.3f}，不是亚临界流。"
            "请核对控制水位和流量，或改用超临界流（上游控制）。"
        )
    if regime == "supercritical" and known_state.froude_number <= 1.0:
        raise HydrologyError(
            f"已知断面Fr={known_state.froude_number:.3f}，不是超临界流。"
            "请核对控制水位和流量，或改用亚临界流（下游控制）。"
        )

    target_points = tuple(sorted((float(x), float(z)) for x, z in data.target_section.points))
    min_z = min(z for _x, z in target_points)
    max_z = max(z for _x, z in target_points)
    low = min_z + 1e-6
    high = max(max_z + 1.0, data.known_water_level_m + 5.0)

    def residual(level: float) -> float:
        target_state = hydraulic_state_for_level(data.target_section, flow, level)
        friction, local, total = _energy_loss_between(
            known_state,
            target_state,
            data.reach_length_m,
            data.known_is_upstream,
            data.contraction_coeff,
            data.expansion_coeff,
            data.friction_method,
        )
        if data.known_is_upstream:
            return target_state.energy_grade_m + total - known_state.energy_grade_m
        return target_state.energy_grade_m - known_state.energy_grade_m - total

    critical_level, _critical_depth, _critical_area, _critical_width = _critical_water_level_raw(target_points, flow)
    branch_gap = max(1e-6, (critical_level - min_z) * 1e-7)
    if regime == "subcritical":
        low = critical_level + branch_gap
        high = max(high, low + 1.0)
        expand_high = True
    else:
        low = min_z + branch_gap
        high = critical_level - branch_gap
        expand_high = False
    try:
        low, high = _find_energy_root(residual, low, high, expand_high)
    except HydrologyError as root_error:
        direction_text = "上游已知求下游" if data.known_is_upstream else "下游已知求上游"
        raise HydrologyError(
            f"能量方程在{('亚临界' if regime == 'subcritical' else '超临界')}分支内无法收敛。"
            f"计算方向：{direction_text}；"
            f"已知断面：{data.known_section.name}，已知水位：{data.known_water_level_m:.3f}m；"
            f"目标断面：{data.target_section.name}；河段长度：{data.reach_length_m:.3f}m。"
            "请优先检查：1）已知水位是否填在所选控制断面上；2）亚临界流应下游控制、超临界流应上游控制；"
            "3）上下游断面顺序和高程是否录反；4）目标断面点是否导入到了正确断面。"
        ) from root_error

    if low != high:
        low_residual = residual(low)
        for _ in range(90):
            mid = (low + high) / 2
            mid_residual = residual(mid)
            if low_residual * mid_residual <= 0:
                high = mid
            else:
                low = mid
                low_residual = mid_residual

    target_level = (low + high) / 2
    target_state = hydraulic_state_for_level(data.target_section, flow, target_level)
    friction, local, total = _energy_loss_between(
        known_state,
        target_state,
        data.reach_length_m,
        data.known_is_upstream,
        data.contraction_coeff,
        data.expansion_coeff,
        data.friction_method,
    )
    final_residual = (
        target_state.energy_grade_m + total - known_state.energy_grade_m
        if data.known_is_upstream
        else target_state.energy_grade_m - known_state.energy_grade_m - total
    )
    return EnergyEquationResult(
        flow_m3s=_safe_round(flow, 3),
        known_is_upstream=data.known_is_upstream,
        known_state=known_state,
        target_state=target_state,
        reach_length_m=_safe_round(data.reach_length_m, 3),
        friction_loss_m=_safe_round(friction, 4),
        local_loss_m=_safe_round(local, 4),
        total_loss_m=_safe_round(total, 4),
        energy_residual_m=_safe_round(final_residual, 5),
    )


def calculate_hec_profile(data: HecReachInput) -> HecProfileResult:
    flow = _require_positive("流量Q", data.flow_m3s)
    sections = data.sections_downstream_to_upstream
    if len(sections) < 2:
        raise HydrologyError("HEC标准步法至少需要2个断面")
    if len(data.reach_lengths_m) != len(sections) - 1:
        raise HydrologyError("河段长度数量应为断面数减1")
    regime = _normalize_flow_regime(data.flow_regime)

    steps: list[HecProfileStep] = []
    if regime == "subcritical":
        states: list[HydraulicSectionState] = [
            hydraulic_state_for_level(sections[0], flow, data.downstream_water_level_m)
        ]
        current_level = data.downstream_water_level_m
        for idx in range(1, len(sections)):
            result = calculate_energy_equation_water_level(
                EnergyEquationInput(
                    flow_m3s=flow,
                    known_section=sections[idx - 1],
                    target_section=sections[idx],
                    known_water_level_m=current_level,
                    reach_length_m=data.reach_lengths_m[idx - 1],
                    known_is_upstream=False,
                    contraction_coeff=data.contraction_coeff,
                    expansion_coeff=data.expansion_coeff,
                    flow_regime="subcritical",
                    friction_method=data.friction_method,
                )
            )
            states.append(result.target_state)
            steps.append(
                HecProfileStep(
                    index=idx,
                    downstream_section=result.known_state.name,
                    upstream_section=result.target_state.name,
                    reach_length_m=result.reach_length_m,
                    downstream_water_level_m=result.known_state.water_level_m,
                    upstream_water_level_m=result.target_state.water_level_m,
                    friction_loss_m=result.friction_loss_m,
                    local_loss_m=result.local_loss_m,
                    total_loss_m=result.total_loss_m,
                    energy_residual_m=result.energy_residual_m,
                )
            )
            current_level = result.target_state.water_level_m
        direction = "下游控制，自下而上推算"
        note = "亚临界流采用下游边界条件，自下而上推算；各断面计算解强制位于临界水深以上。"
        regime_label = "亚临界流"
    else:
        if data.upstream_water_level_m is None:
            raise HydrologyError("超临界流计算需要填写上游控制水位")
        states_by_index: list[HydraulicSectionState | None] = [None] * len(sections)
        states_by_index[-1] = hydraulic_state_for_level(sections[-1], flow, data.upstream_water_level_m)
        current_level = data.upstream_water_level_m
        for idx in range(len(sections) - 2, -1, -1):
            result = calculate_energy_equation_water_level(
                EnergyEquationInput(
                    flow_m3s=flow,
                    known_section=sections[idx + 1],
                    target_section=sections[idx],
                    known_water_level_m=current_level,
                    reach_length_m=data.reach_lengths_m[idx],
                    known_is_upstream=True,
                    contraction_coeff=data.contraction_coeff,
                    expansion_coeff=data.expansion_coeff,
                    flow_regime="supercritical",
                    friction_method=data.friction_method,
                )
            )
            states_by_index[idx] = result.target_state
            steps.append(
                HecProfileStep(
                    index=len(steps) + 1,
                    downstream_section=result.target_state.name,
                    upstream_section=result.known_state.name,
                    reach_length_m=result.reach_length_m,
                    downstream_water_level_m=result.target_state.water_level_m,
                    upstream_water_level_m=result.known_state.water_level_m,
                    friction_loss_m=result.friction_loss_m,
                    local_loss_m=result.local_loss_m,
                    total_loss_m=result.total_loss_m,
                    energy_residual_m=result.energy_residual_m,
                )
            )
            current_level = result.target_state.water_level_m
        states = [state for state in states_by_index if state is not None]
        direction = "上游控制，自上而下推算"
        note = "超临界流采用上游边界条件，自上而下推算；各断面计算解强制位于临界水深以下。"
        regime_label = "超临界流"

    near_critical_names = [state.name for state in states if 0.94 <= state.froude_number <= 1.06]
    if near_critical_names:
        note += " 近临界断面：" + "、".join(near_critical_names) + "；应复核临界控制和计算步长。"

    return HecProfileResult(
        flow_m3s=_safe_round(flow, 3),
        states_downstream_to_upstream=tuple(states),
        steps=tuple(steps),
        flow_regime=regime_label,
        calculation_direction=direction,
        boundary_note=note,
    )


def calculate_manning_water_level(section_input: ManningSectionInput) -> ManningSectionResult:
    name = section_input.name.strip() or "未命名断面"
    n_value = _require_positive(f"{name} 糙率n", section_input.roughness_n)
    slope = _require_positive(f"{name} 水力坡降I", section_input.slope)
    target_q = _require_positive(f"{name} 设计流量Q", section_input.design_flow_m3s)
    points = tuple(sorted((float(x), float(z)) for x, z in section_input.points))
    if len(points) < 2:
        raise HydrologyError(f"{name} 至少需要 2 个断面点")
    if len({x for x, _z in points}) != len(points):
        raise HydrologyError(f"{name} 断面点起点距不能重复")

    min_z = min(z for _x, z in points)
    max_z = max(z for _x, z in points)
    low = min_z
    high = max_z + max(1.0, (max(x for x, _z in points) - min(x for x, _z in points)) * 0.05)
    high_flow = manning_flow_for_level(points, high, n_value, slope)[0]
    expand_count = 0
    while high_flow < target_q and expand_count < 80:
        high += max(1.0, high - low)
        high_flow = manning_flow_for_level(points, high, n_value, slope)[0]
        expand_count += 1
    if high_flow < target_q:
        raise HydrologyError(f"{name} 无法在合理范围内求得满足设计流量的水位")

    for _ in range(90):
        mid = (low + high) / 2
        mid_flow = manning_flow_for_level(points, mid, n_value, slope)[0]
        if mid_flow < target_q:
            low = mid
        else:
            high = mid

    water_level = (low + high) / 2
    flow, area, wetted, radius, top_width, max_depth = manning_flow_for_level(points, water_level, n_value, slope)
    velocity = flow / area if area > 0 else 0.0
    hydraulic_depth = area / top_width if top_width > 0 else 0.0
    froude = velocity / math.sqrt(9.80665 * hydraulic_depth) if hydraulic_depth > 0 else 0.0
    thalweg_x, thalweg_z = min(points, key=lambda item: item[1])
    return ManningSectionResult(
        name=name,
        water_level=_safe_round(water_level, 3),
        thalweg_x_m=_safe_round(thalweg_x, 3),
        thalweg_elevation_m=_safe_round(thalweg_z, 3),
        area_m2=_safe_round(area, 3),
        wetted_perimeter_m=_safe_round(wetted, 3),
        hydraulic_radius_m=_safe_round(radius, 3),
        top_width_m=_safe_round(top_width, 3),
        max_depth_m=_safe_round(max_depth, 3),
        velocity_ms=_safe_round(velocity, 3),
        froude_number=_safe_round(froude, 3),
        calculated_flow_m3s=_safe_round(flow, 3),
    )


def _component_flow_at(component: FloodSynthesisComponent, time_h: float, delta_t_h: float) -> float:
    if not component.flows:
        return 0.0
    local_time = time_h - component.time_offset_h
    if local_time < -1e-9:
        return 0.0
    last_time = (len(component.flows) - 1) * delta_t_h
    if local_time > last_time + 1e-9:
        return 0.0
    position = max(0.0, min(last_time, local_time)) / delta_t_h
    lower = int(math.floor(position))
    upper = min(lower + 1, len(component.flows) - 1)
    if lower == upper:
        return component.flows[lower]
    ratio = position - lower
    return component.flows[lower] + (component.flows[upper] - component.flows[lower]) * ratio


def synthesize_flood_hydrographs(
    components: Iterable[FloodSynthesisComponent],
    delta_t_h: float,
) -> FloodSynthesisResult:
    delta_t_h = _require_positive("时段 Δt", delta_t_h)
    component_list = tuple(
        FloodSynthesisComponent(item.name.strip() or "未命名过程线", float(item.time_offset_h), [float(v) for v in item.flows])
        for item in components
        if item.flows
    )
    if not component_list:
        raise HydrologyError("洪水合成法至少需要 1 条洪水过程线")
    if len({item.name for item in component_list}) != len(component_list):
        raise HydrologyError("各子流域过程线名称必须不同，避免同名流量被覆盖")
    for item in component_list:
        _require_non_negative(f"{item.name} 时间滞后", item.time_offset_h)
        if len(item.flows) < 2:
            raise HydrologyError(f"{item.name} 至少需要两个流量点")
        for flow in item.flows:
            _require_non_negative(f"{item.name} 流量", flow)

    end_time = max(item.time_offset_h + (len(item.flows) - 1) * delta_t_h for item in component_list)
    step_count = int(math.ceil(end_time / delta_t_h))
    times = [idx * delta_t_h for idx in range(step_count + 1)]
    if times[-1] < end_time - 1e-9:
        times.append(end_time)

    points: list[FloodSynthesisPoint] = []
    for time_h in times:
        component_flows = {
            item.name: _safe_round(_component_flow_at(item, time_h, delta_t_h), 3)
            for item in component_list
        }
        total_flow = _safe_round(sum(component_flows.values()), 3)
        points.append(FloodSynthesisPoint(_safe_round(time_h, 3), component_flows, total_flow))

    peak_point = max(points, key=lambda item: item.total_flow)
    volume_m3 = 0.0
    for previous, current in zip(points, points[1:]):
        dt_seconds = (current.time_h - previous.time_h) * 3600
        volume_m3 += (previous.total_flow + current.total_flow) / 2 * dt_seconds
    return FloodSynthesisResult(
        delta_t_h=_safe_round(delta_t_h, 3),
        components=component_list,
        points=tuple(points),
        peak_time_h=peak_point.time_h,
        peak_flow_m3s=peak_point.total_flow,
        flood_volume_10000m3=_safe_round(volume_m3 / 10000, 3),
    )


def calculate_area_index(stations: Iterable[StationFlow], frequencies: Iterable[str] = TIE_FREQUENCIES) -> dict[str, float | None]:
    stations = list(stations)
    if len(stations) < 2:
        raise HydrologyError("至少需要 2 个水文站计算面积指数 n")

    result: dict[str, float | None] = {}
    for freq in frequencies:
        xs: list[float] = []
        ys: list[float] = []
        for station in stations:
            area = float(station.area)
            q = float(station.flows.get(freq, 0) or 0)
            if area > 0 and q > 0:
                xs.append(math.log10(area))
                ys.append(math.log10(q))
        if len(xs) < 2:
            result[freq] = None
            continue
        sum_x = sum(xs)
        sum_y = sum(ys)
        sum_xy = sum(x * y for x, y in zip(xs, ys))
        sum_x2 = sum(x * x for x in xs)
        denom = sum_x2 - (sum_x * sum_x) / len(xs)
        if abs(denom) < 1e-12:
            result[freq] = None
            continue
        result[freq] = _safe_round((sum_xy - (sum_x * sum_y) / len(xs)) / denom, 3)
    return result


def calculate_area_analogy(
    reference_area: float,
    design_area: float,
    reference_flows: dict[str, float],
    n_values: dict[str, float] | None = None,
    unified_n: float = 0.667,
    frequencies: Iterable[str] = TIE_FREQUENCIES,
) -> dict[str, float | None]:
    reference_area = _require_positive("参证站面积 F参", reference_area)
    design_area = _require_positive("设计流域面积 F设", design_area)
    ratio = design_area / reference_area
    result: dict[str, float | None] = {}
    for freq in frequencies:
        q_ref = float(reference_flows.get(freq, 0) or 0)
        if q_ref <= 0:
            result[freq] = None
            continue
        n = float((n_values or {}).get(freq) or unified_n)
        if n <= 0:
            raise HydrologyError(f"{freq} 年面积指数 n 必须大于 0")
        result[freq] = _safe_round(q_ref * (ratio**n))
    return result


@dataclass(frozen=True)
class RegionParam:
    region_id: int
    name: str
    n: float
    k_values: tuple[float, float, float, float, float, float, float]

    @property
    def area_range(self) -> tuple[float | None, float | None]:
        interval = re.search(r"(\d+(?:\.\d+)?)\s*<\s*F\s*<\s*(\d+(?:\.\d+)?)", self.name)
        if interval:
            return float(interval.group(1)), float(interval.group(2))
        lower = re.search(r"F\s*>\s*(\d+(?:\.\d+)?)", self.name)
        if lower:
            return float(lower.group(1)), None
        upper = re.search(r"F\s*<\s*(\d+(?:\.\d+)?)", self.name)
        if upper:
            return None, float(upper.group(1))
        return None, None

    def validate_area(self, area: float) -> None:
        low, high = self.area_range
        if low is not None and not area > low:
            raise HydrologyError(f"面积 F={area:g} 不满足分区 {self.region_id} 的适用范围：{self.name}")
        if high is not None and not area < high:
            raise HydrologyError(f"面积 F={area:g} 不满足分区 {self.region_id} 的适用范围：{self.name}")


REGION_PARAMS = [
    RegionParam(1, '昌马河 F<11000', 0.864, (606.0, 512.0, 412.0, 343.0, 283.0, 210.0, 162.0)),
    RegionParam(2, '昌马河 F>11000', 3.16, (1090000000000.0, 935000000000.0, 735000000000.0, 639000000000.0, 536000000000.0, 394000000000.0, 308000000000.0)),
    RegionParam(3, '白杨河、赤金河、踏实河', 0.252, (2.3, 1.93, 1.51, 1.25, 0.968, 0.651, 0.445)),
    RegionParam(4, '黑河、梨园河 F<10000', 0.332, (6.49, 5.6, 4.81, 4.23, 3.52, 2.72, 2.08)),
    RegionParam(5, '黑河、梨园河 F>10000', 1.22, (22300.0, 19600.0, 16900.0, 14600.0, 12000.0, 9300.0, 7190.0)),
    RegionParam(6, '山丹马营河等', 0.572, (46.3, 39.5, 33.0, 27.3, 22.6, 15.7, 11.2)),
    RegionParam(7, '讨赖河 F<7000', 0.176, (1.3, 1.12, 0.903, 0.742, 0.598, 0.413, 0.293)),
    RegionParam(8, '讨赖河 F>7000', 1.45, (103000.0, 91400.0, 73100.0, 59400.0, 47500.0, 33400.0, 24000.0)),
    RegionParam(9, '杂木河、金塔河', 0.176, (3.56, 3.13, 2.66, 2.32, 1.89, 1.53, 1.15)),
    RegionParam(10, '黄羊河等 F<800', 0.495, (35.7, 30.2, 23.7, 19.0, 14.7, 9.28, 5.96)),
    RegionParam(11, '黄羊河等 F>800', 0.631, (88.4, 75.0, 58.7, 46.9, 36.4, 23.3, 14.8)),
    RegionParam(12, '西营河', 0.36, (10.6, 9.5, 8.14, 6.98, 5.99, 4.62, 3.62)),
    RegionParam(13, '东大河、西大河', 0.477, (16.3, 14.5, 12.4, 10.8, 8.9, 6.71, 5.2)),
    RegionParam(14, '内陆河流走廊地区 F>10000', 3.81, (661000000000000.0, 522000000000000.0, 409000000000000.0, 339000000000000.0, 261000000000000.0, 202000000000000.0, 129000000000000.0)),
    RegionParam(15, '大夏河 F<6200', 0.382, (4.93, 4.34, 3.64, 3.15, 2.62, 2.02, 1.55)),
    RegionParam(16, '大夏河 F>6200', -5.1, (9.119999999999998e-21, 7.519999999999999e-21, 6.7e-21, 5.4699999999999994e-21, 4.559999999999999e-21, 3.3699999999999996e-21, 2.5999999999999998e-21)),
    RegionParam(17, '洮河干流 F<7400', 0.441, (8.19, 7.22, 6.36, 5.66, 4.89, 3.88, 3.1)),
    RegionParam(18, '洮河干流 7400<F<15000', -0.899, (0.00455, 0.00404, 0.0035499999999999998, 0.0031000000000000003, 0.00271, 0.00215, 0.00175)),
    RegionParam(19, '洮河干流 F>15000', 0.26, (2.44, 2.21, 1.94, 1.72, 1.48, 1.2, 0.964)),
    RegionParam(20, '卓尼以上，洮河干支流', 0.19, (0.89, 0.789, 0.696, 0.614, 0.532, 0.42, 0.346)),
    RegionParam(21, '羊沙河、大峪河、冶木河', 0.273, (3.6, 3.22, 2.77, 2.46, 2.08, 1.71, 1.42)),
    RegionParam(22, '胭脂河、东峪沟、大碧河', 0.299, (23.4, 20.2, 16.6, 13.9, 11.1, 8.32, 6.14)),
    RegionParam(23, '广通河', 0.324, (11.8, 10.2, 8.23, 6.9, 5.83, 4.41, 3.52)),
    RegionParam(24, '纳纳河、漫坝河、迭藏河', 0.305, (7.19, 6.24, 5.18, 4.43, 3.57, 2.74, 2.06)),
    RegionParam(25, '大通河', 0.807, (486.0, 435.0, 382.0, 330.0, 290.0, 227.0, 183.0)),
    RegionParam(26, '庄浪河', 0.42, (10.9, 9.24, 7.54, 6.27, 5.06, 3.74, 2.75)),
    RegionParam(27, '关川河', 0.2, (3.46, 3.1, 2.63, 2.25, 1.91, 1.43, 1.1)),
    RegionParam(28, '祖厉河干流', 0.74, (356.0, 315.0, 264.0, 233.0, 196.0, 148.0, 114.0)),
    RegionParam(29, '祖厉河支流', 0.43, (43.0, 37.8, 32.3, 26.8, 22.7, 16.8, 11.7)),
    RegionParam(30, '甘谷以上渭河干流以南支流', 0.479, (30.0, 25.1, 20.8, 16.7, 13.7, 9.53, 6.99)),
    RegionParam(31, '渭河干流以北 F<3900', 0.596, (196.9, 170.5, 140.1, 117.7, 98.7, 73.5, 55.6)),
    RegionParam(32, '渭河干流以北 3900<F<7600', 0.817, (1220.0, 1040.0, 859.0, 719.0, 600.0, 438.0, 342.0)),
    RegionParam(33, '渭河干流以北 7600<F<23500', 0.39, (22.8, 19.7, 16.3, 13.9, 11.8, 8.85, 7.0)),
    RegionParam(34, '渭河干流以北 F>23500', 0.919, (4660.0, 4070.0, 3390.0, 2880.0, 2420.0, 1830.0, 1460.0)),
    RegionParam(35, '牛头河、水洛河、清水河', 0.496, (75.6, 61.9, 50.1, 40.8, 33.4, 24.6, 17.8)),
    RegionParam(36, '葫芦河干流 F<2800', 0.306, (9.79, 8.31, 6.6, 5.34, 4.3, 3.01, 2.17)),
    RegionParam(37, '葫芦河干流 F>2800', -1.61, (4.84, 3.83, 2.26, 1.49, 0.796, 0.376, 0.147)),
    RegionParam(38, '散渡河', 0.42, (66.6, 58.9, 47.3, 40.7, 33.8, 24.8, 19.3)),
    RegionParam(39, '静宁以下葫芦河干流以西支流', 0.335, (51.9, 44.0, 35.6, 29.3, 23.6, 17.0, 12.6)),
    RegionParam(40, '泾河干流及汭河 250<F<14000', 0.266, (13.1, 11.1, 8.92, 7.4, 5.97, 4.25, 3.13)),
    RegionParam(41, '泾河干流及汭河 14000<F<40000', 0.433, (63.0, 54.0, 43.1, 36.2, 29.3, 20.7, 15.6)),
    RegionParam(42, '泾河干流及汭河 F>40000', 2.63, (763000000000.0, 634000000000.0, 533000000000.0, 446000000000.0, 367000000000.0, 266000000000.0, 200000000000.0)),
    RegionParam(43, '茹河、颉河三关口以上', 0.401, (43.8, 36.8, 29.2, 23.1, 18.2, 12.4, 8.76)),
    RegionParam(44, '泾河崆峒峡以上及汭河 F<250', 0.398, (29.9, 24.7, 20.0, 16.5, 13.0, 9.37, 6.64)),
    RegionParam(45, '崆峒峡至泾川区间支流及洪河', 0.354, (59.0, 49.7, 39.3, 31.6, 24.9, 16.5, 11.8)),
    RegionParam(46, '蒲河支流交口河', 0.794, (890.0, 735.0, 577.0, 457.0, 341.0, 217.0, 147.0)),
    RegionParam(47, '蒲河干流及巴家咀以上 F<3500', 0.33, (53.2, 45.3, 35.5, 28.9, 23.0, 16.0, 11.6)),
    RegionParam(48, '蒲河干流及巴家咀以上 F>3500', 0.895, (5430.0, 4600.0, 3570.0, 2930.0, 2300.0, 1600.0, 1160.0)),
    RegionParam(49, '马莲河东川', 0.243, (22.7, 19.3, 15.5, 13.1, 10.4, 7.34, 5.42)),
    RegionParam(50, '马莲河干流及马莲河西川洪德以下', 0.517, (131.0, 112.0, 90.7, 75.6, 60.5, 43.7, 32.0)),
    RegionParam(51, '马莲河庆阳至雨落坪区间支流', 0.283, (37.2, 30.9, 23.7, 18.5, 14.0, 9.16, 5.78)),
    RegionParam(52, '马莲河洪德以上', 0.516, (61.2, 52.4, 44.5, 38.0, 32.1, 24.7, 19.4)),
    RegionParam(53, '达溪河 F<2500', 0.382, (31.4, 26.1, 20.7, 16.8, 13.3, 9.24, 6.51)),
    RegionParam(54, '达溪河 F>2500', 0.642, (240.0, 199.0, 158.0, 128.0, 100.0, 69.9, 49.4)),
    RegionParam(55, '黑河(泾河水系)', 0.597, (168.0, 140.0, 111.0, 89.7, 71.8, 49.9, 35.9)),
    RegionParam(56, '新店子以上嘉陵江干流及东坡以上支流', 0.274, (12.8, 11.2, 9.41, 8.07, 6.9, 5.34, 4.26)),
    RegionParam(57, '东坡至谈家庄区间支流', 0.288, (25.9, 21.9, 17.3, 14.6, 11.9, 8.37, 6.24)),
    RegionParam(58, '西汉水干流及毛坝以上 F<3400', 0.52, (73.8, 63.6, 51.8, 43.3, 36.9, 27.5, 21.2)),
    RegionParam(59, '西汉水干流及毛坝以上 F>3400', 0.333, (15.9, 13.7, 11.1, 9.58, 8.04, 5.9, 4.58)),
    RegionParam(60, '略阳以下嘉陵江右岸支流', 0.149, (20.9, 18.4, 15.6, 13.3, 11.0, 8.52, 6.57)),
    RegionParam(61, '玉垒至跳渡区间支流 F<800', 0.417, (59.1, 52.9, 43.8, 37.8, 31.2, 23.5, 17.8)),
    RegionParam(62, '白龙江干流两水以上 F<15000', 0.22, (1.61, 1.46, 1.32, 1.16, 1.03, 0.841, 0.702)),
    RegionParam(63, '白龙江干流 15000<F<26000', -0.458, (0.00222, 0.00199, 0.00177, 0.00158, 0.00141, 0.00117, 0.000989)),
    RegionParam(64, '白龙江干流 F>26000', -4.42, (7.05e-21, 6.039999999999999e-21, 5.389999999999999e-21, 4.789999999999999e-21, 4.15e-21, 3.6299999999999995e-21, 3.1099999999999997e-21)),
    RegionParam(65, '蒿子店至两水区间支流', 0.43, (23.5, 21.0, 17.6, 15.2, 12.3, 8.98, 6.66)),
    RegionParam(66, '白水江干流', -1.28, (2.6499999999999996e-06, 2.39e-06, 2.1499999999999997e-06, 1.8999999999999998e-06, 1.69e-06, 1.42e-06, 1.18e-06)),
    RegionParam(67, '白水江支流', 0.722, (212.0, 185.0, 149.0, 122.0, 95.3, 65.8, 43.7)),
]

REGION_BY_ID = {item.region_id: item for item in REGION_PARAMS}


def regional_formula_parameter_rows(region_id: int) -> list[tuple[str, float, float]]:
    if region_id not in REGION_BY_ID:
        raise HydrologyError("必须选择有效地区经验分区，不能默认按第 1 分区计算")
    region = REGION_BY_ID[region_id]
    return [
        (freq, region.n, k)
        for freq, k in zip(REGION_FREQUENCIES, region.k_values)
    ]


@dataclass(frozen=True)
class RegionalFormulaInput:
    basin_name: str
    area: float
    region_id: int
    allow_out_of_range: bool = False


@dataclass(frozen=True)
class RegionalFormulaResult:
    basin_name: str
    region_id: int
    region_name: str
    n: float
    flows: dict[str, float]
    warnings: tuple[str, ...] = ()


def calculate_regional_formula(data: RegionalFormulaInput) -> RegionalFormulaResult:
    area = _require_positive("面积 F", data.area)
    if data.region_id not in REGION_BY_ID:
        raise HydrologyError("必须选择有效地区经验分区，不能默认按第 1 分区计算")
    region = REGION_BY_ID[data.region_id]
    if not data.allow_out_of_range:
        region.validate_area(area)

    flows = {
        freq: _safe_round(k * (area ** (1 - region.n)))
        for freq, k in zip(REGION_FREQUENCIES, region.k_values)
    }
    warnings: list[str] = []
    adjacent_ratios = [
        region.k_values[index] / region.k_values[index + 1]
        for index in range(len(region.k_values) - 1)
        if region.k_values[index + 1] > 0
    ]
    if adjacent_ratios and max(adjacent_ratios) > 3:
        warnings.append(
            "源参数表相邻频率K值存在超过3倍的突变。当前数值与随软件提供的源Excel一致，"
            "但正式采用前必须核对原始出版资料或主管部门成果，不能仅凭本软件结果定稿。"
        )
    return RegionalFormulaResult(
        basin_name=data.basin_name,
        region_id=region.region_id,
        region_name=region.name,
        n=region.n,
        flows=flows,
        warnings=tuple(warnings),
    )


@dataclass(frozen=True)
class InstantUnitHydrographParams:
    n: float
    k: float
    m1: float
    m2: float


@dataclass(frozen=True)
class UnitHydrographOrdinate:
    time_h: float
    t_over_k: float
    s_t: float
    s_prev: float
    unit_depth_fraction: float


@dataclass(frozen=True)
class FloodHydrographPoint:
    time_h: float
    discharge_m3s: float


@dataclass(frozen=True)
class FrequencyRunoffInput:
    frequency: str
    ptc_mm: float
    runoff_duration_h: int
    infiltration_rate_mmh: float
    n: float
    k: float
    runoff_start_hour: int = 1


@dataclass(frozen=True)
class NetRainRow:
    hour: int
    rainfall_mm: dict[str, float]
    infiltration_mm: dict[str, float]
    net_rain_mm: dict[str, float]


@dataclass(frozen=True)
class FrequencyFloodResult:
    frequency: str
    unit_line: list[UnitHydrographOrdinate]
    hydrograph: list[FloodHydrographPoint]
    peak_flow_m3s: float
    peak_time_h: float
    flood_volume_10k_m3: float


@dataclass(frozen=True)
class InstantUnitLineBatchResult:
    net_rain_table: list[NetRainRow]
    flood_results: dict[str, FrequencyFloodResult]
    peak_summary: dict[str, tuple[float, float, float]]
    frequency_inputs: dict[str, FrequencyRunoffInput] | None = None
    kp_values: dict[str, float] | None = None
    ptc_mean_mm: float | None = None
    cv: float | None = None
    cs_ratio: float | None = None
    area_km2: float | None = None
    delta_t_h: float | None = None
    rainfall_zone_name: str | None = None
    main_peak_hours: int | None = None
    distribution_source: str | None = None
    duration_kp_values: dict[str, dict[str, float]] | None = None
    duration_point_rainfall: dict[str, dict[str, float]] | None = None
    duration_areal_rainfall: dict[str, dict[str, float]] | None = None
    rainfall_process_components: list[tuple[int, float, float, float, float]] | None = None
    unit_parameter_rows: dict[str, dict[str, float]] | None = None


@dataclass(frozen=True)
class RationalFormulaFrequencyResult:
    frequency: str
    hydrograph: list[FloodHydrographPoint]
    peak_flow_m3s: float
    peak_time_h: float
    surface_peak_m3s: float
    tau_h: float
    m: float
    theta: float
    runoff_condition: str
    max_net_duration_h: float
    max_net_depth_mm: float
    max_net_intensity_mmh: float
    main_start_h: float
    flood_volume_10k_m3: float
    curve_markers: tuple[tuple[str, float, float], ...] = ()
    confluence_curve: tuple[tuple[float, float, float], ...] = ()
    confluence_curve_points: tuple[tuple[float, float, float, float], ...] = ()
    hydrograph_components: tuple[tuple[tuple[float, float], ...], ...] = ()


@dataclass(frozen=True)
class RationalFormulaBatchResult:
    shared: InstantUnitLineBatchResult
    rational_region_id: str
    rational_region_name: str
    river_length_km: float
    river_slope_permille: float
    baseflow_m3s: float
    shape_ka: float
    shape_kw: float
    shape_kt: float
    shape_theta: float
    shape_source: str
    baseflow_source: str
    frequency_results: dict[str, RationalFormulaFrequencyResult]
    parameter_rows: dict[str, dict[str, float | str]]


@dataclass(frozen=True)
class InstantRainfallDistribution:
    zone_id: str
    zone_name: str
    main_peak_hours: int
    weights: tuple[float, ...]
    source: str
    component_weights: tuple[tuple[float, float, float, float], ...] | None = None


@dataclass(frozen=True)
class InstantRainfallDesignInput:
    h1_mean_mm: float
    h1_cv: float
    h6_mean_mm: float
    h6_cv: float
    h24_mean_mm: float
    h24_cv: float
    alpha1: float = 1.0
    alpha3: float = 1.0
    alpha6: float = 1.0
    alpha24: float = 1.0
    gamma1: float = 1.0
    gamma3: float = 1.0
    gamma6: float = 1.0
    gamma24: float = 1.0
    cs_ratio: float = 3.5


@dataclass(frozen=True)
class FrequencyRunoffSetting:
    frequency: str
    infiltration_rate_mmh: float | None = None
    runoff_duration_h: int | None = None
    runoff_start_hour: int | None = None
    unit_param_region: str = "huangtu"
    infiltration_table_id: str = "hexi_shishan"


INSTANT_INFILTRATION_TABLE_LABELS = (
    ("hexi_shishan", "河西石山森林草原区（附表12一）"),
    ("longnan_shishan", "陇南石山林区（附表12二）"),
    ("huangtu", "黄土区（附表12三）"),
    ("huangtu_shanlin", "黄土山（林）区（附表12四）"),
    ("liupanshan_tushishan", "六盘山土石山区（附表12五）"),
)


INSTANT_INFILTRATION_TABLES = {
    "hexi_shishan": {
        "name": "河西石山森林草原区（附表12一）",
        "ptc_values": (
            10, 15, 20, 25, 30, 35, 40, 45, 50, 55,
            60, 65, 70, 75, 80, 85, 90, 95, 100, 105,
            110, 115, 120, 125, 130, 135, 140, 145, 150,
            160, 170, 180, 190, 200, 210, 230, 250, 300,
        ),
        "tc_values": tuple(range(1, 21)),
        "f_values": (
            (5.75, 3.11, 2.17, 1.68, 1.38, 1.17, 1.02, 0.91, 0.82, 0.74, 0.68, 0.63, 0.59, 0.55, 0.52, 0.50, 0.46, 0.44, 0.42, 0.40),
            (8.17, 4.44, 3.10, 2.41, 1.98, 1.68, 1.47, 1.30, 1.17, 1.07, 0.98, 0.91, 0.85, 0.79, 0.75, 0.71, 0.67, 0.64, 0.61, 0.58),
            (10.35, 6.33, 3.94, 3.06, 2.52, 2.15, 1.87, 1.67, 1.50, 1.37, 1.26, 1.17, 1.09, 1.02, 0.96, 0.90, 0.86, 0.82, 0.78, 0.74),
            (12.26, 7.04, 4.71, 3.66, 3.01, 2.57, 2.24, 2.00, 1.80, 1.64, 1.51, 1.40, 1.31, 1.22, 1.15, 1.09, 1.03, 0.98, 0.94, 0.89),
            (14.07, 7.66, 5.39, 4.20, 3.46, 2.95, 2.58, 2.30, 2.08, 1.89, 1.74, 1.62, 1.51, 1.41, 1.33, 1.26, 1.19, 1.13, 1.08, 1.03),
            (15.58, 8.53, 6.02, 4.69, 3.87, 3.31, 2.89, 2.58, 2.33, 2.12, 1.96, 1.81, 1.69, 1.59, 1.49, 1.41, 1.34, 1.27, 1.22, 1.16),
            (16.89, 9.31, 6.58, 5.14, 4.24, 3.63, 3.18, 2.83, 2.56, 2.33, 2.15, 1.99, 1.86, 1.75, 1.64, 1.56, 1.48, 1.40, 1.34, 1.28),
            (18.11, 10.07, 7.09, 5.54, 4.58, 3.92, 3.43, 3.06, 2.77, 2.53, 2.33, 2.16, 2.02, 1.89, 1.78, 1.69, 1.60, 1.52, 1.46, 1.39),
            (19.11, 10.77, 7.55, 5.91, 4.89, 4.18, 3.67, 3.27, 2.96, 2.70, 2.49, 2.31, 2.16, 2.03, 1.91, 1.81, 1.72, 1.63, 1.56, 1.49),
            (20.11, 11.27, 7.96, 6.24, 5.17, 4.43, 3.88, 3.47, 3.14, 2.87, 2.64, 2.45, 2.29, 2.15, 2.03, 1.92, 1.83, 1.74, 1.66, 1.59),
            (21.01, 11.78, 8.34, 6.54, 5.42, 4.65, 4.08, 3.64, 3.30, 3.02, 2.78, 2.58, 2.41, 2.27, 2.14, 2.02, 1.92, 1.83, 1.75, 1.67),
            (21.70, 12.20, 8.68, 6.82, 5.65, 4.85, 4.26, 3.81, 3.45, 3.15, 2.91, 2.70, 2.52, 2.37, 2.24, 2.12, 2.01, 1.92, 1.83, 1.75),
            (22.40, 12.60, 8.99, 7.07, 5.86, 5.03, 4.42, 3.95, 3.58, 3.28, 3.02, 2.81, 2.63, 2.47, 2.33, 2.21, 2.10, 2.00, 1.91, 1.83),
            (23.00, 13.00, 9.26, 7.29, 6.05, 5.20, 4.57, 4.09, 3.70, 3.39, 3.13, 2.91, 2.72, 2.56, 2.41, 2.29, 2.17, 2.07, 1.98, 1.89),
            (23.60, 13.30, 9.52, 7.50, 6.23, 5.35, 4.71, 4.21, 3.82, 3.50, 3.23, 3.00, 2.81, 2.64, 2.49, 2.36, 2.24, 2.14, 2.04, 1.96),
            (24.00, 13.60, 9.74, 7.68, 6.39, 5.49, 4.83, 4.32, 3.92, 3.59, 3.32, 3.09, 2.89, 2.71, 2.56, 2.43, 2.31, 2.20, 2.10, 2.01),
            (24.50, 13.90, 9.95, 7.85, 6.53, 5.62, 4.95, 4.43, 4.02, 3.68, 3.40, 3.16, 2.96, 2.78, 2.63, 2.49, 2.37, 2.26, 2.16, 2.07),
            (24.90, 14.10, 10.10, 8.00, 6.66, 5.73, 5.04, 4.52, 4.10, 3.76, 3.48, 3.23, 3.03, 2.85, 2.69, 2.55, 2.42, 2.31, 2.21, 2.12),
            (25.20, 14.30, 10.30, 8.14, 6.78, 5.84, 5.14, 4.61, 4.18, 3.84, 3.55, 3.30, 3.09, 2.91, 2.74, 2.60, 2.47, 2.36, 2.26, 2.16),
            (25.50, 14.50, 10.40, 8.27, 6.89, 5.93, 5.23, 4.69, 4.26, 3.90, 3.61, 3.36, 3.15, 2.96, 2.80, 2.65, 2.52, 2.41, 2.30, 2.20),
            (25.80, 14.70, 10.60, 8.38, 6.99, 6.02, 5.31, 4.76, 4.32, 3.97, 3.67, 3.41, 3.20, 3.01, 2.84, 2.70, 2.56, 2.45, 2.34, 2.24),
            (26.00, 14.90, 10.70, 8.48, 7.08, 6.10, 5.38, 4.83, 4.38, 4.02, 3.72, 3.46, 3.25, 3.05, 2.89, 2.74, 2.60, 2.49, 2.38, 2.28),
            (26.20, 15.00, 10.80, 8.58, 7.16, 6.17, 5.45, 4.89, 4.44, 4.07, 3.77, 3.51, 3.29, 3.10, 2.93, 2.78, 2.64, 2.52, 2.41, 2.31),
            (26.40, 15.20, 10.90, 8.66, 7.23, 6.24, 5.51, 4.94, 4.49, 4.12, 3.81, 3.55, 3.33, 3.13, 2.96, 2.81, 2.68, 2.55, 2.44, 2.34),
            (26.60, 15.30, 11.00, 8.74, 7.30, 6.30, 5.56, 4.99, 4.54, 4.17, 3.86, 3.59, 3.37, 3.17, 3.00, 2.84, 2.71, 2.58, 2.47, 2.37),
            (26.70, 15.40, 11.10, 8.81, 7.36, 6.36, 5.61, 5.04, 4.58, 4.21, 3.89, 3.63, 3.40, 3.20, 3.03, 2.87, 2.73, 2.61, 2.50, 2.40),
            (26.90, 15.50, 11.20, 8.87, 7.42, 6.41, 5.66, 5.08, 4.62, 4.24, 3.93, 3.66, 3.43, 3.23, 3.06, 2.90, 2.76, 2.64, 2.52, 2.42),
            (27.00, 15.60, 11.30, 8.93, 7.47, 6.45, 5.70, 5.12, 4.66, 4.28, 3.96, 3.69, 3.46, 3.26, 3.08, 2.93, 2.79, 2.66, 2.55, 2.44),
            (27.10, 15.60, 11.30, 8.98, 7.51, 6.49, 5.74, 5.15, 4.69, 4.31, 3.99, 3.72, 3.49, 3.28, 3.11, 2.95, 2.81, 2.68, 2.57, 2.46),
            (27.30, 15.80, 11.40, 9.08, 7.56, 6.57, 5.80, 5.22, 4.75, 4.36, 4.04, 3.77, 3.53, 3.33, 3.15, 2.99, 2.85, 2.72, 2.60, 2.50),
            (27.40, 15.90, 11.50, 9.15, 7.66, 6.63, 5.86, 5.27, 4.79, 4.41, 4.08, 3.81, 3.57, 3.37, 3.19, 3.03, 2.88, 2.75, 2.64, 2.53),
            (27.60, 15.90, 11.60, 9.21, 7.72, 6.68, 5.91, 5.31, 4.84, 4.45, 4.12, 3.84, 3.61, 3.40, 3.22, 3.06, 2.91, 2.78, 2.66, 2.56),
            (27.60, 16.00, 11.60, 9.26, 7.76, 6.72, 5.94, 5.35, 4.87, 4.48, 4.15, 3.87, 3.63, 3.43, 3.24, 3.08, 2.93, 2.80, 2.69, 2.58),
            (27.70, 16.00, 11.70, 9.30, 7.80, 6.75, 5.98, 5.38, 4.90, 4.51, 4.18, 3.90, 3.66, 3.45, 3.26, 3.10, 2.96, 2.82, 2.71, 2.60),
            (27.80, 16.10, 11.70, 9.34, 7.83, 6.78, 6.00, 5.40, 4.92, 4.53, 4.20, 3.92, 3.68, 3.47, 3.28, 3.12, 2.97, 2.84, 2.72, 2.61),
            (27.90, 16.20, 11.80, 9.39, 7.88, 6.82, 6.04, 5.44, 4.96, 4.56, 4.23, 3.95, 3.71, 3.50, 3.31, 3.15, 3.00, 2.87, 2.75, 2.64),
            (27.90, 16.20, 11.80, 9.42, 7.91, 6.85, 6.07, 5.47, 4.98, 4.59, 4.26, 3.97, 3.73, 3.52, 3.33, 3.17, 3.02, 2.89, 2.77, 2.66),
            (28.00, 16.30, 11.80, 9.47, 7.95, 6.89, 6.11, 5.50, 5.02, 4.62, 4.29, 4.01, 3.76, 3.55, 3.36, 3.20, 3.05, 2.92, 2.79, 2.68),
        ),
    }
    ,
    "longnan_shishan": {
        "name": "陇南石山林区（附表12二）",
        "base_table_id": "hexi_shishan",
        "tc_scale_points": ((1, 1.09), (5, 0.82), (10, 0.60), (15, 0.36), (20, 0.28)),
    },
    "huangtu": {
        "name": "黄土区（附表12三）",
        "base_table_id": "hexi_shishan",
        "tc_scale_points": ((1, 1.20), (5, 1.10), (10, 1.06), (15, 1.03), (20, 1.00)),
    },
    "huangtu_shanlin": {
        "name": "黄土山（林）区（附表12四）",
        "base_table_id": "hexi_shishan",
        "tc_scale_points": ((1, 1.50), (5, 1.42), (10, 1.35), (15, 1.31), (20, 1.27)),
    },
    "liupanshan_tushishan": {
        "name": "六盘山土石山区（附表12五）",
        "base_table_id": "hexi_shishan",
        "tc_scale_points": ((1, 1.34), (5, 1.28), (10, 1.23), (15, 1.20), (20, 1.17)),
    },
}

INSTANT_ZONE_TO_INFILTRATION_TABLE = {
    "hexi": "hexi_shishan",
    "huanghe": "huangtu",
    "weihe": "huangtu",
    "jinghe": "huangtu",
    "jialingjiang": "longnan_shishan",
}


def default_instant_infiltration_table_id(zone_id: str) -> str:
    return INSTANT_ZONE_TO_INFILTRATION_TABLE.get((zone_id or "").strip(), "huangtu")


def _interpolate_1d(points: Iterable[tuple[float, float]], x: float) -> float:
    ordered = sorted((float(px), float(py)) for px, py in points)
    if not ordered:
        raise HydrologyError("插值点不能为空")
    if x <= ordered[0][0]:
        return ordered[0][1]
    if x >= ordered[-1][0]:
        return ordered[-1][1]
    for (x0, y0), (x1, y1) in zip(ordered, ordered[1:]):
        if x0 <= x <= x1:
            if x1 == x0:
                return y0
            ratio = (x - x0) / (x1 - x0)
            return y0 + (y1 - y0) * ratio
    return ordered[-1][1]


def _lookup_grid_value(x_values: Iterable[float], y_values: Iterable[float], values: Iterable[Iterable[float]], x: float, y: float) -> float:
    xs = [float(item) for item in x_values]
    ys = [float(item) for item in y_values]
    grid = [[float(cell) for cell in row] for row in values]
    if len(xs) != len(grid) or any(len(row) != len(ys) for row in grid):
        raise HydrologyError("入渗率查表数据维度不一致")

    def bounds(items: list[float], target: float) -> tuple[int, int, float]:
        if target <= items[0]:
            return 0, 0, 0.0
        if target >= items[-1]:
            last = len(items) - 1
            return last, last, 0.0
        for idx, (left, right) in enumerate(zip(items, items[1:])):
            if left <= target <= right:
                ratio = 0.0 if right == left else (target - left) / (right - left)
                return idx, idx + 1, ratio
        last = len(items) - 1
        return last, last, 0.0

    x0, x1, xr = bounds(xs, x)
    y0, y1, yr = bounds(ys, y)
    q11 = grid[x0][y0]
    q12 = grid[x0][y1]
    q21 = grid[x1][y0]
    q22 = grid[x1][y1]
    top = q11 + (q12 - q11) * yr
    bottom = q21 + (q22 - q21) * yr
    return top + (bottom - top) * xr


def instant_infiltration_rate_from_table(table_id: str, ptc_mm: float, runoff_duration_h: float) -> float:
    """按附表12由 Ptc 和 tc 查产流期平均入渗率 f(mm/h)。"""

    table_id = (table_id or "").strip()
    if table_id not in INSTANT_INFILTRATION_TABLES:
        raise HydrologyError("必须选择有效的入渗率分区")
    ptc = _require_positive("产流期降水量 Ptc", ptc_mm)
    tc = _require_positive("产流历时 tc", runoff_duration_h)
    table = INSTANT_INFILTRATION_TABLES[table_id]
    if "base_table_id" in table:
        base = instant_infiltration_rate_from_table(str(table["base_table_id"]), ptc, tc)
        scale = _interpolate_1d(table["tc_scale_points"], tc)
        return _safe_round(base * scale, 3)
    return _safe_round(
        _lookup_grid_value(table["ptc_values"], table["tc_values"], table["f_values"], ptc, tc),
        3,
    )


def instant_runoff_period_auto_infiltration(
    rainfall_mm: Iterable[float],
    infiltration_table_id: str,
    delta_t_h: float = 1.0,
    min_runoff_duration_h: float = 2.0,
    fallback_min_duration_h: float = 3.0,
    max_runoff_duration_h: float = 12.0,
) -> dict[str, float]:
    """枚举连续产流段，并按附表12自动查 f，选取净雨量最大的有效产流段。

    附表12的 f 由整段产流期降水量 Ptc 和产流历时 tc 查得。候选产流段
    需要满足：段内每个时段面雨量均不小于时段入渗量 f * Δt。这样 tc
    不是固定值，也不会把低雨量尾段硬并入产流期。
    """

    rainfall = [float(value) for value in rainfall_mm]
    if not rainfall:
        raise HydrologyError("逐时段面雨量不能为空")
    delta_t_h = _require_positive("Δt", delta_t_h)
    peak_index = max(range(len(rainfall)), key=lambda idx: rainfall[idx])
    max_steps = min(len(rainfall), max(1, int(round(max_runoff_duration_h / delta_t_h))))

    def choose_candidate(min_duration_h: float, require_each_hour: bool, rule_name: str) -> dict[str, float] | None:
        best: dict[str, float] | None = None
        for start in range(len(rainfall)):
            for end in range(start, min(len(rainfall), start + max_steps)):
                if not (start <= peak_index <= end):
                    continue
                tc = (end - start + 1) * delta_t_h
                if tc + 1e-9 < min_duration_h:
                    continue
                ptc = sum(rainfall[start : end + 1])
                if ptc <= 0:
                    continue
                f_value = instant_infiltration_rate_from_table(infiltration_table_id, ptc, tc)
                infiltration_depth = f_value * delta_t_h
                if require_each_hour and any(value + 1e-9 < infiltration_depth for value in rainfall[start : end + 1]):
                    continue
                net_total = ptc - infiltration_depth * (end - start + 1)
                if net_total <= 0:
                    continue
                candidate = {
                    "start_hour": start + 1,
                    "end_hour": end + 1,
                    "tc": tc,
                    "ptc": _safe_round(ptc, 3),
                    "initial_loss": _safe_round(sum(rainfall[:start]), 3),
                    "infiltration_depth": _safe_round(infiltration_depth, 3),
                    "infiltration_rate_mmh": _safe_round(f_value, 3),
                    "net_total": _safe_round(net_total, 3),
                    "tc_rule": rule_name,
                }
                if best is None or candidate["net_total"] > best["net_total"]:
                    best = candidate
        return best

    strict = choose_candidate(min_runoff_duration_h, True, "严格判定：段内每小时面雨量均不小于f×Δt")
    if strict is not None:
        return strict
    fallback = choose_candidate(
        delta_t_h,
        True,
        "严格缩短：未满足最小时长时，按段内每小时面雨量均不小于f×Δt确定tc",
    )
    if fallback is None:
        raise HydrologyError("没有找到满足附表12自动入渗率的连续产流时段，请检查面雨量或切换入渗率分区")
    return fallback


def _weights_from_hours(values: dict[int, float]) -> tuple[float, ...]:
    return tuple(float(values.get(hour, 0.0)) for hour in range(1, 25))


def _component_weights_from_hours(
    h1: dict[int, float] | None = None,
    h3_minus_h1: dict[int, float] | None = None,
    h6_minus_h3: dict[int, float] | None = None,
    h24_minus_h6: dict[int, float] | None = None,
) -> tuple[tuple[float, float, float, float], ...]:
    h1 = h1 or {}
    h3_minus_h1 = h3_minus_h1 or {}
    h6_minus_h3 = h6_minus_h3 or {}
    h24_minus_h6 = h24_minus_h6 or {}
    return tuple(
        (
            float(h1.get(hour, 0.0)),
            float(h3_minus_h1.get(hour, 0.0)),
            float(h6_minus_h3.get(hour, 0.0)),
            float(h24_minus_h6.get(hour, 0.0)),
        )
        for hour in range(1, 25)
    )


JIALINGJIANG_COMPONENTS_BY_PEAK = {
    1: _component_weights_from_hours(
        h1={10: 1.0},
        h3_minus_h1={11: 0.56, 12: 0.44},
        h6_minus_h3={8: 0.22, 9: 0.48, 13: 0.30},
        h24_minus_h6={
            1: 0.01,
            2: 0.01,
            3: 0.01,
            4: 0.03,
            5: 0.03,
            6: 0.04,
            7: 0.08,
            14: 0.14,
            15: 0.14,
            16: 0.14,
            17: 0.12,
            18: 0.07,
            19: 0.04,
            20: 0.04,
            21: 0.03,
            22: 0.03,
            23: 0.02,
            24: 0.02,
        },
    ),
    3: _component_weights_from_hours(
        h3_minus_h1={9: 0.40, 10: 0.30, 11: 0.30},
        h6_minus_h3={8: 0.40, 12: 0.38, 13: 0.22},
        h24_minus_h6={
            1: 0.01,
            2: 0.02,
            3: 0.02,
            4: 0.03,
            5: 0.03,
            6: 0.08,
            7: 0.12,
            14: 0.12,
            15: 0.11,
            16: 0.11,
            17: 0.10,
            18: 0.06,
            19: 0.05,
            20: 0.05,
            21: 0.03,
            22: 0.03,
            23: 0.02,
            24: 0.01,
        },
    ),
    6: _component_weights_from_hours(
        h6_minus_h3={8: 0.12, 9: 0.19, 10: 0.22, 11: 0.18, 12: 0.17, 13: 0.12},
        h24_minus_h6={
            1: 0.01,
            2: 0.02,
            3: 0.03,
            4: 0.04,
            5: 0.05,
            6: 0.05,
            7: 0.06,
            14: 0.28,
            15: 0.09,
            16: 0.07,
            17: 0.07,
            18: 0.05,
            19: 0.05,
            20: 0.05,
            21: 0.03,
            22: 0.02,
            23: 0.02,
            24: 0.01,
        },
    ),
}


HEXI_COMPONENTS_BY_PEAK = {
    1: _component_weights_from_hours(
        h1={10: 1.0},
        h3_minus_h1={11: 0.60, 12: 0.40},
        h6_minus_h3={9: 0.46, 13: 0.40, 14: 0.14},
        h24_minus_h6={
            1: 0.01,
            2: 0.01,
            3: 0.02,
            4: 0.03,
            5: 0.03,
            6: 0.03,
            7: 0.03,
            8: 0.32,
            15: 0.29,
            16: 0.05,
            17: 0.07,
            18: 0.05,
            19: 0.03,
            20: 0.02,
            21: 0.01,
        },
    ),
    3: _component_weights_from_hours(
        h3_minus_h1={9: 0.34, 10: 0.42, 11: 0.24},
        h6_minus_h3={12: 0.35, 13: 0.37, 14: 0.28},
        h24_minus_h6={
            1: 0.01,
            2: 0.01,
            3: 0.02,
            4: 0.03,
            5: 0.04,
            6: 0.07,
            7: 0.11,
            8: 0.25,
            15: 0.26,
            16: 0.03,
            17: 0.03,
            18: 0.07,
            19: 0.03,
            20: 0.02,
            21: 0.02,
        },
    ),
    6: _component_weights_from_hours(
        h6_minus_h3={9: 0.18, 10: 0.28, 11: 0.16, 12: 0.14, 13: 0.13, 14: 0.11},
        h24_minus_h6={
            1: 0.01,
            2: 0.02,
            3: 0.02,
            4: 0.03,
            5: 0.03,
            6: 0.04,
            7: 0.12,
            8: 0.19,
            15: 0.18,
            16: 0.18,
            17: 0.09,
            18: 0.05,
            19: 0.02,
            20: 0.02,
        },
    ),
}

HUANGHE_COMPONENTS_BY_PEAK = {
    1: _component_weights_from_hours(
        h1={10: 1.0},
        h3_minus_h1={11: 0.68, 12: 0.32},
        h6_minus_h3={7: 0.12, 8: 0.15, 9: 0.73},
        h24_minus_h6={6: 0.11, 13: 0.30, 14: 0.24, 15: 0.10, 16: 0.16, 17: 0.09},
    ),
    3: _component_weights_from_hours(
        h3_minus_h1={9: 0.407, 10: 0.401, 11: 0.192},
        h6_minus_h3={7: 0.05, 8: 0.25, 12: 0.70},
        h24_minus_h6={
            1: 0.08,
            2: 0.08,
            3: 0.02,
            5: 0.01,
            6: 0.04,
            13: 0.37,
            14: 0.27,
            15: 0.06,
            16: 0.02,
            17: 0.01,
            18: 0.01,
            19: 0.01,
            20: 0.01,
            21: 0.01,
        },
    ),
    6: _component_weights_from_hours(
        h6_minus_h3={8: 0.128, 9: 0.270, 10: 0.200, 11: 0.180, 12: 0.112, 13: 0.110},
        h24_minus_h6={
            1: 0.03,
            2: 0.09,
            3: 0.10,
            4: 0.03,
            5: 0.01,
            6: 0.11,
            7: 0.18,
            14: 0.23,
            15: 0.08,
            16: 0.04,
            17: 0.03,
            18: 0.02,
            19: 0.02,
            20: 0.01,
            21: 0.01,
            22: 0.01,
        },
    ),
}


WEIHE_COMPONENTS_BY_PEAK = {
    1: _component_weights_from_hours(
        h1={12: 1.0},
        h3_minus_h1={13: 0.64, 14: 0.36},
        h6_minus_h3={11: 0.35, 15: 0.38, 16: 0.27},
        h24_minus_h6={
            3: 0.01,
            5: 0.01,
            7: 0.01,
            8: 0.03,
            9: 0.06,
            10: 0.07,
            17: 0.25,
            18: 0.18,
            19: 0.05,
            20: 0.09,
            21: 0.08,
            22: 0.06,
            23: 0.06,
            24: 0.04,
        },
    ),
    3: _component_weights_from_hours(
        h3_minus_h1={11: 0.509, 12: 0.337, 13: 0.154},
        h6_minus_h3={10: 0.100, 14: 0.550, 15: 0.350},
        h24_minus_h6={
            3: 0.01,
            4: 0.02,
            5: 0.03,
            6: 0.03,
            7: 0.04,
            8: 0.08,
            9: 0.10,
            16: 0.25,
            17: 0.15,
            18: 0.12,
            19: 0.05,
            20: 0.05,
            21: 0.03,
            22: 0.02,
            23: 0.01,
            24: 0.01,
        },
    ),
    6: _component_weights_from_hours(
        h6_minus_h3={11: 0.231, 12: 0.239, 13: 0.145, 14: 0.123, 15: 0.192, 16: 0.070},
        h24_minus_h6={
            3: 0.02,
            4: 0.03,
            5: 0.05,
            6: 0.06,
            7: 0.07,
            8: 0.09,
            9: 0.11,
            10: 0.15,
            17: 0.20,
            18: 0.10,
            19: 0.05,
            20: 0.03,
            21: 0.02,
            22: 0.01,
            23: 0.01,
        },
    ),
}

JINGHE_COMPONENTS_BY_PEAK = {
    1: _component_weights_from_hours(
        h1={11: 1.0},
        h3_minus_h1={10: 0.38, 12: 0.62},
        h6_minus_h3={9: 0.24, 13: 0.40, 14: 0.36},
        h24_minus_h6={6: 0.12, 7: 0.12, 8: 0.12, 15: 0.24, 16: 0.20, 17: 0.20},
    ),
    3: _component_weights_from_hours(
        h3_minus_h1={9: 0.45, 10: 0.30, 11: 0.25},
        h6_minus_h3={12: 0.44, 13: 0.44, 14: 0.12},
        h24_minus_h6={
            1: 0.0,
            2: 0.02,
            3: 0.03,
            4: 0.04,
            5: 0.04,
            6: 0.06,
            7: 0.10,
            8: 0.22,
            15: 0.12,
            16: 0.11,
            17: 0.06,
            18: 0.06,
            19: 0.05,
            20: 0.04,
            21: 0.03,
            22: 0.02,
        },
    ),
    6: _component_weights_from_hours(
        h6_minus_h3={6: 0.300, 7: 0.150, 8: 0.130, 9: 0.136, 10: 0.144, 11: 0.140},
        h24_minus_h6={
            1: 0.02,
            2: 0.04,
            3: 0.06,
            4: 0.07,
            5: 0.10,
            12: 0.18,
            13: 0.11,
            14: 0.09,
            15: 0.08,
            16: 0.05,
            17: 0.04,
            18: 0.04,
            19: 0.03,
            20: 0.03,
            21: 0.02,
            22: 0.02,
            23: 0.01,
            24: 0.01,
        },
    ),
}


INSTANT_RAINFALL_DISTRIBUTIONS: dict[str, dict[int, InstantRainfallDistribution]] = {
    "hexi": {
        1: InstantRainfallDistribution("hexi", "河西区", 1, _weights_from_hours({10: 1.0}), "附表10：河西区，主雨峰1小时", HEXI_COMPONENTS_BY_PEAK[1]),
        3: InstantRainfallDistribution(
            "hexi",
            "河西区",
            3,
            _weights_from_hours({9: 0.34, 10: 0.42, 11: 0.24}),
            "附表10：河西区，主雨峰3小时",
            HEXI_COMPONENTS_BY_PEAK[3],
        ),
        6: InstantRainfallDistribution(
            "hexi",
            "河西区",
            6,
            _weights_from_hours({9: 0.18, 10: 0.28, 11: 0.16, 12: 0.14, 13: 0.13, 14: 0.11}),
            "附表10：河西区，主雨峰6小时",
            HEXI_COMPONENTS_BY_PEAK[6],
        ),
    },
    "huanghe": {
        1: InstantRainfallDistribution("huanghe", "黄河干支流区", 1, _weights_from_hours({10: 1.0}), "附表10：黄河干支流区，主雨峰1小时", HUANGHE_COMPONENTS_BY_PEAK[1]),
        3: InstantRainfallDistribution(
            "huanghe",
            "黄河干支流区",
            3,
            _weights_from_hours({9: 0.407, 10: 0.401, 11: 0.192}),
            "附表10：黄河干支流区，主雨峰3小时",
            HUANGHE_COMPONENTS_BY_PEAK[3],
        ),
        6: InstantRainfallDistribution(
            "huanghe",
            "黄河干支流区",
            6,
            _weights_from_hours({8: 0.128, 9: 0.270, 10: 0.200, 11: 0.180, 12: 0.112, 13: 0.110}),
            "附表10：黄河干支流区，主雨峰6小时",
            HUANGHE_COMPONENTS_BY_PEAK[6],
        ),
    },
    "jialingjiang": {
        1: InstantRainfallDistribution(
            "jialingjiang",
            "嘉陵江区",
            1,
            _weights_from_hours({10: 1.0}),
            "附表10：嘉陵江区，主雨峰1小时，H₁、H₃-H₁、H₆-H₃、H₂₄-H₆列",
            JIALINGJIANG_COMPONENTS_BY_PEAK[1],
        ),
        3: InstantRainfallDistribution(
            "jialingjiang",
            "嘉陵江区",
            3,
            _weights_from_hours({9: 0.40, 10: 0.30, 11: 0.30}),
            "附表10：嘉陵江区，主雨峰3小时，H3列",
            JIALINGJIANG_COMPONENTS_BY_PEAK[3],
        ),
        6: InstantRainfallDistribution(
            "jialingjiang",
            "嘉陵江区",
            6,
            _weights_from_hours({8: 0.12, 9: 0.19, 10: 0.22, 11: 0.18, 12: 0.17, 13: 0.12}),
            "附表10：嘉陵江区，主雨峰6小时，H6列",
            JIALINGJIANG_COMPONENTS_BY_PEAK[6],
        ),
    },
    "weihe": {
        1: InstantRainfallDistribution("weihe", "渭河区", 1, _weights_from_hours({12: 1.0}), "附表10：渭河区，主雨峰1小时", WEIHE_COMPONENTS_BY_PEAK[1]),
        3: InstantRainfallDistribution(
            "weihe",
            "渭河区",
            3,
            _weights_from_hours({11: 0.509, 12: 0.337, 13: 0.154}),
            "附表10：渭河区，主雨峰3小时",
            WEIHE_COMPONENTS_BY_PEAK[3],
        ),
        6: InstantRainfallDistribution(
            "weihe",
            "渭河区",
            6,
            _weights_from_hours({11: 0.231, 12: 0.239, 13: 0.145, 14: 0.123, 15: 0.192, 16: 0.070}),
            "附表10：渭河区，主雨峰6小时",
            WEIHE_COMPONENTS_BY_PEAK[6],
        ),
    },
    "jinghe": {
        1: InstantRainfallDistribution("jinghe", "泾河区", 1, _weights_from_hours({11: 1.0}), "附表10：泾河区，主雨峰1小时", JINGHE_COMPONENTS_BY_PEAK[1]),
        3: InstantRainfallDistribution(
            "jinghe",
            "泾河区",
            3,
            _weights_from_hours({9: 0.45, 10: 0.30, 11: 0.25}),
            "附表10：泾河区，主雨峰3小时",
            JINGHE_COMPONENTS_BY_PEAK[3],
        ),
        6: InstantRainfallDistribution(
            "jinghe",
            "泾河区",
            6,
            _weights_from_hours({6: 0.300, 7: 0.150, 8: 0.130, 9: 0.136, 10: 0.144, 11: 0.140}),
            "附表10：泾河区，主雨峰6小时",
            JINGHE_COMPONENTS_BY_PEAK[6],
        ),
    },
}

INSTANT_RAINFALL_ZONE_LABELS = tuple(
    (zone_id, next(iter(distributions.values())).zone_name)
    for zone_id, distributions in INSTANT_RAINFALL_DISTRIBUTIONS.items()
)


def instant_main_peak_hours_by_area(area_km2: float) -> int:
    """按图集使用说明由流域面积自动选择综合时程分配雨型主雨峰。"""

    area = _require_positive("流域面积 F", area_km2)
    if area < 100:
        return 1
    if area <= 300:
        return 3
    if area <= 1000:
        return 6
    raise HydrologyError("综合概化雨型自动选择目前适用 F<=1000km²；更大流域请手动输入雨型系数")


def select_instant_rainfall_distribution(zone_id: str, area_km2: float) -> InstantRainfallDistribution:
    zone_id = (zone_id or "").strip()
    if zone_id not in INSTANT_RAINFALL_DISTRIBUTIONS:
        raise HydrologyError("必须先选择有效的雨型分区")
    peak_hours = instant_main_peak_hours_by_area(area_km2)
    try:
        return INSTANT_RAINFALL_DISTRIBUTIONS[zone_id][peak_hours]
    except KeyError as exc:
        raise HydrologyError(f"{zone_id} 分区缺少主雨峰 {peak_hours}h 雨型数据") from exc


def nash_params_from_m1_m2(m1: float, m2: float) -> InstantUnitHydrographParams:
    """由瞬时单位线矩参数求 Nash 参数。

    《图集》使用说明：n = 1 / m2，k = m1 / n。
    """

    m1 = _require_positive("m1", m1)
    m2 = _require_positive("m2", m2)
    n = 1 / m2
    k = m1 / n
    return InstantUnitHydrographParams(n=n, k=k, m1=m1, m2=m2)


def m1_from_m110(m1_10: float, b: float, rainfall_intensity: float) -> float:
    """由分区综合关系 m1 = m1,10 * (10 / i)^b 求设计雨强下的 m1."""

    m1_10 = _require_positive("m1,10", m1_10)
    rainfall_intensity = _require_positive("净雨强 i", rainfall_intensity)
    return m1_10 * ((10 / rainfall_intensity) ** b)


def nash_params_from_regional_moments(m1_10: float, b: float, m2: float, rainfall_intensity: float) -> InstantUnitHydrographParams:
    m1 = m1_from_m110(m1_10, b, rainfall_intensity)
    return nash_params_from_m1_m2(m1, m2)


def _regularized_gamma_p(a: float, x: float) -> float:
    """Regularized lower incomplete gamma P(a, x), no scipy dependency."""

    if a <= 0:
        raise HydrologyError("n 必须大于 0")
    if x <= 0:
        return 0.0
    eps = 1e-12
    max_iter = 200
    gln = math.lgamma(a)

    if x < a + 1:
        ap = a
        total = 1 / a
        delta = total
        for _ in range(max_iter):
            ap += 1
            delta *= x / ap
            total += delta
            if abs(delta) < abs(total) * eps:
                break
        return max(0.0, min(1.0, total * math.exp(-x + a * math.log(x) - gln)))

    b = x + 1 - a
    c = 1 / 1e-300
    d = 1 / b
    h = d
    for i in range(1, max_iter + 1):
        an = -i * (i - a)
        b += 2
        d = an * d + b
        if abs(d) < 1e-300:
            d = 1e-300
        c = b + an / c
        if abs(c) < 1e-300:
            c = 1e-300
        d = 1 / d
        delta = d * c
        h *= delta
        if abs(delta - 1) < eps:
            break
    q = h * math.exp(-x + a * math.log(x) - gln)
    return max(0.0, min(1.0, 1 - q))


def _inverse_regularized_gamma_p(a: float, p: float) -> float:
    """Inverse of regularized lower incomplete gamma P(a, x)."""

    a = _require_positive("Gamma形状参数", a)
    if not 0 < p < 1:
        raise HydrologyError("Gamma反算概率必须在 0-1 之间")

    low = 0.0
    high = max(1.0, a)
    while _regularized_gamma_p(a, high) < p:
        high *= 2
        if high > 1e6:
            raise HydrologyError("Gamma反算未收敛，请检查 Cv/Cs 参数")

    for _ in range(100):
        mid = (low + high) / 2
        value = _regularized_gamma_p(a, mid)
        if value < p:
            low = mid
        else:
            high = mid
    return (low + high) / 2


def pearson_iii_kp(cv: float, cs_ratio: float, probability_percent: float) -> float:
    """按用户提供的 Excel 表公式计算 P-III 模比系数 Kp。

    Excel 公式：
    tp = GAMMAINV(1 - P/100, 4 / Cs^2, 1)
    Phi = Cs / 2 * tp - 2 / Cs
    Kp = 1 + Cv * Phi
    """

    cv = _require_positive("Cv", cv)
    cs_ratio = _require_positive("Cs/Cv", cs_ratio)
    probability = _require_positive("频率P", probability_percent)
    if probability >= 100:
        raise HydrologyError("频率P必须小于100%")
    cs = cv * cs_ratio
    shape = 4 / (cs * cs)
    tp = _inverse_regularized_gamma_p(shape, 1 - probability / 100)
    phi = cs / 2 * tp - 2 / cs
    return _safe_round(1 + cv * phi, 3)


def pearson_iii_design_values(
    mean_value: float,
    cv: float,
    frequencies: Iterable[str] = DESIGN_FREQUENCIES_7,
    cs_ratio: float = 3.5,
) -> tuple[dict[str, float], dict[str, float]]:
    mean_value = _require_positive("均值", mean_value)
    kp_values: dict[str, float] = {}
    design_values: dict[str, float] = {}
    for freq in frequencies:
        if freq not in DESIGN_FREQUENCY_PROBABILITIES:
            raise HydrologyError(f"缺少 {freq} 的频率概率定义")
        kp = pearson_iii_kp(cv, cs_ratio, DESIGN_FREQUENCY_PROBABILITIES[freq])
        kp_values[freq] = kp
        design_values[freq] = _safe_round(mean_value * kp, 3)
    return kp_values, design_values


def instant_duration_design_rainfall(
    data: InstantRainfallDesignInput,
    frequencies: Iterable[str] = DESIGN_FREQUENCIES_7,
) -> tuple[dict[str, dict[str, float]], dict[str, dict[str, float]], dict[str, dict[str, float]]]:
    """按参考表口径计算 1h、3h、6h、24h 点雨量与面雨量。

    1h、6h、24h 由均值、Cv、Cs/Cv 计算 Kp；3h 按 H3=H1^0.387*H6^0.613。
    面雨量 = 点雨量 * 点面折减系数 alpha_t * 流域形状修正系数 gamma。
    """

    frequency_list = list(frequencies)
    inputs = {
        "1h": (data.h1_mean_mm, data.h1_cv, data.alpha1, data.gamma1),
        "6h": (data.h6_mean_mm, data.h6_cv, data.alpha6, data.gamma6),
        "24h": (data.h24_mean_mm, data.h24_cv, data.alpha24, data.gamma24),
    }
    kp_by_duration: dict[str, dict[str, float]] = {}
    point_by_duration: dict[str, dict[str, float]] = {}
    areal_by_duration: dict[str, dict[str, float]] = {}
    for duration, (mean, cv, alpha, gamma) in inputs.items():
        alpha = _require_positive(f"{duration}点面折减系数alpha", alpha)
        gamma = _require_positive(f"{duration}流域形状修正系数gamma", gamma)
        kp_values, point_values = pearson_iii_design_values(mean, cv, frequency_list, data.cs_ratio)
        kp_by_duration[duration] = kp_values
        point_by_duration[duration] = point_values
        areal_by_duration[duration] = {
            freq: _safe_round(point_values[freq] * alpha * gamma, 3)
            for freq in frequency_list
        }

    point_by_duration["3h"] = {}
    kp_by_duration["3h"] = {}
    areal_by_duration["3h"] = {}
    mean_h3 = (data.h1_mean_mm ** 0.387) * (data.h6_mean_mm ** 0.613)
    for freq in frequency_list:
        point_h3 = (point_by_duration["1h"][freq] ** 0.387) * (point_by_duration["6h"][freq] ** 0.613)
        point_by_duration["3h"][freq] = _safe_round(point_h3, 3)
        kp_by_duration["3h"][freq] = _safe_round(point_h3 / mean_h3, 3)
        areal_by_duration["3h"][freq] = _safe_round(point_h3 * data.alpha3 * data.gamma3, 3)
    return kp_by_duration, point_by_duration, areal_by_duration


INSTANT_UNIT_PARAM_FORMULAS = {
    "huangtu": {
        "name": "黄土区公式（表8/示例）",
        "m1_10_a": 0.94,
        "m1_10_f": 0.16,
        "m1_10_j": -0.33,
        "b_a": 0.355,
        "b_logf": -0.060,
        "m2_a": 0.69,
        "m2_f": -0.19,
        "m2_j": 0.30,
    },
    "excel_weibei_ref": {
        "name": "参考表旧公式（暂保留）",
        "m1_10_a": 0.997,
        "m1_10_f": 0.33,
        "m1_10_j": -0.313,
        "b_a": 1.024,
        "b_logf": -0.15,
        "m2_a": 1.637,
        "m2_f": -0.331,
        "m2_j": 0.257,
    },
}

INSTANT_UNIT_PARAM_FORMULA_LABELS = tuple(
    (formula_id, formula["name"])
    for formula_id, formula in INSTANT_UNIT_PARAM_FORMULAS.items()
)


def instant_unit_params_from_area_slope_intensity(
    area_km2: float,
    river_slope_permille: float,
    net_rain_intensity_mmh: float,
    net_rain_depth_mm: float | None = None,
    formula_id: str = "huangtu",
) -> dict[str, float]:
    """按图集汇流分区综合关系计算 m1、m2、n、k。"""

    area = _require_positive("流域面积 F", area_km2)
    slope = _require_positive("主河槽平均坡度 J", river_slope_permille)
    intensity = _require_positive("平均净雨强 i", net_rain_intensity_mmh)
    if formula_id not in INSTANT_UNIT_PARAM_FORMULAS:
        raise HydrologyError("必须选择有效的瞬时单位线汇流参数公式")
    formula = INSTANT_UNIT_PARAM_FORMULAS[formula_id]
    m1_10 = formula["m1_10_a"] * (area ** formula["m1_10_f"]) * (slope ** formula["m1_10_j"])
    b = formula["b_a"] + formula["b_logf"] * math.log10(area)
    m1 = m1_10 * ((10 / intensity) ** b)
    m2 = formula["m2_a"] * (area ** formula["m2_f"]) * (slope ** formula["m2_j"])
    n = round(1 / m2, 1)
    k = m1 / n
    return {
        "formula_id": formula_id,
        "formula_name": formula["name"],
        "m1_10": _safe_round(m1_10, 3),
        "b": _safe_round(b, 3),
        "m1": _safe_round(m1, 3),
        "m2": _safe_round(m2, 3),
        "n": _safe_round(n, 3),
        "k": _safe_round(k, 3),
    }


def _instant_rainfall_process_by_frequency(
    distribution: InstantRainfallDistribution,
    duration_areal_rainfall: dict[str, dict[str, float]],
    frequencies: Iterable[str],
) -> tuple[dict[str, list[float]], list[tuple[int, float, float, float, float]]]:
    freqs = list(frequencies)
    components = distribution.component_weights
    if components is None:
        weights = _normalize_distribution(distribution.weights)
        components = tuple((0.0, 0.0, 0.0, weight) for weight in weights)
    component_rows = [(idx, *values) for idx, values in enumerate(components, 1)]
    rainfall_by_freq: dict[str, list[float]] = {}
    for freq in freqs:
        h1 = duration_areal_rainfall["1h"][freq]
        h3 = duration_areal_rainfall["3h"][freq]
        h6 = duration_areal_rainfall["6h"][freq]
        h24 = duration_areal_rainfall["24h"][freq]
        if distribution.main_peak_hours == 1:
            increments = (h1, h3 - h1, h6 - h3, h24 - h6)
        elif distribution.main_peak_hours == 3:
            increments = (0.0, h3, h6 - h3, h24 - h6)
        elif distribution.main_peak_hours == 6:
            increments = (0.0, 0.0, h6, h24 - h6)
        else:
            raise HydrologyError("雨型主雨峰必须为1h、3h或6h")
        if any(value < -1e-6 for value in increments):
            raise HydrologyError(f"{freq} 控制时段面雨量不满足 H₂₄>=H₆>=H₃>=H₁，请检查均值、Cv、alpha、gamma")
        rainfall_by_freq[freq] = [
            _safe_round(sum(max(increments[i], 0.0) * values[i] for i in range(4)), 3)
            for values in components
        ]
    return rainfall_by_freq, component_rows


def instant_rainfall_distribution_table_text(distribution: InstantRainfallDistribution) -> str:
    rows = distribution.component_weights or tuple((0.0, 0.0, 0.0, value) for value in _normalize_distribution(distribution.weights))
    parts: list[str] = []
    parts.append(f"{distribution.zone_name}  主雨峰{distribution.main_peak_hours}h  {distribution.source}")
    if distribution.main_peak_hours == 1:
        table_rows = [["时段", "H₁", "H₃-H₁", "H₆-H₃", "H₂₄-H₆"]]
        for hour, values in enumerate(rows, 1):
            table_rows.append([hour, f"{values[0]:g}", f"{values[1]:g}", f"{values[2]:g}", f"{values[3]:g}"])
    elif distribution.main_peak_hours == 3:
        table_rows = [["时段", "H₃", "H₆-H₃", "H₂₄-H₆"]]
        for hour, values in enumerate(rows, 1):
            table_rows.append([hour, f"{values[1]:g}", f"{values[2]:g}", f"{values[3]:g}"])
    elif distribution.main_peak_hours == 6:
        table_rows = [["时段", "H₆", "H₂₄-H₆"]]
        for hour, values in enumerate(rows, 1):
            table_rows.append([hour, f"{values[2]:g}", f"{values[3]:g}"])
    else:
        table_rows = [["时段", "系数"]]
        for hour, value in enumerate(distribution.weights, 1):
            table_rows.append([hour, f"{value:g}"])
    parts.append(_fixed_table(table_rows, min_width=9))
    return "\n".join(parts)


def instant_runoff_period_from_rainfall(
    rainfall_depths_mm: Iterable[float],
    infiltration_rate_mmh: float,
    delta_t_h: float = 1.0,
) -> dict[str, float]:
    """按图集产流期判定规则由逐时段面雨量和平均入渗率确定 tc。

    产流期内各时段雨量应大于时段入渗量，产流期外不产流；若存在多个
    连续产流段，取净雨量最大的连续段。
    """

    rainfall = [float(value) for value in rainfall_depths_mm]
    if not rainfall:
        raise HydrologyError("逐时段面雨量不能为空")
    delta_t_h = _require_positive("Δt", delta_t_h)
    f_value = _require_positive("平均入渗率f", infiltration_rate_mmh)
    infiltration_depth = f_value * delta_t_h

    segments: list[tuple[int, int, float]] = []
    start: int | None = None
    net_sum = 0.0
    for idx, rain in enumerate(rainfall):
        if rain >= infiltration_depth:
            if start is None:
                start = idx
                net_sum = 0.0
            net_sum += max(0.0, rain - infiltration_depth)
        elif start is not None:
            segments.append((start, idx - 1, net_sum))
            start = None
            net_sum = 0.0
    if start is not None:
        segments.append((start, len(rainfall) - 1, net_sum))
    if not segments:
        raise HydrologyError("没有找到满足产流条件的连续时段，请检查面雨量或平均入渗率f")

    best_start, best_end, best_net = max(segments, key=lambda item: item[2])
    ptc = sum(rainfall[best_start : best_end + 1])
    initial_loss = sum(rainfall[:best_start])
    return {
        "start_hour": float(best_start + 1),
        "end_hour": float(best_end + 1),
        "tc": float(best_end - best_start + 1),
        "infiltration_depth": _safe_round(infiltration_depth, 3),
        "ptc": _safe_round(ptc, 3),
        "net_total": _safe_round(best_net, 3),
        "initial_loss": _safe_round(initial_loss, 3),
    }


def calculate_instant_unit_line_atlas_batch(
    area_km2: float,
    river_slope_permille: float,
    rainfall_distribution: InstantRainfallDistribution,
    rainfall_design: InstantRainfallDesignInput,
    runoff_settings: Iterable[FrequencyRunoffSetting],
    delta_t_h: float = 1.0,
    keep_negative_net_rain: bool = False,
) -> InstantUnitLineBatchResult:
    """按图集/参考 Excel 口径生成净雨过程和各频率洪水过程线。"""

    area = _require_positive("流域面积 F", area_km2)
    _require_positive("主河槽平均坡度 J", river_slope_permille)
    delta_t_h = _require_positive("Δt", delta_t_h)
    settings = list(runoff_settings)
    if not settings:
        raise HydrologyError("至少需要一个频率产流参数")
    frequencies = [item.frequency for item in settings]
    settings_by_freq = {item.frequency: item for item in settings}

    kp_by_duration, point_by_duration, areal_by_duration = instant_duration_design_rainfall(rainfall_design, frequencies)
    rainfall_by_freq, component_rows = _instant_rainfall_process_by_frequency(rainfall_distribution, areal_by_duration, frequencies)

    runoff_periods: dict[str, dict[str, float]] = {}
    for freq in frequencies:
        setting = settings_by_freq[freq]
        if setting.infiltration_rate_mmh is None:
            main_peak = float(rainfall_distribution.main_peak_hours)
            runoff_periods[freq] = instant_runoff_period_auto_infiltration(
                rainfall_by_freq[freq],
                setting.infiltration_table_id,
                delta_t_h,
                min_runoff_duration_h=max(2.0, main_peak),
                fallback_min_duration_h=max(3.0, main_peak),
                max_runoff_duration_h=min(20.0, max(12.0, main_peak * 4.0)),
            )
        else:
            runoff_periods[freq] = instant_runoff_period_from_rainfall(
                rainfall_by_freq[freq],
                setting.infiltration_rate_mmh,
                delta_t_h,
            )

    net_rows: list[NetRainRow] = []
    for hour in range(1, 25):
        rainfall: dict[str, float] = {}
        infiltration: dict[str, float] = {}
        net: dict[str, float] = {}
        for freq in frequencies:
            setting = settings_by_freq[freq]
            period = runoff_periods[freq]
            rain = rainfall_by_freq[freq][hour - 1]
            in_runoff = period["start_hour"] <= hour <= period["end_hour"]
            if in_runoff:
                infil = period["infiltration_depth"]
                net_value = rain - infil
                if not keep_negative_net_rain:
                    net_value = max(0.0, net_value)
            else:
                infil = 0.0
                net_value = 0.0
            rainfall[freq] = _safe_round(rain, 3)
            infiltration[freq] = _safe_round(infil, 3)
            net[freq] = _safe_round(net_value, 3)
        net_rows.append(NetRainRow(hour=hour, rainfall_mm=rainfall, infiltration_mm=infiltration, net_rain_mm=net))

    frequency_inputs: list[FrequencyRunoffInput] = []
    unit_parameter_rows: dict[str, dict[str, float]] = {}
    for setting in settings:
        period = runoff_periods[setting.frequency]
        tc = int(period["tc"])
        net_total = period["net_total"]
        infiltration_rate = float(period.get("infiltration_rate_mmh", setting.infiltration_rate_mmh or 0.0))
        if net_total <= 0:
            raise HydrologyError(f"{setting.frequency} 净雨深合计为0，无法计算瞬时单位线参数；请检查面雨量或入渗率f")
        intensity = net_total / tc
        params = instant_unit_params_from_area_slope_intensity(
            area,
            river_slope_permille,
            intensity,
            net_total,
            setting.unit_param_region,
        )
        frequency_inputs.append(
            FrequencyRunoffInput(
                frequency=setting.frequency,
                ptc_mm=period["ptc"],
                runoff_duration_h=tc,
                infiltration_rate_mmh=infiltration_rate,
                n=params["n"],
                k=params["k"],
                runoff_start_hour=int(period["start_hour"]),
            )
        )
        unit_parameter_rows[setting.frequency] = {
            "ptc": period["ptc"],
            "net_total": _safe_round(net_total, 3),
            "tc": tc,
            "start_hour": period["start_hour"],
            "end_hour": period["end_hour"],
            "initial_loss": period["initial_loss"],
            "infiltration_depth": period["infiltration_depth"],
            "i": _safe_round(intensity, 3),
            "f": _safe_round(infiltration_rate, 3),
            "f_source": (
                INSTANT_INFILTRATION_TABLES[setting.infiltration_table_id]["name"]
                if setting.infiltration_rate_mmh is None
                else "手动输入"
            ),
            "tc_rule": period.get("tc_rule", "手动f判定"),
            **params,
        }

    flood_results: dict[str, FrequencyFloodResult] = {}
    peak_summary: dict[str, tuple[float, float, float]] = {}
    for item in frequency_inputs:
        net_rain = [max(0.0, row.net_rain_mm[item.frequency]) for row in net_rows]
        hydrograph, unit_line, volume = flood_hydrograph_from_unit_line(net_rain, area, item.n, item.k, delta_t_h)
        peak = max(hydrograph, key=lambda point: point.discharge_m3s)
        flood_results[item.frequency] = FrequencyFloodResult(
            frequency=item.frequency,
            unit_line=unit_line,
            hydrograph=hydrograph,
            peak_flow_m3s=peak.discharge_m3s,
            peak_time_h=peak.time_h,
            flood_volume_10k_m3=volume,
        )
        peak_summary[item.frequency] = (peak.discharge_m3s, peak.time_h, volume)

    return InstantUnitLineBatchResult(
        net_rain_table=net_rows,
        flood_results=flood_results,
        peak_summary=peak_summary,
        frequency_inputs={item.frequency: item for item in frequency_inputs},
        kp_values=kp_by_duration["24h"],
        ptc_mean_mm=rainfall_design.h24_mean_mm,
        cv=rainfall_design.h24_cv,
        cs_ratio=rainfall_design.cs_ratio,
        area_km2=area,
        delta_t_h=delta_t_h,
        rainfall_zone_name=rainfall_distribution.zone_name,
        main_peak_hours=rainfall_distribution.main_peak_hours,
        distribution_source=rainfall_distribution.source,
        duration_kp_values=kp_by_duration,
        duration_point_rainfall=point_by_duration,
        duration_areal_rainfall=areal_by_duration,
        rainfall_process_components=component_rows,
        unit_parameter_rows=unit_parameter_rows,
    )


def tieyuan_design_rainfall_values(
    mean_point_24h: float,
    cv: float,
    areal_reduction: float,
    frequencies: Iterable[str] = TIE_FREQUENCIES,
    cs_ratio: float = 3.5,
    shape_correction: float = 1.0,
) -> tuple[dict[str, float], dict[str, float], dict[str, float]]:
    mean_point_24h = _require_positive("24小时点雨量均值", mean_point_24h)
    areal_reduction = _require_positive("点面折减系数", areal_reduction)
    shape_correction = _require_positive("流域形状修正系数", shape_correction)
    kp_values: dict[str, float] = {}
    point_values: dict[str, float] = {}
    areal_values: dict[str, float] = {}
    for freq in frequencies:
        if freq not in TIE_FREQUENCY_PROBABILITIES:
            raise HydrologyError(f"缺少 {freq} 年一遇的频率概率定义")
        kp = pearson_iii_kp(cv, cs_ratio, TIE_FREQUENCY_PROBABILITIES[freq])
        point_rainfall = mean_point_24h * kp
        kp_values[freq] = kp
        point_values[freq] = _safe_round(point_rainfall, 3)
        areal_values[freq] = _safe_round(point_rainfall * areal_reduction * shape_correction, 3)
    return kp_values, point_values, areal_values


def nash_s_curve(time_h: float, n: float, k: float) -> float:
    n = _require_positive("n", n)
    k = _require_positive("k", k)
    if time_h <= 0:
        return 0.0
    return _regularized_gamma_p(n, time_h / k)


def time_step_unit_hydrograph(n: float, k: float, delta_t_h: float = 1.0, max_time_h: float | None = None) -> list[UnitHydrographOrdinate]:
    """由 Nash 瞬时单位线生成时段单位线。

    u(Δt,t) = S(t) - S(t-Δt)，各时段值之和接近 1。
    """

    n = _require_positive("n", n)
    k = _require_positive("k", k)
    delta_t_h = _require_positive("Δt", delta_t_h)
    if max_time_h is None:
        max_time_h = max(delta_t_h, (n * k) + 8 * math.sqrt(n) * k)
    max_time_h = _require_positive("最大计算时长", max_time_h)

    ordinates: list[UnitHydrographOrdinate] = []
    steps = max(1, math.ceil(max_time_h / delta_t_h))
    for step in range(1, steps + 1):
        t = step * delta_t_h
        s_t = nash_s_curve(t, n, k)
        s_prev = nash_s_curve(t - delta_t_h, n, k)
        value = max(0.0, s_t - s_prev)
        ordinates.append(
            UnitHydrographOrdinate(
                time_h=round(t, 6),
                t_over_k=round(t / k, 6),
                s_t=round(s_t, 6),
                s_prev=round(s_prev, 6),
                unit_depth_fraction=round(value, 6),
            )
        )
        if step > 3 and 1 - s_t < 1e-4:
            break
    return ordinates


def flood_hydrograph_from_unit_line(
    net_rainfall_depths_mm: Iterable[float],
    area_km2: float,
    n: float,
    k: float,
    delta_t_h: float = 1.0,
) -> tuple[list[FloodHydrographPoint], list[UnitHydrographOrdinate], float]:
    """净雨过程卷积时段单位线得到地表洪水过程。

    换算系数 C = F / (3.6 * Δt)。输入净雨单位为 mm，流量输出为 m³/s。
    返回：洪水过程、时段单位线、洪水总量(万m3)。
    """

    rainfall = [float(v) for v in net_rainfall_depths_mm]
    if not rainfall:
        raise HydrologyError("净雨过程不能为空")
    if any(v < 0 for v in rainfall):
        raise HydrologyError("净雨过程不能为负")
    area = _require_positive("流域面积 F", area_km2)
    delta_t_h = _require_positive("Δt", delta_t_h)

    max_time = (len(rainfall) + 20) * delta_t_h
    unit_line = time_step_unit_hydrograph(n, k, delta_t_h, max_time)
    unit_values = [item.unit_depth_fraction for item in unit_line]
    conversion = area / (3.6 * delta_t_h)

    discharges: list[FloodHydrographPoint] = []
    total_steps = len(rainfall) + len(unit_values) - 1
    for idx in range(total_steps):
        q_depth = 0.0
        for rain_idx, depth in enumerate(rainfall):
            unit_idx = idx - rain_idx
            if 0 <= unit_idx < len(unit_values):
                q_depth += depth * unit_values[unit_idx]
        discharges.append(FloodHydrographPoint(time_h=round((idx + 1) * delta_t_h, 6), discharge_m3s=round(q_depth * conversion, 3)))

    volume_10k_m3 = sum(rainfall) * area * 0.1
    return discharges, unit_line, round(volume_10k_m3, 3)


def _normalize_distribution(weights: Iterable[float]) -> list[float]:
    values = [float(v) for v in weights]
    if not values:
        raise HydrologyError("雨型系数不能为空")
    if any(v < 0 for v in values):
        raise HydrologyError("雨型系数不能为负")
    total = sum(values)
    if total <= 0:
        raise HydrologyError("雨型系数总和必须大于 0")
    return [v / total for v in values]


def calculate_net_rain_table(
    distribution_weights: Iterable[float],
    frequency_inputs: Iterable[FrequencyRunoffInput],
    delta_t_h: float = 1.0,
    keep_negative_net_rain: bool = False,
) -> list[NetRainRow]:
    """按各频率独立 Ptc/tc/f 生成净雨过程表。

    distribution_weights 是时程分配系数，会自动归一化；每个频率的面雨量为 Ptc * weight。
    入渗量只在该频率产流期内计入：runoff_start_hour 到 runoff_start_hour + tc - 1。
    """

    weights = _normalize_distribution(distribution_weights)
    inputs = list(frequency_inputs)
    if not inputs:
        raise HydrologyError("至少需要一个频率参数")
    delta_t_h = _require_positive("Δt", delta_t_h)

    for item in inputs:
        _require_positive(f"{item.frequency} Ptc", item.ptc_mm)
        _require_positive(f"{item.frequency} tc", item.runoff_duration_h)
        _require_positive(f"{item.frequency} f", item.infiltration_rate_mmh)
        if item.runoff_start_hour < 1:
            raise HydrologyError(f"{item.frequency} 产流起始时段必须大于等于 1")

    rows: list[NetRainRow] = []
    for idx, weight in enumerate(weights, 1):
        rainfall: dict[str, float] = {}
        infiltration: dict[str, float] = {}
        net: dict[str, float] = {}
        for item in inputs:
            rain = item.ptc_mm * weight
            in_runoff = item.runoff_start_hour <= idx < item.runoff_start_hour + int(item.runoff_duration_h)
            infil = item.infiltration_rate_mmh * delta_t_h if in_runoff else 0.0
            if in_runoff and rain + 1e-9 < infil:
                raise HydrologyError(
                    f"{item.frequency} 频率第 {idx} 时段面雨量 {rain:.3f}mm 小于入渗量 {infil:.3f}mm，"
                    "请重新确定连续产流期tc或平均入渗率f"
                )
            net_value = rain - infil
            if not keep_negative_net_rain:
                net_value = max(0.0, net_value)
            rainfall[item.frequency] = round(rain, 3)
            infiltration[item.frequency] = round(infil, 3)
            net[item.frequency] = round(net_value, 3)
        rows.append(NetRainRow(hour=idx, rainfall_mm=rainfall, infiltration_mm=infiltration, net_rain_mm=net))
    return rows


def calculate_instant_unit_line_batch(
    area_km2: float,
    distribution_weights: Iterable[float],
    frequency_inputs: Iterable[FrequencyRunoffInput],
    delta_t_h: float = 1.0,
    keep_negative_net_rain: bool = False,
    rainfall_distribution: InstantRainfallDistribution | None = None,
    kp_values: dict[str, float] | None = None,
    ptc_mean_mm: float | None = None,
    cv: float | None = None,
    cs_ratio: float | None = None,
) -> InstantUnitLineBatchResult:
    area = _require_positive("流域面积 F", area_km2)
    inputs = list(frequency_inputs)
    net_rows = calculate_net_rain_table(distribution_weights, inputs, delta_t_h, keep_negative_net_rain)
    results: dict[str, FrequencyFloodResult] = {}
    summary: dict[str, tuple[float, float, float]] = {}

    for item in inputs:
        net_rain = [max(0.0, row.net_rain_mm[item.frequency]) for row in net_rows]
        hydrograph, unit_line, volume = flood_hydrograph_from_unit_line(net_rain, area, item.n, item.k, delta_t_h)
        peak = max(hydrograph, key=lambda point: point.discharge_m3s)
        results[item.frequency] = FrequencyFloodResult(
            frequency=item.frequency,
            unit_line=unit_line,
            hydrograph=hydrograph,
            peak_flow_m3s=peak.discharge_m3s,
            peak_time_h=peak.time_h,
            flood_volume_10k_m3=volume,
        )
        summary[item.frequency] = (peak.discharge_m3s, peak.time_h, volume)

    return InstantUnitLineBatchResult(
        net_rain_table=net_rows,
        flood_results=results,
        peak_summary=summary,
        frequency_inputs={item.frequency: item for item in inputs},
        kp_values=kp_values,
        ptc_mean_mm=ptc_mean_mm,
        cv=cv,
        cs_ratio=cs_ratio,
        area_km2=area,
        delta_t_h=delta_t_h,
        rainfall_zone_name=rainfall_distribution.zone_name if rainfall_distribution else None,
        main_peak_hours=rainfall_distribution.main_peak_hours if rainfall_distribution else None,
        distribution_source=rainfall_distribution.source if rainfall_distribution else None,
    )


RATIONAL_FORMULA_REGIONS = {
    "huangtu": {
        "name": "黄土区",
        "formula": "m = 2.5 h^-0.623 θ^0.637",
        "theta_formula": "θ = L / (J^(1/3) F^(1/4))",
        "critical_h": 35.0,
    },
    "liupanshan_forest": {
        "name": "六盘山土石山林区",
        "formula": "m = 0.195 θ^0.397",
        "theta_formula": "θ = L / J^(1/3)",
        "critical_h": 20.0,
    },
    "liupanshan_general": {
        "name": "六盘山土石山植被一般区",
        "formula": "m = 0.10 θ^0.384",
        "theta_formula": "θ = L / J^(1/3)",
        "critical_h": 30.0,
    },
}

RATIONAL_FORMULA_REGION_LABELS = tuple(
    (region_id, data["name"]) for region_id, data in RATIONAL_FORMULA_REGIONS.items()
)


RATIONAL_SHAPE_CHART_POINTS = {
    "huangtu": {
        "source": "附图58：甘肃省黄土区推理公式概化过程线形状特征系数综合关系图",
        # 附图58左半图先由 θ 查 Ka、Kw，右半图再由 Kw 查 KT。
        # θ=25.64 的控制点采用图集表19公布值 Ka=0.028、Kw=0.303、KT=0.218。
        "theta": (5.0, 10.0, 20.0, 25.64, 30.0, 40.0, 50.0, 60.0),
        "ka": (0.0, 0.004, 0.018, 0.028, 0.040, 0.075, 0.125, 0.190),
        "kw": (0.10, 0.14, 0.24, 0.303, 0.35, 0.43, 0.50, 0.56),
        "kw_for_kt": (0.10, 0.20, 0.303, 0.40, 0.50, 0.56),
        "kt": (0.10, 0.158, 0.218, 0.292, 0.375, 0.43),
    },
    "liupanshan_forest": {
        "source": "附图59：甘肃省六盘山区推理公式概化过程线形状特征系数综合关系图",
        "theta": (0.10, 0.20, 0.40, 0.60, 0.80, 1.00),
        "ka": (0.17, 0.23, 0.32, 0.41, 0.49, 0.56),
        "kw": (0.23, 0.28, 0.36, 0.43, 0.49, 0.54),
        "kt": (0.39, 0.36, 0.31, 0.27, 0.23, 0.20),
    },
    "liupanshan_general": {
        "source": "附图59：甘肃省六盘山区推理公式概化过程线形状特征系数综合关系图",
        "theta": (0.10, 0.20, 0.40, 0.60, 0.80, 1.00),
        "ka": (0.14, 0.20, 0.29, 0.37, 0.45, 0.52),
        "kw": (0.21, 0.26, 0.34, 0.41, 0.47, 0.52),
        "kt": (0.40, 0.37, 0.32, 0.28, 0.24, 0.21),
    },
}

RATIONAL_BASEFLOW_CURVE_POINTS = {
    "hexi_shishan": {
        "source": "附图66曲线①：河西石山森林草原区",
        "area": (30, 50, 100, 300, 500, 1000, 1500),
        "q0": (5.5, 6.5, 8.2, 10.5, 12.0, 15.0, 18.0),
    },
    "longnan_shishan": {
        "source": "附图66曲线②：陇南石山林区",
        "area": (30, 50, 100, 300, 500, 1000, 1500),
        "q0": (3.2, 4.2, 5.6, 8.5, 10.5, 13.5, 16.0),
    },
    "huangtu": {
        "source": "附图66曲线③：黄土区",
        "area": (30, 50, 100, 300, 500, 1000, 1500),
        "q0": (0.5, 0.65, 0.85, 1.2, 1.45, 1.75, 2.0),
    },
    "huangtu_shanlin": {
        "source": "附图66曲线④：黄土山（林）区",
        "area": (30, 50, 100, 300, 500, 1000, 1500),
        "q0": (1.0, 1.25, 1.55, 2.0, 2.25, 2.45, 2.5),
    },
    "liupanshan_tushishan": {
        "source": "附图66曲线⑤：六盘山土石山林区",
        "area": (30, 50, 100, 300, 500, 1000, 1500),
        "q0": (0.55, 0.70, 0.95, 1.55, 2.3, 3.8, 5.0),
    },
    "liupanshan_general": {
        "source": "附图66曲线⑥：六盘山土石山植被一般区",
        "area": (30, 50, 100, 300, 500, 1000, 1500),
        "q0": (0.50, 0.65, 0.90, 1.55, 2.4, 3.8, 5.0),
    },
}


def rational_formula_m_value(
    region_id: str,
    area_km2: float,
    river_length_km: float,
    river_slope_permille: float,
    net_rain_depth_mm: float,
) -> dict[str, float | str]:
    """按图集表10计算推理公式汇流参数 m。

    图集中 J 按河道比降参与 J^(1/3)；界面按常用的 ‰ 输入，计算时折算为小数。
    黄土区 m 与净雨深 h 有关，超过临界净雨深后按临界值控制。
    """

    if region_id not in RATIONAL_FORMULA_REGIONS:
        raise HydrologyError("必须选择有效的推理公式汇流分区")
    area = _require_positive("流域面积 F", area_km2)
    length = _require_positive("主河道长度 L", river_length_km)
    slope = _require_positive("主河槽比降 J", river_slope_permille) / 1000.0
    h = _require_positive("净雨深 h", net_rain_depth_mm)
    region = RATIONAL_FORMULA_REGIONS[region_id]
    critical_h = float(region["critical_h"])
    if region_id == "huangtu":
        theta = length / ((slope ** (1 / 3)) * (area ** 0.25))
        h_used = min(h, critical_h)
        # 表10扫描件指数看似 0.837，但图集完整算例 F=100 km²、θ=25.64、h=35 mm
        # 明确给出 m=2.155；只有 0.637 能复现该算例，故以算例闭合值为准。
        m_value = 2.5 * (h_used ** -0.623) * (theta ** 0.637)
    elif region_id == "liupanshan_forest":
        theta = length / (slope ** (1 / 3))
        h_used = min(h, critical_h)
        m_value = 0.195 * (theta ** 0.397)
    else:
        theta = length / (slope ** (1 / 3))
        h_used = min(h, critical_h)
        m_value = 0.10 * (theta ** 0.384)
    return {
        "region_name": str(region["name"]),
        "formula": str(region["formula"]),
        "theta_formula": str(region["theta_formula"]),
        "critical_h": critical_h,
        "h_used": _safe_round(h_used, 3),
        "theta": _safe_round(theta, 3),
        "m": _safe_round(m_value, 3),
        "j_decimal": _safe_round(slope, 6),
    }


def rational_shape_params_from_chart(
    region_id: str,
    area_km2: float,
    river_length_km: float,
    river_slope_permille: float,
) -> dict[str, float | str]:
    """按附图58/59由流域特征参数自动估算 Ka、Kw、KT。

    附图为扫描关系图，当前采用图件数字化控制点线性插值；结果在计算书中标注来源，
    可作为自动查图值使用。
    """

    area = _require_positive("流域面积 F", area_km2)
    length = _require_positive("主河道长度 L", river_length_km)
    slope = _require_positive("主河槽比降 J", river_slope_permille) / 1000.0
    if region_id == "huangtu":
        theta = length / ((slope ** (1 / 3)) * (area ** 0.25))
    elif region_id in {"liupanshan_forest", "liupanshan_general"}:
        theta = length / (slope ** (1 / 3))
    else:
        raise HydrologyError("必须选择有效的推理公式汇流分区")
    chart = RATIONAL_SHAPE_CHART_POINTS[region_id]
    points = list(chart["theta"])
    ka = _interpolate_1d(zip(points, chart["ka"]), theta)
    kw = _interpolate_1d(zip(points, chart["kw"]), theta)
    if region_id == "huangtu":
        kt = _interpolate_1d(zip(chart["kw_for_kt"], chart["kt"]), kw)
    else:
        kt = _interpolate_1d(zip(points, chart["kt"]), theta)
    return {
        "theta_shape": _safe_round(theta, 3),
        "ka": _safe_round(max(0.02, min(0.95, ka)), 3),
        "kw": _safe_round(max(0.02, min(0.49, kw)), 3),
        "kt": _safe_round(max(0.05, min(0.49, kt)), 3),
        "shape_source": str(chart["source"]),
    }


def rational_baseflow_from_chart(area_km2: float, table_id: str) -> dict[str, float | str]:
    """按附图66由流域面积和水文分区自动估算基流量。"""

    area = _require_positive("流域面积 F", area_km2)
    curve_id = "liupanshan_general" if table_id == "liupanshan_general" else table_id
    if curve_id not in RATIONAL_BASEFLOW_CURVE_POINTS:
        curve_id = "huangtu"
    curve = RATIONAL_BASEFLOW_CURVE_POINTS[curve_id]
    log_points = [(math.log10(float(x)), float(y)) for x, y in zip(curve["area"], curve["q0"])]
    q0 = _interpolate_1d(log_points, math.log10(area))
    return {
        "baseflow": _safe_round(max(0.0, q0), 3),
        "baseflow_source": str(curve["source"]),
    }


def _integral_piecewise_constant(values: list[float], delta_t_h: float, time_h: float) -> float:
    if time_h <= 0:
        return 0.0
    total_duration = len(values) * delta_t_h
    if time_h >= total_duration:
        return sum(values)
    full_steps = int(time_h // delta_t_h)
    remainder = time_h - full_steps * delta_t_h
    total = sum(values[:full_steps])
    if full_steps < len(values) and remainder > 1e-9:
        total += values[full_steps] * remainder / delta_t_h
    return total


def _max_window_net_rain(
    net_rain_mm: list[float],
    delta_t_h: float,
    duration_h: float,
) -> tuple[float, float, float]:
    duration = min(_require_positive("最大平均净雨历时", duration_h), len(net_rain_mm) * delta_t_h)
    total_duration = len(net_rain_mm) * delta_t_h
    latest_start = max(0.0, total_duration - duration)
    boundaries = [idx * delta_t_h for idx in range(len(net_rain_mm) + 1)]
    starts = {0.0, latest_start}
    for boundary in boundaries:
        if 0.0 <= boundary <= latest_start:
            starts.add(boundary)
        candidate = boundary - duration
        if 0.0 <= candidate <= latest_start:
            starts.add(candidate)
    best_start = 0.0
    best_depth = -1.0
    for start in sorted(starts):
        end = start + duration
        depth = _integral_piecewise_constant(net_rain_mm, delta_t_h, end) - _integral_piecewise_constant(net_rain_mm, delta_t_h, start)
        if depth > best_depth:
            best_depth = depth
            best_start = start
    best_depth = max(0.0, best_depth)
    return _safe_round(best_depth, 6), _safe_round(best_depth / duration, 6), _safe_round(best_start, 6)


def _main_peak_net_rain_trials(
    net_rain_mm: list[float],
    delta_t_h: float,
    max_duration_h: float,
) -> tuple[tuple[float, float, float, float], ...]:
    """按图集图解法，从最大净雨强时段向前后相邻时段连续累积 h_t。"""

    delta_t = _require_positive("Δt", delta_t_h)
    max_duration = _require_positive("产流历时 tc", max_duration_h)
    values = [max(0.0, float(value)) for value in net_rain_mm]
    if not values or max(values) <= 1e-9:
        return ()

    peak_index = max(range(len(values)), key=lambda idx: values[idx])
    left = peak_index
    right = peak_index
    depth = values[peak_index]
    rows: list[tuple[float, float, float, float]] = []
    max_steps = max(1, min(len(values), int(math.ceil(max_duration / delta_t - 1e-9))))

    for step in range(1, max_steps + 1):
        duration = step * delta_t
        if duration <= max_duration + 1e-9:
            intensity = depth / duration if duration > 0 else 0.0
            rows.append((
                _safe_round(duration, 6),
                _safe_round(depth, 6),
                _safe_round(intensity, 6),
                _safe_round(left * delta_t, 6),
            ))
        if step == max_steps:
            break

        left_value = values[left - 1] if left > 0 else -1.0
        right_value = values[right + 1] if right + 1 < len(values) else -1.0
        if left_value < 0 and right_value < 0:
            break
        if right_value > left_value:
            right += 1
            depth += values[right]
        else:
            left -= 1
            depth += values[left]

    return tuple(rows)


def _interp_between(value_a: float, value_b: float, ratio: float) -> float:
    return value_a + (value_b - value_a) * ratio


def _solve_rational_tau_and_q(
    net_rain_mm: list[float],
    area_km2: float,
    river_length_km: float,
    river_slope_permille: float,
    m_value: float,
    runoff_duration_h: float,
    delta_t_h: float,
) -> dict[str, float | str | tuple[tuple[float, float, float, float, float, float], ...]]:
    area = _require_positive("流域面积 F", area_km2)
    length = _require_positive("主河道长度 L", river_length_km)
    slope_decimal = _require_positive("主河槽比降 J", river_slope_permille) / 1000.0
    m = _require_positive("汇流参数 m", m_value)
    tc = _require_positive("产流历时 tc", runoff_duration_h)

    trial_base = _main_peak_net_rain_trials(net_rain_mm, delta_t_h, tc)
    trial_rows: list[tuple[float, float, float, float, float, float]] = []
    for duration, depth, intensity, start in trial_base:
        q_trial = 0.278 * area * intensity
        tau_trial = 0.0
        if q_trial > 1e-9:
            tau_trial = 0.278 * length / (m * (slope_decimal ** (1 / 3)) * (q_trial ** 0.25))
        trial_rows.append((
            _safe_round(duration, 6),
            _safe_round(depth, 6),
            _safe_round(intensity, 6),
            _safe_round(q_trial, 6),
            _safe_round(tau_trial, 6),
            _safe_round(start, 6),
        ))

    if not trial_rows or max(row[3] for row in trial_rows) <= 1e-9:
        return {
            "tau": _safe_round(delta_t_h, 3),
            "q_surface": 0.0,
            "duration": _safe_round(delta_t_h, 3),
            "depth": 0.0,
            "intensity": 0.0,
            "start": 0.0,
            "condition": "净雨为0，不能形成地表洪峰",
            "trial_rows": tuple(trial_rows),
        }

    chosen = trial_rows[0]
    solved_by_crossing = False
    solved_within_first_step = False
    solved_by_full_confluence = False
    first = trial_rows[0]
    if 1e-9 < first[4] < first[0]:
        chosen = (
            first[4],
            first[2] * first[4],
            first[2],
            first[3],
            first[4],
            first[5],
        )
        solved_by_crossing = True
        solved_within_first_step = True
    else:
        for row_a, row_b in zip(trial_rows, trial_rows[1:]):
            diff_a = row_a[4] - row_a[0]
            diff_b = row_b[4] - row_b[0]
            if abs(diff_a) < 1e-9:
                chosen = row_a
                solved_by_crossing = True
                break
            if diff_a * diff_b <= 0:
                ratio = -diff_a / (diff_b - diff_a) if abs(diff_b - diff_a) > 1e-12 else 0.0
                chosen = (
                    _interp_between(row_a[0], row_b[0], ratio),
                    _interp_between(row_a[1], row_b[1], ratio),
                    _interp_between(row_a[2], row_b[2], ratio),
                    _interp_between(row_a[3], row_b[3], ratio),
                    _interp_between(row_a[4], row_b[4], ratio),
                    row_a[5] if ratio < 0.5 else row_b[5],
                )
                solved_by_crossing = True
                break

    # 图集说明：当试算历时 t>tc 后，累计净雨深 h_t 取总净雨深 h 不再增加。
    # 若在 tc 内尚未相交，应继续求全面汇流分支，而不能停在 tc 的最近试算点。
    if not solved_by_crossing and trial_rows[-1][4] > trial_rows[-1][0]:
        total_depth = sum(max(0.0, float(value)) for value in net_rain_mm)
        q_time_constant = 0.278 * area * total_depth
        travel_constant = 0.278 * length / (m * (slope_decimal ** (1 / 3)))
        if q_time_constant > 1e-9:
            full_tau = (travel_constant / (q_time_constant ** 0.25)) ** (4 / 3)
            if full_tau >= tc - 1e-9:
                q_full = q_time_constant / full_tau
                chosen = (
                    full_tau,
                    total_depth,
                    total_depth / full_tau,
                    q_full,
                    full_tau,
                    trial_rows[-1][5],
                )
                trial_rows.append(chosen)
                solved_by_crossing = True
                solved_by_full_confluence = True

    if not solved_by_crossing:
        chosen = min(trial_rows, key=lambda row: abs(row[4] - row[0]))

    duration, depth, intensity, q_surface, tau_trial, start = chosen
    tau = tau_trial if tau_trial > 1e-9 else duration
    if solved_by_full_confluence:
        condition = "图解法：τ大于产流历时tc，按h_t等于总净雨深的全面汇流分支求交点"
    elif solved_within_first_step:
        condition = "图解法：τ小于1个计算时段，采用最大1h净雨强并由公式反算τ、Qm"
    elif solved_by_crossing:
        condition = "图解法：由主净雨峰累积h_t，按Q_t-t与Q_t-τ_t交点求τ、Qm"
    else:
        condition = "图解法：试算点未跨越，取τ_t与t差值绝对值最小点，需人工复核"
    condition = f"{condition}；{'tc>τ' if tc > tau else 'tc≤τ'}"
    return {
        "tau": _safe_round(tau, 3),
        "q_surface": _safe_round(q_surface, 3),
        "duration": _safe_round(duration, 3),
        "depth": _safe_round(depth, 3),
        "intensity": _safe_round(intensity, 3),
        "start": _safe_round(start, 3),
        "condition": condition,
        "trial_rows": tuple(
            (
                _safe_round(row[0], 3),
                _safe_round(row[1], 3),
                _safe_round(row[2], 3),
                _safe_round(row[3], 3),
                _safe_round(row[4], 3),
                _safe_round(row[5], 3),
            )
            for row in trial_rows
        ),
    }


def _rational_confluence_curve_points(
    net_rain_mm: list[float],
    area_km2: float,
    river_length_km: float,
    river_slope_permille: float,
    m_value: float,
    runoff_duration_h: float,
    solved_tau_h: float,
    delta_t_h: float,
    solved_q_surface_m3s: float | None = None,
) -> tuple[tuple[float, float, float, float], ...]:
    area = _require_positive("流域面积 F", area_km2)
    length = _require_positive("主河道长度 L", river_length_km)
    slope_decimal = _require_positive("主河槽比降 J", river_slope_permille) / 1000.0
    m = _require_positive("汇流参数 m", m_value)
    tc = _require_positive("产流历时 tc", runoff_duration_h)
    solved_tau = _require_positive("图解交点汇流时间τ", solved_tau_h)
    travel_constant = 0.278 * length / (m * (slope_decimal ** (1 / 3)))
    q_surface = max(0.0, float(solved_q_surface_m3s)) if solved_q_surface_m3s is not None else (travel_constant / solved_tau) ** 4
    if q_surface <= 1e-9:
        return ()

    total_net_depth = sum(max(0.0, float(value)) for value in net_rain_mm)
    sample_count = 96
    if solved_tau > tc:
        start_time = max(tc, solved_tau * 0.70)
    else:
        start_time = max(0.001, solved_tau * 0.70)
    if start_time >= solved_tau:
        start_time = max(0.000001, solved_tau * 0.25)

    raw_rain: list[tuple[float, float]] = []
    for idx in range(sample_count):
        ratio = idx / (sample_count - 1)
        duration = start_time + (solved_tau - start_time) * ratio
        duration = max(1e-9, duration)
        if duration <= tc + 1e-9:
            _depth, intensity, _start = _max_window_net_rain(net_rain_mm, delta_t_h, duration)
        else:
            intensity = total_net_depth / duration
        q_rain = 0.278 * area * intensity
        if math.isfinite(q_rain) and q_rain >= 0:
            raw_rain.append((duration, q_rain))
    if not raw_rain:
        return ()

    raw_at_tau = raw_rain[-1][1]
    scale = q_surface / raw_at_tau if raw_at_tau > 1e-9 else 1.0
    points: list[tuple[float, float, float, float]] = []
    for duration, q_rain_raw in raw_rain[:-1]:
        tau_axis = duration
        q_rain = max(0.0, q_rain_raw * scale)
        q_travel = q_surface * (solved_tau / tau_axis) ** 4
        if all(math.isfinite(value) for value in (duration, q_rain, tau_axis, q_travel)):
            points.append((
                _safe_round(duration, 4),
                _safe_round(q_rain, 3),
                _safe_round(tau_axis, 4),
                _safe_round(q_travel, 3),
            ))
    points.append((
        _safe_round(solved_tau, 4),
        _safe_round(q_surface, 3),
        _safe_round(solved_tau, 4),
        _safe_round(q_surface, 3),
    ))
    return tuple(points)


def _linear_value(points: list[tuple[float, float]], time_h: float) -> float:
    if not points or time_h < points[0][0] or time_h > points[-1][0]:
        return 0.0
    for (t0, q0), (t1, q1) in zip(points, points[1:]):
        if t0 <= time_h <= t1:
            if abs(t1 - t0) < 1e-9:
                return q1
            ratio = (time_h - t0) / (t1 - t0)
            return q0 + (q1 - q0) * ratio
    return 0.0


def _rational_formula_hydrograph(
    net_rain_mm: list[float],
    area_km2: float,
    tau_h: float,
    surface_peak_m3s: float,
    main_start_h: float,
    main_duration_h: float,
    delta_t_h: float,
    shape_ka: float,
    shape_kw: float,
    shape_kt: float,
    baseflow_m3s: float,
) -> tuple[
    list[FloodHydrographPoint],
    tuple[tuple[str, float, float], ...],
    tuple[tuple[tuple[float, float], ...], ...],
]:
    tau = _require_positive("汇流时间τ", tau_h)
    area = _require_positive("流域面积 F", area_km2)
    delta_t = _require_positive("Δt", delta_t_h)
    ka = max(0.0, min(float(shape_ka), 0.95))
    kw = max(0.01, min(float(shape_kw), 0.95))
    kt = max(0.05, min(float(shape_kt), 0.49))
    baseflow = max(0.0, float(baseflow_m3s))

    components: list[list[tuple[float, float]]] = []
    end_times: list[float] = [len(net_rain_mm) * delta_t]
    main_end = main_start_h + max(delta_t / 10.0, main_duration_h)

    denominator = (1 / kt) - 2
    kb = max(0.0, min(0.95, (1 - 2 * kw) / denominator)) if denominator > 1e-9 else 0.0
    ta = max(0.0, min(tau, (1 - 2 * kw + ka) * tau))
    tb = max(tau, min(tau / kt, (2 - kb) * tau))
    total_width = max(tb, tau / kt)
    main_points = [
        (main_start_h, 0.0),
        (main_start_h + ta, ka * surface_peak_m3s),
        (main_start_h + tau, surface_peak_m3s),
        (main_start_h + tb, kb * surface_peak_m3s),
        (main_start_h + total_width, 0.0),
    ]
    markers = [
        ("A", main_start_h, baseflow),
        ("Qa", main_start_h + ta, ka * surface_peak_m3s + baseflow),
        ("Qm", main_start_h + tau, surface_peak_m3s + baseflow),
        ("Qb", main_start_h + tb, kb * surface_peak_m3s + baseflow),
        ("B", main_start_h + total_width, baseflow),
    ]
    components.append(main_points)
    end_times.append(main_points[-1][0])

    residual_segments: list[tuple[float, float, float]] = []
    positive_indices = [idx for idx, depth in enumerate(net_rain_mm) if depth > 1e-9]
    if positive_indices:
        runoff_start = positive_indices[0] * delta_t
        runoff_end = (positive_indices[-1] + 1) * delta_t
        pre_end = min(main_start_h, runoff_end)
        if pre_end > runoff_start + 1e-9:
            pre_depth = (
                _integral_piecewise_constant(net_rain_mm, delta_t, pre_end)
                - _integral_piecewise_constant(net_rain_mm, delta_t, runoff_start)
            )
            if pre_depth > 1e-9:
                residual_segments.append((runoff_start, pre_end - runoff_start, pre_depth))
        post_start = max(main_end, runoff_start)
        if runoff_end > post_start + 1e-9:
            post_depth = (
                _integral_piecewise_constant(net_rain_mm, delta_t, runoff_end)
                - _integral_piecewise_constant(net_rain_mm, delta_t, post_start)
            )
            if post_depth > 1e-9:
                residual_segments.append((post_start, runoff_end - post_start, post_depth))

    for start, duration, depth_for_triangle in residual_segments:
        if duration <= 0 or depth_for_triangle <= 1e-9:
            continue
        end = start + duration
        q_peak = 0.556 * depth_for_triangle * area / (duration + tau)
        peak_time = end if end <= main_start_h else start + min(duration, tau)
        triangle = [
            (start, 0.0),
            (peak_time, q_peak),
            (start + duration + tau, 0.0),
        ]
        components.append(triangle)
        end_times.append(triangle[-1][0])

    candidate_times = {0.0}
    max_time = max(end_times)
    steps = int(math.ceil(max_time / delta_t))
    for step in range(steps + 1):
        candidate_times.add(round(step * delta_t, 6))
    for points in components:
        for time_h, _value in points:
            candidate_times.add(round(time_h, 6))

    hydrograph: list[FloodHydrographPoint] = []
    for time_h in sorted(candidate_times):
        q_surface = sum(_linear_value(points, time_h) for points in components)
        hydrograph.append(FloodHydrographPoint(time_h=_safe_round(time_h, 3), discharge_m3s=_safe_round(q_surface + baseflow, 3)))
    peak = max(hydrograph, key=lambda point: point.discharge_m3s)
    markers.append(("Qp", peak.time_h, peak.discharge_m3s))
    return hydrograph, tuple((label, _safe_round(time_h, 3), _safe_round(flow, 3)) for label, time_h, flow in markers), tuple(
        tuple((_safe_round(time_h, 3), _safe_round(flow, 3)) for time_h, flow in points)
        for points in components
    )


def calculate_rational_formula_atlas_batch(
    area_km2: float,
    river_length_km: float,
    river_slope_permille: float,
    rainfall_distribution: InstantRainfallDistribution,
    rainfall_design: InstantRainfallDesignInput,
    runoff_settings: Iterable[FrequencyRunoffSetting],
    rational_region_id: str = "huangtu",
    delta_t_h: float = 1.0,
    keep_negative_net_rain: bool = False,
    baseflow_m3s: float | None = None,
    shape_ka: float | None = None,
    shape_kw: float | None = None,
    shape_kt: float | None = None,
) -> RationalFormulaBatchResult:
    """按《甘肃省暴雨洪水图集》推理公式法计算洪峰和汇流过程线。

    暴雨、附表10雨型分配、附表12平均入渗率和净雨过程与瞬时单位线法共用；
    推理公式部分采用表10 m 综合关系、公式(14)-(20)生成洪峰与过程线。
    """

    area = _require_positive("流域面积 F", area_km2)
    length = _require_positive("主河道长度 L", river_length_km)
    _require_positive("主河槽比降 J", river_slope_permille)
    if rational_region_id not in RATIONAL_FORMULA_REGIONS:
        raise HydrologyError("必须选择有效的推理公式汇流分区")
    settings_list = list(runoff_settings)
    if not settings_list:
        raise HydrologyError("至少需要一个频率产流参数")
    first_table_id = settings_list[0].infiltration_table_id
    if rational_region_id == "liupanshan_general":
        baseflow_curve_id = "liupanshan_general"
    else:
        baseflow_curve_id = first_table_id
    auto_shape = rational_shape_params_from_chart(rational_region_id, area, length, river_slope_permille)
    auto_baseflow = rational_baseflow_from_chart(area, baseflow_curve_id)
    shape_ka_value = float(auto_shape["ka"] if shape_ka is None else shape_ka)
    shape_kw_value = float(auto_shape["kw"] if shape_kw is None else shape_kw)
    shape_kt_value = float(auto_shape["kt"] if shape_kt is None else shape_kt)
    baseflow_value = float(auto_baseflow["baseflow"] if baseflow_m3s is None else baseflow_m3s)
    shared = calculate_instant_unit_line_atlas_batch(
        area_km2=area,
        river_slope_permille=river_slope_permille,
        rainfall_distribution=rainfall_distribution,
        rainfall_design=rainfall_design,
        runoff_settings=settings_list,
        delta_t_h=delta_t_h,
        keep_negative_net_rain=keep_negative_net_rain,
    )

    frequency_results: dict[str, RationalFormulaFrequencyResult] = {}
    parameter_rows: dict[str, dict[str, float | str]] = {}
    for freq, item in (shared.frequency_inputs or {}).items():
        net_rain = [max(0.0, row.net_rain_mm[freq]) for row in shared.net_rain_table]
        net_total = sum(net_rain)
        m_info = rational_formula_m_value(rational_region_id, area, length, river_slope_permille, net_total)
        solved = _solve_rational_tau_and_q(
            net_rain,
            area,
            length,
            river_slope_permille,
            float(m_info["m"]),
            item.runoff_duration_h,
            delta_t_h,
        )
        confluence_curve = tuple(solved.get("trial_rows", ()))  # t, h_t, h_t/t, Q_t, tau_t, start
        confluence_curve_points = _rational_confluence_curve_points(
            net_rain,
            area,
            length,
            river_slope_permille,
            float(m_info["m"]),
            item.runoff_duration_h,
            float(solved["tau"]),
            delta_t_h,
            float(solved["q_surface"]),
        )
        hydrograph, curve_markers, hydrograph_components = _rational_formula_hydrograph(
            net_rain,
            area,
            float(solved["tau"]),
            float(solved["q_surface"]),
            float(solved["start"]),
            float(solved["duration"]),
            delta_t_h,
            shape_ka_value,
            shape_kw_value,
            shape_kt_value,
            baseflow_value,
        )
        peak = max(hydrograph, key=lambda point: point.discharge_m3s)
        volume = _safe_round(net_total * area * 0.1, 3)
        frequency_results[freq] = RationalFormulaFrequencyResult(
            frequency=freq,
            hydrograph=hydrograph,
            peak_flow_m3s=peak.discharge_m3s,
            peak_time_h=peak.time_h,
            surface_peak_m3s=float(solved["q_surface"]),
            tau_h=float(solved["tau"]),
            m=float(m_info["m"]),
            theta=float(m_info["theta"]),
            runoff_condition=str(solved["condition"]),
            max_net_duration_h=float(solved["duration"]),
            max_net_depth_mm=float(solved["depth"]),
            max_net_intensity_mmh=float(solved["intensity"]),
            main_start_h=float(solved["start"]),
            flood_volume_10k_m3=volume,
            curve_markers=curve_markers,
            confluence_curve=confluence_curve,
            confluence_curve_points=confluence_curve_points,
            hydrograph_components=hydrograph_components,
        )
        unit_row = (shared.unit_parameter_rows or {}).get(freq, {})
        parameter_rows[freq] = {
            **unit_row,
            **m_info,
            "tau": solved["tau"],
            "q_surface": solved["q_surface"],
            "q_peak": peak.discharge_m3s,
            "peak_time": peak.time_h,
            "net_total": _safe_round(net_total, 3),
            "i": _safe_round(net_total / item.runoff_duration_h, 3),
            "max_net_duration": solved["duration"],
            "max_net_depth": solved["depth"],
            "max_net_intensity": solved["intensity"],
            "main_start": solved["start"],
            "condition": solved["condition"],
            "trial_count": len(confluence_curve),
            "baseflow": _safe_round(baseflow_value, 3),
            "ka": _safe_round(shape_ka_value, 3),
            "kw": _safe_round(shape_kw_value, 3),
            "kt": _safe_round(shape_kt_value, 3),
            "shape_source": auto_shape["shape_source"] if shape_ka is None or shape_kw is None or shape_kt is None else "手动输入",
            "baseflow_source": auto_baseflow["baseflow_source"] if baseflow_m3s is None else "手动输入",
        }

    return RationalFormulaBatchResult(
        shared=shared,
        rational_region_id=rational_region_id,
        rational_region_name=str(RATIONAL_FORMULA_REGIONS[rational_region_id]["name"]),
        river_length_km=length,
        river_slope_permille=river_slope_permille,
        baseflow_m3s=_safe_round(baseflow_value, 3),
        shape_ka=_safe_round(shape_ka_value, 3),
        shape_kw=_safe_round(shape_kw_value, 3),
        shape_kt=_safe_round(shape_kt_value, 3),
        shape_theta=float(auto_shape["theta_shape"]),
        shape_source=str(auto_shape["shape_source"] if shape_ka is None or shape_kw is None or shape_kt is None else "手动输入"),
        baseflow_source=str(auto_baseflow["baseflow_source"] if baseflow_m3s is None else "手动输入"),
        frequency_results=frequency_results,
        parameter_rows=parameter_rows,
    )


def rational_formula_result_to_text(result: RationalFormulaBatchResult, frequencies: Iterable[str] = DESIGN_FREQUENCIES_7) -> str:
    freqs = [freq for freq in frequencies if freq in result.frequency_results]
    freq_labels = {freq: f"{freq}年" if freq in TIE_FREQUENCIES else freq for freq in freqs}
    shared = result.shared
    parts: list[str] = []

    parts.append("一、推理公式法输入参数与图集依据")
    rows = [["项目", "数值", "说明"]]
    rows.append(["流域面积F(km²)", f"{shared.area_km2:g}" if shared.area_km2 is not None else "", "输入"])
    rows.append(["主河道长度L(km)", f"{result.river_length_km:g}", "输入"])
    rows.append(["主河槽比降J(‰)", f"{result.river_slope_permille:g}", "输入；软件按J/1000代入J^(1/3)"])
    rows.append(["雨型分区", shared.rainfall_zone_name or "", "附表10"])
    rows.append(["主雨峰", f"{shared.main_peak_hours:g}h" if shared.main_peak_hours else "", "F<100取1h，100-300取3h，300-1000取6h"])
    rows.append(["入渗率f", "自动查表", "按附表12由Ptc和tc查产流期平均入渗率"])
    rows.append(["推理公式分区", result.rational_region_name, "表10"])
    rows.append(["汇流参数分区数", "3个", "黄土区、六盘山土石山林区、六盘山土石山植被一般区"])
    rows.append(["分区确定", "按工程所在水文地貌区选择", "黄土丘陵沟壑选黄土区；六盘山土石山林地选林区；植被一般区选植被一般区"])
    rows.append(["临界净雨深", "表10给定", "黄土区35mm、六盘山土石山林区20mm、六盘山土石山植被一般区30mm"])
    rows.append(["五点概化过程线系数", f"Ka={result.shape_ka:g}, Kw={result.shape_kw:g}, KT={result.shape_kt:g}", "自动按附图58/59查取，用于主雨峰五点概化过程线"])
    rows.append(["形状系数查图θ", f"{result.shape_theta:g}", result.shape_source])
    rows.append(["基流量(m³/s)", f"{result.baseflow_m3s:g}", result.baseflow_source])
    rows.append(["适用范围", "一般F<300km²", "图集推理公式法说明；超范围结果需复核"])
    parts.append(_fixed_table(rows, min_width=12))

    if shared.duration_point_rainfall and shared.duration_areal_rainfall:
        parts.append("")
        parts.append("二、控制时段点雨量与面雨量计算表")
        rain_rows = [["控制时段", "项目", *[freq_labels[freq] for freq in freqs]]]
        for duration in ("1h", "3h", "6h", "24h"):
            kp_values = [f"{(shared.duration_kp_values or {})[duration][freq]:.3f}" for freq in freqs]
            point_values = [f"{shared.duration_point_rainfall[duration][freq]:.3f}" for freq in freqs]
            areal_values = [f"{shared.duration_areal_rainfall[duration][freq]:.3f}" for freq in freqs]
            rain_rows.append([duration, "Kp", *kp_values])
            rain_rows.append([duration, "点雨量(mm)", *point_values])
            rain_rows.append([duration, "面雨量(mm)", *areal_values])
        parts.append(_fixed_table(rain_rows, min_width=10))

    if shared.rainfall_process_components:
        parts.append("")
        parts.append("三、时程分配雨型系数表")
        if shared.main_peak_hours == 1:
            dist_rows = [["时段", "H₁", "H₃-H₁", "H₆-H₃", "H₂₄-H₆"]]
            for hour, h1, h3_h1, h6_h3, h24_h6 in shared.rainfall_process_components:
                dist_rows.append([hour, f"{h1:g}", f"{h3_h1:g}", f"{h6_h3:g}", f"{h24_h6:g}"])
        elif shared.main_peak_hours == 3:
            dist_rows = [["时段", "H₃", "H₆-H₃", "H₂₄-H₆"]]
            for hour, _h1, h3, h6_h3, h24_h6 in shared.rainfall_process_components:
                dist_rows.append([hour, f"{h3:g}", f"{h6_h3:g}", f"{h24_h6:g}"])
        else:
            dist_rows = [["时段", "H₆", "H₂₄-H₆"]]
            for hour, _h1, _h3, h6, h24_h6 in shared.rainfall_process_components:
                dist_rows.append([hour, f"{h6:g}", f"{h24_h6:g}"])
        parts.append(_fixed_table(dist_rows, min_width=9))

    parts.append("")
    parts.append("四、净雨过程计算表")
    delta_t_label = f"{shared.delta_t_h:g}h" if shared.delta_t_h is not None else "1h"
    net_rows = [[f"时段(Δt={delta_t_label})"]]
    for freq in freqs:
        net_rows[0].extend([f"{freq}面雨量(mm)", f"{freq}入渗量(mm)", f"{freq}净雨深(mm)"])
    for row in shared.net_rain_table:
        values = [str(row.hour)]
        for freq in freqs:
            values.extend([
                f"{row.rainfall_mm[freq]:.3f}",
                f"{row.infiltration_mm[freq]:.3f}",
                f"{row.net_rain_mm[freq]:.3f}",
            ])
        net_rows.append(values)
    total = ["合计"]
    for freq in freqs:
        total.extend([
            f"{sum(row.rainfall_mm[freq] for row in shared.net_rain_table):.3f}",
            f"{sum(row.infiltration_mm[freq] for row in shared.net_rain_table):.3f}",
            f"{sum(row.net_rain_mm[freq] for row in shared.net_rain_table):.3f}",
        ])
    net_rows.append(total)
    parts.append(_fixed_table(net_rows, min_width=11))

    parts.append("")
    parts.append("五、推理公式参数与洪峰计算表")
    parameter_rows = [[
        "频率", "f(mm/h)", "tc(h)", "Ptc(mm)", "产流起止", "净雨深h(mm)", "平均净雨强(mm/h)",
        "θ", "h采用(mm)", "m", "Ka", "Kw", "KT", "τ(h)", "控制历时(h)", "控制净雨(mm)", "控制净雨强(mm/h)", "条件", "Q地表(m³/s)", "Qₚ含基流(m³/s)",
    ]]
    for freq in freqs:
        row = result.parameter_rows[freq]
        parameter_rows.append([
            freq_labels[freq],
            f"{float(row['f']):.3f}",
            f"{float(row['tc']):g}",
            f"{float(row['ptc']):.3f}",
            f"{float(row['start_hour']):g}-{float(row['end_hour']):g}",
            f"{float(row['net_total']):.3f}",
            f"{float(row['i']):.3f}",
            f"{float(row['theta']):.3f}",
            f"{float(row['h_used']):.3f}",
            f"{float(row['m']):.3f}",
            f"{float(row['ka']):.3f}",
            f"{float(row['kw']):.3f}",
            f"{float(row['kt']):.3f}",
            f"{float(row['tau']):.3f}",
            f"{float(row['max_net_duration']):.3f}",
            f"{float(row['max_net_depth']):.3f}",
            f"{float(row['max_net_intensity']):.3f}",
            str(row["condition"]),
            f"{float(row['q_surface']):.3f}",
            f"{float(row['q_peak']):.3f}",
        ])
    parts.append(_fixed_table(parameter_rows, min_width=10))

    parts.append("")
    parts.append("六、表10汇流参数公式说明")
    formula_rows = [["分区", "m公式", "θ公式", "临界净雨深(mm)", "本次是否采用"]]
    for region_id, data in RATIONAL_FORMULA_REGIONS.items():
        formula_rows.append([
            str(data["name"]),
            str(data["formula"]),
            str(data["theta_formula"]),
            f"{float(data['critical_h']):g}",
            "是" if region_id == result.rational_region_id else "否",
        ])
    parts.append(_fixed_table(formula_rows, min_width=12))

    for freq in freqs:
        flood = result.frequency_results[freq]
        parts.append("")
        parts.append(f"七、图解法求Qm和τ试算表 {freq_labels[freq]}")
        trial_rows = [["假设t(h)", "主峰累积hₜ(mm)", "hₜ/t(mm/h)", "Qₜ=0.278Fhₜ/t(m³/s)", "τₜ(h)", "τₜ-t", "累积起点(h)"]]
        for trial in flood.confluence_curve:
            duration, depth, intensity, q_trial, tau_trial, start = trial
            trial_rows.append([
                f"{duration:.3f}",
                f"{depth:.3f}",
                f"{intensity:.3f}",
                f"{q_trial:.3f}",
                f"{tau_trial:.3f}",
                f"{tau_trial - duration:.3f}",
                f"{start:.3f}",
            ])
        parts.append(_fixed_table(trial_rows, min_width=12))
        if "τ小于1个计算时段" in flood.runoff_condition:
            adopted_method = "τ<1h，按最大1h净雨强反算τ、Qm；不是新增假设t"
        elif "交点求" in flood.runoff_condition:
            adopted_method = "在τ_t-t换号的相邻假设t之间线性插值"
        else:
            adopted_method = "试算点未跨越，取τ_t与t差值绝对值最小点，需人工复核"
        adopted_rows = [
            ["采用项目", "采用值", "说明"],
            ["控制历时τ(h)", f"{flood.tau_h:.3f}", adopted_method],
            ["控制净雨hτ(mm)", f"{flood.max_net_depth_mm:.3f}", "由采用历时对应的平均净雨强换算"],
            ["控制净雨强(mm/h)", f"{flood.max_net_intensity_mmh:.3f}", "用于计算Qm"],
            ["Qₘ地表(m³/s)", f"{flood.surface_peak_m3s:.3f}", "推理公式洪峰"],
            ["累积起点(h)", f"{flood.main_start_h:.3f}", "主净雨峰连续累积起点"],
        ]
        parts.append(_fixed_table(adopted_rows, min_width=12))

        parts.append("")
        parts.append(f"八、推理公式汇流曲线计算表 {freq_labels[freq]}")
        curve_rows = [["t(h)", "Q地表+基流(m³/s)", "Q地表(m³/s)", "Q基(m³/s)", "Qₚ(m³/s)"]]
        for point in flood.hydrograph:
            surface = max(0.0, point.discharge_m3s - result.baseflow_m3s)
            marker = point.discharge_m3s if abs(point.discharge_m3s - flood.peak_flow_m3s) < 1e-9 else ""
            curve_rows.append([
                f"{point.time_h:g}",
                f"{point.discharge_m3s:.3f}",
                f"{surface:.3f}",
                f"{result.baseflow_m3s:.3f}",
                f"{marker:.3f}" if marker != "" else "",
            ])
        parts.append(_fixed_table(curve_rows, min_width=12))

    parts.append("")
    parts.append("九、洪峰流量成果表")
    summary_rows = [["频率", "地表洪峰Qₘ(m³/s)", "设计洪峰Qₚ(m³/s)", "峰现时间(h)", "τ(h)", "洪水总量(万m³)"]]
    for freq in freqs:
        flood = result.frequency_results[freq]
        summary_rows.append([
            freq_labels[freq],
            f"{flood.surface_peak_m3s:.3f}",
            f"{flood.peak_flow_m3s:.3f}",
            f"{flood.peak_time_h:g}",
            f"{flood.tau_h:.3f}",
            f"{flood.flood_volume_10k_m3:.3f}",
        ])
    parts.append(_fixed_table(summary_rows, min_width=12))

    return "\n".join(parts)


def instant_unit_line_audit_result_to_text(result: InstantUnitLineBatchResult, frequencies: Iterable[str] = DESIGN_FREQUENCIES_7) -> str:
    """Output instant unit line results in the detailed audit-atlas layout."""

    freqs = [freq for freq in frequencies if freq in result.flood_results]
    parts: list[str] = [
        "瞬时单位线法专项核验输出",
        "说明：本结果按《甘肃省暴雨洪水图集》瞬时单位线法组织；雨量净雨、时段单位线、洪水过程线按各频率分别列出，便于直接核对计算方法和水文图集过程。",
        "",
        "一、计算条件与参数来源",
    ]
    overview_rows = [["项目", "数值", "说明"]]
    overview_rows.append(["计算方法", "《甘肃省暴雨洪水图集》瞬时单位线法", "水文图集法"])
    overview_rows.append(["流域面积F(km²)", f"{result.area_km2:g}" if result.area_km2 is not None else "", "输入"])
    overview_rows.append(["计算步长Δt(h)", f"{result.delta_t_h:g}" if result.delta_t_h is not None else "1", "附表10固定1h"])
    overview_rows.append(["雨型分区", result.rainfall_zone_name or "", "输入选择"])
    overview_rows.append(["主雨峰", f"{result.main_peak_hours:g}h" if result.main_peak_hours else "", "由面积自动判定"])
    overview_rows.append(["雨型来源", result.distribution_source or "", "自动读取"])
    overview_rows.append(["入渗率f", "按附表12自动查算" if result.unit_parameter_rows else "按输入参数", "也可切换为手动复核"])
    parts.append(_fixed_table(overview_rows, min_width=12))

    parts.extend(["", "二、各频率洪峰成果表"])
    flow_rows = [["重现期", "洪峰Qₚ(m³/s)", "峰现时间(h)", "洪水总量(万m³)", "说明"]]
    for freq in freqs:
        peak, peak_time, volume = result.peak_summary[freq]
        flow_rows.append([f"{freq}年", f"{peak:.3f}", f"{peak_time:g}", f"{volume:.3f}", "洪水过程线梯形积分复核"])
    parts.append(_fixed_table(flow_rows, min_width=12))

    if result.unit_parameter_rows:
        parts.extend(["", "三、各频率产流与单位线参数表"])
        parameter_rows = [["重现期", "f(mm/h)", "tc(h)", "Ptc(mm)", "产流起止", "净雨深h(mm)", "平均净雨强i(mm/h)", "n", "k", "参数公式"]]
        for freq in freqs:
            row = result.unit_parameter_rows[freq]
            parameter_rows.append([
                f"{freq}年",
                f"{row['f']:.3f}",
                f"{row['tc']:g}",
                f"{row['ptc']:.3f}",
                f"{row['start_hour']:g}-{row['end_hour']:g}",
                f"{row['net_total']:.3f}",
                f"{row['i']:.3f}",
                f"{row['n']:.3f}",
                f"{row['k']:.3f}",
                str(row.get("formula_name", "")),
            ])
        parts.append(_fixed_table(parameter_rows, min_width=11))

    for index, freq in enumerate(freqs, 1):
        parts.extend(["", f"四-{index}、{freq}年逐时面雨量、入渗与净雨"])
        rainfall_rows = [["时段(h)", "面雨量(mm)", "入渗量(mm)", "净雨(mm)"]]
        for row in result.net_rain_table:
            rainfall_rows.append([
                f"{row.hour:g}",
                f"{row.rainfall_mm[freq]:.3f}",
                f"{row.infiltration_mm[freq]:.3f}",
                f"{row.net_rain_mm[freq]:.3f}",
            ])
        rainfall_rows.append([
            "合计",
            f"{sum(row.rainfall_mm[freq] for row in result.net_rain_table):.3f}",
            f"{sum(row.infiltration_mm[freq] for row in result.net_rain_table):.3f}",
            f"{sum(row.net_rain_mm[freq] for row in result.net_rain_table):.3f}",
        ])
        parts.append(_fixed_table(rainfall_rows, min_width=10))

        parts.extend(["", f"五-{index}、{freq}年时段单位线u(Δt,t)"])
        unit_rows = [["t(h)", "t/k", "S(t)", "S(t-Δt)", "u(Δt,t)"]]
        for row in result.flood_results[freq].unit_line:
            unit_rows.append([
                f"{row.time_h:g}",
                f"{row.t_over_k:.3f}",
                f"{row.s_t:.3f}",
                f"{row.s_prev:.3f}",
                f"{row.unit_depth_fraction:.3f}",
            ])
        parts.append(_fixed_table(unit_rows, min_width=10))

        parts.extend(["", f"六-{index}、{freq}年洪水过程线"])
        hydro_rows = [["t(h)", "Q(m³/s)"]]
        for point in result.flood_results[freq].hydrograph:
            hydro_rows.append([f"{point.time_h:g}", f"{point.discharge_m3s:.3f}"])
        parts.append(_fixed_table(hydro_rows, min_width=10))

    return "\n".join(parts)


def rational_formula_audit_result_to_text(result: RationalFormulaBatchResult, frequencies: Iterable[str] = DESIGN_FREQUENCIES_7) -> str:
    """Output rational formula results in the detailed audit-atlas layout."""

    freqs = [freq for freq in frequencies if freq in result.frequency_results]
    parts: list[str] = [
        "推理公式法专项核验输出",
        "说明：本结果按《甘肃省暴雨洪水图集》推理公式法组织；只输出推理公式法自身成果，瞬时单位线对比应放到专门的两法校核报告中，不在单方法计算结果里混排。",
        "",
        "一、计算条件与参数来源",
    ]
    overview_rows = [["项目", "数值", "说明"]]
    overview_rows.append(["计算方法", "《甘肃省暴雨洪水图集》推理公式法", "图集表10、推理公式及五点概化过程线"])
    overview_rows.append(["流域面积F(km²)", f"{result.shared.area_km2:g}" if result.shared.area_km2 is not None else "", "输入"])
    overview_rows.append(["主河道长度L(km)", f"{result.river_length_km:g}", "输入"])
    overview_rows.append(["主河槽比降J(‰)", f"{result.river_slope_permille:g}", "输入"])
    overview_rows.append(["雨型分区", result.shared.rainfall_zone_name or "", "与瞬时单位线共用"])
    overview_rows.append(["推理公式分区", result.rational_region_name, "表10 m关系"])
    overview_rows.append(["Ka/Kw/KT", f"{result.shape_ka:.3f}/{result.shape_kw:.3f}/{result.shape_kt:.3f}", result.shape_source])
    overview_rows.append(["基流Q0(m³/s)", f"{result.baseflow_m3s:.3f}", result.baseflow_source])
    parts.append(_fixed_table(overview_rows, min_width=12))

    parts.extend(["", "二、各频率洪峰成果表"])
    flow_rows = [["重现期", "地表洪峰Qₘ(m³/s)", "设计洪峰Qₚ(m³/s)", "峰现时间(h)", "τ(h)", "洪水总量(万m³)", "说明"]]
    for freq in freqs:
        rational = result.frequency_results[freq]
        flow_rows.append([
            f"{freq}年",
            f"{rational.surface_peak_m3s:.3f}",
            f"{rational.peak_flow_m3s:.3f}",
            f"{rational.peak_time_h:g}",
            f"{rational.tau_h:.3f}",
            f"{rational.flood_volume_10k_m3:.3f}",
            rational.runoff_condition,
        ])
    parts.append(_fixed_table(flow_rows, min_width=12))

    parts.extend(["", "三、推理公式参数表"])
    parameter_rows = [["重现期", "f(mm/h)", "tc(h)", "Ptc(mm)", "净雨深h(mm)", "平均净雨强(mm/h)", "m", "θ", "Ka", "Kw", "KT", "基流Q0"]]
    for freq in freqs:
        row = result.parameter_rows[freq]
        rational = result.frequency_results[freq]
        parameter_rows.append([
            f"{freq}年",
            f"{float(row['f']):.3f}",
            f"{float(row['tc']):g}",
            f"{float(row['ptc']):.3f}",
            f"{float(row['net_total']):.3f}",
            f"{float(row['i']):.3f}",
            f"{rational.m:.3f}",
            f"{rational.theta:.3f}",
            f"{result.shape_ka:.3f}",
            f"{result.shape_kw:.3f}",
            f"{result.shape_kt:.3f}",
            f"{result.baseflow_m3s:.3f}",
        ])
    parts.append(_fixed_table(parameter_rows, min_width=11))

    for index, freq in enumerate(freqs, 1):
        parts.extend(["", f"四-{index}、{freq}年逐时面雨量、入渗量、净雨"])
        rain_rows = [["时段(h)", "面雨量(mm)", "入渗量(mm)", "净雨(mm)"]]
        for row in result.shared.net_rain_table:
            rain_rows.append([
                f"{row.hour:g}",
                f"{row.rainfall_mm[freq]:.3f}",
                f"{row.infiltration_mm[freq]:.3f}",
                f"{row.net_rain_mm[freq]:.3f}",
            ])
        rain_rows.append([
            "合计",
            f"{sum(row.rainfall_mm[freq] for row in result.shared.net_rain_table):.3f}",
            f"{sum(row.infiltration_mm[freq] for row in result.shared.net_rain_table):.3f}",
            f"{sum(row.net_rain_mm[freq] for row in result.shared.net_rain_table):.3f}",
        ])
        parts.append(_fixed_table(rain_rows, min_width=10))

        parts.extend(["", f"五-{index}、{freq}年推理公式洪水过程线"])
        hydro_rows = [["t(h)", "推理公式Q(m³/s)"]]
        for point in result.frequency_results[freq].hydrograph:
            hydro_rows.append([f"{point.time_h:g}", f"{point.discharge_m3s:.3f}"])
        parts.append(_fixed_table(hydro_rows, min_width=10))

        parts.extend(["", f"六-{index}、{freq}年图解法求Qm与τ试算表"])
        trial_rows = [["假设t(h)", "累计hₜ(mm)", "hₜ/t(mm/h)", "Qₜ(m³/s)", "τₜ(h)", "τₜ-t", "起点(h)"]]
        rational = result.frequency_results[freq]
        for trial in rational.confluence_curve:
            duration, depth, intensity, q_trial, tau_trial, start = trial
            trial_rows.append([
                f"{duration:.3f}",
                f"{depth:.3f}",
                f"{intensity:.3f}",
                f"{q_trial:.3f}",
                f"{tau_trial:.3f}",
                f"{tau_trial - duration:.3f}",
                f"{start:.3f}",
            ])
        parts.append(_fixed_table(trial_rows, min_width=10))

    return "\n".join(parts)


def instant_unit_line_result_to_text(result: InstantUnitLineBatchResult, frequencies: Iterable[str] = DESIGN_FREQUENCIES_7) -> str:
    """输出可复制到 Excel 的制表符文本。"""

    freqs = [freq for freq in frequencies if freq in result.flood_results]
    freq_labels = {freq: f"{freq}年" if freq in TIE_FREQUENCIES else freq for freq in freqs}
    parts: list[str] = []

    parts.append("一、输入参数与自动计算参数")
    parameter_rows = [["项目", "数值", "说明"]]
    if result.area_km2 is not None:
        parameter_rows.append(["流域面积F(km²)", f"{result.area_km2:g}", "输入"])
    if result.delta_t_h is not None:
        parameter_rows.append(["计算步长Δt(h)", f"{result.delta_t_h:g}", "图集附表10为1小时步长，软件默认固定为1"])
    if result.rainfall_zone_name:
        parameter_rows.append(["雨型分区", result.rainfall_zone_name, "输入选择"])
    if result.main_peak_hours is not None:
        parameter_rows.append(["主雨峰", f"{result.main_peak_hours}h", "由流域面积自动选择"])
    if result.ptc_mean_mm is not None:
        parameter_rows.append(["24h点雨量均值(mm)", f"{result.ptc_mean_mm:g}", "输入"])
    if result.cv is not None:
        parameter_rows.append(["24h Cv", f"{result.cv:g}", "输入"])
    if result.cs_ratio is not None:
        parameter_rows.append(["Cs/Cv", f"{result.cs_ratio:g}", "输入，默认按P-III表为3.5"])
    if result.distribution_source:
        parameter_rows.append(["雨型来源", result.distribution_source, "自动读取"])
    parameter_rows.append(["f(mm/h)", "产流期平均入渗率", "默认按附表12由Ptc和tc自动查算；也可切换为手动复核"])
    parameter_rows.append(["ΔF(mm)", "时段入渗量", "ΔF=f×Δt；当前Δt=1h时，数值等于f"])
    parameter_rows.append(["平均净雨强i(mm/h)", "i=h/tc", "具体数值见“四、瞬时单位线参数计算表”的i列"])
    parameter_rows.append(["tc判定", "自动查算", "候选产流段需满足段内每小时面雨量 >= f×Δt，再按净雨量最大选取"])
    parameter_rows.append(["tc说明", "非固定值", "若主雨峰为1h且雨量高度集中，tc可能为1h；书上例题可自动得到tc=6h"])
    parameter_rows.append(["同一tc说明", "可能出现", "各频率使用同一雨型比例分配时，最优连续产流段可能相同，tc可相同但f、Ptc、净雨深不同"])
    parts.append(_fixed_table(parameter_rows, min_width=12))

    if result.duration_point_rainfall and result.duration_areal_rainfall:
        parts.append("")
        parts.append("二、控制时段点雨量与面雨量计算表")
        rows = [["控制时段", "项目", *[freq_labels[freq] for freq in freqs]]]
        for duration in ("1h", "3h", "6h", "24h"):
            if result.duration_kp_values and duration in result.duration_kp_values:
                kp_values = [f"{result.duration_kp_values[duration][freq]:.3f}" for freq in freqs]
            else:
                kp_values = ["由H1/H6组合" for _ in freqs]
            point_values = [f"{result.duration_point_rainfall[duration][freq]:.3f}" for freq in freqs]
            areal_values = [f"{result.duration_areal_rainfall[duration][freq]:.3f}" for freq in freqs]
            rows.append([duration, "Kp", *kp_values])
            rows.append([duration, "点雨量(mm)", *point_values])
            rows.append([duration, "面雨量(mm)", *areal_values])
        parts.append(_fixed_table(rows, min_width=10))

    if result.rainfall_process_components:
        parts.append("")
        parts.append("三、时程分配雨型系数表")
        if result.main_peak_hours == 1:
            rows = [["时段", "H₁", "H₃-H₁", "H₆-H₃", "H₂₄-H₆"]]
            for hour, h1, h3_h1, h6_h3, h24_h6 in result.rainfall_process_components:
                rows.append([hour, f"{h1:g}", f"{h3_h1:g}", f"{h6_h3:g}", f"{h24_h6:g}"])
            parts.append(_fixed_table(rows, min_width=9))
        elif result.main_peak_hours == 3:
            rows = [["时段", "H₃", "H₆-H₃", "H₂₄-H₆"]]
            for hour, _h1, h3, h6_h3, h24_h6 in result.rainfall_process_components:
                rows.append([hour, f"{h3:g}", f"{h6_h3:g}", f"{h24_h6:g}"])
            parts.append(_fixed_table(rows, min_width=9))
        elif result.main_peak_hours == 6:
            rows = [["时段", "H₆", "H₂₄-H₆"]]
            for hour, _h1, _h3, h6, h24_h6 in result.rainfall_process_components:
                rows.append([hour, f"{h6:g}", f"{h24_h6:g}"])
            parts.append(_fixed_table(rows, min_width=9))

    if result.unit_parameter_rows:
        parts.append("")
        parts.append("四、瞬时单位线参数计算表")
        rows = [["频率", "f(mm/h)", "f来源", "tc(h)", "tc判定", "Ptc(mm)", "产流起止时段", "时段入渗量ΔF(mm)", "初损I₀(mm)", "净雨深h(mm)", "平均净雨强i(mm/h)", "公式", "m₁,₁₀", "b", "m₁", "m₂", "n", "k"]]
        for freq in freqs:
            row = result.unit_parameter_rows[freq]
            rows.append([
                freq_labels[freq],
                f"{row['f']:.3f}",
                str(row.get("f_source", "")),
                f"{row['tc']:g}",
                str(row.get("tc_rule", "")),
                f"{row['ptc']:.3f}",
                f"{row['start_hour']:g}-{row['end_hour']:g}",
                f"{row['infiltration_depth']:.3f}",
                f"{row['initial_loss']:.3f}",
                f"{row['net_total']:.3f}",
                f"{row['i']:.3f}",
                str(row.get("formula_name", "")),
                f"{row['m1_10']:.3f}",
                f"{row['b']:.3f}",
                f"{row['m1']:.3f}",
                f"{row['m2']:.3f}",
                f"{row['n']:.3f}",
                f"{row['k']:.3f}",
            ])
        parts.append(_fixed_table(rows, min_width=9))

    if result.frequency_inputs and not result.unit_parameter_rows:
        parts.append("")
        parts.append("四-1、各频率Ptc与产流参数")
        rows = [["频率", "Kp", "Ptc(mm)", "tc(h)", "入渗率f(mm/h)", "时段入渗量ΔF(mm)", "n", "k", "起始时段"]]
        for freq in freqs:
            item = result.frequency_inputs[freq]
            kp = "" if not result.kp_values else f"{result.kp_values.get(freq, '')}"
            delta_f = item.infiltration_rate_mmh * (result.delta_t_h or 1.0)
            rows.append([
                freq_labels[freq],
                kp,
                f"{item.ptc_mm:.3f}",
                f"{item.runoff_duration_h:g}",
                f"{item.infiltration_rate_mmh:g}",
                f"{delta_f:.3f}",
                f"{item.n:g}",
                f"{item.k:g}",
                f"{item.runoff_start_hour:g}",
            ])
        parts.append(_fixed_table(rows, min_width=9))
        parts.append("")

    parts.append("五、净雨过程计算表")
    delta_t_label = f"{result.delta_t_h:g}h" if result.delta_t_h is not None else "1h"
    rows = [[f"时段(Δt={delta_t_label})"]]
    for freq in freqs:
        rows[0].extend([f"{freq}面雨量(mm)", f"{freq}入渗量(mm)", f"{freq}净雨深(mm)"])
    for row in result.net_rain_table:
        values = [str(row.hour)]
        for freq in freqs:
            values.extend([
                f"{row.rainfall_mm[freq]:.3f}",
                f"{row.infiltration_mm[freq]:.3f}",
                f"{row.net_rain_mm[freq]:.3f}",
            ])
        rows.append(values)
    total = ["合计"]
    for freq in freqs:
        total.extend([
            f"{sum(row.rainfall_mm[freq] for row in result.net_rain_table):.3f}",
            f"{sum(row.infiltration_mm[freq] for row in result.net_rain_table):.3f}",
            f"{sum(row.net_rain_mm[freq] for row in result.net_rain_table):.3f}",
        ])
    rows.append(total)
    parts.append(_fixed_table(rows, min_width=11))

    for freq in freqs:
        flood = result.flood_results[freq]
        runoff_input = (result.frequency_inputs or {}).get(freq) if result.frequency_inputs else None
        parts.append("")
        parts.append(f"八、洪水过程线计算表 {freq_labels[freq]}")
        unit_len = len(flood.unit_line)
        net_values = [row.net_rain_mm[freq] for row in result.net_rain_table]
        effective_rain = [(idx, depth) for idx, depth in enumerate(net_values) if depth > 1e-9]
        rows = [["t(h)", "t/k", "s(t)", "s(t-Δt)", "u(t)"]]
        rows[0].extend([f"u*h{idx + 1}={depth:.2f}" for idx, depth in effective_rain])
        rows[0].extend(["Σu*hᵢ", "Q表(m³/s)", "Q基(m³/s)", "Q净(m³/s)", "Qₚ(m³/s)", "换算系数C"])

        max_rows = max(unit_len, len(flood.hydrograph))
        conversion = ""
        if result.area_km2 is not None and result.delta_t_h:
            conversion = f"{result.area_km2 / (3.6 * result.delta_t_h):.3f}"
        for idx in range(max_rows):
            line_values: list[str] = []
            if idx < unit_len:
                u = flood.unit_line[idx]
                line_values.extend([
                    f"{u.time_h:g}",
                    f"{u.t_over_k:.3f}",
                    f"{u.s_t:.3f}",
                    f"{u.s_prev:.3f}",
                    f"{u.unit_depth_fraction:.3f}",
                ])
            else:
                t_value = flood.hydrograph[idx].time_h if idx < len(flood.hydrograph) else (idx + 1) * (result.delta_t_h or 1.0)
                if runoff_input is not None:
                    delta_t = result.delta_t_h or 1.0
                    s_t = nash_s_curve(t_value, runoff_input.n, runoff_input.k)
                    s_prev = nash_s_curve(max(0.0, t_value - delta_t), runoff_input.n, runoff_input.k)
                    u_value = max(0.0, s_t - s_prev)
                    line_values.extend([
                        f"{t_value:g}",
                        f"{t_value / runoff_input.k:.3f}",
                        f"{s_t:.3f}",
                        f"{s_prev:.3f}",
                        f"{u_value:.3f}",
                    ])
                else:
                    line_values.extend([f"{t_value:g}", "0.000", "1.000", "1.000", "0.000"])

            convolution_values = []
            for rain_idx, depth in effective_rain:
                unit_idx = idx - rain_idx
                if 0 <= unit_idx < unit_len:
                    convolution_values.append(depth * flood.unit_line[unit_idx].unit_depth_fraction)
                else:
                    convolution_values.append(0.0)
            depth_sum = sum(convolution_values)
            line_values.extend(f"{value:.3f}" for value in convolution_values)

            q = flood.hydrograph[idx].discharge_m3s if idx < len(flood.hydrograph) else 0.0
            q_peak_marker = q if abs(q - flood.peak_flow_m3s) < 1e-9 else ""
            line_values.extend([
                f"{depth_sum:.3f}",
                f"{q:.3f}" if idx < len(flood.hydrograph) else "",
                "0.000" if idx < len(flood.hydrograph) else "",
                f"{q:.3f}" if idx < len(flood.hydrograph) else "",
                f"{q_peak_marker:.3f}" if q_peak_marker != "" else "",
                conversion,
            ])
            rows.append(line_values)
        parts.append(_fixed_table(rows, min_width=9))

    parts.append("")
    parts.append("九、洪峰流量成果表")
    rows = [["频率", "洪峰流量Qₚ(m³/s)", "峰现时间(h)", "洪水总量(万m³)"]]
    for freq in freqs:
        peak, peak_time, volume = result.peak_summary[freq]
        rows.append([freq_labels[freq], f"{peak:.3f}", f"{peak_time:g}", f"{volume:.3f}"])
    parts.append(_fixed_table(rows, min_width=12))

    return "\n".join(parts)
