考虑了高塔情况。

main.py

@@ -3,6 +3,35 @@ import ezdxf
 import numpy as np
 
 
+class Draw:
+    def __init__(self):
+        self._doc = ezdxf.new(dxfversion="R2010")
+        self._doc.layers.add("EGM", color=2)
+
+    def draw(self, i_curt, u_ph, h_gav, h_cav, dgc):
+        doc = self._doc
+        msp = doc.modelspace()
+        rs = rs_fun(i_curt)
+        rc = rc_fun(i_curt, u_ph)
+        rg = rg_fun(i_curt, h_cav)
+        msp.add_circle((0, h_gav), rs)
+        msp.add_line((0, 0), (0, h_gav))  # 地线
+        msp.add_circle((dgc, h_cav), rc)
+        msp.add_line((dgc, 0), (dgc, h_cav))  # 导线
+        msp.add_line((0, h_gav), (dgc, h_cav))
+        msp.add_line((0, rg), (200, rg))
+        # 计算圆交点
+        circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
+        msp.add_line((0, h_gav), circle_intersection)  # 地线
+        msp.add_line((dgc, h_cav), circle_intersection)  # 导线
+        circle_line_section = solve_circle_line_intersection(rc, rg, h_cav, dgc)
+        msp.add_line((0, 0), circle_line_section)  # 导线和圆的交点
+
+    def save(self):
+        doc = self._doc
+        doc.saveas("egm.dxf")
+
+
 # 圆交点
 def solve_circle_intersection(rs, rc, hgav, hcav, dgc):
     # x = Symbol('x', real=True)
@@ -45,26 +74,6 @@ def solve_circle_line_intersection(rc, rg, hcav, dgc):
     return [r, rg]
 
 
-def draw(rs, rc, rg, h_gav, h_cav, dgc):
-    doc = ezdxf.new(dxfversion="R2010")
-    doc.layers.add("EGM", color=2)
-    msp = doc.modelspace()
-    msp.add_circle((0, h_gav), rs)
-    msp.add_line((0, 0), (0, h_gav))  # 地线
-    msp.add_circle((dgc, h_cav), rc)
-    msp.add_line((dgc, 0), (dgc, h_cav))  # 导线
-    msp.add_line((0, h_gav), (dgc, h_cav))
-    msp.add_line((0, rg), (200, rg))
-    # 计算圆交点
-    circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
-    msp.add_line((0, h_gav), circle_intersection)  # 地线
-    msp.add_line((dgc, h_cav), circle_intersection)  # 导线
-    circle_line_section = solve_circle_line_intersection(rc, rg, h_cav, dgc)
-    msp.add_line((0, 0), circle_line_section)  # 导线和圆的交点
-    doc.saveas("egm.dxf")
-    solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
-
-
 def min_i(string_len, u_ph):
     u_50 = 530 * string_len + 35
     z_0 = 300  # 雷电波阻抗
@@ -74,7 +83,7 @@ def min_i(string_len, u_ph):
 
 
 def thunder_density(i):  # l雷电流幅值密度函数
-    r = -10 ** (-i / 44) * math.log(10) * (-1 / 44)
+    r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44)
     return r
 
 
@@ -113,10 +122,14 @@ def intersection_angel(dgc, h_gav, h_cav, i_curt, u_ph):  # 暴露弧的角度
     np_circle_line_intersection = np.array(circle_line_intersection)
     theta1_line = np_circle_line_intersection - np.array([dgc, h_cav])
     theta1 = math.atan(theta1_line[1] / theta1_line[0])
+    if theta1 < 0:
+        # print(f"θ_1角度为负数{theta1:.4f},人为设置为0")
+        theta1 = 0
     return np.array([theta1, theta2])
 
 
-def bd_area(i_curt, u_ph, theta1, theta2):  # 暴露弧的投影面积
+def bd_area(i_curt, u_ph, dgc, h_gav, h_cav):  # 暴露弧的投影面积
+    theta1, theta2 = intersection_angel(dgc, h_gav, h_cav, i_curt, u_ph)
     rc = rc_fun(i_curt, u_ph)
     # 暂时不考虑雷电入射角的影响
     r = (math.cos(theta1) - math.cos(theta2)) * rc
