from __future__ import annotations

from dataclasses import dataclass
import math


class ChannelHydraulicError(ValueError):
    pass


def _positive(value: float, name: str) -> float:
    if value <= 0:
        raise ChannelHydraulicError(f"{name}必须大于0")
    return value


def _non_negative(value: float, name: str) -> float:
    if value < 0:
        raise ChannelHydraulicError(f"{name}不能小于0")
    return value


def _manning_chezy(n: float, radius: float) -> float:
    return (1.0 / n) * radius ** (1.0 / 6.0)


@dataclass(frozen=True)
class OpenChannelInput:
    bottom_width: float
    water_depth: float
    side_slope: float
    roughness: float
    bed_slope: float
    enlarge_coefficient: float = 1.15


@dataclass(frozen=True)
class OpenChannelResult:
    area: float
    wetted_perimeter: float
    hydraulic_radius: float
    chezy: float
    velocity: float
    flow: float
    enlarged_flow: float
    enlarged_depth: float
    enlarged_area: float
    freeboard: float
    berm_width: float
    bend_radius: float
    formula: str


def calculate_open_channel(data: OpenChannelInput) -> OpenChannelResult:
    b = _positive(data.bottom_width, "渠底宽b")
    h = _positive(data.water_depth, "水深h")
    m = _non_negative(data.side_slope, "边坡系数m")
    n = _positive(data.roughness, "糙率n")
    i = _positive(data.bed_slope, "渠底比降i")
    coeff = _positive(data.enlarge_coefficient, "加大流量系数")

    area = (b + m * h) * h
    wetted = b + 2.0 * h * math.sqrt(1.0 + m * m)
    radius = area / wetted
    chezy = _manning_chezy(n, radius)
    velocity = chezy * math.sqrt(radius * i)
    flow = area * velocity
    enlarged_flow = flow * coeff
    enlarged_area = area * coeff
    enlarged_depth = solve_trapezoid_depth_by_area(enlarged_area, b, m)
    freeboard = 0.25 * enlarged_depth + 0.2
    berm_width = enlarged_depth + 0.3
    bend_radius = 11.0 * velocity * velocity * math.sqrt(area) + 12.0
    return OpenChannelResult(
        area=area,
        wetted_perimeter=wetted,
        hydraulic_radius=radius,
        chezy=chezy,
        velocity=velocity,
        flow=flow,
        enlarged_flow=enlarged_flow,
        enlarged_depth=enlarged_depth,
        enlarged_area=enlarged_area,
        freeboard=freeboard,
        berm_width=berm_width,
        bend_radius=bend_radius,
        formula="Q=A*C*(R*i)^0.5；A=(b+m*h)*h；x=b+2*h*(1+m^2)^0.5；R=A/x；C=R^(1/6)/n",
    )


def open_channel_flow(bottom_width: float, water_depth: float, side_slope: float, roughness: float, bed_slope: float) -> float:
    b = _positive(bottom_width, "渠底宽b")
    h = _positive(water_depth, "水深h")
    m = _non_negative(side_slope, "边坡系数m")
    n = _positive(roughness, "糙率n")
    i = _positive(bed_slope, "渠底比降i")
    area = (b + m * h) * h
    wetted = b + 2.0 * h * math.sqrt(1.0 + m * m)
    radius = area / wetted
    return area * _manning_chezy(n, radius) * math.sqrt(radius * i)


def solve_open_channel_depth(
    target_flow: float,
    bottom_width: float,
    side_slope: float,
    roughness: float,
    bed_slope: float,
) -> float:
    q = _positive(target_flow, "目标流量Q")
    high = 0.1
    while open_channel_flow(bottom_width, high, side_slope, roughness, bed_slope) < q:
        high *= 2.0
        if high > 10000:
            raise ChannelHydraulicError("未能求得水深，请检查流量、比降或糙率")
    low = 0.0
    for _ in range(90):
        mid = (low + high) / 2.0
        if open_channel_flow(bottom_width, mid, side_slope, roughness, bed_slope) < q:
            low = mid
        else:
            high = mid
    return (low + high) / 2.0


def solve_trapezoid_depth_by_area(target_area: float, bottom_width: float, side_slope: float) -> float:
    area = _positive(target_area, "目标过水面积A")
    b = _positive(bottom_width, "渠底宽b")
    m = _non_negative(side_slope, "边坡系数m")
    if m == 0:
        return area / b
    return (-b + math.sqrt(b * b + 4.0 * m * area)) / (2.0 * m)


@dataclass(frozen=True)
class CriticalDepthInput:
    flow: float
    bottom_width: float
    side_slope: float = 0.0
    alpha: float = 1.0
    gravity: float = 9.81


@dataclass(frozen=True)
class CriticalDepthResult:
    critical_depth: float
    unit_discharge: float
    top_width: float
    area: float
    criterion: float
    formula: str


