1.每一次计算都重新画画。

版本号 Version: 1.2.0.2

.gitignore

@@ -0,0 +1,13 @@
+*.dxf
+build
+__pycache__
+CSharp
+.idea
+dist
+*.spec
+*.dwg
+历史
+.venv
+*.toml
+launch.json
+settings.json

Makefile

@@ -0,0 +1,14 @@
+target: dist build
+	create-version-file metadata.yml --outfile build/file_version_info.txt
+	pyinstaller -F main.py --version-file build/file_version_info.txt -n Lightening
+	
+dist:
+	mkdir dist
+
+build:
+	mkdir build
+	
+.PHONY:clean
+
+clean:
+	rm -fr build

animation.py

@@ -0,0 +1,94 @@
+import matplotlib.pyplot as plt
+from functools import wraps
+import numpy as np
+
+
+class Animation:
+    def __init__(self) -> None:
+        fig, ax = plt.subplots()
+        self._fig = fig
+        self._ax = ax
+        self._ticks = 0
+        self._disable = False
+        self.init_fig()
+        pass
+
+    @staticmethod
+    def switch_decorator(func):
+        @wraps(func)
+        def not_run(cls, *args, **kwargs):
+            # print("not run")
+            pass
+
+        @wraps(func)
+        def wrapTheFunction(cls, *args, **kwargs):
+            if not cls._disable:
+                # print("desc")
+                return func(cls, *args, **kwargs)
+            return not_run(cls, *args, **kwargs)
+
+        return wrapTheFunction
+
+    def disable(self, _disable):
+        self._disable = _disable
+
+    @switch_decorator
+    def init_fig(self):
+        ax = self._ax
+        ax.set_aspect(1)
+        ax.set_xlim([-500, 500])
+        ax.set_ylim([-500, 500])
+
+    @switch_decorator
+    def show(self):
+        self._fig.show()
+
+    @switch_decorator
+    def add_rg_line(self, line_func):
+        ax = self._ax
+        x = np.linspace(0, 300)
+        y = line_func(x)
+        ax.plot(x, y)
+
+    @switch_decorator
+    def add_rs(self, rs, rs_x, rs_y):
+        ax = self._ax
+        ax.add_artist(plt.Circle((rs_x, rs_y), rs, fill=False))
+
+    @switch_decorator
+    def add_rc(self, rc, rc_x, rc_y):
+        ax = self._ax
+        ax.add_artist(plt.Circle((rc_x, rc_y), rc, fill=False))
+
+    # 增加暴露弧范围
+    @switch_decorator
+    def add_expose_area(
+        self,
+        rc_x,
+        rc_y,
+        intersection_x1,
+        intersection_y1,
+        intersection_x2,
+        intersection_y2,
+    ):
+        ax = self._ax
+        ax.plot([rc_x, intersection_x1], [rc_y, intersection_y1], color="red")
+        ax.plot([rc_x, intersection_x2], [rc_y, intersection_y2], color="red")
+        pass
+
+    @switch_decorator
+    def clear(self):
+        ax = self._ax
+        ax.cla()
+
+    @switch_decorator
+    def pause(self):
+        ax = self._ax
+        self._ticks += 1
+        ticks = self._ticks
+        ax.set_title(f"{ticks}")
+        plt.pause(0.02)
+        self.clear()
+        self.init_fig()
+
+    pass

article-第一版论文的参数.toml

@@ -0,0 +1,18 @@
+title = "绕击跳闸率计算文件"
+[parameter]
+rated_voltage=750 #额定电压等级
+h_c_sag = 14.43  # 导线弧垂
+h_g_sag = 11.67  # 地线弧垂
+h_whole = 100  # 杆塔全高
+insulator_c_len = 6.8  # 导线串子绝缘长度
+string_c_len = 9.2  # 导线串长
+string_g_len = 0.5  # 地线串长
+h_arm = [100,80]  # 导、地线挂点垂直距离
+gc_x = [17.9,17]  # 导、地线水平坐标,计算的是中相
+ground_angels = [0]  # 地面倾角，向下为正，单位°
+altitude = 1000  # 海拔，单位米
+max_i = 200  # 最大尝试电流，单位kA
+td=20#雷暴日
+
+[optional]
+voltage_n=3 #计算时电压分成多少份

article.toml

@@ -0,0 +1,23 @@
+title = "绕击跳闸率计算文件"
+[parameter]
+rated_voltage = 750  # 额定电压等级
+h_c_sag = 14.43  # 导线弧垂
+h_g_sag = 11.67  # 地线弧垂
+insulator_c_len = 7.0  # 导线串子绝缘长度
+string_c_len = 9.2  # 导线串长
+string_g_len = 0.5  # 地线串长
+h_arm = [100, 80]  # 导、地线挂点垂直距离,计算的是中相
+gc_x = [17.9, 17]  # 导、地线水平坐标,计算的是中相
+ground_angels = [0]  # 地面倾角，向下为正，单位°
+altitude = 1500  # 海拔，单位米
+td = 20  # 雷暴日
+[advance]
+# ng=29.6 #地闪密度 ！！注意！！ 如果地闪密度大于0，则不会通过雷暴日计算地闪密度。填-1则忽略该项数据。
+# Ip_a=19.896 #概率密度曲线系数a ！！注意！！ 如果该值大于0，则不会通过雷暴日计算雷电流概率。填-1则忽略该项数据。
+# Ip_b=3.012 #概率密度曲线系数b  ！！注意！！ 如果该值大于0，则不会通过雷暴日计算雷电流概率。填-1则忽略该项数据。
+ng = -1  # 地闪密度 ！！注意！！ 如果地闪密度大于0，则不会通过雷暴日计算地闪密度。填-1则忽略该项数据。
+Ip_a = -1  # 概率密度曲线系数a ！！注意！！ 如果该值大于0，则不会通过雷暴日计算雷电流概率。填-1则忽略该项数据。
+Ip_b = -1  # 概率密度曲线系数b  ！！注意！！ 如果该值大于0，则不会通过雷暴日计算雷电流概率。填-1则忽略该项数据。
+[optional]
+voltage_n = 3  # 计算时电压分成多少份
+max_i = 200  # 最大尝试电流，单位kA

core.py

