from __future__ import annotations

from dataclasses import dataclass
import math

from .core import HydrologyError


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


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


def _manning_rectangular_flow(width: float, depth: float, roughness: float, slope: float) -> float:
    area = width * depth
    wetted = width + 2.0 * depth
    radius = area / wetted
    return area * radius ** (1.0 / 6.0) / roughness * math.sqrt(radius * slope)


def _manning_trapezoid_flow(width: float, depth: float, side_slope: float, roughness: float, slope: float) -> float:
    area = (width + side_slope * depth) * depth
    wetted = width + 2.0 * depth * math.sqrt(1.0 + side_slope**2)
    radius = area / wetted
    return area * radius ** (1.0 / 6.0) / roughness * math.sqrt(radius * slope)


def _solve_depth_for_flow(flow_func, target_flow: float) -> float:
    q = _positive(target_flow, "目标流量Q")
    low = 1e-9
    high = 0.1
    while flow_func(high) < q:
        high *= 2.0
        if high > 10000:
            raise HydrologyError("未能求得水深，请检查流量、糙率或比降")
    for _ in range(100):
        mid = (low + high) / 2.0
        if flow_func(mid) < q:
            low = mid
        else:
            high = mid
    return (low + high) / 2.0


@dataclass(frozen=True)
class SluiceHydraulicInput:
    flow_m3s: float
    upstream_head_m: float
    downstream_head_m: float
    gate_design_head_m: float
    freeboard_m: float
    discharge_coefficient: float = 0.40
    contraction_coefficient: float = 1.0
    submergence_coefficient: float = 1.0
    opening_count: float = 1.0
    chamber_coefficient: float = 1.0
    groove_depth_ratio: float = 1.4
    gravity_m_s2: float = 9.81


@dataclass(frozen=True)
class SluiceHydraulicResult:
    total_clear_width_m: float
    chamber_length_m: float
    groove_width_cm: float
    groove_depth_cm: float
    formula: str


def calculate_sluice_hydraulic(data: SluiceHydraulicInput) -> SluiceHydraulicResult:
    q = _positive(data.flow_m3s, "过闸流量Q")
    h0 = _positive(data.upstream_head_m, "上游总水头H0")
    h_design = _positive(data.gate_design_head_m, "闸孔设计水头H")
    downstream = _non_negative(data.downstream_head_m, "下游水头")
    freeboard = _non_negative(data.freeboard_m, "闸室安全加高a")
    m = _positive(data.discharge_coefficient, "流量系数m")
    eps = _positive(data.contraction_coefficient, "侧收缩系数epsilon")
    sigma = _positive(data.submergence_coefficient, "淹没系数sigma")
    n = _positive(data.opening_count, "孔数n")
    k = _positive(data.chamber_coefficient, "闸室长度系数K")
    g = _positive(data.gravity_m_s2, "重力加速度g")
    width = q / (sigma * eps * m * n * math.sqrt(2.0 * g) * h0 ** 1.5)
    chamber_length = (h0 + freeboard) * (1.0 + 0.1 * abs(h0 - downstream)) * k
    groove_width = 3.5 * h_design * math.sqrt(width) + 5.0
    groove_depth = groove_width / _positive(data.groove_depth_ratio, "槽深折减比")
    return SluiceHydraulicResult(
        width,
        chamber_length,
        groove_width,
        groove_depth,
        "Q=sigma*epsilon*m*n*b*sqrt(2g)*H0^(3/2); L=(H0+a)(1+0.1*DeltaH)K; W=3.5Hsqrt(b0)+5; D=W/k",
    )


@dataclass(frozen=True)
class SpillwayChuteInput:
    flow_m3s: float
    bottom_width_m: float
    bed_slope: float
    roughness: float
    side_slope: float = 0.0
    aeration_coefficient: float = 1.0
    gravity_m_s2: float = 9.81


@dataclass(frozen=True)
class SpillwayChuteResult:
    normal_depth_m: float
    velocity_m_s: float
    aerated_depth_m: float
    freeboard_m: float
    sidewall_height_m: float
    formula: str