def calculate_critical_depth(data: CriticalDepthInput) -> CriticalDepthResult:
    q_total = _positive(data.flow, "流量Q")
    b = _positive(data.bottom_width, "渠底宽b")
    m = _non_negative(data.side_slope, "边坡系数m")
    alpha = _positive(data.alpha, "动能修正系数alpha")
    g = _positive(data.gravity, "重力加速度g")
    if m == 0:
        unit_q = q_total / b
        depth = (alpha * unit_q * unit_q / g) ** (1.0 / 3.0)
        area = b * depth
        top_width = b
        criterion = alpha * q_total * q_total * top_width / (g * area ** 3)
        formula = "矩形：h=(alpha*q^2/g)^(1/3)，q=Q/b"
    else:
        depth = _solve_critical_trapezoid_by_book_formula(q_total, b, m, g)
        unit_q = q_total / (b + 0.8 * m * depth)
        area = (b + m * depth) * depth
        top_width = b + 2.0 * m * depth
        left = (unit_q * unit_q / g) ** (1.0 / 3.0)
        right = ((1.0 + m * depth / b) * depth) / ((1.0 + 2.0 * m * depth / b) ** (1.0 / 3.0))
        criterion = left / right
        formula = "梯形：q=Q/(b+0.8*m*h)；(q^2/g)^(1/3)=[(1+m*h/b)*h]/(1+2*m*h/b)^(1/3)"
    return CriticalDepthResult(
        critical_depth=depth,
        unit_discharge=unit_q,
        top_width=top_width,
        area=area,
        criterion=criterion,
        formula=formula,
    )


def _solve_critical_trapezoid(flow: float, bottom_width: float, side_slope: float, alpha: float, gravity: float) -> float:
    def criterion(depth: float) -> float:
        area = (bottom_width + side_slope * depth) * depth
        top_width = bottom_width + 2.0 * side_slope * depth
        return alpha * flow * flow * top_width / (gravity * area ** 3)

    low = 1e-9
    high = 0.1
    while criterion(high) > 1.0:
        high *= 2.0
        if high > 10000:
            raise ChannelHydraulicError("未能求得临界水深，请检查流量或断面参数")
    for _ in range(90):
        mid = (low + high) / 2.0
        if criterion(mid) > 1.0:
            low = mid
        else:
            high = mid
    return (low + high) / 2.0


def _solve_critical_trapezoid_by_book_formula(flow: float, bottom_width: float, side_slope: float, gravity: float) -> float:
    def residual(depth: float) -> float:
        unit_q = flow / (bottom_width + 0.8 * side_slope * depth)
        left = (unit_q * unit_q / gravity) ** (1.0 / 3.0)
        right = ((1.0 + side_slope * depth / bottom_width) * depth) / (
            (1.0 + 2.0 * side_slope * depth / bottom_width) ** (1.0 / 3.0)
        )
        return left - right

    low = 1e-9
    high = 0.1
    while residual(high) > 0:
        high *= 2.0
        if high > 10000:
            raise ChannelHydraulicError("未能求得临界水深，请检查流量或断面参数")
    for _ in range(90):
        mid = (low + high) / 2.0
        if residual(mid) > 0:
            low = mid
        else:
            high = mid
    return (low + high) / 2.0


@dataclass(frozen=True)
class BestHydraulicSectionInput:
    flow: float
    roughness: float
    bed_slope: float
    side_slope: float = 0.0


@dataclass(frozen=True)
class BestHydraulicSectionResult:
    section_type: str
    water_depth: float
    bottom_width: float
    width_depth_ratio: float
    area: float
    wetted_perimeter: float
    hydraulic_radius: float
    chezy: float
    velocity: float
    check_flow: float
    formula: str


def calculate_best_hydraulic_section(data: BestHydraulicSectionInput) -> BestHydraulicSectionResult:
    q = _positive(data.flow, "设计流量Q")
    n = _positive(data.roughness, "糙率n")
    i = _positive(data.bed_slope, "渠底比降i")
    m = _non_negative(data.side_slope, "边坡系数m")
    if m == 0:
        depth = ((q * n * (2.0 ** (2.0 / 3.0))) / (2.0 * math.sqrt(i))) ** (3.0 / 8.0)
        bottom_width = 2.0 * depth
        section_type = "矩形"
        formula = "矩形水力最佳断面：b/h=2；Q=A*C*(R*i)^0.5"
    else:
        m1 = 2.0 * math.sqrt(1.0 + m * m)
        depth = 1.189 * (n * q / ((m1 - m) * math.sqrt(i))) ** (3.0 / 8.0)
        bottom_width = 2.0 * depth * (math.sqrt(1.0 + m * m) - m)
        section_type = "梯形"
        formula = "梯形水力最佳断面：h0=1.189*[nQ/((m1-m)*i^0.5)]^(3/8)，b0=2h0*((1+m^2)^0.5-m)"
    area = (bottom_width + m * depth) * depth
    wetted = bottom_width + 2.0 * depth * math.sqrt(1.0 + m * m)
    radius = area / wetted
    chezy = _manning_chezy(n, radius)
    velocity = chezy * math.sqrt(radius * i)
    check_flow = area * velocity
    return BestHydraulicSectionResult(
        section_type=section_type,
        water_depth=depth,
        bottom_width=bottom_width,
        width_depth_ratio=bottom_width / depth,
        area=area,
        wetted_perimeter=wetted,
        hydraulic_radius=radius,
        chezy=chezy,
        velocity=velocity,
        check_flow=check_flow,
        formula=formula,
    )