@@ -0,0 +1,492 @@
+import math
+import ezdxf
+import numpy as np
+from typing import List
+
+gCAD = None
+gMSP = None
+gCount = 1
+
+
+class Parameter:
+    h_g_sag: float  # 地线弧垂
+    h_c_sag: float  # 导线弧垂
+    voltage_n: int  # 工作电压分成多少份来计算
+    td: int  # 雷暴日
+    insulator_c_len: float  # 串子绝缘长度
+    string_c_len: float
+    string_g_len: float
+    gc_x: List[float]  # 导、地线水平坐标
+    ground_angels: List[float]  # 地面倾角，向下为正
+    h_arm: float  # 导、地线垂直坐标
+    altitude: int  # 海拔，单位米
+    max_i: float  # 最大尝试电流，单位kA
+    rated_voltage: float  # 额定电压
+    ng: float  # 地闪密度  次/(每平方公里·每年)
+    Ip_a: float  # 概率密度曲线系数a
+    Ip_b: float  # 概率密度曲线系数b
+
+
+para = Parameter()
+
+
+def rg_line_function_factory(_rg, ground_angel):  # 返回一个地面捕雷线的直线方程
+    y_d = _rg / math.cos(ground_angel)  # y轴上的截距
+    # 利用公式y-y0=k(x-x0) 得到直线公式
+    y0 = y_d
+    x0 = 0
+    k = math.tan(math.pi - ground_angel)
+
+    def f(x):
+        return y0 + k * (x - x0)
+
+    return f
+
+
+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,
+        rs_x,
+        rs_y,
+        rc_x,
+        rc_y,
+        rg_x,
+        rg_y,
+        rg_type,
+        ground_angel,
+        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, 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":
+            ground_angel_func = rg_line_function_factory(rg, ground_angel)
+            msp.add_line(
+                (0, ground_angel_func(0)),
+                (2000, ground_angel_func(2000)),
+                dxfattribs={"color": color},
+            )
+            circle_line_section = solve_circle_line_intersection(
+                rc, rc_x, rc_y, ground_angel_func
+            )
+            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
+            )
+            if rg_rc_intersection:
+                msp.add_line(
+                    (rc_x, rc_y), rg_rc_intersection, dxfattribs={"color": color}
+                )  # 圆和圆的交点
+        # 计算圆交点
+
+        # msp.add_line((dgc, h_cav), circle_intersection)  # 导线
+
+    def save(self):
+        doc = self._doc
+        doc.saveas("egm.dxf")
+
+    def save_as(self, file_name):
+        doc = self._doc
+        doc.saveas(file_name)
+
+
+# 圆交点
+def solve_circle_intersection(
+    radius1,
+    radius2,
+    center_x1,
+    center_y1,
+    center_x2,
+    center_y2,
+):
+    # 用牛顿法求解
+    x = radius2 + center_x2  # 初始值
+    y = radius2 + center_y2  # 初始值
+    # TODO 考虑出现2个解的情况
+    for _ in range(0, 10):
+        A = [
+            [-2 * (x - center_x1), -2 * (y - center_y1)],
+            [-2 * (x - center_x2), -2 * (y - center_y2)],
+        ]
+        b = [
+            (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]
+        y += X_set[1]
+        if np.all(np.abs(X_set) < 1e-5):
+            return [x, y]
+    return []
+
+
+# 圆与捕雷线交点
+def solve_circle_line_intersection(
+    radius, center_x, center_y, ground_surface_func
+):  # 返回交点的x和y坐标
+    x0 = 0
+    y0 = ground_surface_func(x0)
+    x1 = 1
+    y1 = ground_surface_func(x1)
+    k = (y1 - y0) / (x1 - x0)
+    distance = distance_point_line(center_x, center_y, x0, y0, k)  # 捕雷线到暴露圆中点的距离
+    if distance > radius:
+        return []
+    else:
+        # r = (radius ** 2 - (rg - center_y) ** 2) ** 0.5 + center_x
+        a = center_x
+        b = center_y
+        c = y0
+        d = x0
+        bb = -2 * a + 2 * c * k - 2 * d * (k**2) - 2 * b * k
+        aa = 1 + k**2
+        rr = radius
+        cc = (
+            a**2
+            + c**2
+            - 2 * c * k * d
+            + (k**2) * (d**2)
+            - 2 * b * (c - k * d)
+            + b**2
+            - rr**2
+        )
+        _x = (-bb + (bb**2 - 4 * aa * cc) ** 0.5) / 2 / aa
+        _y = ground_surface_func(_x)
+        # 验算结果
+        equ = (center_x - _x) ** 2 + (center_y - _y) ** 2 - radius**2
+        assert abs(equ) < 1e-5
+        return [_x, _y]
+
+
+def min_i(string_len, u_ph):
+    # 海拔修正
+    altitude = para.altitude
+    if altitude > 1000:
+        k_a = math.exp((altitude - 1000) / 8150)  # 气隙海拔修正
+    else:
+        k_a = 1
+    u_50 = 1 / k_a * (530 * string_len + 35)  # 50045 上附录的公式，实际应该用负极性电压的公式
+    # u_50 = 1 / k_a * (533 * string_len + 132)  # 串放电路径 1000m海拔
+    # u_50 = 1 / k_a * (477 * string_len + 99)  # 串放电路径 2000m海拔
+    # u_50 = 615 * string_len  # 导线对塔身放电 1000m海拔
+    # u_50= 263.32647401+533.90081562*string_len
+    z_0 = 300  # 雷电波阻抗
+    z_c = 251  # 导线波阻抗
+    # 新版大手册公式 3-277
+    r = (u_50 + 2 * z_0 / (2 * z_0 + z_c) * u_ph) * (2 * z_0 + z_c) / (z_0 * z_c)
+    # r = 2 * (u_50 - u_ph) / z_c
+    return r
+
+
+def thunder_density(i, td, ip_a, ip_b):  # 雷电流幅值密度函数
+    # td = para.td
+    r = None
+    # ip_a = para.Ip_a
+    # ip_b = para.Ip_b
+    if ip_a > 0 and ip_b > 0:
+        r = -(
+            -ip_b / ip_a / ((1 + (i / ip_a) ** ip_b) ** 2) * ((i / ip_a) ** (ip_b - 1))
+        )
+        return r
+    else:
+        if td == 20:
+            r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44)  # 雷暴日20d
+            return r
+        if td == 40:
+            r = -(10 ** (-i / 88)) * math.log(10) * (-1 / 88)  # 雷暴日40d
+            return r
+        raise Exception("检查雷电参数!")
+
+
+def angel_density(angle):  # 入射角密度函数 angle单位是弧度
+    r = 0.75 * abs((np.cos(angle - math.pi / 2) ** 3))  # 新版大手册公式3-275
+    return r
+
+
+def rs_fun(i):
+    r = 10 * (i**0.65)  # 新版大手册公式3-271
+    return r
+
+
+def rc_fun(i, u_ph):
+    r = 1.63 * ((5.015 * (i**0.578) - 0.001 * u_ph * 1) ** 1.125)  # 新版大手册公式3-272
+    return r
+
+
+# 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)  # 新版大手册公式3-273
+        else:
+            rg = 5.5 * (i_curt**0.65)  # 新版大手册公式3-273
+    elif typ == "c":  # 此时返回的是圆半径
+        rg = rc_fun(i_curt, u_ph)
+    return rg
+
+
+def intersection_angle(
+    rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_angel, rg_type
+):  # 暴露弧的角度
+    rs = rs_fun(i_curt)
+    rc = rc_fun(i_curt, u_ph)
+    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
+    rg_line_func = rg_line_function_factory(rg, ground_angel)
+    if rg_type == "g":
+        circle_line_or_rg_intersection = solve_circle_line_intersection(
+            rc, rc_x, rc_y, rg_line_func
+        )  # 暴露圆和补雷线的交点
+    if rg_type == "c":
+        circle_line_or_rg_intersection = solve_circle_intersection(
+            rg, rc, rg_x, rg_y, rc_x, rc_y
+        )  # 两圆的交点
+    # TODO 应该是不存在落到地面线以下的情况，先把以下注释掉
+
+    # if circle_line_or_rg_intersection:
+    #     (
+    #         circle_line_or_rg_intersection_x,
+    #         circle_line_or_rg_intersection_y,
+    #     ) = circle_line_or_rg_intersection
+    #     if (
+    #         ground_surface(rg, 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
+    #         )
+    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)
+    if not circle_line_or_rg_intersection:
+        if rc_y - rc > rg_line_func(rc_x):  # 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])
+
+
+# 点到直线的距离
+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
+
+
+def func_calculus_pw(theta, max_w):
+    w_fineness = 0.01
+    segments = int(max_w / w_fineness)
+    if segments < 2:  # 最大最小太小，没有可以积分的
+        return 0
+    w_samples, d_w = np.linspace(0, max_w, segments, retstep=True)
+    # 童中宇 750KV信洛Ⅰ线雷电防护性能研究  公式 3-10
+    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
+
+
+def calculus_bd(theta, rc, rs, rg, rc_x, rc_y, rs_x, rs_y):  # 对θ进行积分
+    max_w = 0
+    # 求暴露弧上一点的切线
+    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(rs_x, rs_y, line_x, line_y, k)
+    if d_to_rs < rs:  # 相交
+        # 要用过这一点到保护弧的切线
+        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:
+            max_w = math.atan(new_k) + math.pi
+        # TODO to be removed
+        # global gCount
+        # gCount = gCount+1
+        # if gCount % 100 == 0:
+        # 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, new_k * (-500 - line_x) + line_y),
+        #     (500, new_k * (500 - line_x) + line_y),
+        # )
+        # gCAD.save()
+        # pass
+    else:
+        max_w = theta + math.pi / 2  # 入射角
+        # TODO to be removed
+        if gCount % 200 == 0:
+            # # 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),
+            # )
+            # gCAD.save()
+            pass
+    # 童中宇 750KV信洛Ⅰ线雷电防护性能研究  公式 3-10
+    r = rc / math.cos(theta) * func_calculus_pw(theta, max_w)
+    return r
+
+
+def bd_area(
+    i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, ground_angel, rg_type
+):  # 暴露弧的投影面积
+    theta1, theta2 = intersection_angle(
+        rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_angel, rg_type
+    )  # θ角度
+    theta_fineness = 0.01
+    rc = rc_fun(i_curt, u_ph)
+    rs = rs_fun(i_curt)
+    rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
+    theta_segments = int((theta2 - theta1) / theta_fineness)
+    if theta_segments < 2:
+        return 0
+    theta_sample, d_theta = np.linspace(theta1, theta2, theta_segments, 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, 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 += (
+    #         (
+    #             calculus_bd(calculus_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
+    #     )
+    return r_bd
+
+
+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不存在
+        k_candidate = [99999999, 99999999]
+    else:
+        for ind, k_cdi in enumerate(list(k_candidate)):
+            k = k_candidate[ind]
+            k_candidate[ind] = None
+            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)
+
+                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-3:
+                    dd = distance_point_line(center_x, center_y, line_x, line_y, k)
+                    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:
+        #     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:
+            k_angle.append(math.pi + math.atan(kk))
+    # 返回相对x轴最大的角度k
+    return np.array(k_candidate)[np.max(k_angle) == k_angle].tolist()[-1]
+
+
+def func_ng(td):  # 地闪密度，通过雷暴日计算
+    if para.ng > 0:
+        r = para.ng
+    else:
+        r = 0.023 * (td**1.3)
+    return r
+
+
+# 圆和地面线的交点，只去正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
+
+
+# u_ph是相电压
+# insulator_c_len绝缘子闪络距离
+def arc_possibility(rated_voltage, insulator_c_len):  # 建弧率
+    # 50064 中附录给的公式
+    # TODO 需要区分交直流
+    # _e = rated_voltage / (3**0.5) / insulator_c_len #交流
+    _e = abs(rated_voltage) / (1) / insulator_c_len  # 直流
+    r = (4.5 * (_e**0.75) - 14) * 1e-2
+    return r

main.py

@@ -1,333 +1,370 @@
 import math
-import ezdxf
+import os.path
-import numpy as np
+import sys
-
+import time
-gCAD = None
+import tomli
-gMSP = None
+from loguru import logger
+from core import *
+import timeit
+from animation import Animation
 
 
-class Draw:
+# 打印参数
-    def __init__(self):
+def parameter_display(para_dis: Parameter):
-        self._doc = ezdxf.new(dxfversion="R2010")
+    logger.info(f"额定电压 kV {para_dis.rated_voltage}")
-        self._doc.layers.add("EGM", color=2)
+    logger.info(f"导线弧垂 m {para_dis.h_c_sag}")
-        global gCAD
+    logger.info(f"地线弧垂 m {para_dis.h_g_sag}")
-        gCAD = self
+    logger.info(f"全塔高 m {para_dis.h_arm[0]}")
-
+    logger.info(f"串绝缘距离 m {para_dis.insulator_c_len}")
-    def draw(self, i_curt, u_ph, h_gav, h_cav, dgc, color):
+    logger.info(f"导线串长 m {para_dis.string_c_len}")
-        doc = self._doc
+    logger.info(f"地线串长 m {para_dis.string_g_len}")
-        msp = doc.modelspace()
+    logger.info(f"挂点垂直坐标 m {para_dis.h_arm}")
-        global gMSP
+    logger.info(f"挂点水平坐标 m {para_dis.gc_x}")
-        gMSP = msp
+    logger.info(f"地面倾角 ° {[an * 180 / math.pi for an in para_dis.ground_angels]}")
-        rs = rs_fun(i_curt)
+    logger.info(f"海拔高度 m {para_dis.altitude}")
-        rc = rc_fun(i_curt, u_ph)
+    if para_dis.ng > 0:
-        rg = rg_fun(i_curt, h_cav)
+        logger.info("不采用雷暴日计算地闪密度和雷电流密度")
-        msp.add_circle((0, h_gav), rs, dxfattribs={"color": color})
+        logger.info(f"地闪密度 次/(每平方公里·每年) {para_dis.ng}")
-        msp.add_line((0, 0), (0, h_gav))  # 地线
+        logger.info(f"概率密度曲线系数a {para_dis.Ip_a}")
-        msp.add_circle((dgc, h_cav), rc, dxfattribs={"color": color})
+        logger.info(f"概率密度曲线系数b {para_dis.Ip_b}")
-        msp.add_line((dgc, 0), (dgc, h_cav))  # 导线
+        pass
-        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 = rc  # 初始值
-    y = rc  # 初始值
-    for bar in range(0, 10):
-        A = [[-2 * x, -2 * (y - hgav)], [-2 * (x - dgc), -2 * (y - hcav)]]
-        b = [
-            x ** 2 + (y - hgav) ** 2 - rs ** 2,
-            (x - dgc) ** 2 + (y - hcav) ** 2 - rc ** 2,
-        ]
-        X_set = np.linalg.solve(A, b)
-        x += X_set[0]
-        y += X_set[1]
-        if np.all(np.abs(X_set) < 1e-5):
-            return [x, y]
-    return []
-
-
-# 圆与地面线交点
-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 min_i(string_len, u_ph):
-    u_50 = 530 * string_len + 35
-    z_0 = 300  # 雷电波阻抗
-    z_c = 251  # 导线波阻抗
-    r = (u_50 + 2 * z_0 / (2 * z_0 + z_c) * u_ph) * (2 * z_0 + z_c) / (z_0 * z_c)
-    return r
-
-
-def thunder_density(i):  # l雷电流幅值密度函数
-    r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44)
-    return r
-
-
-def angel_density(angle):  # 入射角密度函数 angle单位是弧度
-    r = 0.75 * (math.cos(angle - math.pi / 2) ** 3)
-    return r
-
-
-def rs_fun(i):
-    r = 10 * (i ** 0.65)
-    return r
-
-
-def rc_fun(i, u_ph):
-    r = 1.63 * ((5.015 * (i ** 0.578) - 0.001 * u_ph) ** 1.125)
-    return r
-
-
-def rg_fun(i, 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)
+        logger.info(f"雷暴日 d {para_dis.td}")
-    return rg
 
 
-def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph):  # 暴露弧的角度
+def read_parameter(toml_file_path):
-    rs = rs_fun(i_curt)
+    with open(toml_file_path, "rb") as toml_fs:
-    rc = rc_fun(i_curt, u_ph)
+        toml_dict = tomli.load(toml_fs)
-    rg = rg_fun(i_curt, h_cav)
+        toml_parameter = toml_dict["parameter"]
-    circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)  # 两圆的交点
+        para.h_g_sag = toml_parameter["h_g_sag"]  # 地线弧垂
-    circle_line_intersection = solve_circle_line_intersection(
+        para.h_c_sag = toml_parameter["h_c_sag"]  # 导线弧垂
-        rc, rg, h_cav, dgc
+        # para.h_whole = toml_parameter["h_whole"]  # 杆塔全高
-    )  # 暴露圆和补雷线的交点
+        para.td = toml_parameter["td"]  # 雷暴日
-    np_circle_intersection = np.array(circle_intersection)
+        para.insulator_c_len = toml_parameter["insulator_c_len"]  # 串子绝缘长度
-    if not circle_intersection:
+        para.string_c_len = toml_parameter["string_c_len"]
-        abc = 123
+        para.string_g_len = toml_parameter["string_g_len"]
-    theta2_line = np_circle_intersection - np.array([dgc, h_cav])
+        para.gc_x = toml_parameter["gc_x"]  # 导、地线水平坐标
-    theta2 = math.atan(theta2_line[1] / theta2_line[0])
+        para.ground_angels = [
-    np_circle_line_intersection = np.array(circle_line_intersection)
+            angel / 180 * math.pi for angel in toml_parameter["ground_angels"]
-    theta1_line = np_circle_line_intersection - np.array([dgc, h_cav])
+        ]  # 地面倾角，向下为正
-    theta1 = math.atan(theta1_line[1] / theta1_line[0])
+        para.h_arm = toml_parameter["h_arm"]
-    # 考虑雷电入射角度，所以theta1可以小于0，即计算从侧面击中的雷
+        para.altitude = toml_parameter["altitude"]
-    # if theta1 < 0:
+        para.rated_voltage = toml_parameter["rated_voltage"]
-    #     # print(f"θ_1角度为负数{theta1:.4f},人为设置为0")
+        toml_advance = toml_dict["advance"]
-    #     theta1 = 0
+        para.ng = toml_advance["ng"]  # 地闪密度
-    return np.array([theta1, theta2])
+        para.Ip_a = toml_advance["Ip_a"]  # 概率密度曲线系数a
-
+        para.Ip_b = toml_advance["Ip_b"]  # 概率密度曲线系数b
-
+        toml_optional = toml_dict["optional"]
-def distance_point_line(point_x, point_y, line_x, line_y, k):
+        para.voltage_n = toml_optional["voltage_n"]  # 工作电压分成多少份来计算
-    d = abs(k * point_x - point_y - k * line_x + line_y) / ((k ** 2 + 1) ** 0.5)
+        para.max_i = toml_optional["max_i"]
-    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():
-    for u_bar in range(1):
+    if len(sys.argv) < 2:
-        u_ph = math.sqrt(2) * 750 * math.cos(2 * math.pi / 6 * 0) / 1.732  # 运行相电压
+        toml_file_path = r"内自500kV-ZCK上相.toml"
-        h_cav = 90  # 导线对地平均高
+    else:
-        h_gav = h_cav + 9.5 + 2.7
+        toml_file_path = sys.argv[1]
-        dgc = 2.9  # 导地线水平距离
+    if not os.path.exists(toml_file_path):
-        # 迭代法计算最大电流
+        logger.info(f"无法找到数据文件{toml_file_path}，程序退出。")
-        i_max = 0
+        sys.exit(0)
-        _min_i = 20  # 尝试的最小电流
+    logger.info(f"读取文件{toml_file_path}")
-        _max_i = 300  # 尝试的最大电流
+    read_parameter(toml_file_path)
-        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)
+    parameter_display(para)
-            # if not np.isreal(rs):
+    h_whole = para.h_arm[0]  # 挂点高
-            #     continue
+    string_g_len = para.string_g_len
-            rc = rc_fun(i_bar, u_ph)
+    string_c_len = para.string_c_len
-            # if not np.isreal(rc):
+    h_g_sag = para.h_g_sag
-            #     continue
+    h_c_sag = para.h_c_sag
-            rg = rg_fun(i_bar, h_cav)
+    gc_x = para.gc_x
-            # if not np.isreal(rg):
+    h_arm = para.h_arm
-            #     continue
+    gc_y = [
-            circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
+        h_whole - string_g_len - h_g_sag * 2 / 3,  # 地线对地平均高
-            if not circle_intersection:  # if circle_intersection is []
+    ]
-                continue
+    if len(h_arm) > 1:
-            circle_line_intersection = solve_circle_line_intersection(
+        for hoo in h_arm[1:]:
-                rc, rg, h_cav, dgc
+            gc_y.append(hoo - string_c_len - h_c_sag * 2 / 3)  # 导线平均高
-            )
+    if len(gc_y) > 2:  # 双回路
-            min_distance_intersection = (
+        phase_n = 3  # 边相导线数量
-                np.sum(
+    else:
-                    (np.array(circle_intersection) - np.array(circle_line_intersection))
+        phase_n = 1
-                    ** 2
+    # 地闪密度 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
+    td = para.td
+    ng = func_ng(td)
+    avr_n_sf = 0  # 考虑电压的影响计算的跳闸率
+    ground_angels = para.ground_angels
+    # 初始化动画
+    animate = Animation()
+    animate.disable(False)
+    # animate.show()
+    for ground_angel in ground_angels:
+        logger.info(f"地面倾角{ground_angel/math.pi*180:.3f}°")
+        rg_type = None
+        rg_x = None
+        rg_y = None
+        voltage_n = para.voltage_n
+        n_sf_phases = np.zeros((phase_n, voltage_n))  # 存储每一相的跳闸率
+        if np.any(np.array(gc_y) < 0):
+            logger.info("导线可能掉地面下了，程序退出。")
+            return 0
+        for phase_conductor_foo in range(phase_n):
+            exposed_curve_shielded = False
+            rs_x = gc_x[phase_conductor_foo]
+            rs_y = gc_y[phase_conductor_foo]
+            rc_x = gc_x[phase_conductor_foo + 1]
+            rc_y = gc_y[phase_conductor_foo + 1]
+            if phase_n == 1:
+                rg_type = "g"
+            if phase_n > 1:  # 多回路
+                if phase_conductor_foo < 2:
+                    rg_type = "c"  # 捕捉弧由下面一相导线的击距代替
+                    rg_x = gc_x[phase_conductor_foo + 2]
+                    rg_y = gc_y[phase_conductor_foo + 2]
+                else:
+                    rg_type = "g"
+            # TODO 保护角公式可能有问题，后面改
+            shield_angle_at_avg_height = (
+                math.atan(
+                    (rc_x - rs_x)
+                    / (
+                        (h_arm[0] - string_g_len - h_arm[phase_conductor_foo + 1])
+                        + string_c_len
+                    )
                 )
-                ** 0.5
+                * 180
-            )  # 计算两圆交点和地面直线交点的最小距离
+                / math.pi
-            i_max = i_bar
+            )  # 挂点处保护角
-            if min_distance_intersection < 0.1:
+            logger.info(f"挂点处保护角{shield_angle_at_avg_height:.3f}°")
-                break
+            logger.debug(f"最低相防护标识{rg_type}")
-        i_min = min_i(6.78, u_ph / 1.732)
+            rated_voltage = para.rated_voltage
-        cad = Draw()
+            for u_bar in range(voltage_n):  # 计算不同工作电压下的跳闸率
-        cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2)
+                # TODO 需要区分交、直流
-        cad.draw(i_max, u_ph, h_gav, h_cav, dgc, 6)
+                # u_ph = (
-        cad.save()
+                #     math.sqrt(2)
-        if abs(i_max - _max_i) < 1e-5:
+                #     * rated_voltage
-            print("无法找到最大电流，可能是杆塔较高。")
+                #     * math.cos(2 * math.pi / voltage_n * u_bar)
-            i_max = 300
+                #     / 1.732
-            print(f"最大电流设置为自然界最大电流{i_max}kA")
+                # )  # 运行相电压
-        print(f"最大电流为{i_max:.2f}")
+                u_ph = rated_voltage / 1.732
-        print(f"最小电流为{i_min:.2f}")
+                logger.info(f"计算第{phase_conductor_foo + 1}相，电压为{u_ph:.2f}kV")
-        curt_fineness = 0.1  # 电流积分细度
+                # 迭代法计算最大电流
-        if i_min > i_max or abs(i_min - i_max) < curt_fineness:
+                i_max = 0
-            print("最大电流小于最小电流，没有暴露弧，程序结束。")
+                insulator_c_len = para.insulator_c_len
-            return
+                # i_min = min_i(insulator_c_len, u_ph / 1.732)
-        # 开始积分
+                # TODO 需要考虑交、直流
-        curt_segment_n = int((i_max - i_min) / curt_fineness)  # 分成多少份
+                i_min = min_i(insulator_c_len, u_ph)
-        calculus = 0
+                _min_i = i_min  # 尝试的最小电流
-        i_curt_samples, d_curt = np.linspace(
+                _max_i = para.max_i  # 尝试的最大电流
-            i_min, i_max, curt_segment_n + 1, retstep=True
+                # 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) / 1)
+                ):  # 雷电流
+                    logger.info(f"尝试计算电流为{i_bar:.2f}")
+                    rs = rs_fun(i_bar)
+                    animate.add_rs(rs, rs_x, rs_y)
+                    rc = rc_fun(i_bar, u_ph)
+                    animate.add_rc(rc, rc_x, rc_y)
+                    rg = rg_fun(i_bar, rc_y, u_ph, typ=rg_type)
+                    rg_line_func = None
+                    if rg_type == "g":
+                        rg_line_func = rg_line_function_factory(rg, ground_angel)
+                        animate.add_rg_line(rg_line_func)
+                    rs_rc_circle_intersection = solve_circle_intersection(
+                        rs, rc, rs_x, rs_y, rc_x, rc_y
+                    )
+                    i_max = i_bar
+                    if not rs_rc_circle_intersection:  # if circle_intersection is []
+                        logger.debug("保护弧和暴露弧无交点，检查设置参数。")
+                        continue
+                    circle_rc_or_rg_line_intersection = None
+                    if rg_type == "g":
+                        circle_rc_or_rg_line_intersection = (
+                            solve_circle_line_intersection(rc, rc_x, rc_y, rg_line_func)
+                        )
+                    elif rg_type == "c":
+                        circle_rc_or_rg_line_intersection = solve_circle_intersection(
+                            rg, rc, rg_x, rg_y, rc_x, rc_y
+                        )
+                    if not circle_rc_or_rg_line_intersection:
+                        # 暴露弧和捕捉弧无交点
+                        if rg_type == "g":
+                            if rg_line_func(rc_x) > rc_y:
+                                i_min = i_bar  # 用于后面判断最小和最大电流是否相等，相等意味着暴露弧一直被屏蔽
+                                logger.info(
+                                    f"捕捉面在暴露弧之上，设置最小电流为{i_min:.2f}"
+                                )
+                            else:
+                                logger.info("暴露弧和地面捕捉弧无交点，检查设置参数。")
+                            continue
+                        else:
+                            logger.info(
+                                "上面的导地线无法保护下面的导地线，检查设置参数。"
+                            )
+                            continue
+                    animate.add_expose_area(
+                        rc_x,
+                        rc_y,
+                        *rs_rc_circle_intersection,
+                        *circle_rc_or_rg_line_intersection,
+                    )
+                    cad = Draw()
+                    cad.draw(
+                        i_min,
+                        u_ph,
+                        rs_x,
+                        rs_y,
+                        rc_x,
+                        rc_y,
+                        rg_x,
+                        rg_y,
+                        rg_type,
+                        ground_angel,
+                        2,
+                    )  # 最小电流时
+                    cad.draw(
+                        i_max,
+                        u_ph,
+                        rs_x,
+                        rs_y,
+                        rc_x,
+                        rc_y,
+                        rg_x,
+                        rg_y,
+                        rg_type,
+                        ground_angel,
+                        6,
+                    )  # 最大电流时
+                    cad.save_as(f"egm{phase_conductor_foo + 1}.dxf")
+                    min_distance_intersection = (
+                        np.sum(
+                            (
+                                np.array(rs_rc_circle_intersection)
+                                - np.array(circle_rc_or_rg_line_intersection)
+                            )
+                            ** 2
+                        )
+                        ** 0.5
+                    )  # 计算两圆交点和地面直线交点的最小距离
+                    if min_distance_intersection < 0.1:
+                        break  # 已经找到了最大电流
+                    # 判断是否以完全被保护
+                    if (
+                        rs_rc_circle_intersection[1]
+                        < circle_rc_or_rg_line_intersection[1]
+                    ):
+                        circle_rs_line_or_rg_intersection = None
+                        if rg_type == "g":
+                            circle_rs_line_or_rg_intersection = (
+                                solve_circle_line_intersection(
+                                    rs, rs_x, rs_y, rg_line_func
+                                )  # 保护弧和捕雷弧的交点
+                            )
+                        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:
+                            logger.info(f"电流为{i_bar}kV时，暴露弧已经完全被屏蔽")
+                            exposed_curve_shielded = True
+                            break
+                    animate.pause()
+                # 判断是否导线已经被完全保护
+                if abs(i_max - _max_i) < 1e-5:
+                    logger.info("无法找到最大电流，可能是杆塔较高。")
+                    logger.info(f"最大电流设置为自然界最大电流{i_max}kA")
+                logger.info(f"最大电流为{i_max:.2f}")
+                logger.info(f"最小电流为{i_min:.2f}")
+                # if exposed_curve_shielded:
+                #     logger.info("暴露弧已经完全被屏蔽,不会跳闸。")
+                #     continue
+                curt_fineness = 0.1  # 电流积分细度
+                if i_min > i_max or abs(i_min - i_max) < curt_fineness:
+                    logger.info("最大电流小于等于最小电流，没有暴露弧。")
+                    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)
+                td = para.td
+                ip_a = para.Ip_a
+                ip_b = para.Ip_b
+                bd_area_vec_result = bd_area_vec(
+                    i_curt_samples,
+                    u_ph,
+                    rc_x,
+                    rc_y,
+                    rs_x,
+                    rs_y,
+                    rg_x,
+                    rg_y,
+                    ground_angel,
+                    rg_type,
+                )
+                thunder_density_result = thunder_density(
+                    i_curt_samples, td, ip_a, ip_b
+                )  # 雷电流幅值密度函数
+                cal_bd_np = bd_area_vec_result * thunder_density_result
+                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
+                rated_voltage = para.rated_voltage
+                n_sf = (
+                    2
+                    * ng
+                    / 10
+                    * calculus
+                    * arc_possibility(rated_voltage, insulator_c_len)
+                )
+                avr_n_sf += n_sf / voltage_n
+                n_sf_phases[phase_conductor_foo][u_bar] = n_sf
+                logger.info(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.16f}次/(100km·a)")
+        logger.info(f"线路跳闸率是{avr_n_sf:.16f}次/(100km·a)")
+        logger.info(
+            f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=16)}次/(100km·a)"
         )