def calculate_spillway_chute(data: SpillwayChuteInput) -> SpillwayChuteResult:
    q = _positive(data.flow_m3s, "泄槽流量Q")
    b = _positive(data.bottom_width_m, "泄槽底宽B")
    i = _positive(data.bed_slope, "泄槽底坡i")
    n = _positive(data.roughness, "糙率n")
    m = _non_negative(data.side_slope, "边坡系数m")
    if m == 0:
        depth = _solve_depth_for_flow(lambda h: _manning_rectangular_flow(b, h, n, i), q)
        area = b * depth
    else:
        depth = _solve_depth_for_flow(lambda h: _manning_trapezoid_flow(b, h, m, n, i), q)
        area = (b + m * depth) * depth
    velocity = q / area
    aerated = (1.0 + _non_negative(data.aeration_coefficient, "掺气系数zeta") * velocity / 100.0) * depth
    freeboard = 0.61 + 0.0372 * velocity * depth ** (1.0 / 3.0)
    return SpillwayChuteResult(
        depth,
        velocity,
        aerated,
        freeboard,
        aerated + freeboard,
        "Q=A*C*sqrt(R*i); hb=(1+zeta*v/100)h; Delta=0.61+0.0372*v*h^(1/3); H=hb+Delta",
    )


@dataclass(frozen=True)
class StillingBasinInput:
    flow_m3s: float
    basin_width_m: float
    upstream_depth_m: float
    approach_velocity_m_s: float
    drop_height_m: float
    downstream_depth_m: float
    slope_ratio: float = 5.0
    velocity_coefficient: float = 0.95
    safety_submergence: float = 1.10
    slab_coefficient: float = 0.141
    apron_coefficient: float = 10.33
    allowed_velocity_m_s: float = 3.0
    gravity_m_s2: float = 9.81


@dataclass(frozen=True)
class StillingBasinResult:
    unit_discharge_m2s: float
    first_conjugate_depth_m: float
    second_conjugate_depth_m: float
    basin_depth_m: float
    basin_length_m: float
    basin_width_m: float
    sidewall_height_m: float
    slab_thickness_m: float
    apron_length_m: float
    trench_depth_m: float
    trench_length_m: float
    formula: str


def _solve_first_conjugate_depth(unit_q: float, total_energy: float, phi: float, gravity: float) -> float:
    c = unit_q**2 / (2.0 * gravity * phi**2)

    def f(depth: float) -> float:
        return depth**3 - total_energy * depth**2 + c

    roots: list[float] = []
    prev = 1e-8
    fp = f(prev)
    step = max(total_energy / 20000.0, 0.0005)
    x = step
    while x <= total_energy * 1.5 + 1.0:
        fx = f(x)
        if fp * fx < 0:
            lo, hi = prev, x
            for _ in range(80):
                mid = (lo + hi) / 2.0
                if f(lo) * f(mid) <= 0:
                    hi = mid
                else:
                    lo = mid
            roots.append((lo + hi) / 2.0)
        prev, fp = x, fx
        x += step
    if not roots:
        raise HydrologyError("未能求得跃前水深，请检查消力池参数")
    return min(roots)