@dataclass(frozen=True)
class IntakeAngleInput:
    river_velocity: float
    canal_velocity: float


@dataclass(frozen=True)
class IntakeAngleResult:
    velocity_ratio: float
    angle_degrees: float
    formula: str


def calculate_intake_angle(data: IntakeAngleInput) -> IntakeAngleResult:
    v0 = _positive(data.river_velocity, "河道流速v0")
    v1 = _positive(data.canal_velocity, "渠道流速v1")
    ratio = v0 / v1
    if ratio > 1.0:
        raise ChannelHydraulicError("书中公式要求v0/v1不大于1；通常应满足渠道流速大于河道流速")
    return IntakeAngleResult(
        velocity_ratio=ratio,
        angle_degrees=math.degrees(math.acos(ratio)),
        formula="δ=arccos(v0/v1)",
    )


@dataclass(frozen=True)
class NonScourVelocityInput:
    bottom_width: float
    water_depth: float
    side_slope: float
    roughness: float
    bed_slope: float
    enlarge_coefficient: float
    soil_coefficient: float
    block_length: float
    block_width: float
    block_height: float
    rock_lining_water_depth: float | None = None


@dataclass(frozen=True)
class NonScourVelocityResult:
    design_depth: float
    design_area: float
    design_wetted_perimeter: float
    design_hydraulic_radius: float
    rock_lining_water_depth: float
    soil_velocity: float
    equivalent_block_diameter: float
    rock_lining_velocity: float
    formula: str


def calculate_non_scour_velocity(data: NonScourVelocityInput) -> NonScourVelocityResult:
    open_result = calculate_open_channel(
        OpenChannelInput(
            bottom_width=data.bottom_width,
            water_depth=data.water_depth,
            side_slope=data.side_slope,
            roughness=data.roughness,
            bed_slope=data.bed_slope,
            enlarge_coefficient=data.enlarge_coefficient,
        )
    )
    b = _positive(data.bottom_width, "渠底宽b")
    design_h = _positive(data.water_depth, "水深h")
    m = _non_negative(data.side_slope, "边坡系数m")
    a_coeff = _positive(data.soil_coefficient, "不冲流速系数A")
    length = _positive(data.block_length, "块石长度a")
    width = _positive(data.block_width, "块石宽度b")
    height = _positive(data.block_height, "块石高度c")

    design_area = open_result.area
    design_wetted = open_result.wetted_perimeter
    design_radius = open_result.hydraulic_radius
    rock_h = data.rock_lining_water_depth
    if rock_h is None:
        rock_h = open_result.enlarged_depth
    rock_h = _positive(rock_h, "砌石护面渠道加大流量水深h")
    area = (b + m * rock_h) * rock_h
    wetted = b + 2.0 * rock_h * math.sqrt(1.0 + m * m)
    radius = area / wetted
    soil_velocity = a_coeff * design_radius ** 0.4
    equivalent_diameter = 1.24 * (length * width * height) ** (1.0 / 3.0)
    rock_velocity = (
        2.7
        * (height / math.sqrt(length * width)) ** (-0.5)
        * equivalent_diameter ** (5.0 / 14.0)
        * rock_h ** (1.0 / 7.0)
    )
    return NonScourVelocityResult(
        design_depth=design_h,
        design_area=design_area,
        design_wetted_perimeter=design_wetted,
        design_hydraulic_radius=design_radius,
        rock_lining_water_depth=rock_h,
        soil_velocity=soil_velocity,
        equivalent_block_diameter=equivalent_diameter,
        rock_lining_velocity=rock_velocity,
        formula="土渠：vp=A*R^0.4；砌石护面：vf=2.7*(c/sqrt(ab))^(-1/2)*d^(5/14)*h^(1/7)，d=1.24*(abc)^(1/3)",
    )


@dataclass(frozen=True)
class NonSiltingVelocityInput:
    coefficient: float
    minimum_flow: float


@dataclass(frozen=True)
class NonSiltingVelocityResult:
    velocity: float
    formula: str


