"""
铁塔计算模块 (Tower Calculation Module)
依据《架空输电线路杆塔结构设计技术规定》(DL/T 5154-2012)
和《建筑结构荷载规范》(GB 50009-2012) 进行铁塔结构计算。
"""

import math

MATERIAL_PROPS = {
    "Q235": {"fy": 235, "f": 215, "E": 206000, "name": "Q235碳素结构钢"},
    "Q345": {"fy": 345, "f": 310, "E": 206000, "name": "Q345低合金高强度钢"},
    "Q390": {"fy": 390, "f": 350, "E": 206000, "name": "Q390低合金高强度钢"},
    "Q420": {"fy": 420, "f": 375, "E": 206000, "name": "Q420低合金高强度钢"},
}

WIND_LEVEL_TABLE = {
    8: 20.7, 9: 24.4, 10: 28.3, 11: 32.6,
    12: 37.0, 13: 41.4, 14: 46.1, 15: 50.9,
    16: 56.0, 17: 61.2,
}

VOLTAGE_PARAMS = {
    "10kV":  {"typical_height": 12, "base_width_ratio": 0.12, "top_width": 0.6},
    "35kV":  {"typical_height": 18, "base_width_ratio": 0.13, "top_width": 0.8},
    "110kV": {"typical_height": 30, "base_width_ratio": 0.18, "top_width": 1.2},
    "220kV": {"typical_height": 45, "base_width_ratio": 0.20, "top_width": 1.8},
    "500kV": {"typical_height": 60, "base_width_ratio": 0.22, "top_width": 2.5},
}

WIRE_TABLE = {
    "LGJ-120/20": {"area": 134.85, "weight": 466.8, "diameter": 15.07},
    "LGJ-185/25": {"area": 211.29, "weight": 706.0, "diameter": 18.90},
    "LGJ-240/30": {"area": 275.96, "weight": 921.0, "diameter": 21.60},
    "LGJ-300/25": {"area": 325.21, "weight": 1058.0, "diameter": 23.40},
    "LGJ-400/35": {"area": 425.24, "weight": 1349.0, "diameter": 26.82},
    "LGJ-500/45": {"area": 531.68, "weight": 1688.0, "diameter": 30.00},
}