def calculate_stilling_basin(data: StillingBasinInput) -> StillingBasinResult:
    q = _positive(data.flow_m3s, "流量Q")
    b = _positive(data.basin_width_m, "消力池宽度B")
    unit_q = q / b
    g = _positive(data.gravity_m_s2, "重力加速度g")
    t0 = (
        _positive(data.upstream_depth_m, "上游水深H")
        + _positive(data.approach_velocity_m_s, "行近流速v") ** 2 / (2.0 * g)
        + _positive(data.drop_height_m, "跌坎高度P")
    )
    phi = _positive(data.velocity_coefficient, "流速系数phi")
    h1 = _solve_first_conjugate_depth(unit_q, t0, phi, g)
    h2 = 0.5 * h1 * (math.sqrt(1.0 + 8.0 * unit_q**2 / (g * h1**3)) - 1.0)
    ht = _positive(data.downstream_depth_m, "下游水深ht")
    delta_z = unit_q**2 / (2.0 * g * phi**2 * ht**2) - unit_q**2 / (2.0 * g * h2**2)
    basin_depth = _positive(data.safety_submergence, "淹没安全系数sigma0") * h2 - ht - delta_z
    hydraulic_jump_length = 6.9 * (h2 - h1)
    basin_length = 0.8 * hydraulic_jump_length
    approach_length = _positive(data.slope_ratio, "斜坡坡比Z") * basin_depth
    sidewall_height = h2 + 0.5
    slab = _positive(data.slab_coefficient, "底板厚度系数k") * math.sqrt(unit_q * math.sqrt(data.drop_height_m))
    apron = _positive(data.apron_coefficient, "海漫长度系数ks") * math.sqrt(unit_q * math.sqrt(data.drop_height_m))
    trench = max(0.0, 1.1 * unit_q / _positive(data.allowed_velocity_m_s, "允许流速[v]") - ht)
    return StillingBasinResult(
        unit_q,
        h1,
        h2,
        basin_depth,
        basin_length,
        b + 2.0 * approach_length * math.tan(math.radians(6.0)),
        sidewall_height,
        slab,
        apron,
        trench,
        2.0 * trench,
        "h2/h1=0.5*(sqrt(1+8q^2/(g*h1^3))-1); h1^3-T0*h1^2+q^2/(2gphi^2)=0; d=sigma0*h2-ht-DeltaZ; L=0.8*6.9(h2-h1)",
    )


@dataclass(frozen=True)
class SteepSlopeInput:
    flow_m3s: float
    entrance_head_m: float
    added_head_m: float
    chute_slope: float
    roughness: float
    downstream_depth_m: float
    drop_height_m: float
    entrance_coefficient: float = 0.385
    contraction_coefficient: float = 1.0
    entrance_side_slope: float = 0.0


@dataclass(frozen=True)
class SteepSlopeResult:
    entrance_bottom_width_m: float
    transition_length_m: float
    chute_depth_m: float
    chute_velocity_m_s: float
    aerated_depth_m: float
    basin_depth_m: float
    basin_length_m: float
    basin_width_m: float
    slab_thickness_m: float
    sidewall_height_m: float
    formula: str


def calculate_steep_slope(data: SteepSlopeInput) -> SteepSlopeResult:
    q = _positive(data.flow_m3s, "流量Q")
    h0 = _positive(data.entrance_head_m, "进口水头h0")
    m0 = _positive(data.entrance_coefficient, "进口流量系数m")
    eps = _positive(data.contraction_coefficient, "侧收缩系数epsilon")
    nc = _non_negative(data.entrance_side_slope, "进口边坡系数nc")
    # Rectangular inlet is nc=0; for trapezoid the solved value is the bottom width bc.
    width = q / (eps * m0 * math.sqrt(2.0 * 9.81) * h0 ** 1.5) - 0.8 * nc * h0
    transition = 3.0 * _positive(data.added_head_m, "加大水深h加")
    chute = calculate_spillway_chute(
        SpillwayChuteInput(q, width, _positive(data.chute_slope, "陡坡坡降i"), _positive(data.roughness, "糙率n"))
    )
    basin = calculate_stilling_basin(
        StillingBasinInput(
            flow_m3s=q,
            basin_width_m=width,
            upstream_depth_m=h0,
            approach_velocity_m_s=q / (width * h0),
            drop_height_m=_positive(data.drop_height_m, "落差P"),
            downstream_depth_m=_positive(data.downstream_depth_m, "下游水深ht"),
        )
    )
    return SteepSlopeResult(
        width,
        transition,
        chute.normal_depth_m,
        chute.velocity_m_s,
        chute.aerated_depth_m,
        basin.basin_depth_m,
        basin.basin_length_m,
        basin.basin_width_m,
        basin.slab_thickness_m,
        basin.sidewall_height_m,
        "入口Q=epsilon*(bc+0.8nc*hc)*m*sqrt(2g)*h0^(3/2); L渐=3h加; 陡坡按曼宁公式，消力池按水跃公式联算",
    )


@dataclass(frozen=True)
class VerticalDropInput:
    flow_m3s: float
    downstream_channel_width_m: float
    entrance_head_m: float
    added_head_m: float
    drop_height_m: float
    entrance_coefficient: float = 0.350
    contraction_coefficient: float = 1.0