def calculate_non_silting_velocity(data: NonSiltingVelocityInput) -> NonSiltingVelocityResult:
    c0 = _positive(data.coefficient, "不淤流速系数C0")
    q = _positive(data.minimum_flow, "渠道最小流量Q")
    return NonSiltingVelocityResult(
        velocity=c0 * math.sqrt(q),
        formula="vcd=C0*sqrt(Q)",
    )


@dataclass(frozen=True)
class PracticalEconomicSectionInput:
    flow: float
    roughness: float
    bed_slope: float
    side_slope: float
    area_ratio: float


@dataclass(frozen=True)
class PracticalEconomicSectionResult:
    best_hydraulic_depth: float
    depth_ratio: float
    water_depth: float
    bottom_width: float
    width_depth_ratio: float
    area: float
    wetted_perimeter: float
    hydraulic_radius: float
    velocity: float
    check_flow: float
    formula: str


def calculate_practical_economic_section(data: PracticalEconomicSectionInput) -> PracticalEconomicSectionResult:
    q = _positive(data.flow, "设计流量Q")
    n = _positive(data.roughness, "糙率n")
    i = _positive(data.bed_slope, "渠底比降i")
    m = _non_negative(data.side_slope, "边坡系数m")
    area_ratio = _positive(data.area_ratio, "面积偏离系数a")
    if area_ratio < 1.0:
        raise ChannelHydraulicError("最佳实用经济断面要求面积偏离系数a不小于1")

    best = calculate_best_hydraulic_section(
        BestHydraulicSectionInput(flow=q, roughness=n, bed_slope=i, side_slope=m)
    )
    if abs(area_ratio - 1.0) < 1e-12:
        depth_ratio = 1.0
    else:
        discriminant = area_ratio ** 5 - area_ratio
        if discriminant < 0:
            raise ChannelHydraulicError("面积偏离系数无法求得实用经济断面，请检查a值")
        # The shallow practical section is the common economic alternative in the book examples.
        depth_ratio = area_ratio ** 2.5 - math.sqrt(discriminant)
    depth = best.water_depth * depth_ratio
    beta = (
        area_ratio
        / (depth_ratio * depth_ratio)
        * (2.0 * math.sqrt(1.0 + m * m) - m)
        - m
    )
    if beta <= 1.0:
        raise ChannelHydraulicError("求得b/h不大于1，不满足书中最佳实用经济断面提示")
    bottom_width = beta * depth
    area = (bottom_width + m * depth) * depth
    wetted = bottom_width + 2.0 * depth * math.sqrt(1.0 + m * m)
    radius = area / wetted
    chezy = _manning_chezy(n, radius)
    velocity = chezy * math.sqrt(radius * i)
    check_flow = area * velocity
    return PracticalEconomicSectionResult(
        best_hydraulic_depth=best.water_depth,
        depth_ratio=depth_ratio,
        water_depth=depth,
        bottom_width=bottom_width,
        width_depth_ratio=beta,
        area=area,
        wetted_perimeter=wetted,
        hydraulic_radius=radius,
        velocity=velocity,
        check_flow=check_flow,
        formula="hm按水力最佳断面求得；(h/hm)^2-2*a^(5/2)*(h/hm)+a=0；b/h=a/(h/hm)^2*(2*sqrt(1+m^2)-m)-m",
    )


@dataclass(frozen=True)
class LinedSeepageInput:
    channel_length: float
    lining_coefficient: float
    bottom_width: float
    water_depth: float
    side_slope: float
    roughness: float
    bed_slope: float


@dataclass(frozen=True)
class LinedSeepageResult:
    area: float
    wetted_perimeter: float
    hydraulic_radius: float
    chezy: float
    velocity: float
    seepage_loss: float
    formula: str


def calculate_lined_seepage_loss(data: LinedSeepageInput) -> LinedSeepageResult:
    length = _positive(data.channel_length, "渠道长度l")
    c1 = _positive(data.lining_coefficient, "衬砌系数c1")
    b = _positive(data.bottom_width, "渠底宽b")
    h = _positive(data.water_depth, "水深h")
    m = _non_negative(data.side_slope, "边坡系数m")
    n = _positive(data.roughness, "糙率n")
    i = _positive(data.bed_slope, "渠底比降i")

    area = (b + m * h) * h
    wetted = b + 2.0 * h * math.sqrt(1.0 + m * m)
    radius = area / wetted
    chezy = _manning_chezy(n, radius)
    velocity = chezy * math.sqrt(radius * i)
    loss = 0.45 * c1 * (wetted * length) / (4.0e6 + 3650.0 * math.sqrt(velocity)) * h ** (1.0 / 3.0)
    return LinedSeepageResult(
        area=area,
        wetted_perimeter=wetted,
        hydraulic_radius=radius,
        chezy=chezy,
        velocity=velocity,
        seepage_loss=loss,
        formula="S=0.45*c1*x*l/(4*10^6+3650*sqrt(v))*h^(1/3)",
    )


