1.处理了90°，k不存在的情况。

2.重命名了一些函数。
3.增加了一些不绕击的判断。

main.py

@@ -4,6 +4,7 @@ import numpy as np
 
 gCAD = None
 gMSP = None
+gCount = 1
 
 
 class Draw:
@@ -59,9 +60,9 @@ def solve_circle_intersection(rs, rc, hgav, hcav, dgc):
 
 
 # 圆与地面线交点
-def solve_circle_line_intersection(rc, rg, h_cav, dgc):
+def solve_circle_line_intersection(radius, rg, center_y, center_x):
     # TODO: 需要考虑地面捕雷线与暴露弧完全没交点的情况
-    r = (rc ** 2 - (rg - h_cav) ** 2) ** 0.5 + dgc
+    r = (radius ** 2 - (rg - center_y) ** 2) ** 0.5 + center_x
     return [r, rg]
 
 
@@ -93,11 +94,11 @@ def rc_fun(i, u_ph):
     return r
 
 
-def rg_fun(i, h_cav):
+def rg_fun(i_curt, h_cav):
     if h_cav < 40:
         rg = (3.6 + 1.7 ** math.log(43 - h_cav)) ** 0.65
     else:
-        rg = 5.5 * (i ** 0.65)
+        rg = 5.5 * (i_curt ** 0.65)
     return rg
 
 
@@ -148,17 +149,17 @@ def fun_calculus_pw(theta, max_w):
     return r_pw
 
 
-def cal_bd(theta, rc, rs, rg, dgc, h_cav, h_gav):
+def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav):  # 对θ进行积分
     # 求暴露弧上一点的切线
     line_x = math.cos(theta) * rc + dgc
     line_y = math.sin(theta) * rc + h_cav
-    max_w = theta + math.pi / 2  # 入射角
+    k = math.tan(theta + math.pi / 2)  # 入射角
-    k = math.tan(max_w)
     # 求保护弧到直线的距离，判断是否相交
     d_to_rs = distance_point_line(0, h_gav, 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)
+
         if not new_k:
             a = 12
             tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k)
@@ -169,18 +170,21 @@ def cal_bd(theta, rc, rs, rg, dgc, h_cav, h_gav):
         if max_w < 0:
             abc = 123
             tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k)
+        global gCount
+        gCount += 1
+        if gCount % 1000 == 0:
             # intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph)
-            # gMSP.add_circle((0, h_gav), rs)
+            gMSP.add_circle((0, h_gav), rs)
-            # gMSP.add_circle((dgc, h_cav), rc)
+            gMSP.add_circle((dgc, h_cav), rc)
-            # gMSP.add_line((dgc, h_cav), (line_x, line_y))
+            gMSP.add_line((dgc, h_cav), (line_x, line_y))
-            # gMSP.add_line(
+            gMSP.add_line(
-            #     (-500, k * (-500 - line_x) + line_y), (500, k * (500 - line_x) + line_y)
+                (-500, new_k * (-500 - line_x) + line_y),
-            # )
+                (500, new_k * (500 - line_x) + line_y),
-            # gMSP.add_line((0, rg), (1000, rg))
+            )
-            # gCAD.save()
+            gMSP.add_line((0, rg), (1000, rg))
-            tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k)
+            gCAD.save()
-
+    else:
-        pass
+        max_w = theta + math.pi / 2  # 入射角
     r = rc / math.cos(theta) * fun_calculus_pw(theta, max_w)
     return r
 
@@ -195,11 +199,11 @@ def bd_area(i_curt, u_ph, dgc, h_gav, h_cav):  # 暴露弧的投影面积
     theta_sample, d_theta = np.linspace(
         theta1, theta2, int((theta2 - theta1) / theta_fineness) + 1, retstep=True
     )