-        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
-            )
-        n_sf = 2 * 2.7 / 10 * calculus  # 调整率
-        print(f"跳闸率是{n_sf:.6}")
 
-    # draw(rs, rc, rg, h_gav, h_cav, dgc)
+
+def speed():
+    a = 0
+    for bar in range(100000000):
+        a += bar
 
 
 if __name__ == "__main__":
-    thunder_density(2)
+    logger.remove()
+    logger.add(sys.stderr, level="DEBUG")
     egm()
-    print("Finished.")
+    # run_time = timeit.timeit("egm()", globals=globals(), number=1)
+    # logger.info(f"运行时间:{run_time:.2f}s")
+    # input('enter any key to exit.')
+    logger.info("Finished.")

metadata.yml

@@ -0,0 +1,7 @@
+Version: 1.2.0.2
+#CompanyName: My Imaginary Company
+#FileDescription: Simple App
+#InternalName: Simple App
+#LegalCopyright: © My Imaginary Company. All rights reserved.
+#OriginalFilename: SimpleApp.exe
+ProductName: Lightening

parameter_test.py

@@ -0,0 +1,75 @@
+import datetime
+
+import pymongo
+from easyprocess import EasyProcess
+import numpy as np
+import re
+
+if __name__ == "__main__":
+    result = {}
+    utc_now = datetime.datetime.utcnow()
+    print(f"当前utc时间{utc_now}")
+    client = pymongo.MongoClient(
+        "mongodb+srv://dmy:uqXLtjkJRKzWYnSpK8jF@cluster0.x0ged.mongodb.net/db?retryWrites=true&w=majority"
+    )
+    db = client["db"]
+    result_collection = db["result_collection"]
+    with open("article_parameter.toml", encoding="utf-8") as toml_file:
+        toml_string = toml_file.read()
+    for altitude in np.arange(1500, 4000, 500):
+        for height in np.arange(100, 160, 10):
+            for ground_arm_length in np.arange(17, 25, 0.01):
+                replaced_toml_string = toml_string
+                replaced_toml_string = replaced_toml_string.replace(
+                    "ARM", f"{ground_arm_length:.3f}"
+                )
+                replaced_toml_string = replaced_toml_string.replace(
+                    "HEIGHT", str(height)
+                )
+                replaced_toml_string = replaced_toml_string.replace(
+                    "MID", str(height - 20)
+                )
+                replaced_toml_string = replaced_toml_string.replace(
+                    "ALTITUDE", str(altitude)
+                )
+                with open(
+                    "abc_parameter.toml", "w", encoding="utf-8"
+                ) as parameter_toml_file:
+                    parameter_toml_file.write(replaced_toml_string)
+                print(f"计算{altitude}m海拔，{height}m全塔高,{ground_arm_length:.3f}m横担长的跳闸率")
+                console_output = (
+                    EasyProcess(
+                        [
+                            "python",
+                            "d:/code/EGM/main.py",
+                            "d:/code/EGM/abc_parameter.toml",
+                        ]
+                    )
+                    .call()
+                    .stderr
+                )
+                tripout_rate = re.findall("不同相跳闸率是\[(.*)\]", console_output)[0]
+                print(f"跳闸率是{tripout_rate}")
+                if float(tripout_rate) < 0.14:
+                    print(f"最短地线横担是{ground_arm_length:.3f}m")
+                    if altitude not in result:
+                        result[altitude] = {}
+                    if height not in result[altitude]:
+                        result[altitude][height] = {}
+                    result[altitude][height] = ground_arm_length
+                    record_id = result_collection.insert_one(
+                        {
+                            "cate": "最短地线横担长度",
+                            "utc_time": utc_now,
+                            "altitude": float(altitude),
+                            "height": float(height),
+                            "地线横担长度": float(ground_arm_length),
+                            "跳闸率": tripout_rate,
+                        }
+                    )
+                    print(f"id{record_id}")
+                    break
+    for altitude in np.arange(1500, 4000, 500):
+        for height in np.arange(100, 160, 10):
+            ground_arm_length = result[altitude][height]
+            print(f"{altitude}m海拔，{height}m全塔高,需要最短{ground_arm_length:.3f}m横担")