@dataclass(frozen=True)
class UnlinedSeepageInput:
    channel_length_km: float
    permeability_coefficient: float
    net_flow: float


@dataclass(frozen=True)
class UnlinedSeepageResult:
    loss_per_km: float
    total_loss: float
    formula: str


def calculate_unlined_seepage_loss(data: UnlinedSeepageInput) -> UnlinedSeepageResult:
    length = _positive(data.channel_length_km, "渠道长度L")
    kp = _positive(data.permeability_coefficient, "渗透系数kp")
    q = _positive(data.net_flow, "渠道净流量Q净")
    loss_per_km = 0.063 * kp * math.sqrt(q)
    return UnlinedSeepageResult(
        loss_per_km=loss_per_km,
        total_loss=loss_per_km * length,
        formula="S=0.063*kp*sqrt(Q净)，总损失=S*L",
    )


@dataclass(frozen=True)
class UChannelInput:
    water_width: float
    water_depth: float
    roughness: float
    bed_slope: float


@dataclass(frozen=True)
class UChannelResult:
    water_width: float
    water_depth: float
    semicircle_radius: float
    area: float
    wetted_perimeter: float
    hydraulic_radius: float
    chezy: float
    velocity: float
    flow: float
    formula: str


def calculate_u_channel(data: UChannelInput) -> UChannelResult:
    width = _positive(data.water_width, "U形渠水面宽B")
    depth = _positive(data.water_depth, "U形渠水深h")
    n = _positive(data.roughness, "糙率n")
    i = _positive(data.bed_slope, "渠底比降i")
    if depth <= width / 2.0:
        raise ChannelHydraulicError("U形渠道水深应大于半圆半径B/2")
    area = math.pi / 2.0 * (width / 2.0) ** 2 + width * (depth - width / 2.0)
    wetted = math.pi * width / 2.0 + 2.0 * (depth - width / 2.0)
    radius = area / wetted
    chezy = _manning_chezy(n, radius)
    velocity = chezy * math.sqrt(radius * i)
    flow = area * velocity
    return UChannelResult(
        water_width=width,
        water_depth=depth,
        semicircle_radius=width / 2.0,
        area=area,
        wetted_perimeter=wetted,
        hydraulic_radius=radius,
        chezy=chezy,
        velocity=velocity,
        flow=flow,
        formula="ω=π/2*(B/2)^2+B*(h-B/2)；x=πB/2+2*(h-B/2)；Q=ω*C*sqrt(R*i)",
    )


def solve_u_channel_by_width_depth_ratio(
    target_flow: float,
    width_depth_ratio: float,
    roughness: float,
    bed_slope: float,
) -> UChannelResult:
    q = _positive(target_flow, "设计流量Q")
    ratio = _positive(width_depth_ratio, "水深与水面宽比h/B")
    if ratio <= 0.5:
        raise ChannelHydraulicError("U形渠道h/B必须大于0.5")
    n = _positive(roughness, "糙率n")
    i = _positive(bed_slope, "渠底比降i")
    low = 1e-9
    high = 0.1
    while calculate_u_channel(UChannelInput(high, ratio * high, n, i)).flow < q:
        high *= 2.0
        if high > 10000:
            raise ChannelHydraulicError("未能求得U形渠道尺寸，请检查流量、糙率或比降")
    for _ in range(90):
        mid = (low + high) / 2.0
        if calculate_u_channel(UChannelInput(mid, ratio * mid, n, i)).flow < q:
            low = mid
        else:
            high = mid
    width = (low + high) / 2.0
    return calculate_u_channel(UChannelInput(width, ratio * width, n, i))


@dataclass(frozen=True)
class RebarGroupInput:
    name: str
    diameter_mm: float
    unit_weight_kg_per_m: float
    structure_length_m: float
    protection_total_thickness_m: float
    spacing_m: float
    straight_length_m: float
    bends_count: int
    segment_count: float = 1.0
    layer_count: float = 1.0


@dataclass(frozen=True)
class RebarGroupResult:
    name: str
    hook_length_m: float
    single_bar_length_m: float
    root_count: float
    total_length_m: float
    net_weight_kg: float


@dataclass(frozen=True)
class ReinforcementInput:
    groups: tuple[RebarGroupInput, ...]
    concrete_volume_m3: float
    loss_rate_percent: float


@dataclass(frozen=True)
class ReinforcementResult:
    groups: tuple[RebarGroupResult, ...]
    total_net_weight_kg: float
    total_gross_weight_kg: float
    total_gross_weight_t: float
    net_ratio_kg_per_m3: float
    gross_ratio_kg_per_m3: float
    formula: str