@dataclass(frozen=True)
class VerticalDropResult:
    entrance_bottom_width_m: float
    transition_length_m: float
    unit_discharge_m2s: float
    drop_flow_parameter: float
    basin_length_m: float
    sill_height_m: float
    formula: str


def calculate_vertical_drop(data: VerticalDropInput) -> VerticalDropResult:
    q = _positive(data.flow_m3s, "流量Q")
    h0 = _positive(data.entrance_head_m, "进口水头h0")
    width = q / (_positive(data.contraction_coefficient, "侧收缩系数epsilon") * _positive(data.entrance_coefficient, "流量系数m") * math.sqrt(2.0 * 9.81) * h0 ** 1.5)
    unit_q = q / width
    p = _positive(data.drop_height_m, "跌差P")
    d_param = unit_q**2 / (9.81 * p**3)
    hc0 = 1.66 * d_param**0.27 * p
    hc = 0.54 * d_param**0.425 * p
    ld = 4.3 * d_param**0.27 * p
    lj = 1.9 * hc0 - hc
    length = ld + 0.8 * lj
    sill = max(0.0, hc0 - _positive(data.downstream_channel_width_m, "下游渠道底宽") * 0.0 - (h0 - 0.76))
    return VerticalDropResult(
        width,
        3.0 * _positive(data.added_head_m, "加大水深h加"),
        unit_q,
        d_param,
        length,
        sill,
        "D=q^2/(gP^3); ld=4.3D^0.27P; hc0=1.66D^0.27P; hc=0.54D^0.425P; ls=ld+0.8(1.9hc0-hc)",
    )


@dataclass(frozen=True)
class AqueductInput:
    flow_m3s: float
    trough_length_m: float
    bed_slope: float
    roughness: float
    upstream_channel_width_m: float
    upstream_channel_depth_m: float
    upstream_velocity_m_s: float
    upstream_bed_elevation_m: float
    trough_bottom_width_m: float
    inlet_loss_coefficient: float = 0.20
    outlet_recovery_ratio: float = 1.0 / 3.0
    inlet_transition_factor: float = 2.0
    outlet_transition_factor: float = 3.0


@dataclass(frozen=True)
class AqueductResult:
    trough_depth_m: float
    trough_velocity_m_s: float
    inlet_water_drop_m: float
    trough_drop_m: float
    outlet_recovery_m: float
    total_head_loss_m: float
    inlet_transition_length_m: float
    outlet_transition_length_m: float
    trough_start_bed_m: float
    trough_end_bed_m: float
    downstream_bed_m: float
    freeboard_cm: float
    formula: str


def calculate_aqueduct(data: AqueductInput) -> AqueductResult:
    q = _positive(data.flow_m3s, "渡槽流量Q")
    b = _positive(data.trough_bottom_width_m, "槽身底宽B")
    depth = _solve_depth_for_flow(
        lambda h: _manning_rectangular_flow(b, h, _positive(data.roughness, "糙率n"), _positive(data.bed_slope, "槽底比降i")),
        q,
    )
    velocity = q / (b * depth)
    v1 = _positive(data.upstream_velocity_m_s, "上游渠道流速v1")
    inlet_drop = (1.0 + _non_negative(data.inlet_loss_coefficient, "进口损失系数zeta")) * max(velocity**2 - v1**2, 0.0) / (2.0 * 9.81)
    trough_drop = _positive(data.trough_length_m, "槽身长度L") * data.bed_slope
    outlet_recovery = _positive(data.outlet_recovery_ratio, "出口水头恢复比例") * inlet_drop
    total = inlet_drop + trough_drop - outlet_recovery
    width_delta = max(_positive(data.upstream_channel_width_m, "上游渠底宽") - b, 0.0)
    start_bed = _positive(data.upstream_bed_elevation_m, "上游渠底高程") + _positive(data.upstream_channel_depth_m, "上游渠水深") - inlet_drop - depth
    end_bed = start_bed - trough_drop
    downstream_bed = data.upstream_bed_elevation_m - total
    freeboard_cm = depth / 12.0 * 100.0 + 5.0
    return AqueductResult(
        depth,
        velocity,
        inlet_drop,
        trough_drop,
        outlet_recovery,
        total,
        data.inlet_transition_factor * width_delta,
        data.outlet_transition_factor * width_delta,
        start_bed,
        end_bed,
        downstream_bed,
        freeboard_cm,
        "DeltaZ1=(1+zeta)(v2^2-v1^2)/(2g); DeltaZ2=iL; DeltaZ3=DeltaZ1/3; DeltaZ=DeltaZ1+DeltaZ2-DeltaZ3",
    )


