完善了双回路EGM模型代码。

core.py

@@ -15,26 +15,44 @@ class Draw:
         global gCAD
         gCAD = self
 
-    def draw(self, i_curt, u_ph, h_gav, h_cav, dgc, color):
+    def draw(self, i_curt, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, color):
         doc = self._doc
         msp = doc.modelspace()
         global gMSP
         gMSP = msp
         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, dxfattribs={"color": color})
-        msp.add_line((0, 0), (0, h_gav))  # 地线
-        msp.add_circle((dgc, h_cav), rc, dxfattribs={"color": color})
-        msp.add_line((dgc, 0), (dgc, h_cav))  # 导线
-        msp.add_line((0, h_gav), (dgc, h_cav))
-        msp.add_line((0, rg), (2000, rg), dxfattribs={"color": color})
+        rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
+        msp.add_circle((rs_x, rs_y), rs, dxfattribs={"color": color})
+        msp.add_line((0, 0), (rs_x, rs_y))  # 地线
+        msp.add_circle((rc_x, rc_y), rc, dxfattribs={"color": color + 2})
+        msp.add_line((rc_x, 0), (rc_x, rc_y))  # 导线
+        msp.add_line((rs_x, rs_y), (rc_x, rc_y))
+        # 角度线
+        circle_intersection = solve_circle_intersection(rs, rc, rs_x, rs_y, rc_x, rc_y)
+        msp.add_line(
+            (rc_x, rc_y), circle_intersection, dxfattribs={"color": color}
+        )  # 地线
+        if rg_type == "g":
+            msp.add_line((0, rg), (2000, rg), dxfattribs={"color": color})
+            circle_line_section = solve_circle_line_intersection(rc, rg, rc_x, rc_y)
+            if not circle_line_section:
+                pass
+            else:
+                msp.add_line(
+                    (rc_x, rc_y), circle_line_section, dxfattribs={"color": color}
+                )  # 导线和圆的交点
+        if rg_type == "c":
+            msp.add_circle((rg_x, rg_y), rg, dxfattribs={"color": color})
+            rg_rc_intersection = solve_circle_intersection(
+                rg, rc, rg_x, rg_y, rc_x, rc_y
+            )
+            msp.add_line(
+                (rc_x, rc_y), rg_rc_intersection, dxfattribs={"color": color}
+            )  # 圆和圆的交点
         # 计算圆交点
-        # 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
@@ -42,15 +60,26 @@ class Draw:
 
 
 # 圆交点
-def solve_circle_intersection(rs, rc, h_gav, h_cav, dgc):
+def solve_circle_intersection(
+    radius1,
+    radius2,
+    center_x1,
+    center_y1,
+    center_x2,
+    center_y2,
+):
     # 用牛顿法求解
-    x = rc  # 初始值
-    y = rc  # 初始值
+    x = radius2  # 初始值
+    y = radius2  # 初始值
+    # TODO 考虑出现2个解的情况
     for bar in range(0, 10):
-        A = [[-2 * x, -2 * (y - h_gav)], [-2 * (x - dgc), -2 * (y - h_cav)]]
+        A = [
+            [-2 * (x - center_x1), -2 * (y - center_y1)],
+            [-2 * (x - center_x2), -2 * (y - center_y2)],
+        ]
         b = [
-            x ** 2 + (y - h_gav) ** 2 - rs ** 2,
-            (x - dgc) ** 2 + (y - h_cav) ** 2 - rc ** 2,
+            (x - center_x1) ** 2 + (y - center_y1) ** 2 - radius1 ** 2,
+            (x - center_x2) ** 2 + (y - center_y2) ** 2 - radius2 ** 2,
         ]
         X_set = np.linalg.solve(A, b)
         x += X_set[0]
@@ -84,7 +113,7 @@ def thunder_density(i):  # l雷电流幅值密度函数
 
 
 def angel_density(angle):  # 入射角密度函数 angle单位是弧度
-    r = 0.75 * (np.cos(angle - math.pi / 2) ** 3)
+    r = 0.75 * abs((np.cos(angle - math.pi / 2) ** 3))
     return r
 
 
@@ -95,35 +124,62 @@ def rs_fun(i):
 
 def rc_fun(i, u_ph):
     r = 1.63 * ((5.015 * (i ** 0.578) - 0.001 * u_ph) ** 1.125)