def calculate_reinforcement(data: ReinforcementInput) -> ReinforcementResult:
    if not data.groups:
        raise ChannelHydraulicError("至少需要输入一组钢筋")
    concrete_volume = _positive(data.concrete_volume_m3, "混凝土体积")
    loss_rate = _non_negative(data.loss_rate_percent, "钢筋损耗率")
    results: list[RebarGroupResult] = []
    total_net = 0.0
    for group in data.groups:
        diameter_m = _positive(group.diameter_mm, f"{group.name}直径") / 1000.0
        unit_weight = _positive(group.unit_weight_kg_per_m, f"{group.name}每米重量")
        structure_length = _positive(group.structure_length_m, f"{group.name}结构长度")
        cover = _non_negative(group.protection_total_thickness_m, f"{group.name}保护层总厚")
        spacing = _positive(group.spacing_m, f"{group.name}钢筋间距")
        straight = _positive(group.straight_length_m, f"{group.name}直段长度")
        if group.bends_count < 0:
            raise ChannelHydraulicError(f"{group.name}弯钩端数不能小于0")
        segment_count = _positive(group.segment_count, f"{group.name}分段数")
        layer_count = _positive(group.layer_count, f"{group.name}层数")
        effective_length = structure_length - cover
        if effective_length <= 0:
            raise ChannelHydraulicError(f"{group.name}结构长度扣保护层后必须大于0")
        hook_length = 6.25 * diameter_m * group.bends_count
        single_bar_length = straight + hook_length
        root_count = (effective_length / spacing + 1.0) * segment_count * layer_count
        total_length = root_count * single_bar_length
        net_weight = total_length * unit_weight
        total_net += net_weight
        results.append(
            RebarGroupResult(
                name=group.name,
                hook_length_m=hook_length,
                single_bar_length_m=single_bar_length,
                root_count=root_count,
                total_length_m=total_length,
                net_weight_kg=net_weight,
            )
        )
    gross = total_net * (1.0 + loss_rate / 100.0)
    return ReinforcementResult(
        groups=tuple(results),
        total_net_weight_kg=total_net,
        total_gross_weight_kg=gross,
        total_gross_weight_t=gross / 1000.0,
        net_ratio_kg_per_m3=total_net / concrete_volume,
        gross_ratio_kg_per_m3=gross / concrete_volume,
        formula="根数=(结构长度-保护层总厚)/间距+1；弯钩长=6.25d；总长=根数*(直段长+弯钩长)；总重=净重*(1+损耗率)",
    )


@dataclass(frozen=True)
class ChannelSidewallStabilityInput:
    soil_friction_angle_deg: float
    soil_unit_weight: float
    wall_height: float
    groundwater_depth: float = 0.0
    water_unit_weight: float = 9.81
    channel_inner_water_depth: float = 0.0
    structure_factor: float = 1.2
    importance_factor: float = 1.0
    design_condition_factor: float = 1.0
    concrete_fc: float = 9.6
    wall_thickness_mm: float = 200.0
    cover_thickness_mm: float = 30.0
    steel_diameter_mm: float = 12.0
    steel_fy: float = 300.0
    steel_spacing_mm: float = 200.0
    steel_layers: int = 1


@dataclass(frozen=True)
class ChannelSidewallStabilityResult:
    earth_pressure_coefficient: float
    earth_pressure: float
    groundwater_pressure: float
    inner_water_pressure: float
    total_pressure: float
    bending_moment: float
    effective_height_mm: float
    resistance_coefficient: float
    compression_zone_height_ratio: float
    required_steel_area_mm2: float
    provided_steel_area_mm2: float
    steel_ratio_percent: float
    stable: bool
    formula: str