@dataclass(frozen=True)
class InvertedSiphonInput:
    flow_m3s: float
    design_velocity_m_s: float
    pipe_length_m: float
    headloss_coefficient_sum: float
    friction_factor: float
    upstream_bed_m: float
    upstream_depth_m: float
    downstream_depth_m: float
    allowable_stress_kg_cm2: float
    design_head_m: float


@dataclass(frozen=True)
class InvertedSiphonResult:
    pipe_diameter_m: float
    total_headloss_m: float
    wall_thickness_mm: float
    minimum_wall_thickness_mm: float
    downstream_bed_m: float
    inlet_basin_length_m: float
    inlet_basin_width_m: float
    outlet_basin_length_m: float
    outlet_basin_width_m: float
    bend_radius_m: float
    formula: str


def calculate_inverted_siphon(data: InvertedSiphonInput) -> InvertedSiphonResult:
    q = _positive(data.flow_m3s, "倒虹吸流量Q")
    v = _positive(data.design_velocity_m_s, "管内设计流速v")
    diameter = math.sqrt(4.0 * q / (math.pi * v))
    headloss = (_non_negative(data.headloss_coefficient_sum, "局部损失系数和") + _positive(data.friction_factor, "沿程阻力系数lambda") * _positive(data.pipe_length_m, "管长L") / diameter) * v**2 / (2.0 * 9.81)
    thickness_mm = 1000.0 * _positive(data.design_head_m, "设计水头Hp") * diameter / (2.0 * _positive(data.allowable_stress_kg_cm2, "允许应力"))
    min_thickness = max(6.0, diameter * 1000.0 / 130.0)
    downstream_bed = data.upstream_bed_m + data.upstream_depth_m - data.downstream_depth_m - headloss
    return InvertedSiphonResult(
        diameter,
        headloss,
        thickness_mm,
        min_thickness,
        downstream_bed,
        10.0,
        1.5,
        12.0,
        2.0,
        4.0 * diameter,
        "D=sqrt(4Q/(pi*v)); hw=(sum(zeta)+lambda*L/D)*v^2/(2g); delta=10*Hp*D/(2*[sigma]); Hd=Hu+hu-hd-hw",
    )


@dataclass(frozen=True)
class BridgeHydraulicInput:
    flow_m3s: float
    upstream_normal_depth_m: float
    safety_freeboard_m: float
    water_level_difference_m: float
    shape_width_coefficient: float
    pier_coefficient: float
    bridge_deck_thickness_m: float = 0.0
    gravity_m_s2: float = 9.81


@dataclass(frozen=True)
class BridgeHydraulicResult:
    flow_area_m2: float
    total_clear_width_m: float
    design_clear_height_m: float
    single_same_section_flow_m3s: float
    single_different_section_flow_m3s: float
    multi_same_section_flow_m3s: float
    multi_different_section_flow_m3s: float
    formula: str


def calculate_bridge_hydraulic(data: BridgeHydraulicInput) -> BridgeHydraulicResult:
    q = _positive(data.flow_m3s, "过桥流量Q")
    h0 = _positive(data.upstream_normal_depth_m, "上游渠道正常水深h0")
    dz = _positive(data.water_level_difference_m, "桥上游壅高水头DeltaZ")
    mu = _positive(data.shape_width_coefficient, "进口出口流量系数mu1") * _positive(data.pier_coefficient, "桥墩形状流量系数mu2")
    omega = q / (mu * math.sqrt(2.0 * _positive(data.gravity_m_s2, "重力加速度g") * dz))
    width = omega / h0
    height = h0 + _non_negative(data.safety_freeboard_m, "安全超高") + _non_negative(data.bridge_deck_thickness_m, "桥面厚度")
    single_same = mu * omega * math.sqrt(2.0 * data.gravity_m_s2 * dz)
    return BridgeHydraulicResult(
        omega,
        width,
        height,
        single_same,
        q,
        single_same * 1.111,
        q,
        "Q=mu*omega*sqrt(2g*DeltaZ0); omega=Bh0; H桥=h0+安全超高+桥面厚",
    )