plot.py

@@ -0,0 +1,100 @@
+import matplotlib
+from plot_data import *
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mticker
+
+matplotlib.use("Qt5Agg")
+# 解决中文乱码
+plt.rcParams["font.sans-serif"] = ["simsun"]
+plt.rcParams["font.family"] = "sans-serif"
+# plt.rcParams["font.weight"] = "bold"
+# 解决负号无法显示的问题
+plt.rcParams["axes.unicode_minus"] = False
+plt.rcParams["savefig.dpi"] = 1200  # 图片像素
+# plt.savefig("port.png", dpi=600, bbox_inches="tight")
+fontsize = 12
+################################################
+witdh_of_bar=0.3
+color=plt.cm.BuPu(np.linspace(152/255, 251/255,152/255))
+percent1 = data_150m塔高_不同地线保护角[:, 1] / data_150m塔高_不同地线保护角[:, 0]
+# percent1 = data_66m串长_不同塔高[:, 1] / data_66m串长_不同塔高[:, 0]
+# percent2 = data_68m串长_不同塔高[:, 1] / data_68m串长_不同塔高[:, 0]
+fig, ax = plt.subplots()
+x = np.arange(len(category_names_150m塔高_不同地线保护角))  # the label locations
+p1 = ax.bar(category_names_150m塔高_不同地线保护角, percent1, witdh_of_bar, label="绕击/反击跳闸率比值",color=color,hatch='-')
+# p1 = ax.bar(x - 0.3 / 2, percent1, 0.3, label="6.6m绝缘距离")
+# p2 = ax.bar(x + 0.3 / 2, percent2, 0.3, label="6.8m绝缘距离")
+ax.xaxis.set_major_locator(mticker.FixedLocator(x))
+ax.set_xticklabels(category_names_150m塔高_不同地线保护角)
+ax.set_ylabel("比值", fontsize=fontsize)
+ax.set_xlabel("地线保护角(°)", fontsize=fontsize)
+# ax.set_xlabel("接地电阻(Ω)", fontsize=fontsize)
+plt.xticks(fontsize=fontsize)
+plt.yticks(fontsize=fontsize)
+ax.bar_label(p1, padding=0, fontsize=fontsize)
+# ax.bar_label(p2, padding=0, fontsize=fontsize)
+ax.legend(fontsize=fontsize)
+
+fig.tight_layout()
+plt.show()
+
+
+# results = {
+#     "100m": 100 * data[0, :] / np.sum(data[0, :]),
+#     "110m": data[1, :] / np.sum(data[1, :]),
+#     "120m": data[2, :] / np.sum(data[2, :]),
+#     "130m": data[3, :] / np.sum(data[3, :]),
+#     "140m": data[4, :] / np.sum(data[4, :]),
+#     "150m": data[5, :] / np.sum(data[5, :]),
+# }
+
+
+# def survey(results, category_names):
+#     """
+#     Parameters
+#     ----------
+#     results : dict
+#         A mapping from question labels to a list of answers per category.
+#         It is assumed all lists contain the same number of entries and that
+#         it matches the length of *category_names*.
+#     category_names : list of str
+#         The category labels.
+#     """
+#     labels = list(results.keys())
+#     data = np.array(list(results.values()))
+#     data_cum = data.cumsum(axis=1)
+#     category_colors = plt.get_cmap("RdYlGn")(np.linspace(0.15, 0.85, data.shape[1]))
+#
+#     fig, ax = plt.subplots(figsize=(9.2, 5))
+#     ax.invert_yaxis()
+#     ax.xaxis.set_visible(False)
+#     ax.set_xlim(0, np.sum(data, axis=1).max())
+#
+#     for i, (colname, color) in enumerate(zip(category_names, category_colors)):
+#         widths = data[:, i]
+#         starts = data_cum[:, i] - widths
+#         rects = ax.barh(
+#             labels, widths, left=starts, height=0.5, label=colname, color=color
+#         )
+#
+#         r, g, b, _ = color
+#         text_color = "white" if r * g * b < 0.5 else "darkgrey"
+#         ax.bar_label(rects, label_type="center", color=text_color)
+#     ax.legend(
+#         ncol=len(category_names),
+#         bbox_to_anchor=(0, 1),
+#         loc="lower left",
+#         fontsize="small",
+#     )
+#
+#     return fig, ax
+
+# percent=data/np.sum(data,axis=1)[:,None]*100
+# percent = data[:, 1] / data[:, 0]
+# plt.bar(category_names, percent, 0.3, label="黑")
+# # plt.bar(category_names, percent[:,0], 0.2, label="r")
+#
+# # plt.bar(category_names, [0.014094 / 100, 0.025094 / 100], 0.2, label="h")
+# plt.legend()
+# # survey(results, category_names)
+# plt.show()