def calculate_channel_sidewall_stability(data: ChannelSidewallStabilityInput) -> ChannelSidewallStabilityResult:
    phi = _non_negative(data.soil_friction_angle_deg, "填土有效内摩擦角φ")
    gamma = _positive(data.soil_unit_weight, "填土容重γ")
    height = _positive(data.wall_height, "渠道侧墙高度H")
    groundwater_depth = _non_negative(data.groundwater_depth, "墙后水深Hw")
    water_gamma = _positive(data.water_unit_weight, "水容重γw")
    inner_depth = _non_negative(data.channel_inner_water_depth, "渠道内水深H内")
    yd = _positive(data.structure_factor, "结构系数γd")
    y0 = _positive(data.importance_factor, "结构重要性系数γ0")
    psi = _positive(data.design_condition_factor, "设计状况系数ψ")
    fc = _positive(data.concrete_fc, "混凝土轴心抗压强度fc")
    wall_thickness = _positive(data.wall_thickness_mm, "墙厚h")
    cover = _non_negative(data.cover_thickness_mm, "保护层厚度c")
    diameter = _positive(data.steel_diameter_mm, "钢筋直径d")
    fy = _positive(data.steel_fy, "钢筋抗拉强度fy")
    spacing = _positive(data.steel_spacing_mm, "钢筋间距")
    layers = int(_positive(float(data.steel_layers), "钢筋层数"))

    k = 1.0 - math.sin(math.radians(phi))
    earth_pressure = 0.5 * k * gamma * height * height
    groundwater_pressure = 0.5 * water_gamma * groundwater_depth * groundwater_depth
    inner_water_pressure = 0.5 * water_gamma * inner_depth * inner_depth
    total_pressure = earth_pressure + groundwater_pressure - inner_water_pressure
    if total_pressure < 0:
        total_pressure = 0.0
    moment = total_pressure * height / 3.0
    effective_height = wall_thickness - (cover + diameter / 2.0)
    if effective_height <= 0:
        raise ChannelHydraulicError("墙厚扣保护层和半径后必须大于0")
    combined_factor = yd * y0 * psi
    resistance_coefficient = combined_factor * moment * 1_000_000.0 / (fc * 1000.0 * effective_height * effective_height)
    if resistance_coefficient >= 0.5:
        raise ChannelHydraulicError("截面抵抗矩系数过大，需增大墙厚或调整配筋")
    epsilon = 1.0 - math.sqrt(1.0 - 2.0 * resistance_coefficient)
    required_area = fc * epsilon * 1000.0 * effective_height / fy
    provided_area = math.pi * (diameter / 2.0) ** 2 * (1000.0 / spacing) * layers
    steel_ratio = required_area / (1000.0 * effective_height) * 100.0
    stable = provided_area >= required_area and steel_ratio >= 0.15 and epsilon < 0.544
    return ChannelSidewallStabilityResult(
        earth_pressure_coefficient=k,
        earth_pressure=earth_pressure,
        groundwater_pressure=groundwater_pressure,
        inner_water_pressure=inner_water_pressure,
        total_pressure=total_pressure,
        bending_moment=moment,
        effective_height_mm=effective_height,
        resistance_coefficient=resistance_coefficient,
        compression_zone_height_ratio=epsilon,
        required_steel_area_mm2=required_area,
        provided_steel_area_mm2=provided_area,
        steel_ratio_percent=steel_ratio,
        stable=stable,
        formula="E=0.5*K*γ*H^2；K=1-sinφ；M=ΣE*H/3；as=γM/(fc*b*h0^2)；As=fc*ε*b*h0/fy",
    )


@dataclass(frozen=True)
class CompoundChannelInput:
    total_depth: float
    main_depth: float
    main_bottom_width: float
    main_side_slope: float
    floodplain_side_slope: float
    left_floodplain_width: float
    right_floodplain_width: float
    main_roughness: float
    floodplain_roughness: float
    bed_slope: float


@dataclass(frozen=True)
class CompoundChannelResult:
    total_area: float
    total_wetted_perimeter: float
    composite_roughness: float
    integral_flow: float
    main_area: float
    floodplain_area: float
    main_flow: float
    floodplain_flow: float
    divided_flow: float
    formula: str


def calculate_compound_channel(data: CompoundChannelInput) -> CompoundChannelResult:
    h = _positive(data.total_depth, "渠底至水面水深h")
    h1 = _positive(data.main_depth, "渠底至滩面水深h1")
    if h <= h1:
        raise ChannelHydraulicError("复式断面总水深h必须大于深槽水深h1")
    b1 = _positive(data.main_bottom_width, "深槽底宽b1")
    m1 = _non_negative(data.main_side_slope, "深槽边坡m1")
    m2 = _non_negative(data.floodplain_side_slope, "滩槽边坡m2")
    b2l = _non_negative(data.left_floodplain_width, "左滩底宽")
    b2r = _non_negative(data.right_floodplain_width, "右滩底宽")
    n1 = _positive(data.main_roughness, "深槽糙率n1")
    n2 = _positive(data.floodplain_roughness, "滩槽糙率n2")
    i = _positive(data.bed_slope, "渠底比降i")

    over_depth = h - h1
    total_area = (b1 + m1 * h1) * h1 + ((b1 + 2.0 * m1 * h1) + (b2l + b2r) + m2 * over_depth) * over_depth
    total_wetted = b1 + 2.0 * h1 * math.sqrt(1.0 + m1 * m1) + b2l + b2r + 2.0 * over_depth * math.sqrt(1.0 + m2 * m2)
    main_area = (b1 + m1 * h1) * h1 + (b1 + 2.0 * m1 * h1) * over_depth
    floodplain_area = (b2l + b2r + m2 * over_depth) * over_depth
    main_wetted = b1 + 2.0 * h1 * math.sqrt(1.0 + m1 * m1)
    floodplain_wetted = b2l + b2r + 2.0 * over_depth * math.sqrt(1.0 + m2 * m2)
    composite_n = math.sqrt((main_wetted * n1 * n1 + floodplain_wetted * n2 * n2) / (main_wetted + floodplain_wetted))
    total_radius = total_area / total_wetted
    integral_flow = total_area * (total_radius ** (1.0 / 6.0) / composite_n) * math.sqrt(total_radius * i)
    r1 = main_area / main_wetted
    r2 = floodplain_area / floodplain_wetted
    q1 = main_area * (r1 ** (1.0 / 6.0) / n1) * math.sqrt(r1 * i)
    q2 = floodplain_area * (r2 ** (1.0 / 6.0) / n2) * math.sqrt(r2 * i)
    return CompoundChannelResult(
        total_area=total_area,
        total_wetted_perimeter=total_wetted,
        composite_roughness=composite_n,
        integral_flow=integral_flow,
        main_area=main_area,
        floodplain_area=floodplain_area,
        main_flow=q1,
        floodplain_flow=q2,
        divided_flow=q1 + q2,
        formula="整体：ω、x按复式断面整体求，n综=sqrt((x1*n1^2+x2*n2^2)/(x1+x2))；分块：Q=Q1+Q2",
    )