@dataclass(frozen=True)
class TunnelHydraulicInput:
    mode: str
    flow_m3s: float
    upstream_head_m: float
    coefficient: float
    submergence_coefficient: float = 1.0
    height_width_ratio: float = 1.5
    arch_angle_deg: float = 120.0
    length_m: float = 300.0
    bed_slope: float = 0.001
    roughness: float = 0.015
    pressure_section_area_m2: float | None = None
    pressure_head_correction_m: float = 0.0
    gravity_m_s2: float = 9.81


@dataclass(frozen=True)
class TunnelHydraulicResult:
    bottom_width_m: float
    tunnel_height_m: float
    arch_radius_m: float
    circular_diameter_m: float
    short_tunnel_flow_m3s: float
    long_tunnel_flow_m3s: float
    formula: str


def calculate_tunnel_hydraulic(data: TunnelHydraulicInput) -> TunnelHydraulicResult:
    mode = data.mode.strip()
    q = _positive(data.flow_m3s, "过洞流量Q")
    h0 = _positive(data.upstream_head_m, "上游总水头H0")
    g = _positive(data.gravity_m_s2, "重力加速度g")
    if mode == "无压城门洞形":
        m = _positive(data.coefficient, "流量系数m")
        sigma = _positive(data.submergence_coefficient, "淹没系数sigma")
        width = q / (m * sigma * math.sqrt(2.0 * g) * h0 ** 1.5)
        height = _positive(data.height_width_ratio, "洞高宽比") * width
        radius = 0.5 * width / math.sin(math.radians(_positive(data.arch_angle_deg, "拱中心角") / 2.0))
        return TunnelHydraulicResult(width, height, radius, 0.0, q, q, "Q=m*sigma*b*sqrt(2g)*H0^(3/2); r=0.5b/sin(alpha/2)")
    if mode == "无压圆形":
        d = 1.57 * (_positive(data.roughness, "糙率n") * q / math.sqrt(_positive(data.bed_slope, "洞底比降i"))) ** (3.0 / 8.0)
        return TunnelHydraulicResult(0.0, 0.0, 0.0, d, q, q, "d=1.57*(nQ/sqrt(i))^(3/8)")
    if mode == "有压矩形":
        area = data.pressure_section_area_m2
        if area is None:
            width = math.sqrt(q / (_positive(data.coefficient, "流量系数mu") * math.sqrt(2.0 * g * (h0 - data.pressure_head_correction_m))) / data.height_width_ratio)
            area = width * data.height_width_ratio * width
        else:
            width = math.sqrt(area / data.height_width_ratio)
        height = data.height_width_ratio * width
        pressure_flow = data.coefficient * area * math.sqrt(2.0 * g * (h0 - data.pressure_head_correction_m))
        return TunnelHydraulicResult(width, height, 0.0, 0.0, pressure_flow, pressure_flow, "Q=mu*omega*sqrt(2g(T0-hp))")
    raise HydrologyError("隧洞类型必须为：无压城门洞形、无压圆形、有压矩形")


@dataclass(frozen=True)
class CulvertHydraulicInput:
    flow_m3s: float
    upstream_head_m: float
    coefficient: float
    submergence_coefficient: float = 0.98
    height_width_ratio: float = 1.5
    mode: str = "无压箱涵"
    culvert_depth_m: float = 1.0
    correction_coefficient: float = 0.748
    length_m: float = 200.0
    bed_slope: float = 0.001
    gravity_m_s2: float = 9.81


@dataclass(frozen=True)
class CulvertHydraulicResult:
    bottom_width_m: float
    culvert_height_m: float
    inner_depth_m: float
    free_flow_m3s: float
    submerged_flow_m3s: float
    formula: str


