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

@@ -1,6 +1,7 @@
 import math
 import ezdxf
 import numpy as np
+from typing import List
 
 gCAD = None
 gMSP = None
@@ -10,17 +11,20 @@ gCount = 1
 class Parameter:
     h_g_sag: float  # 地线弧垂
     h_c_sag: float  # 导线弧垂
-    h_whole: float  # 杆塔全高
     voltage_n: int  # 工作电压分成多少份来计算
     td: int  # 雷暴日
     insulator_c_len: float  # 串子绝缘长度
     string_c_len: float
     string_g_len: float
-    gc_x: [float]  # 导、地线水平坐标
-    ground_angels: [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()
@@ -128,14 +132,14 @@ def solve_circle_intersection(
     x = radius2 + center_x2  # 初始值
     y = radius2 + center_y2  # 初始值
     # TODO 考虑出现2个解的情况
-    for bar in range(0, 10):
+    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 - 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]
@@ -163,22 +167,22 @@ def solve_circle_line_intersection(
         b = center_y
         c = y0
         d = x0
-        bb = -2 * a + 2 * c * k - 2 * d * (k ** 2) - 2 * b * k
-        aa = 1 + k ** 2
+        bb = -2 * a + 2 * c * k - 2 * d * (k**2) - 2 * b * k
+        aa = 1 + k**2
         rr = radius
         cc = (
-            a ** 2
-            + c ** 2
+            a**2
+            + c**2
             - 2 * c * k * d
-            + (k ** 2) * (d ** 2)
+            + (k**2) * (d**2)
             - 2 * b * (c - k * d)
-            + b ** 2
-            - rr ** 2
+            + b**2
+            - rr**2
         )
-        _x = (-bb + (bb ** 2 - 4 * aa * cc) ** 0.5) / 2 / aa
+        _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
+        equ = (center_x - _x) ** 2 + (center_y - _y) ** 2 - radius**2
         assert abs(equ) < 1e-5
         return [_x, _y]
 
@@ -186,23 +190,41 @@ def solve_circle_line_intersection(
 def min_i(string_len, u_ph):
     # 海拔修正
     altitude = para.altitude
-    k_a = math.exp(altitude / 8150)  # 气隙海拔修正
+    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):  # l雷电流幅值密度函数
-    td = para.td
+def thunder_density(i, td, ip_a, ip_b):  # 雷电流幅值密度函数
+    # td = para.td
     r = None
-    if td == 20:
-        r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44)  # 雷暴日20d
-    if td == 40:
-        r = -(10 ** (-i / 88)) * math.log(10) * (-1 / 88)  # 雷暴日40d
-    return r
+    # 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单位是弧度
@@ -211,12 +233,12 @@ def angel_density(angle):  # 入射角密度函数 angle单位是弧度
 
 
 def rs_fun(i):
-    r = 10 * (i ** 0.65)  # 新版大手册公式3-271
+    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.125)  # 新版大手册公式3-272
+    r = 1.63 * ((5.015 * (i**0.578) - 0.001 * u_ph * 1) ** 1.125)  # 新版大手册公式3-272
     return r
 
 
@@ -225,9 +247,9 @@ 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
+            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
+            rg = 5.5 * (i_curt**0.65)  # 新版大手册公式3-273
     elif typ == "c":  # 此时返回的是圆半径
         rg = rc_fun(i_curt, u_ph)
     return rg
@@ -284,7 +306,7 @@ def intersection_angle(
 
 # 点到直线的距离
 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)
+    d = abs(k * point_x - point_y - k * line_x + line_y) / ((k**2 + 1) ** 0.5)
     return d
 
 
@@ -394,14 +416,14 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
             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)
+                    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
+                    - 2 * (radius**2) * k
                 )
                 if abs(d_fk) < 1e-5 and abs(line_x - center_x - radius) < 1e-5:
                     # k不存在，角度为90°，k取一个很大的正数
@@ -433,8 +455,12 @@ def tangent_line_k(line_x, line_y, center_x, center_y, radius, init_k=10.0):
     return np.array(k_candidate)[np.max(k_angle) == k_angle].tolist()[-1]
 
 
-def func_ng(td):  # 地闪密度
-    return 0.023 * (td ** 1.3)
+def func_ng(td):  # 地闪密度，通过雷暴日计算
+    if para.ng > 0:
+        r = para.ng
+    else:
+        r = 0.023 * (td**1.3)
+    return r
 
 
 # 圆和地面线的交点，只去正x轴上的。