def solve_compound_total_depth_for_flow(
    target_flow: float,
    data: CompoundChannelInput,
    method: str = "integral",
) -> float:
    q = _positive(target_flow, "目标流量Q")
    if method not in {"integral", "divided"}:
        raise ChannelHydraulicError("复式断面求解方法必须为integral或divided")
    low = data.main_depth + 1e-9
    high = max(data.total_depth, data.main_depth + 0.1)
    def flow_at(depth: float) -> float:
        result = calculate_compound_channel(
            CompoundChannelInput(
                total_depth=depth,
                main_depth=data.main_depth,
                main_bottom_width=data.main_bottom_width,
                main_side_slope=data.main_side_slope,
                floodplain_side_slope=data.floodplain_side_slope,
                left_floodplain_width=data.left_floodplain_width,
                right_floodplain_width=data.right_floodplain_width,
                main_roughness=data.main_roughness,
                floodplain_roughness=data.floodplain_roughness,
                bed_slope=data.bed_slope,
            )
        )
        return result.integral_flow if method == "integral" else result.divided_flow
    while flow_at(high) < q:
        high *= 2.0
        if high > 10000:
            raise ChannelHydraulicError("未能求得复式断面水深，请检查流量或断面参数")
    for _ in range(90):
        mid = (low + high) / 2.0
        if flow_at(mid) < q:
            low = mid
        else:
            high = mid
    return (low + high) / 2.0


@dataclass(frozen=True)
class LongitudinalPointInput:
    station_no: int
    ground_elevation: float
    distance_from_previous: float = 0.0


@dataclass(frozen=True)
class LongitudinalDesignInput:
    bed_slope: float
    design_depth: float
    minimum_depth: float
    enlarged_depth: float
    freeboard: float
    first_bed_elevation: float
    points: tuple[LongitudinalPointInput, ...]


@dataclass(frozen=True)
class LongitudinalPointResult:
    station_no: int
    ground_elevation: float
    bed_elevation: float
    minimum_water_level: float
    design_water_level: float
    enlarged_water_level: float
    top_elevation: float
    excavation_depth: float
    fill_height: float


@dataclass(frozen=True)
class LongitudinalDesignResult:
    points: tuple[LongitudinalPointResult, ...]
    formula: str


def calculate_channel_longitudinal_design(data: LongitudinalDesignInput) -> LongitudinalDesignResult:
    slope = _positive(data.bed_slope, "渠底比降i")
    design_depth = _positive(data.design_depth, "设计水深")
    minimum_depth = _non_negative(data.minimum_depth, "最小水深")
    enlarged_depth = _positive(data.enlarged_depth, "加大水深")
    freeboard = _non_negative(data.freeboard, "超高")
    bed = data.first_bed_elevation
    if not data.points:
        raise ChannelHydraulicError("至少需要一个桩号点")
    results: list[LongitudinalPointResult] = []
    for index, point in enumerate(data.points):
        distance = 0.0 if index == 0 else _non_negative(point.distance_from_previous, "相邻桩距")
        if index > 0:
            bed -= distance * slope
        min_level = bed + minimum_depth
        design_level = bed + design_depth
        enlarged_level = bed + enlarged_depth
        top_level = enlarged_level + freeboard
        excavation = max(0.0, point.ground_elevation - bed)
        fill = max(0.0, top_level - point.ground_elevation)
        results.append(
            LongitudinalPointResult(
                station_no=point.station_no,
                ground_elevation=point.ground_elevation,
                bed_elevation=bed,
                minimum_water_level=min_level,
                design_water_level=design_level,
                enlarged_water_level=enlarged_level,
                top_elevation=top_level,
                excavation_depth=excavation,
                fill_height=fill,
            )
        )
    return LongitudinalDesignResult(
        points=tuple(results),
        formula="设计渠底高程=设计水位-设计水深；最低/设计/加大水位=渠底高程+相应水深；渠顶高程=加大水位+超高；下一桩渠底=上一桩渠底-桩距*i",
    )