def calculate_culvert_hydraulic(data: CulvertHydraulicInput) -> CulvertHydraulicResult:
    q = _positive(data.flow_m3s, "过涵流量Q")
    h0 = _positive(data.upstream_head_m, "涵洞上游总水头H0")
    g = _positive(data.gravity_m_s2, "重力加速度g")
    mode = data.mode.strip()
    if mode == "无压箱涵":
        width = q / (_positive(data.coefficient, "流量系数m") * _positive(data.submergence_coefficient, "淹没系数sigma") * math.sqrt(2.0 * g) * h0 ** 1.5)
        height = _positive(data.height_width_ratio, "涵洞高宽比") * width
        return CulvertHydraulicResult(width, height, 0.0, q, q * data.submergence_coefficient, "Q=m*sigma*b*sqrt(2g)*H0^(3/2)")
    if mode == "有压箱涵":
        depth = _positive(data.culvert_depth_m, "涵内水深h")
        width = q / (_positive(data.coefficient, "流量系数mu") * depth * math.sqrt(2.0 * g * (h0 + _positive(data.bed_slope, "涵底坡降i") * _positive(data.length_m, "涵长l") - _positive(data.correction_coefficient, "修正收缩水深系数n") * _positive(data.height_width_ratio, "非淹没时涵洞高度a"))))
        return CulvertHydraulicResult(width, data.height_width_ratio, depth, q, q, "Q=mu*omega*sqrt(2g(H0+il-na))")
    raise HydrologyError("涵洞类型必须为：无压箱涵、有压箱涵")


@dataclass(frozen=True)
class SideWeirInput:
    mode: str
    flow_m3s: float
    weir_width_m: float | None
    discharge_coefficient: float
    approach_velocity_m_s: float = 0.0
    approach_depth_m: float = 1.0
    diversion_angle_deg: float = 90.0
    upstream_weir_head_m: float = 1.0
    downstream_weir_head_m: float = 1.0
    normal_weir_coefficient: float = 0.428
    gravity_m_s2: float = 9.81


@dataclass(frozen=True)
class SideWeirResult:
    mean_weir_head_m: float
    side_weir_width_m: float
    side_weir_flow_m3s: float
    formula: str


def calculate_side_weir(data: SideWeirInput) -> SideWeirResult:
    q = _positive(data.flow_m3s, "侧堰流量Q")
    g = _positive(data.gravity_m_s2, "重力加速度g")
    mode = data.mode.strip()
    if mode == "锐角分水":
        coeff = _positive(data.discharge_coefficient, "流量系数m")
        v1 = _non_negative(data.approach_velocity_m_s, "侧堰首端平均流速v1")
        h1 = _positive(data.approach_depth_m, "侧堰首端水深h1")
        angle = math.radians(_positive(data.diversion_angle_deg, "分洪或引水角度alpha"))
        head = _positive(data.upstream_weir_head_m, "侧堰首端堰顶水头H1")
        factor = 1.0 - v1 / math.sqrt(g * h1) * math.sin(angle)
        width = data.weir_width_m or q / (coeff * factor * math.sqrt(2.0 * g) * head ** 1.5)
        flow = coeff * factor * width * math.sqrt(2.0 * g) * head ** 1.5
        return SideWeirResult(head, width, flow, "Q=m(1-v1/sqrt(g*h1)*sin(alpha))*b*sqrt(2g)*H1^(3/2)")
    if mode == "直角分水":
        head = (_positive(data.upstream_weir_head_m, "侧堰首端堰顶水头H1") + _positive(data.downstream_weir_head_m, "侧堰末端堰顶水头H2")) / 2.0
        coeff = _positive(data.normal_weir_coefficient, "正堰流量系数m")
        c = _positive(data.discharge_coefficient, "侧堰流量系数C")
        width = data.weir_width_m or q / (c * math.sqrt(2.0 * g) * head ** 1.5)
        flow = c * width * math.sqrt(2.0 * g) * head ** 1.5
        return SideWeirResult(head, width, flow, "Q=C*b*sqrt(2g)*H^(3/2); H=(H1+H2)/2")
    raise HydrologyError("侧堰类型必须为：锐角分水、直角分水")
