考虑了切线的k为负的情况

main.py

@@ -2,18 +2,48 @@ import math
 import ezdxf
 import numpy as np
 
+gCAD = None
+gMSP = None
+
+
+class Draw:
+    def __init__(self):
+        self._doc = ezdxf.new(dxfversion="R2010")
+        self._doc.layers.add("EGM", color=2)
+        global gCAD
+        gCAD = self
+
+    def draw(self, i_curt, u_ph, h_gav, h_cav, dgc, 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})
+        # 计算圆交点
+        # 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)
-    # y = Symbol('y', real=True)
-    # equ = [
-    #     x ** 2 + (y - hgav) ** 2 - rs ** 2,
-    #     (x - dgc) ** 2 + (y - hcav) ** 2 - rc ** 2,
-    # ]
     # 用牛顿法求解
-    x = 300
-    y = 300
+    x = rc  # 初始值
+    y = rc  # 初始值
     for bar in range(0, 10):
         A = [[-2 * x, -2 * (y - hgav)], [-2 * (x - dgc), -2 * (y - hcav)]]
         b = [
@@ -26,45 +56,15 @@ def solve_circle_intersection(rs, rc, hgav, hcav, dgc):
         if np.all(np.abs(X_set) < 1e-5):
             return [x, y]
     return []
-    # list_set = list(X_set)
-    # solve_set = nonlinsolve(equ, [x, y])
-    # print(ask(Q.real(solve_set)))
-    # list_set = list(solve_set)
-    # pprint(list_set)
-    # if not np.all(np.isreal(list_set)):
-    #     return []
-    # for value in list_set:
-    #     if value[0] > 0 and value[1] > 1:
-    #         return value
-    # return []
 
 
 # 圆与地面线交点
-def solve_circle_line_intersection(rc, rg, hcav, dgc):
-    r = (rc ** 2 - (rg - hcav) ** 2) ** 0.5 + dgc
+def solve_circle_line_intersection(rc, rg, h_cav, dgc):
+    # TODO: 需要考虑地面捕雷线与暴露弧完全没交点的情况
+    r = (rc ** 2 - (rg - h_cav) ** 2) ** 0.5 + 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,12 +74,12 @@ 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
 
 
 def angel_density(angle):  # 入射角密度函数 angle单位是弧度
-    r = 0.75 * (math.cos(angle) ** 3)
+    r = 0.75 * (math.cos(angle - math.pi / 2) ** 3)
     return r
 
 
@@ -101,94 +101,228 @@ def rg_fun(i, h_cav):
     return rg
 
 
-def intersection_angel(dgc, h_gav, h_cav, i_curt, u_ph):  # 暴露弧的角度
+def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph):  # 暴露弧的角度
     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, h_cav, dgc)
+    circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)  # 两圆的交点
+    circle_line_intersection = solve_circle_line_intersection(
+        rc, rg, h_cav, dgc
+    )  # 暴露圆和补雷线的交点
     np_circle_intersection = np.array(circle_intersection)
+    if not circle_intersection:
+        abc = 123
     theta2_line = np_circle_intersection - np.array([dgc, h_cav])
     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])
     theta1 = math.atan(theta1_line[1] / theta1_line[0])
+    # 考虑雷电入射角度，所以theta1可以小于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):  # 暴露弧的投影面积
-    rc = rc_fun(i_curt, u_ph)
-    # 暂时不考虑雷电入射角的影响
-    r = (math.cos(theta1) - math.cos(theta2)) * rc
+def distance_point_line(point_x, point_y, line_x, line_y, k):
+    d = abs(k * point_x - point_y - k * line_x + line_y) / ((k ** 2 + 1) ** 0.5)
+    return d
+
+
+def fun_calculus_pw(theta, max_w):
+    w_fineness = 0.01
+    r_pw = 0
+    if int(max_w / w_fineness) < 0:
+        abc = 123
+        pass
+    w_samples, d_w = np.linspace(0, max_w, int(max_w / w_fineness) + 1, retstep=True)
+    for cal_w in w_samples:
+        r_pw += (
+            (
+                abs(angel_density(cal_w)) * math.sin(theta - cal_w + math.pi)
+                + abs(angel_density(cal_w + d_w))
+                * math.sin(theta - cal_w + math.pi - d_w)
+            )
+            / 2
+        ) * d_w
+    return r_pw
+
+
+def cal_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(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)
+        if new_k >= 0:
+            max_w = math.atan(new_k)  # 用于保护弧相切的角度
+        elif new_k < 0:
+            max_w = math.atan(new_k) + math.pi
+        if max_w < 0:
+            abc = 123
+            tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k)
+            # intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph)
+            # gMSP.add_circle((0, h_gav), rs)
+            # gMSP.add_circle((dgc, h_cav), rc)
+            # gMSP.add_line((dgc, h_cav), (line_x, line_y))
+            # gMSP.add_line(
+            #     (-500, k * (-500 - line_x) + line_y), (500, k * (500 - line_x) + line_y)
+            # )
+            # gMSP.add_line((0, rg), (1000, rg))
+            # gCAD.save()
+            tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k)
+
+        pass
+    r = rc / math.cos(theta) * fun_calculus_pw(theta, max_w)
     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)