-    # r=14.7*(i**0.42)
     return r
 
 
-def rg_fun(i_curt, h_cav):
-    if h_cav < 40:
-        rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65)
-    else:
-        rg = 5.5 * (i_curt ** 0.65)
+# typ 如果是g，代表捕雷线公式，c代表暴露弧公式
+def rg_fun(i_curt, h_cav, u_ph, typ="g"):
+    rg = None
+    if typ == "g":
+        if h_cav < 40:
+            rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65)
+        else:
+            rg = 5.5 * (i_curt ** 0.65)
+    elif typ == "c":  # 此时返回的是圆半径
+        rg = rc_fun(i_curt, u_ph)
     return rg
 
 
-def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph):  # 暴露弧的角度
+def intersection_angle(
+    rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_surface, rg_type
+):  # 暴露弧的角度
     rs = rs_fun(i_curt)
     rc = rc_fun(i_curt, u_ph)
-    rg = rg_fun(i_curt, h_cav)
-    circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)  # 两圆的交点
-    circle_line_intersection = solve_circle_line_intersection(
-        rc, rg, dgc, h_cav
-    )  # 暴露圆和补雷线的交点
+    rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
+    circle_intersection = solve_circle_intersection(
+        rs, rc, rs_x, rs_y, rc_x, rc_y
+    )  # 两圆的交点
+    circle_line_or_rg_intersection = None
+    if rg_type == "g":
+        circle_line_or_rg_intersection = solve_circle_line_intersection(
+            rc, rg, rc_x, rc_y
+        )  # 暴露圆和补雷线的交点
+    if rg_type == "c":
+        circle_line_or_rg_intersection = solve_circle_intersection(
+            rg, rc, rg_x, rg_y, rc_x, rc_y
+        )  # 两圆的交点
+    (
+        circle_line_or_rg_intersection_x,
+        circle_line_or_rg_intersection_y,
+    ) = circle_line_or_rg_intersection
+    if (
+        ground_surface(circle_line_or_rg_intersection_x)
+        > circle_line_or_rg_intersection_y
+    ):  # 交点在地面线以下，就可以不积分
+        # 找到暴露弧和地面线的交点
+        circle_line_or_rg_intersection = circle_ground_surface_intersection(
+            rc, rc_x, rc_y, ground_surface
+        )
     np_circle_intersection = np.array(circle_intersection)
-    theta2_line = np_circle_intersection - np.array([dgc, h_cav])
+    theta2_line = np_circle_intersection - np.array([rc_x, rc_y])
     theta2 = math.atan(theta2_line[1] / theta2_line[0])
-    np_circle_line_intersection = np.array(circle_line_intersection)
-    theta1_line = np_circle_line_intersection - np.array([dgc, h_cav])
+    np_circle_line_or_rg_intersection = np.array(circle_line_or_rg_intersection)
+    theta1_line = np_circle_line_or_rg_intersection - np.array([rc_x, rc_y])
     theta1 = math.atan(theta1_line[1] / theta1_line[0])
     return np.array([theta1, theta2])
 
 
+# 点到直线的距离
 def distance_point_line(point_x, point_y, line_x, line_y, k) -> float:
     d = abs(k * point_x - point_y - k * line_x + line_y) / ((k ** 2 + 1) ** 0.5)
     return d
@@ -137,17 +193,17 @@ def func_calculus_pw(theta, max_w):
     return r_pw
 
 
-def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav):  # 对θ进行积分
+def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y):  # 对θ进行积分
     max_w = 0
     # 求暴露弧上一点的切线
-    line_x = math.cos(theta) * rc + dgc
-    line_y = math.sin(theta) * rc + h_cav
+    line_x = math.cos(theta) * rc + rc_x
+    line_y = math.sin(theta) * rc + rc_y
     k = math.tan(theta + math.pi / 2)  # 入射角
     # 求保护弧到直线的距离，判断是否相交
-    d_to_rs = distance_point_line(0, h_gav, line_x, line_y, k)
+    d_to_rs = distance_point_line(rs_x, rs_y, line_x, line_y, k)
     if d_to_rs < rs:  # 相交
         # 要用过直线上一点到暴露弧的切线