-    for cal_theta in theta_sample[:-1]:
+    for calculus_theta in theta_sample[:-1]:
         r_bd += (
             (
-                cal_bd(cal_theta, rc, rs, rg, dgc, h_cav, h_gav)
+                calculus_bd(calculus_theta, rc, rs, rg, dgc, h_cav, h_gav)
-                + cal_bd(cal_theta + d_theta, rc, rs, rg, dgc, h_cav, h_gav)
+                + calculus_bd(calculus_theta + d_theta, rc, rs, rg, dgc, h_cav, h_gav)
             )
             / 2
             * d_theta
@@ -217,31 +221,44 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
     # 牛顿法求解k
     # f(k)=(k*x1-y1-k*x0+y0)**2-R**2*(k**2+1) x1,y1是圆心
     # TODO:需要检验k值不存在的情况
-    k_candidate = [-100, 100]
-    for ind, k_cdi in enumerate(list(k_candidate)):
-        k = k_candidate[ind]
-        k_candidate[ind] = None
-        for bar in range(0, 30):
-            fk = (k * center_x - center_y - k * line_x + line_y) ** 2 - (
-                radius ** 2
-            ) * (k ** 2 + 1)
 
-            d_fk = (
+    k_candidate = [-100, 100]
-                2
+    if abs(center_y - line_y) < 1e-4 and abs(line_x - center_x - radius) < 1e-4:
-                * (k * center_x - center_y - k * line_x + line_y)
+        # k不存在
-                * (center_x - line_x)
+        k_candidate = [99999999, 99999999]
-                - 2 * (radius ** 2) * k
+    else:
-            )
+        for ind, k_cdi in enumerate(list(k_candidate)):
-            d_k = -fk / d_fk
+            k = k_candidate[ind]
-            k += d_k
+            k_candidate[ind] = None
-            if abs(d_k) < 1e-5:
+            for bar in range(0, 30):
-                dd = distance_point_line(center_x, center_y, line_x, line_y, k)
+                fk = (k * center_x - center_y - k * line_x + line_y) ** 2 - (
-                if abs(dd - radius) < 1e-5:
+                    radius ** 2
-                    k_candidate[ind] = k
+                ) * (k ** 2 + 1)
+
+                d_fk = (
+                    2
+                    * (k * center_x - center_y - k * line_x + line_y)
+                    * (center_x - line_x)
+                    - 2 * (radius ** 2) * k
+                )
+                if abs(d_fk) < 1e-5 and abs(line_x - center_x - radius) < 1e-5:
+                    # k不存在，角度为90°，k取一个很大的正数
+                    k_candidate[ind] = 99999999999999
                     break
+                d_k = -fk / d_fk
+                k += d_k
+                if abs(d_k) < 1e-4:
+                    dd = distance_point_line(center_x, center_y, line_x, line_y, k)
+                    if abs(dd - radius) < 1e-5:
+                        k_candidate[ind] = k
+                        break
     # 把k转化成相应的角度，从x开始，逆时针为正
     k_angle = []
     for kk in k_candidate:
+        if kk == 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:
@@ -253,9 +270,9 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
 def egm():
     for u_bar in range(1):
         u_ph = math.sqrt(2) * 750 * math.cos(2 * math.pi / 6 * 0) / 1.732  # 运行相电压
-        h_cav = 90  # 导线对地平均高
+        h_gav = 30
-        h_gav = h_cav + 9.5 + 2.7
+        h_cav = h_gav - 9.5 - 2.7 - 5  # 导线对地平均高
-        dgc = 2.9  # 导地线水平距离
+        dgc = -2.9  # 导地线水平距离
         # 迭代法计算最大电流
         i_max = 0
         _min_i = 20  # 尝试的最小电流
@@ -263,23 +280,20 @@ def egm():
         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)
-            # if not np.isreal(rs):
-            #     continue
             rc = rc_fun(i_bar, u_ph)
-            # if not np.isreal(rc):
-            #     continue
             rg = rg_fun(i_bar, h_cav)
-            # if not np.isreal(rg):
-            #     continue
             circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
             if not circle_intersection:  # if circle_intersection is []
                 continue
-            circle_line_intersection = solve_circle_line_intersection(
+            circle_rc_line_intersection = solve_circle_line_intersection(
                 rc, rg, h_cav, dgc
             )
             min_distance_intersection = (
                 np.sum(
-                    (np.array(circle_intersection) - np.array(circle_line_intersection))
+                    (
+                        np.array(circle_intersection)
+                        - np.array(circle_rc_line_intersection)
+                    )
                     ** 2
                 )
                 ** 0.5
@@ -287,14 +301,29 @@ def egm():
             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, h_gav, 0
+                )
+                # 判断与保护弧的交点是否在暴露弧外面
+                distance = (
+                    np.sum(
+                        (np.array(circle_rs_line_intersection) - np.array([dgc, h_cav]))
+                        ** 2
+                    )
+                    ** 0.5
+                )
+                if distance > rc:
+                    print("暴露弧已经完全被屏蔽")
+                    break
         i_min = min_i(6.78, u_ph / 1.732)
         cad = Draw()
         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("无法找到最大电流，可能是杆塔较高。")
-            i_max = 300
             print(f"最大电流设置为自然界最大电流{i_max}kA")
         print(f"最大电流为{i_max:.2f}")
         print(f"最小电流为{i_min:.2f}")
@@ -328,6 +357,6 @@ def egm():
 
 
 if __name__ == "__main__":
-    thunder_density(2)
+    tangent_line_k(1, 0, 0, 0, 1)
     egm()
     print("Finished.")