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

.gitignore

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

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

core.py

@@ -1,6 +1,7 @@
 import math
 import ezdxf
 import numpy as np
+from typing import List
 
 gCAD = None
 gMSP = None
@@ -15,8 +16,8 @@ class Parameter:
     insulator_c_len: float  # 串子绝缘长度
     string_c_len: float
     string_g_len: float
-    gc_x: [float]  # 导、地线水平坐标
+    gc_x: List[float]  # 导、地线水平坐标
-    ground_angels: [float]  # 地面倾角，向下为正
+    ground_angels: List[float]  # 地面倾角，向下为正
     h_arm: float  # 导、地线垂直坐标
     altitude: int  # 海拔，单位米
     max_i: float  # 最大尝试电流，单位kA
@@ -131,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_x1) ** 2 + (y - center_y1) ** 2 - radius1**2,
-            (x - center_x2) ** 2 + (y - center_y2) ** 2 - radius2 ** 2,
+            (x - center_x2) ** 2 + (y - center_y2) ** 2 - radius2**2,
         ]
         X_set = np.linalg.solve(A, b)
         x += X_set[0]
@@ -166,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
+        bb = -2 * a + 2 * c * k - 2 * d * (k**2) - 2 * b * k
-        aa = 1 + k ** 2
+        aa = 1 + k**2
         rr = radius
         cc = (
-            a ** 2
+            a**2
-            + c ** 2
+            + c**2
             - 2 * c * k * d
-            + (k ** 2) * (d ** 2)
+            + (k**2) * (d**2)
             - 2 * b * (c - k * d)
-            + b ** 2
+            + b**2
-            - rr ** 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]
 
@@ -189,12 +190,20 @@ def solve_circle_line_intersection(
 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
 
 
@@ -207,12 +216,15 @@ def thunder_density(i, td, ip_a, ip_b):  # 雷电流幅值密度函数
         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单位是弧度
@@ -221,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
 
 
@@ -235,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
@@ -294,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
 
 
@@ -404,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
+                    radius**2
-                ) * (k ** 2 + 1)
+                ) * (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取一个很大的正数
@@ -447,7 +459,7 @@ def func_ng(td):  # 地闪密度，通过雷暴日计算
     if para.ng > 0:
         r = para.ng
     else:
-        r = 0.023 * (td ** 1.3)
+        r = 0.023 * (td**1.3)
     return r
 
 
@@ -456,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]
@@ -473,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
+    # TODO 需要区分交直流
-    r = (4.5 * (_e ** 0.75) - 14) * 1e-2
+    # _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,10 +1,12 @@
 import math
 import os.path
 import sys
+import time
 import tomli
 from loguru import logger
 from core import *
 import timeit
+from animation import Animation
 
 
 # 打印参数
@@ -59,7 +61,7 @@ def read_parameter(toml_file_path):
 
 def egm():
     if len(sys.argv) < 2:
-        toml_file_path = r"default.toml"
+        toml_file_path = r"内自500kV-ZCK上相.toml"
     else:
         toml_file_path = sys.argv[1]
     if not os.path.exists(toml_file_path):
@@ -70,7 +72,7 @@ def egm():
     #########################################################
     # 以上是需要设置的参数
     parameter_display(para)
-    h_whole = para.h_arm[0]  # 塔全高
+    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
@@ -82,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:
@@ -92,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):
@@ -113,13 +118,13 @@ 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 = (
+            shield_angle_at_avg_height = (
                 math.atan(
                     (rc_x - rs_x)
                     / (
@@ -129,78 +134,113 @@ def egm():
                 )
                 * 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 = (
+                # TODO 需要区分交、直流
-                    math.sqrt(2)
+                # u_ph = (
-                    * 750
+                #     math.sqrt(2)
-                    * math.cos(2 * math.pi / voltage_n * u_bar)
+                #     * rated_voltage
-                    / 1.732
+                #     * 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)
-                    #######
+                        animate.add_rg_line(rg_line_func)
-                    # cccCount += 1
+                    rs_rc_circle_intersection = solve_circle_intersection(
-                    # 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(
                         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(rs_rc_circle_intersection)
-                                - np.array(circle_rc_line_or_rg_intersection)
+                                - np.array(circle_rc_or_rg_line_intersection)
                             )
                             ** 2
                         )
@@ -210,8 +250,8 @@ def egm():
                         break  # 已经找到了最大电流
                     # 判断是否以完全被保护
                     if (
-                        rg_rc_circle_intersection[1]
+                        rs_rc_circle_intersection[1]
-                        < circle_rc_line_or_rg_intersection[1]
+                        < circle_rc_or_rg_line_intersection[1]
                     ):
                         circle_rs_line_or_rg_intersection = None
                         if rg_type == "g":
@@ -241,43 +281,16 @@ def egm():
                             logger.info(f"电流为{i_bar}kV时，暴露弧已经完全被屏蔽")
                             exposed_curve_shielded = True
                             break
-                # if phase_conductor_foo == 2:
+                    animate.pause()
-                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")
                 # 判断是否导线已经被完全保护
                 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:
+                # if exposed_curve_shielded:
-                    logger.info("暴露弧已经完全被屏蔽,不会跳闸。")
+                #     logger.info("暴露弧已经完全被屏蔽,不会跳闸。")
-                    continue
+                #     continue
                 curt_fineness = 0.1  # 电流积分细度
                 if i_min > i_max or abs(i_min - i_max) < curt_fineness:
                     logger.info("最大电流小于等于最小电流，没有暴露弧。")
@@ -291,8 +304,7 @@ def egm():
                 td = para.td
                 ip_a = para.Ip_a
                 ip_b = para.Ip_b
-                cal_bd_np = (
+                bd_area_vec_result = bd_area_vec(
-                    bd_area_vec(
                     i_curt_samples,
                     u_ph,
                     rc_x,
@@ -304,8 +316,10 @@ def egm():
                     ground_angel,
                     rg_type,
                 )
-                    * thunder_density(i_curt_samples, td, ip_a, ip_b)
+                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)
@@ -333,10 +347,10 @@ def egm():
                 )
                 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}")
+                logger.info(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.16f}次/(100km·a)")
-        logger.info(f"跳闸率是{avr_n_sf:.16f}")
+        logger.info(f"线路跳闸率是{avr_n_sf:.16f}次/(100km·a)")
         logger.info(
-            f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=16)}"
+            f"不同相跳闸率是{np.array2string(np.mean(n_sf_phases,axis=1),precision=16)}次/(100km·a)"
         )
 
 
@@ -349,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)
+    egm()
-    print(f"运行时间:{run_time:.2f}s")
+    # run_time = timeit.timeit("egm()", globals=globals(), number=1)
-    print("Finished.")
+    # logger.info(f"运行时间:{run_time:.2f}s")
+    # input('enter any key to exit.')
+    logger.info("Finished.")

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