plot_data.py

@@ -0,0 +1,72 @@
+import numpy as np
+
+category_names_7欧电阻_随塔高变化 = ["100", "110", "120", "130", "140", "150"]
+# 第1列反击 ，第2列绕击
+data_7欧电阻_随塔高变化 = np.array(
+    [
+        [0.000002, 0.019094],
+        [0.000003, 0.043287],
+        [0.000006, 0.073033],
+        [0.000010, 0.103132],
+        [0.000019, 0.130923],
+        [0.000032, 0.155414],
+    ]
+)
+
+category_names_15欧电阻_随塔高变化 = ["100", "110", "120", "130", "140", "150"]
+
+data_15欧电阻_随塔高变化 = np.array(
+    [
+        [0.000039, 0.019094],
+        [0.000064, 0.043287],
+        [0.000098, 0.073033],
+        [0.000170, 0.103132],
+        [0.000287, 0.130923],
+        [0.000440, 0.155414],
+    ]
+)
+
+
+category_names_130m塔高_不同接地电阻 = ["7", "10", "15", "20", "25", "30"]
+
+data_130m塔高_不同接地电阻 = np.array(
+    [
+        [0.000010, 0.103132],
+        [0.000101, 0.103132],
+        [0.000171, 0.103132],
+        [0.000333, 0.103132],
+        [0.000563, 0.103132],
+        [0.000950, 0.103132],
+    ]
+)
+
+category_names_150m塔高_不同地线保护角 = ["-1", "-3", "-6"]
+
+data_150m塔高_不同地线保护角 = np.array(
+    [
+        [0.000440, 0.155414],
+        [0.000433, 0.128227],
+        [0.000429, 0.099819],
+    ]
+)
+
+
+category_names_66m串长_不同塔高 = ["100", "120", "140"]
+
+data_66m串长_不同塔高 = np.array(
+    [
+        [0.000053, 0.023285],
+        [0.000139, 0.083229],
+        [0.000470, 0.145586],
+    ]
+)
+
+category_names_68m串长_不同塔高 = [100, 120, 140]
+
+data_68m串长_不同塔高 = np.array(
+    [
+        [0.000039, 0.019094],
+        [0.000098, 0.073033],
+        [0.000287, 0.130923],
+    ]
+)