@@ -442,12 +468,12 @@ def circle_ground_surface_intersection(radius, center_x, center_y, ground_surfac
     # 最笨的办法，一个个去试
     x_series = np.linspace(0, radius, int(radius / 0.001)) + center_x
     part_to_be_squared = (
-        radius ** 2 - (x_series - center_x) ** 2
+        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
+    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]
@@ -459,6 +485,8 @@ def circle_ground_surface_intersection(radius, center_x, center_y, ground_surfac
 # insulator_c_len绝缘子闪络距离
 def arc_possibility(rated_voltage, insulator_c_len):  # 建弧率
     # 50064 中附录给的公式
-    _e = rated_voltage / (3 ** 0.5) / insulator_c_len
-    r = (4.5 * (_e ** 0.75) - 14) * 1e-2
+    # 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,42 +1,44 @@
 import math
+import os.path
 import sys
+import time
 import tomli
 from loguru import logger
 from core import *
 import timeit
+from animation import Animation
 
 
-def egm():
-    if len(sys.argv) < 2:
-        toml_file_path = "default.toml"
-        # # logger.info('没指定计算文件!程序结束。')
-        # # sys.exit(0)
-        # h_g_sag = 20  # 地线弧垂
-        # h_c_sag = 20  # 导线弧垂
-        # h_whole = 106.1  # 杆塔全高
-        # voltage_n = 3  # 工作电压分成多少份来计算
-        # td = 20  # 雷暴日
-        # insulator_c_len = 6.98  # 串子绝缘长度
-        # string_c_len = 9.2
-        # string_g_len = 0.63
-        # gc_x = [32.2 / 2, 32.2 / 2, 15, 17.0]  # 导、地线水平坐标
-        # ground_angels = [40 / 180 * math.pi]  # 地面倾角，向下为正
-        # h_arm = 34
-        # gc_y = [
-        #     h_whole - string_g_len - h_g_sag * 2 / 3,  # 地线对地平均高
-        #     # h_whole - string_c_len - h_c_sag - 2.7,  # 导线对地平均高
-        #     h_whole - string_c_len - h_c_sag * 2 / 3 - h_arm,  # 导线对地平均高
-        #     # h_whole - string_c_len - h_c_sag - 35.7,  # 导线对地平均高
-        # ]
+# 打印参数
+def parameter_display(para_dis: Parameter):
+    logger.info(f"额定电压 kV {para_dis.rated_voltage}")
+    logger.info(f"导线弧垂 m {para_dis.h_c_sag}")
+    logger.info(f"地线弧垂 m {para_dis.h_g_sag}")
+    logger.info(f"全塔高 m {para_dis.h_arm[0]}")
+    logger.info(f"串绝缘距离 m {para_dis.insulator_c_len}")
+    logger.info(f"导线串长 m {para_dis.string_c_len}")
+    logger.info(f"地线串长 m {para_dis.string_g_len}")
+    logger.info(f"挂点垂直坐标 m {para_dis.h_arm}")
+    logger.info(f"挂点水平坐标 m {para_dis.gc_x}")
+    logger.info(f"地面倾角 ° {[an * 180 / math.pi for an in para_dis.ground_angels]}")
+    logger.info(f"海拔高度 m {para_dis.altitude}")
+    if para_dis.ng > 0:
+        logger.info("不采用雷暴日计算地闪密度和雷电流密度")
+        logger.info(f"地闪密度 次/(每平方公里·每年) {para_dis.ng}")
+        logger.info(f"概率密度曲线系数a {para_dis.Ip_a}")
+        logger.info(f"概率密度曲线系数b {para_dis.Ip_b}")
+        pass
     else:
-        toml_file_path = sys.argv[1]
-    logger.info(f"读取文件{toml_file_path}")
+        logger.info(f"雷暴日 d {para_dis.td}")
+
+
+def read_parameter(toml_file_path):
     with open(toml_file_path, "rb") as toml_fs:
         toml_dict = tomli.load(toml_fs)
         toml_parameter = toml_dict["parameter"]
         para.h_g_sag = toml_parameter["h_g_sag"]  # 地线弧垂
         para.h_c_sag = toml_parameter["h_c_sag"]  # 导线弧垂
-        para.h_whole = toml_parameter["h_whole"]  # 杆塔全高
+        # para.h_whole = toml_parameter["h_whole"]  # 杆塔全高
         para.td = toml_parameter["td"]  # 雷暴日
         para.insulator_c_len = toml_parameter["insulator_c_len"]  # 串子绝缘长度
         para.string_c_len = toml_parameter["string_c_len"]
@@ -47,12 +49,30 @@ def egm():
         ]  # 地面倾角，向下为正
         para.h_arm = toml_parameter["h_arm"]
         para.altitude = toml_parameter["altitude"]
-        para.max_i = toml_parameter["max_i"]
+        para.rated_voltage = toml_parameter["rated_voltage"]
+        toml_advance = toml_dict["advance"]
+        para.ng = toml_advance["ng"]  # 地闪密度
+        para.Ip_a = toml_advance["Ip_a"]  # 概率密度曲线系数a
+        para.Ip_b = toml_advance["Ip_b"]  # 概率密度曲线系数b
         toml_optional = toml_dict["optional"]
         para.voltage_n = toml_optional["voltage_n"]  # 工作电压分成多少份来计算
+        para.max_i = toml_optional["max_i"]
+
+
+def egm():
+    if len(sys.argv) < 2:
+        toml_file_path = r"内自500kV-ZCK上相.toml"
+    else:
+        toml_file_path = sys.argv[1]
+    if not os.path.exists(toml_file_path):
+        logger.info(f"无法找到数据文件{toml_file_path}，程序退出。")
+        sys.exit(0)
+    logger.info(f"读取文件{toml_file_path}")
+    read_parameter(toml_file_path)
     #########################################################
     # 以上是需要设置的参数
-    h_whole = para.h_whole
+    parameter_display(para)
+    h_whole = para.h_arm[0]  # 挂点高
     string_g_len = para.string_g_len
     string_c_len = para.string_c_len
     h_g_sag = para.h_g_sag
@@ -64,7 +84,7 @@ def egm():
     ]
     if len(h_arm) > 1:
         for hoo in h_arm[1:]:
-            gc_y.append(hoo - string_c_len - h_c_sag * 2 / 3)
+            gc_y.append(hoo - string_c_len - h_c_sag * 2 / 3)  # 导线平均高
     if len(gc_y) > 2:  # 双回路
         phase_n = 3  # 边相导线数量
     else:
@@ -74,12 +94,15 @@ def egm():
     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
-        cad = Draw()
         voltage_n = para.voltage_n
         n_sf_phases = np.zeros((phase_n, voltage_n))  # 存储每一相的跳闸率
         if np.any(np.array(gc_y) < 0):
@@ -95,88 +118,129 @@ def egm():
                 rg_type = "g"
             if phase_n > 1:  # 多回路
                 if phase_conductor_foo < 2:
-                    rg_type = "c"  # 捕捉弧有下面一相导线的击距代替
+                    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 = (
-                math.atan((rc_x - rs_x) / ((rs_y - rc_y) + string_c_len))
+            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
+                    )
+                )
                 * 180
                 / math.pi
-            )  # 保护角
-            logger.info(f"保护角{shield_angle:.3f}°")
+            )  # 挂点处保护角
+            logger.info(f"挂点处保护角{shield_angle_at_avg_height:.3f}°")
             logger.debug(f"最低相防护标识{rg_type}")
+            rated_voltage = para.rated_voltage
             for u_bar in range(voltage_n):  # 计算不同工作电压下的跳闸率
-                u_ph = (
-                    math.sqrt(2)
-                    * 750
-                    * math.cos(2 * math.pi / voltage_n * u_bar)
-                    / 1.732
-                )  # 运行相电压
+                # TODO 需要区分交、直流
+                # u_ph = (
+                #     math.sqrt(2)
+                #     * rated_voltage
+                #     * math.cos(2 * math.pi / voltage_n * u_bar)
+                #     / 1.732
+                # )  # 运行相电压
+                u_ph = rated_voltage / 1.732
                 logger.info(f"计算第{phase_conductor_foo + 1}相，电压为{u_ph:.2f}kV")
                 # 迭代法计算最大电流
                 i_max = 0
                 insulator_c_len = para.insulator_c_len
