修复几个边界条件bug

core.py

@@ -158,24 +158,31 @@ def intersection_angle(
         circle_line_or_rg_intersection = solve_circle_intersection(
             rg, rc, rg_x, rg_y, rc_x, rc_y
         )  # 两圆的交点
-    (
+    if circle_line_or_rg_intersection:
-        circle_line_or_rg_intersection_x,
+        (
-        circle_line_or_rg_intersection_y,
+            circle_line_or_rg_intersection_x,
-    ) = circle_line_or_rg_intersection
+            circle_line_or_rg_intersection_y,
-    if (
+        ) = circle_line_or_rg_intersection
-        ground_surface(circle_line_or_rg_intersection_x)
+        if (
-        > circle_line_or_rg_intersection_y
+            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
+            circle_line_or_rg_intersection = circle_ground_surface_intersection(
-        )
+                rc, rc_x, rc_y, ground_surface
+            )
+    theta1 = None
     np_circle_intersection = np.array(circle_intersection)
     theta2_line = np_circle_intersection - np.array([rc_x, rc_y])
     theta2 = math.atan(theta2_line[1] / theta2_line[0])
     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])
+    if not circle_line_or_rg_intersection:
-    theta1 = math.atan(theta1_line[1] / theta1_line[0])
+        if rc_y - rc > rg:  # rg在rc下面
+            # 捕捉线太低了，对高塔无保护，θ_1从-90°开始计算。
+            theta1 = -math.pi / 2
+    else:
+        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])
 
 
@@ -187,7 +194,10 @@ def distance_point_line(point_x, point_y, line_x, line_y, k) -> float:
 
 def func_calculus_pw(theta, max_w):
     w_fineness = 0.01
-    w_samples, d_w = np.linspace(0, max_w, int(max_w / w_fineness), retstep=True)
+    segments = int(max_w / w_fineness)
+    if segments < 2:  # 最大最小太小，没有可以积分的
+        return 0
+    w_samples, d_w = np.linspace(0, max_w, segments, retstep=True)
     cal_w_np = abs(angel_density(w_samples)) * np.sin(theta - (w_samples - math.pi / 2))
     r_pw = np.sum((cal_w_np[:-1] + cal_w_np[1:])) / 2 * d_w
     return r_pw
@@ -282,7 +292,8 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
         for ind, k_cdi in enumerate(list(k_candidate)):
             k = k_candidate[ind]
             k_candidate[ind] = None
-            for bar in range(0, 30):
+            max_iteration = 30
+            for bar in range(0, max_iteration):
                 fk = (k * center_x - center_y - k * line_x + line_y) ** 2 - (
                     radius ** 2
                 ) * (k ** 2 + 1)
@@ -304,13 +315,17 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
                     if abs(dd - radius) < 1:
                         k_candidate[ind] = k
                         break
+                # 解决数值稳定性
+                if bar == max_iteration - 1:
+                    if abs(math.atan(k)) * 180 / math.pi > 89:
+                        k_candidate[ind] = k
     # 把k转化成相应的角度，从x开始，逆时针为正
     k_angle = []
     for kk in k_candidate:
-        # if kk is None:
+        if kk is None:
-        #     abc = 123
+            abc = 123
-        #     # tangent_line_k(line_x, line_y, center_x, center_y, radius)
+            tangent_line_k(line_x, line_y, center_x, center_y, radius)
-        #     pass
+            pass
         if kk >= 0:
             k_angle.append(math.atan(kk))
         if kk < 0:
@@ -327,7 +342,9 @@ def func_ng(td):  # 地闪密度
 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 = (
+        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

main.py

@@ -1,3 +1,5 @@
+import numpy as np
+
 from core import *
 import timeit
 
@@ -5,13 +7,13 @@ import timeit
 def egm():
     h_g_avr_sag = 11.67 * 2 / 3
     h_c_avr_sag = 14.43 * 2 / 3
-    h_whole = 140  # 杆塔全高
+    h_whole = 260  # 杆塔全高
     voltage_n = 3  # 工作电压分成多少份来计算
     td = 20  # 雷暴日
     insulator_c_len = 6.8  # 串子绝缘长度
     string_c_len = 9.2
     string_g_len = 0.5
-    gc_x = [17.9, 16, 15, 16]
+    gc_x = [17.9, 17, 15, 17]
 
     # 以后考虑地形角度,地面线
     def ground_surface(x):
@@ -34,6 +36,7 @@ def egm():
     rg_x = None
     rg_y = None
     cad = Draw()
+    n_sf_phases = np.zeros((phase_n, voltage_n))  # 计算每一相的跳闸率
     for phase_conductor in range(phase_n):
         rs_x = gc_x[phase_conductor]
         rs_y = gc_y[phase_conductor]
@@ -149,7 +152,7 @@ def egm():
                     if distance > rc:
                         print("暴露弧已经完全被屏蔽")
                         break
-            if phase_conductor == 1:
+            if phase_conductor == 2:
                 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()
@@ -205,8 +208,10 @@ def egm():
                 2 * ng / 10 * calculus
             )  # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
             avr_n_sf += n_sf / voltage_n
+            n_sf_phases[phase_conductor][u_bar] = n_sf
             print(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}")
     print(f"跳闸率是{avr_n_sf:.6f}")
+    print(f"不同相跳闸率是{np.mean(n_sf_phases,axis=1)}")
 
 
 def speed():