@@ -129,12 +142,12 @@ def bd_area(i_curt, u_ph, theta1, theta2):  # 暴露弧的投影面积
 
 def egm():
     u_ph = 750 / 1.732  # 运行相电压
-    h_cav = 60  # 导线对地平均高
-    h_gav = h_cav + 9.5 + 7.2
-    dgc = 2
+    h_cav = 160  # 导线对地平均高
+    h_gav = h_cav + 9.5 + 2.2
+    dgc = 2  # 导地线水平距离
     # 迭代法计算最大电流
     i_max = 0
-    _min_i = 30  # 尝试的最小电流
+    _min_i = 20  # 尝试的最小电流
     _max_i = 80  # 尝试的最大电流
     for i_bar in np.linspace(_min_i, _max_i, int((_max_i - _min_i) / 0.01)):  # 雷电流
         print(f"尝试计算电流为{i_bar:.2f}")
@@ -158,27 +171,33 @@ def egm():
             )
             ** 0.5
         )  # 计算两圆交点和地面直线交点的最小距离
-        if min_distance_intersection < 0.01:
-            i_max = i_bar
-            draw(rs, rc, rg, h_gav, h_cav, dgc)
+        i_max = i_bar
+        if min_distance_intersection < 0.1:
             break
-    print(f"最大电流为{i_max:.2f}")
     i_min = min_i(6.78, 750 / 1.732)
+    cad = Draw()
+    cad.draw(i_min, u_ph, h_gav, h_cav, dgc)
+    cad.draw(i_max, u_ph, h_gav, h_cav, dgc)
+    cad.save()
+    if abs(i_max - _max_i) < 1e-5:
+        print("无法找到最大电流，可能是杆塔较高。")
+        i_max = 300
+        print(f"最大电流设置为自然界最大电流{i_max}kA")
+    print(f"最大电流为{i_max:.2f}")
     print(f"最小电流为{i_min:.2f}")
+    if i_min > i_max:
+        print("最大电流小于最小电流，没有暴露弧，程序结束。")
+        return
     # 开始积分
-    curt_fineness = 0.1  # 电流积分细度
-    curt_segment_n = int((i_max - i_min) / curt_fineness)
-    d_curt = (i_max - i_min) / curt_segment_n
+    curt_fineness = 0.001  # 电流积分细度
+    curt_segment_n = int((i_max - i_min) / curt_fineness)  # 分成多少份
     calculus = 0
-    for curt in np.linspace(i_min, i_max, curt_segment_n):
-        cal_thyta_first = intersection_angel(dgc, h_gav, h_cav, curt, u_ph)
-        cal_bd_first = bd_area(curt, u_ph, cal_thyta_first[0], cal_thyta_first[1])
-        cal_thyta_second = intersection_angel(dgc, h_gav, h_cav, curt + d_curt, u_ph)
-        cal_bd_second = bd_area(
-            curt + d_curt, u_ph, cal_thyta_second[0], cal_thyta_second[1]
-        )
-        cal_thunder_density_first = thunder_density(curt)
-        cal_thunder_density_second = thunder_density(curt + d_curt)
+    i_curt_samples, d_curt = np.linspace(i_min, i_max, curt_segment_n + 1, retstep=True)
+    for i_curt in i_curt_samples:
+        cal_bd_first = bd_area(i_curt, u_ph, dgc, h_gav, h_cav)
+        cal_bd_second = bd_area(i_curt + d_curt, u_ph, dgc, h_gav, h_cav)
+        cal_thunder_density_first = thunder_density(i_curt)
+        cal_thunder_density_second = thunder_density(i_curt + d_curt)
         calculus += (
             (
                 cal_bd_first * cal_thunder_density_first
@@ -187,8 +206,8 @@ def egm():
             / 2
             * d_curt
         )
-    n_sf=2*2.7/10*calculus
-    print(f'跳闸率是{n_sf}')
+    n_sf = 2 * 2.7 / 10 * calculus  # 调整率
+    print(f"跳闸率是{n_sf:.6}")
 
     # draw(rs, rc, rg, h_gav, h_cav, dgc)