-                i_min = min_i(insulator_c_len, u_ph / 1.732)
+                # i_min = min_i(insulator_c_len, u_ph / 1.732)
+                # TODO 需要考虑交、直流
+                i_min = min_i(insulator_c_len, u_ph)
                 _min_i = i_min  # 尝试的最小电流
                 _max_i = para.max_i  # 尝试的最大电流
                 # 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)
+                    _min_i, _max_i, int((_max_i - _min_i) / 1)
                 ):  # 雷电流
-                    # logger.info(f"尝试计算电流为{i_bar:.2f}")
+                    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)
-                    #######
-                    # 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(
+                        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 rg_rc_circle_intersection:  # if circle_intersection is []
+                    if not rs_rc_circle_intersection:  # if circle_intersection is []
                         logger.debug("保护弧和暴露弧无交点，检查设置参数。")
                         continue
-                    circle_rc_line_or_rg_intersection = None
+                    circle_rc_or_rg_line_intersection = None
                     if rg_type == "g":
-                        circle_rc_line_or_rg_intersection = (
+                        circle_rc_or_rg_line_intersection = (
                             solve_circle_line_intersection(rc, rc_x, rc_y, rg_line_func)
                         )
                     elif rg_type == "c":
-                        circle_rc_line_or_rg_intersection = solve_circle_intersection(
+                        circle_rc_or_rg_line_intersection = solve_circle_intersection(
                             rg, rc, rg_x, rg_y, rc_x, rc_y
                         )
-                    if not circle_rc_line_or_rg_intersection:
+                    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}")
+                                logger.info(
+                                    f"捕捉面在暴露弧之上，设置最小电流为{i_min:.2f}"
+                                )
                             else:
                                 logger.info("暴露弧和地面捕捉弧无交点，检查设置参数。")
                             continue
                         else:
-                            logger.info("上面的导地线无法保护下面的导地线，检查设置参数。")
+                            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(rg_rc_circle_intersection)
-                                - np.array(circle_rc_line_or_rg_intersection)
+                                np.array(rs_rc_circle_intersection)
+                                - np.array(circle_rc_or_rg_line_intersection)
                             )
                             ** 2
                         )
@@ -186,8 +250,8 @@ def egm():
                         break  # 已经找到了最大电流
                     # 判断是否以完全被保护
                     if (
-                        rg_rc_circle_intersection[1]
-                        < circle_rc_line_or_rg_intersection[1]
+                        rs_rc_circle_intersection[1]
+                        < circle_rc_or_rg_line_intersection[1]
                     ):
                         circle_rs_line_or_rg_intersection = None
                         if rg_type == "g":
@@ -214,46 +278,19 @@ def egm():
                             ** 0.5
                         )
                         if distance > rc:
-                            logger.info("暴露弧已经完全被屏蔽")
+                            logger.info(f"电流为{i_bar}kV时，暴露弧已经完全被屏蔽")
                             exposed_curve_shielded = True
                             break
-                # if phase_conductor_foo == 2:
-                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")
+                    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
+                # 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("最大电流小于等于最小电流，没有暴露弧。")
@@ -264,21 +301,25 @@ def egm():
                     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_angel,
-                        rg_type,
-                    )
-                    * thunder_density(i_curt_samples)
+                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)
@@ -296,13 +337,20 @@ def egm():
                 #     if abs(calculus-0.05812740052770032)<1e-5:
                 #         abc=123
                 #         pass
-                n_sf = 2 * ng / 10 * calculus * arc_possibility(750, insulator_c_len)
+                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:.6}")
-        logger.info(f"跳闸率是{avr_n_sf:.6f}")
+                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=6)}"
+            f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=16)}次/(100km·a)"
         )
 
 
@@ -315,6 +363,8 @@ def speed():
 if __name__ == "__main__":
     logger.remove()
     logger.add(sys.stderr, level="DEBUG")
-    run_time = timeit.timeit("egm()", globals=globals(), number=1)
-    print(f"运行时间:{run_time:.2f}s")
-    print("Finished.")
+    egm()
+    # 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()

server.py

@@ -14,6 +14,7 @@ class CalculationParameter(BaseModel):
     h_whole: float  # 杆塔全高
     gc_x: tuple[float]  # 导、地线水平坐标
     ground_angels: tuple[float]  # 地面倾角，向下为正
+    rated_voltage: float  # 额定电压
 
 
 fastapi_app = FastAPI()

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()