server.py

@@ -0,0 +1,34 @@
+from fastapi import FastAPI
+import uvicorn
+from pydantic import BaseModel
+
+
+class CalculationParameter(BaseModel):
+    loop: str  # single or double 回路数
+    insulator_c_len: float  # 绝缘长度
+    string_c_len: float  # 导线串长
+    string_g_len: float  # 地线串长
+    td: int  # 雷暴日
+    h_g_avr_sag: float  # 地线平均弧垂
+    h_c_avr_sag: float  # 导线平均弧垂
+    h_whole: float  # 杆塔全高
+    gc_x: tuple[float]  # 导、地线水平坐标
+    ground_angels: tuple[float]  # 地面倾角，向下为正
+    rated_voltage: float  # 额定电压
+
+
+fastapi_app = FastAPI()
+
+
+@fastapi_app.post("/calculation")
+async def calculation(cp: CalculationParameter):
+    print(cp)
+    return {"Hello": "World"}
+
+
+def server_start():
+    pass
+
+
+if __name__ == "__main__":
+    uvicorn.run("server:fastapi_app", host="127.0.0.1", port=5000, log_level="info")

untest.py

@@ -0,0 +1,29 @@
+import unittest
+from core import thunder_density
+from main import read_parameter
+import numpy as np
+import sympy
+
+
+class Testing(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        # read_parameter('default.toml')
+        pass
+
+    def test_thunder_density(self):
+        i = sympy.symbols("i")
+        a = np.random.random()
+        b = np.random.random()
+        p = 1-1 / (1 + (i / a) ** b)
+        d_p = sympy.diff(p, i)
+        random_i = np.random.randint(10, 100)
+        v_from_thunder_density = thunder_density(random_i, 0, a, b)
+        v_from_diff = d_p.evalf(subs={i: random_i})
+        self.assertTrue(
+            abs(v_from_thunder_density - v_from_diff) < 1e-5, "与自动微分结果不一致"
+        )  # add assertion here
+
+
+if __name__ == "__main__":
+    unittest.main()