-        new_k = tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k)
+        new_k = tangent_line_k(line_x, line_y, rs_x, rs_y, rs, init_k=k)
         if new_k >= 0:
             max_w = math.atan(new_k)  # 用于保护弧相切的角度
         elif new_k < 0:
@@ -183,20 +239,23 @@ def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav):  # 对θ进行积分
     return r
 
 
-def bd_area(i_curt, u_ph, dgc, h_gav, h_cav):  # 暴露弧的投影面积
-    theta1, theta2 = intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph)  # θ角度
+def bd_area(
+    i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, ground_surface, rg_type
+):  # 暴露弧的投影面积
+    theta1, theta2 = intersection_angle(
+        rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_surface, rg_type
+    )  # θ角度
     theta_fineness = 0.01
     rc = rc_fun(i_curt, u_ph)
     rs = rs_fun(i_curt)
-    rg = rg_fun(i_curt, h_cav)
-    r_bd = 0
+    rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
     theta_sample, d_theta = np.linspace(
         theta1, theta2, int((theta2 - theta1) / theta_fineness), retstep=True
     )
     if len(theta_sample) < 2:
         return 0
     vec_calculus_bd = np.vectorize(calculus_bd)
-    calculus_bd_np = vec_calculus_bd(theta_sample, rc, rs, rg, dgc, h_cav, h_gav)
+    calculus_bd_np = vec_calculus_bd(theta_sample, rc, rs, rg, rc_x, rc_y, rs_x, rs_y)
     r_bd = np.sum(calculus_bd_np[:-1] + calculus_bd_np[1:]) / 2 * d_theta
     # for calculus_theta in theta_sample[:-1]:
     #     r_bd += (
@@ -214,7 +273,7 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
     # 直线方程为 y-y0=k(x-x0)，x0和y0为经过直线的任意一点
     # 牛顿法求解k
     # f(k)=(k*x1-y1-k*x0+y0)**2-R**2*(k**2+1) x1,y1是圆心
-
+    # 已考虑两个解的判别
     k_candidate = [-100, 100]
     if abs(center_y - line_y) < 1 and abs(line_x - center_x - radius) < 1:
         # k不存在
@@ -248,10 +307,10 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
     # 把k转化成相应的角度，从x开始，逆时针为正
     k_angle = []
     for kk in k_candidate:
-        if kk is None:
-            abc = 123
-            # tangent_line_k(line_x, line_y, center_x, center_y, radius)
-            pass
+        # if kk is None:
+        #     abc = 123
+        #     # tangent_line_k(line_x, line_y, center_x, center_y, radius)
+        #     pass
         if kk >= 0:
             k_angle.append(math.atan(kk))
         if kk < 0:
@@ -262,3 +321,19 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
 
 def func_ng(td):  # 地闪密度
     return 0.023 * (td ** 1.3)
+
+
+# 圆和地面线的交点，只去正x轴上的。
+def circle_ground_surface_intersection(radius, center_x, center_y, ground_surface):
+    # 最笨的办法，一个个去试
+    x_series = np.linspace(0, radius, int(radius / 0.001)) + center_x
+    part_to_be_squared = radius ** 2 - (x_series - center_x) ** 2  # 有可能出现-0.00001的数值，只是一个数值稳定问题。
+    part_to_be_squared[
+        (part_to_be_squared < 0) & (abs(part_to_be_squared) < 1e-3)
+    ] = 0  # 强制为0
+    y_series = center_y - part_to_be_squared ** 0.5
+    ground_surface_y = ground_surface(x_series)
+    equal_location = np.abs(ground_surface_y - y_series) < 0.5
+    r_x = x_series[equal_location][0]
+    r_y = ground_surface(r_x)
+    return r_x, r_y

main.py

@@ -1,135 +1,212 @@
-import numpy as np
-
 from core import *
 import timeit
 
 
 def egm():
-    avr_n_sf = 0  # 考虑电压的影响计算的跳闸率
-    voltage_n = 3  # 工作电压分成多少份来计算
-    ng = func_ng(20)
+    h_g_avr_sag = 11.67 * 2 / 3
+    h_c_avr_sag = 14.43 * 2 / 3
     h_whole = 140  # 杆塔全高