class TowerCalculator:
    """铁塔计算器主类"""

    def __init__(self, wind_level, tower_height, tower_type, material,
                 voltage="110kV", wire_type="LGJ-240/30", wire_count=3,
                 terrain_type="B", safety_level=2):
        self.wind_level = wind_level
        self.tower_height = tower_height
        self.tower_type = tower_type
        self.material = material
        self.voltage = voltage
        self.wire_type = wire_type
        self.wire_count = wire_count
        self.terrain_type = terrain_type
        self.safety_level = safety_level
        self.gamma_0 = {1: 1.1, 2: 1.0, 3: 0.9}.get(safety_level, 1.0)
        self.mat = MATERIAL_PROPS.get(material, MATERIAL_PROPS["Q345"])
        self.wire = WIRE_TABLE.get(wire_type, WIRE_TABLE["LGJ-240/30"])
        self.volt_params = VOLTAGE_PARAMS.get(voltage, VOLTAGE_PARAMS["110kV"])
        self.results = {}

    def calc_basic_wind_pressure(self):
        """基本风压 w0 = v^2 / 1600 (kN/m2)"""
        v = WIND_LEVEL_TABLE.get(self.wind_level, 28.3)
        w0 = v ** 2 / 1600.0
        self.results["设计风速"] = v
        self.results["基本风压_w0"] = w0
        return w0

    def calc_wind_pressure_variation(self, z):
        """风压高度变化系数"""
        terrain_coeffs = {
            "A": {"alpha": 0.12, "k": 1.28},
            "B": {"alpha": 0.15, "k": 1.00},
            "C": {"alpha": 0.22, "k": 0.54},
            "D": {"alpha": 0.30, "k": 0.16},
        }
        tc = terrain_coeffs.get(self.terrain_type, terrain_coeffs["B"])
        z = max(z, 1.0)
        mu_z = tc["k"] * (z / 10.0) ** (2 * tc["alpha"])
        return mu_z

    def calc_wind_load_tower(self):
        """计算塔身分段风荷载"""
        w0 = self.calc_basic_wind_pressure()
        n_segments = 10
        dz = self.tower_height / n_segments
        mu_s = {"自立塔": 1.3, "拉线塔": 1.2, "单杆塔": 0.8}.get(self.tower_type, 1.3)
        beta_z = 1.0 + 0.05 * math.sqrt(self.tower_height)

        segments = []
        total_force = 0.0
        total_moment = 0.0
        base_width = self.tower_height * self.volt_params["base_width_ratio"]
        top_width = self.volt_params["top_width"]

        for i in range(n_segments):
            z_bottom = i * dz
            z_top = (i + 1) * dz
            z_mid = (z_bottom + z_top) / 2.0
            mu_z = self.calc_wind_pressure_variation(z_mid)
            width = top_width + (base_width - top_width) * (1 - z_mid / self.tower_height)
            q = beta_z * mu_s * mu_z * w0 * width
            F = q * dz
            M = F * z_mid
            segments.append({"z_bottom": z_bottom, "z_top": z_top, "z_mid": z_mid,
                             "mu_z": mu_z, "width": width, "q_wind": q, "F_wind": F, "M_base": M})
            total_force += F
            total_moment += M

        self.results["塔身分段"] = segments
        self.results["塔身总风荷载"] = total_force
        self.results["塔身风荷载弯矩"] = total_moment
        self.results["风振系数"] = beta_z
        self.results["体型系数"] = mu_s
        self.results["塔底宽度"] = base_width
        self.results["塔顶宽度"] = top_width
        return segments

    def calc_wind_load_wire(self):
        """计算导线、地线风荷载"""
        w0 = self.results.get("基本风压_w0", self.calc_basic_wind_pressure())
        d_wire = self.wire["diameter"] / 1000.0
        mu_s_wire = 1.2
        mu_z_wire = self.calc_wind_pressure_variation(self.tower_height * 0.9)
        span = {"10kV": 80, "35kV": 150, "110kV": 250, "220kV": 350, "500kV": 450}.get(self.voltage, 250)
        alpha_wire = 0.85
        W_wire = alpha_wire * mu_s_wire * mu_z_wire * d_wire * span * w0
        W_wire_total = W_wire * self.wire_count
        d_ground = 0.009
        W_ground = alpha_wire * mu_s_wire * mu_z_wire * d_ground * span * w0
        z_wire = self.tower_height * 0.92
        z_ground = self.tower_height * 0.98
        M_wire = W_wire_total * z_wire
        M_ground = W_ground * z_ground
        self.results["导线风荷载"] = W_wire_total
        self.results["地线风荷载"] = W_ground
        self.results["导线风荷载弯矩"] = M_wire
        self.results["地线风荷载弯矩"] = M_ground
        self.results["计算挡距"] = span
        self.results["导线高度"] = z_wire
        self.results["地线高度"] = z_ground
        return W_wire_total, W_ground

    def calc_internal_forces(self):
        """计算塔身内力（弯矩、剪力、轴力）"""
        span = self.results.get("计算挡距", 250)
        q_wire_weight = self.wire["weight"] * 9.81 / 1e6
        G_wire = q_wire_weight * span * self.wire_count * 2
        G_ground = 3.14159 * (0.009/2)**2 * 7850 * 9.81 * span * 2 / 1000
        if self.tower_type == "自立塔":
            q_tower_weight = 1.5 + 0.08 * self.tower_height
        elif self.tower_type == "拉线塔":
            q_tower_weight = 0.5 + 0.03 * self.tower_height
        else:
            q_tower_weight = 0.8 + 0.04 * self.tower_height

        tower_segments = self.results.get("塔身分段", [])
        n_points = 11
        dz = self.tower_height / (n_points - 1)
        internal_forces = []

        for i in range(n_points):
            z = i * dz
            M = 0.0
            V = 0.0
            for seg in tower_segments:
                if seg["z_mid"] >= z:
                    M += seg["F_wind"] * (seg["z_mid"] - z)
                    V += seg["F_wind"]
            z_wire = self.tower_height * 0.92
            z_ground = self.tower_height * 0.98
            if z_wire >= z:
                M += self.results.get("导线风荷载", 0) * (z_wire - z)
                V += self.results.get("导线风荷载", 0)
            if z_ground >= z:
                M += self.results.get("地线风荷载", 0) * (z_ground - z)
                V += self.results.get("地线风荷载", 0)
            N = q_tower_weight * (self.tower_height - z)
            if z <= z_wire:
                N += G_wire + G_ground
            elif z <= z_ground:
                N += G_ground
            internal_forces.append({"z": z, "M": M, "V": V, "N": N})

        base_force = internal_forces[0]
        self.results["内力分段"] = internal_forces
        self.results["基底弯矩"] = base_force["M"]
        self.results["基底剪力"] = base_force["V"]
        self.results["基底轴力"] = base_force["N"]
        self.results["塔身单位自重"] = q_tower_weight
        self.results["导线重力"] = G_wire
        self.results["地线重力"] = G_ground
        return internal_forces

    def calc_main_member(self):
        """主材截面选择与应力验算"""
        M_max = self.results["基底弯矩"]
        N_max = self.results["基底轴力"]
        base_width = self.results["塔底宽度"]
        N_main = M_max / base_width + N_max / 4.0

        angle_steel = [
            {"name": "L80x8",   "area": 12.30, "i": 2.44},
            {"name": "L90x8",   "area": 13.94, "i": 2.76},
            {"name": "L100x10", "area": 19.26, "i": 3.05},
            {"name": "L110x10", "area": 21.26, "i": 3.38},
            {"name": "L125x12", "area": 28.91, "i": 3.83},
            {"name": "L140x14", "area": 37.57, "i": 4.30},
            {"name": "L160x16", "area": 49.07, "i": 4.89},
            {"name": "L180x18", "area": 61.86, "i": 5.52},
            {"name": "L200x20", "area": 76.54, "i": 6.14},
        ]

        l0 = self.tower_height / 10.0 * 0.8 * 100  # cm (panel length)
        f = self.mat["f"]
        E = self.mat["E"]
        selected = None
        for angle in angle_steel:
            A = angle["area"]  # cm^2
            i = angle["i"]  # cm
            lam = l0 / i  # 长细比
            # 稳定系数 (b类截面简化公式)
            if lam <= 100:
                phi = 1.0 / (1.0 + 0.43 * (lam / 100.0) ** 2)
            else:
                phi = 100.0 / (lam * (1.0 + 0.16 * (lam / 100.0) ** 2))
            # 稳定应力 (N/mm^2)
            sigma = N_main * 1000.0 / (A * 100 * phi)  # kN->N, cm^2->mm^2
            # 强度应力
            sigma_strength = N_main * 1000.0 / (A * 100)
            if sigma <= f and sigma_strength <= f:
                selected = {**angle, "lambda": lam, "phi": phi,
                            "sigma_stability": sigma, "sigma_strength": sigma_strength,
                            "stress_ratio": sigma / f}
                break

        if selected is None:
            selected = {**angle_steel[-1], "lambda": l0 / angle_steel[-1]["i"],
                        "phi": 0.3, "sigma_stability": N_main * 1000.0 / (angle_steel[-1]["area"] * 100 * 0.3),
                        "sigma_strength": N_main * 1000.0 / (angle_steel[-1]["area"] * 100),
                        "stress_ratio": 1.5}

        self.results["主材轴力"] = N_main
        self.results["主材规格"] = selected["name"]
        self.results["主材面积"] = selected["area"]
        self.results["长细比"] = selected["lambda"]
        self.results["稳定系数"] = selected["phi"]
        self.results["稳定应力"] = selected["sigma_stability"]
        self.results["强度应力"] = selected["sigma_strength"]
        self.results["应力比"] = selected["stress_ratio"]
        self.results["主材验算通过"] = selected["stress_ratio"] <= 1.0
        return selected

    def calc_foundation(self):
        """基础底面积计算"""
        N = self.results["基底轴力"]
        M = self.results["基底弯矩"]
        V = self.results["基底剪力"]

        # 假设地基承载力特征值
        fa = 200  # kN/m^2 (中硬土)
        gamma_m = 20.0  # 基础及覆土加权重度 kN/m^3
        d = 2.0  # 基础埋深 m

        # 修正后的地基承载力
        fa_mod = fa + gamma_m * d * (1.0 - 0.5)

        # 迭代计算基础底面积 (同时考虑偏心受压和抗倾覆)
        h_f = 2.0  # 基础高度
        b = math.sqrt(N / fa_mod)
        for _ in range(200):
            A_req = b * b
            W = b ** 3 / 6.0
            p_max = N / A_req + M / W
            G_foundation = 25.0 * b * b * h_f
            M_resist = (N + G_foundation) * b / 2.0 + 0.5 * gamma_m * b * h_f * d * b / 3.0
            k_ot = M_resist / M if M > 0 else 99.0
            if p_max <= fa_mod and k_ot >= 1.6:
                break
            b += 0.05
        p_min = N / A_req - M / W

        # 基础抗倾覆验算
        k_overturn = k_ot

        self.results["地基承载力"] = fa
        self.results["修正承载力"] = fa_mod
        self.results["基础底面积"] = A_req
        self.results["基础边长"] = b
        self.results["基础埋深"] = d
        self.results["基底最大应力"] = p_max
        self.results["基底最小应力"] = p_min
        self.results["抗倾覆安全系数"] = k_overturn
        self.results["基础验算通过"] = p_max <= fa_mod and k_overturn >= 1.6
        return {"b": b, "A": A_req, "p_max": p_max, "p_min": p_min, "k": k_overturn}

    def calc_stability(self):
        """整体稳定验算"""
        M = self.results["基底弯矩"]
        N = self.results["基底轴力"]
        V = self.results["基底剪力"]
        base_width = self.results["塔底宽度"]

        # 抗倾覆稳定系数
        M_resist = (N + 25.0 * base_width * base_width * 2.0) * base_width / 2.0 + 80.0 * base_width
        k_overturn = M_resist / M if M > 0 else 99.0

        # 抗滑移稳定系数
        mu_f = 0.4  # 摩擦系数
        F_resist = mu_f * (N + 25.0 * base_width * base_width * 2.0)
        k_sliding = F_resist / V if V > 0 else 99.0

        self.results["整体抗倾覆系数"] = k_overturn
        self.results["整体抗滑移系数"] = k_sliding
        self.results["整体稳定验算通过"] = k_overturn >= 1.5 and k_sliding >= 1.3
        return {"k_overturn": k_overturn, "k_sliding": k_sliding}

    def run_all(self):
        """执行全部计算"""
        self.calc_wind_load_tower()
        self.calc_wind_load_wire()
        self.calc_internal_forces()
        self.calc_main_member()
        self.calc_foundation()
        self.calc_stability()
        return self.results