+    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
+    theta_sample, d_theta = np.linspace(
+        theta1, theta2, int((theta2 - theta1) / theta_fineness) + 1, retstep=True
+    )
+    for cal_theta in theta_sample[:-1]:
+        r_bd += (
+            (
+                cal_bd(cal_theta, rc, rs, rg, dgc, h_cav, h_gav)
+                + cal_bd(cal_theta + d_theta, rc, rs, rg, dgc, h_cav, h_gav)
+            )
+            / 2
+            * d_theta
+        )
+    return r_bd
+
     # r1=rc*(-math.cos(thyta2)+math.cos(thyta1))
     # 入射角密度函数积分
     # arrival_angle_fineness=0.0001
     # for calculus_arv_angle in np.linspace()
 
 
+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是圆心
+    # 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 = (
+                2
+                * (k * center_x - center_y - k * line_x + line_y)
+                * (center_x - line_x)
+                - 2 * (radius ** 2) * k
+            )
+            d_k = -fk / d_fk
+            k += d_k
+            if abs(d_k) < 1e-5:
+                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 >= 0:
+            k_angle.append(math.atan(kk))
+        if kk < 0:
+            k_angle.append(math.pi + math.atan(kk))
+    # 返回相对x轴最大的角度k
+    return np.array(k_candidate)[np.max(k_angle) == k_angle].tolist()[-1]
+
+
 def egm():
-    u_ph = 750 / 1.732  # 运行相电压
-    h_cav = 60  # 导线对地平均高
-    h_gav = h_cav + 9.5 + 7.2
-    dgc = 2
-    # 迭代法计算最大电流
-    i_max = 0
-    _min_i = 30  # 尝试的最小电流
-    _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}")
-        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(rc, rg, h_cav, dgc)
-        min_distance_intersection = (
-            np.sum(
-                (np.array(circle_intersection) - np.array(circle_line_intersection))
-                ** 2
+    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 = h_cav + 9.5 + 2.7
+        dgc = 2.9  # 导地线水平距离
+        # 迭代法计算最大电流
+        i_max = 0
+        _min_i = 20  # 尝试的最小电流
+        _max_i = 300  # 尝试的最大电流
+        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(
+                rc, rg, h_cav, dgc
             )
-            ** 0.5
-        )  # 计算两圆交点和地面直线交点的最小距离
-        if min_distance_intersection < 0.01:
+            min_distance_intersection = (
+                np.sum(
+                    (np.array(circle_intersection) - np.array(circle_line_intersection))
+                    ** 2
+                )
+                ** 0.5
+            )  # 计算两圆交点和地面直线交点的最小距离
             i_max = i_bar
-            draw(rs, rc, rg, h_gav, h_cav, dgc)
-            break
-    print(f"最大电流为{i_max:.2f}")
-    i_min = min_i(6.78, 750 / 1.732)
-    print(f"最小电流为{i_min:.2f}")
-    # 开始积分
-    curt_fineness = 0.1  # 电流积分细度
-    curt_segment_n = int((i_max - i_min) / curt_fineness)
-    d_curt = (i_max - i_min) / curt_segment_n
-    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]
+            if min_distance_intersection < 0.1:
+                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}")
+        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
         )
-        cal_thunder_density_first = thunder_density(curt)
-        cal_thunder_density_second = thunder_density(curt + d_curt)
-        calculus += (
-            (
-                cal_bd_first * cal_thunder_density_first
-                + cal_bd_second * cal_thunder_density_second
+        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
+                    + cal_bd_second * cal_thunder_density_second
+                )
+                / 2
+                * d_curt
             )
-            / 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)