+    voltage_n = 3  # 工作电压分成多少份来计算
+    td = 20  # 雷暴日
     insulator_c_len = 6.8  # 串子绝缘长度
     string_c_len = 9.2
     string_g_len = 0.5
-    dgc = -0.0  # 导地线水平距离
-    vertical_dgc = 2.7  # 导地线挂点垂直距离
-    h_g_avr_sag = 11.67 * 2 / 3
-    h_c_avr_sag = 14.43 * 2 / 3
-    h_gav = h_whole - string_g_len - h_g_avr_sag  # 地线对地平均高
-    h_cav = h_whole - string_c_len - vertical_dgc - h_c_avr_sag  # 导线对地平均高
-    shield_angle = math.atan(dgc / (vertical_dgc + string_c_len)) * 180 / math.pi
-    print(f"保护角{shield_angle:.3f}°")
-    for u_bar in range(voltage_n):
-        u_ph = (
-            math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732
-        )  # 运行相电压
-        # 迭代法计算最大电流
-        i_max = 0
-        i_min = min_i(insulator_c_len, u_ph / 1.732)
-        _min_i = i_min  # 尝试的最小电流
-        _max_i = 200  # 尝试的最大电流
-        cad = Draw()
-        cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2)
-        for i_bar in np.linspace(_min_i, _max_i, int((_max_i - _min_i) / 0.1)):  # 雷电流
-            # print(f"尝试计算电流为{i_bar:.2f}")
-            rs = rs_fun(i_bar)
-            rc = rc_fun(i_bar, u_ph)
-            rg = rg_fun(i_bar, h_cav)
-            #######
-            # cccCount += 1
-            # if cccCount % 30 == 0:
-            #     import core
-            #
-            #     core.gMSP.add_circle((0, h_gav), rs)
-            #     core.gMSP.add_circle(
-            #         (dgc, h_cav), rc_fun(i_bar, -u_ph), dxfattribs={"color": 4}
-            #     )
-            #     core.gMSP.add_circle((dgc, h_cav), rc)
-            #######
-            circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
-            if not circle_intersection:  # if circle_intersection is []
-                # print("保护弧和暴露弧无交点，检查设置参数。程序退出。")
-                continue
-            circle_rc_line_intersection = solve_circle_line_intersection(
-                rc, rg, dgc, h_cav
-            )
-            if not circle_rc_line_intersection:
-                continue
-            min_distance_intersection = (
-                np.sum(
-                    (
-                        np.array(circle_intersection)
-                        - np.array(circle_rc_line_intersection)
+    gc_x = [17.9, 16, 15, 16]
+
+    # 以后考虑地形角度,地面线
+    def ground_surface(x):
+        return 0
+
+    gc_y = [
+        h_whole - string_g_len - h_g_avr_sag,  # 地线对地平均高
+        h_whole - string_c_len - h_c_avr_sag - 2.7,  # 导线对地平均高
+        h_whole - string_c_len - h_c_avr_sag - 20,  # 导线对地平均高
+        h_whole - string_c_len - h_c_avr_sag - 35.7,  # 导线对地平均高
+    ]
+    if len(gc_y) > 2:  # 双回路
+        phase_n = 3  # 边相导线数量
+    else:
+        phase_n = 1
+    #########################################################
+    rg_type = None
+    # 以上是需要设置的参数
+    avr_n_sf = 0  # 考虑电压的影响计算的跳闸率
+    rg_x = None
+    rg_y = None
+    cad = Draw()
+    for phase_conductor in range(phase_n):
+        rs_x = gc_x[phase_conductor]
+        rs_y = gc_y[phase_conductor]
+        rc_x = gc_x[phase_conductor + 1]
+        rc_y = gc_y[phase_conductor + 1]
+        if phase_n == 1:
+            rg_type = "g"
+        if phase_n > 1:  # 多回路
+            if phase_conductor < 2:
+                rg_type = "c"
+                rg_x = gc_x[phase_conductor + 2]
+                rg_y = gc_y[phase_conductor + 2]
+            else:
+                rg_type = "g"
+        # TODO 保护角公式可能有问题，后面改
+        shield_angle = (
+            math.atan(rc_x / ((rc_y - rs_y) + string_c_len)) * 180 / math.pi
+        )  # 保护角
+        print(f"保护角{shield_angle:.3f}°")
+        print(f"最低相防护标识{rg_type}")
+        ng = func_ng(td)
+        for u_bar in range(voltage_n):
+            u_ph = (
+                math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732
+            )  # 运行相电压
+            print(f"计算第{phase_conductor + 1}相，电压为{u_ph:.2f}kV")
+            # 迭代法计算最大电流
+            i_max = 0
+            i_min = min_i(insulator_c_len, u_ph / 1.732)
+            _min_i = i_min  # 尝试的最小电流
+            _max_i = 200  # 尝试的最大电流
+            # cad.draw(i_min, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 2)
+            for i_bar in np.linspace(
+                _min_i, _max_i, int((_max_i - _min_i) / 0.1)
+            ):  # 雷电流
+                # print(f"尝试计算电流为{i_bar:.2f}")
+                rs = rs_fun(i_bar)
+                rc = rc_fun(i_bar, u_ph)
+                rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type)
+                #######
+                # cccCount += 1
+                # if cccCount % 30 == 0:
+                #     import core
+                #
+                #     core.gMSP.add_circle((0, h_gav), rs)
+                #     core.gMSP.add_circle(
+                #         (dgc, h_cav), rc_fun(i_bar, -u_ph), dxfattribs={"color": 4}
+                #     )
+                #     core.gMSP.add_circle((dgc, h_cav), rc)
+                #######
+                rg_rc_circle_intersection = solve_circle_intersection(
+                    rs, rc, rs_x, rs_y, rc_x, rc_y
+                )
+                i_max = i_bar
+                if not rg_rc_circle_intersection:  # if circle_intersection is []
+                    print("保护弧和暴露弧无交点，检查设置参数。")
+                    continue
+                circle_rc_line_or_rg_intersection = None
+                if rg_type == "g":
+                    circle_rc_line_or_rg_intersection = solve_circle_line_intersection(
+                        rc, rg, rc_x, rc_y
                     )
-                    ** 2
-                )
-                ** 0.5
-            )  # 计算两圆交点和地面直线交点的最小距离
-            i_max = i_bar
-            if min_distance_intersection < 0.1:
-                break
-            if circle_intersection[1] < circle_rc_line_intersection[1]:
-                circle_rs_line_intersection = solve_circle_line_intersection(
-                    rs, rg, 0, h_gav
-                )
-                # 判断与保护弧的交点是否在暴露弧外面
-                distance = (
+                elif rg_type == "c":
+                    circle_rc_line_or_rg_intersection = solve_circle_intersection(
+                        rg, rc, rg_x, rg_y, rc_x, rc_y
+                    )
+                if not circle_rc_line_or_rg_intersection:
+                    # 暴露弧和捕捉弧无交点
+                    if rg_type == "g":
+                        if rg > rc_y:
+                            i_min = i_bar
+                            print(f"捕捉弧在暴露弧之上，设置最小电流为{i_min:.2f}")
+                        else:
+                            print("暴露弧和捕捉弧无交点，检查设置参数。")
+                        continue
+                    else:
+                        print("暴露弧和捕捉弧无交点，检查设置参数。")
+                        continue
+                min_distance_intersection = (
                     np.sum(
-                        (np.array(circle_rs_line_intersection) - np.array([dgc, h_cav]))
+                        (
+                            np.array(rg_rc_circle_intersection)
+                            - np.array(circle_rc_line_or_rg_intersection)
+                        )
                         ** 2
                     )
                     ** 0.5
-                )
-                if distance > rc:
-                    print("暴露弧已经完全被屏蔽")
+                )  # 计算两圆交点和地面直线交点的最小距离
+                if min_distance_intersection < 0.1:
                     break
-        cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2)
-        cad.draw(i_max, u_ph, h_gav, h_cav, dgc, 6)
-        cad.save()
-        # 判断是否导线已经被完全保护
-        if abs(i_max - _max_i) < 1e-5:
-            print("无法找到最大电流，可能是杆塔较高。")
-            print(f"最大电流设置为自然界最大电流{i_max}kA")
-        print(f"最大电流为{i_max:.2f}")
-        print(f"最小电流为{i_min:.2f}")
-        curt_fineness = 0.1  # 电流积分细度
-        if i_min > i_max or abs(i_min - i_max) < curt_fineness:
-            print("最大电流小于最小电流，没有暴露弧，程序结束。")
-            return
-        # 开始积分
-        curt_segment_n = int((i_max - i_min) / curt_fineness)  # 分成多少份
-        calculus = 0
-        i_curt_samples, d_curt = np.linspace(
-            i_min, i_max, curt_segment_n + 1, retstep=True
-        )
-        bd_area_vec = np.vectorize(bd_area)
-        cal_bd_np = bd_area_vec(
-            i_curt_samples, u_ph, dgc, h_gav, h_cav
-        ) * thunder_density(i_curt_samples)
-        calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt
-        # for i_curt in i_curt_samples[:-1]:
-        #     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
-        #             + cal_bd_second * cal_thunder_density_second
-        #         )
-        #         / 2
-        #         * d_curt
-        #     )
-        #     if abs(calculus-0.05812740052770032)<1e-5:
-        #         abc=123
-        #         pass
-        n_sf = (
-            2 * ng / 10 * calculus
-        )  # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
-        avr_n_sf += n_sf / voltage_n
-        print(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}")
-    print(f"跳闸率是{avr_n_sf:.6}")
+                # 判断是否以完全被保护
+                if rg_rc_circle_intersection[1] < circle_rc_line_or_rg_intersection[1]:
+                    circle_rs_line_or_rg_intersection = None
+                    if rg_type == "g":
+                        circle_rs_line_or_rg_intersection = (
+                            solve_circle_line_intersection(rs, rg, rs_x, rs_y)
+                        )
+                    if rg_type == "c":
+                        circle_rs_line_or_rg_intersection = solve_circle_intersection(
+                            rs, rg, rs_x, rs_y, rg_x, rg_y
+                        )
+                    # 判断与保护弧的交点是否在暴露弧外面
+                    distance = (
+                        np.sum(
+                            (
+                                np.array(circle_rs_line_or_rg_intersection)
+                                - np.array([rc_x, rc_y])
+                            )
+                            ** 2
+                        )
+                        ** 0.5
+                    )
+                    if distance > rc:
+                        print("暴露弧已经完全被屏蔽")
+                        break
+            if phase_conductor == 1:
+                cad.draw(i_min, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 2)
+                cad.draw(i_max, u_ph, rs_x, rs_y, rc_x, rc_y, rg_x, rg_y, rg_type, 6)
+                cad.save()
+            # 判断是否导线已经被完全保护
+            if abs(i_max - _max_i) < 1e-5:
+                print("无法找到最大电流，可能是杆塔较高。")
+                print(f"最大电流设置为自然界最大电流{i_max}kA")
+            print(f"最大电流为{i_max:.2f}")
+            print(f"最小电流为{i_min:.2f}")
+            curt_fineness = 0.1  # 电流积分细度
+            if i_min > i_max or abs(i_min - i_max) < curt_fineness:
+                print("最大电流小于最小电流，没有暴露弧。")
+                continue
+            # 开始积分
+            curt_segment_n = int((i_max - i_min) / curt_fineness)  # 分成多少份
+            i_curt_samples, d_curt = np.linspace(
+                i_min, i_max, curt_segment_n + 1, retstep=True
+            )
+            bd_area_vec = np.vectorize(bd_area)
+            cal_bd_np = (
+                bd_area_vec(
+                    i_curt_samples,
+                    u_ph,
+                    rc_x,
+                    rc_y,
+                    rs_x,
+                    rs_y,
+                    rg_x,
+                    rg_y,
+                    ground_surface,
+                    rg_type,
+                )
+                * thunder_density(i_curt_samples)
+            )
+            calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt
+            # for i_curt in i_curt_samples[:-1]:
+            #     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
+            #             + cal_bd_second * cal_thunder_density_second
+            #         )
+            #         / 2
+            #         * d_curt
+            #     )
+            #     if abs(calculus-0.05812740052770032)<1e-5:
+            #         abc=123
+            #         pass
+            n_sf = (
+                2 * ng / 10 * calculus
+            )  # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
+            avr_n_sf += n_sf / voltage_n
+            print(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}")
+    print(f"跳闸率是{avr_n_sf:.6f}")
 
 
 def speed():