1.增加Unittest

2.考虑利用实际地闪密度和雷电流曲线。
3.版本号v1.2.0

core.py

@@ -21,6 +21,9 @@ class Parameter:
     altitude: int  # 海拔，单位米
     max_i: float  # 最大尝试电流，单位kA
     rated_voltage: float  # 额定电压
+    ng: float  # 地闪密度  次/(每平方公里·每年)
+    Ip_a: float  # 概率密度曲线系数a
+    Ip_b: float  # 概率密度曲线系数b
 
 
 para = Parameter()
@@ -186,7 +189,7 @@ def solve_circle_line_intersection(
 def min_i(string_len, u_ph):
     # 海拔修正
     altitude = para.altitude
-    k_a = math.exp((altitude-1000) / 8150)  # 气隙海拔修正
+    k_a = math.exp((altitude - 1000) / 8150)  # 气隙海拔修正
     u_50 = 1 / k_a * (530 * string_len + 35)  # 50045 上附录的公式，实际应该用负极性电压的公式
     z_0 = 300  # 雷电波阻抗
     z_c = 251  # 导线波阻抗
@@ -195,13 +198,20 @@ def min_i(string_len, u_ph):
     return r
 
 
-def thunder_density(i):  # l雷电流幅值密度函数
+def thunder_density(i, td, ip_a, ip_b):  # 雷电流幅值密度函数
-    td = para.td
+    # td = para.td
     r = None
-    if td == 20:
+    # ip_a = para.Ip_a
-        r = -(10 ** (-i / 44)) * math.log(10) * (-1 / 44)  # 雷暴日20d
+    # ip_b = para.Ip_b
-    if td == 40:
+    if ip_a > 0 and ip_b > 0:
-        r = -(10 ** (-i / 88)) * math.log(10) * (-1 / 88)  # 雷暴日40d
+        r = -(
+            -ip_b / ip_a / ((1 + (i / ip_a) ** ip_b) ** 2) * ((i / ip_a) ** (ip_b - 1))
+        )
+    else:
+        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
 
 
@@ -433,7 +443,7 @@ 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):  # 地闪密度
+def func_ng(td):  # 地闪密度，通过雷暴日计算
     return 0.023 * (td ** 1.3)
 
 

main.py

@@ -8,30 +8,29 @@ import timeit
 
 
 # 打印参数
-def parameter_display(para: Parameter):
+def parameter_display(para_dis: Parameter):
-    logger.info(f"额定电压 kV {para.rated_voltage}")
+    logger.info(f"额定电压 kV {para_dis.rated_voltage}")
-    logger.info(f"导线弧垂 m {para.h_c_sag}")
+    logger.info(f"导线弧垂 m {para_dis.h_c_sag}")
-    logger.info(f"地线弧垂 m {para.h_g_sag}")
+    logger.info(f"地线弧垂 m {para_dis.h_g_sag}")
-    logger.info(f"全塔高 m {para.h_arm[0]}")
+    logger.info(f"全塔高 m {para_dis.h_arm[0]}")
-    logger.info(f"串绝缘距离 m {para.insulator_c_len}")
+    logger.info(f"串绝缘距离 m {para_dis.insulator_c_len}")
-    logger.info(f"导线串长 m {para.string_c_len}")
+    logger.info(f"导线串长 m {para_dis.string_c_len}")
-    logger.info(f"地线串长 m {para.string_g_len}")
+    logger.info(f"地线串长 m {para_dis.string_g_len}")
-    logger.info(f"挂点垂直坐标 m {para.h_arm}")
+    logger.info(f"挂点垂直坐标 m {para_dis.h_arm}")
-    logger.info(f"挂点水平坐标 m {para.gc_x}")
+    logger.info(f"挂点水平坐标 m {para_dis.gc_x}")
-    logger.info(f"地面倾角 ° {[an*180/math.pi for an in para.ground_angels]}")
+    logger.info(f"地面倾角 ° {[an * 180 / math.pi for an in para_dis.ground_angels]}")
-    logger.info(f"海拔高度 m {para.altitude}")
+    logger.info(f"海拔高度 m {para_dis.altitude}")
-    logger.info(f"雷暴日 d {para.td}")
+    if para_dis.ng > 0:
-
+        logger.info("不采用雷暴日计算地闪密度和雷电流密度")
-
+        logger.info(f"地闪密度 次/(每平方公里·每年) {para_dis.ng}")
-def egm():
+        logger.info(f"概率密度曲线系数a {para_dis.Ip_a}")
-    if len(sys.argv) < 2:
+        logger.info(f"概率密度曲线系数b {para_dis.Ip_b}")
-        toml_file_path = r"article.toml"
+        pass
     else:
-        toml_file_path = sys.argv[1]
+        logger.info(f"雷暴日 d {para_dis.td}")
-    if not os.path.exists(toml_file_path):
+
-        logger.info(f'无法找到数据文件{toml_file_path}，程序退出。')
+
-        sys.exit(0)
+def read_parameter(toml_file_path):
-    logger.info(f"读取文件{toml_file_path}")
     with open(toml_file_path, "rb") as toml_fs:
         toml_dict = tomli.load(toml_fs)
         toml_parameter = toml_dict["parameter"]
@@ -49,9 +48,25 @@ def egm():
         para.h_arm = toml_parameter["h_arm"]
         para.altitude = toml_parameter["altitude"]
         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"default.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)
     #########################################################
     # 以上是需要设置的参数
     parameter_display(para)
@@ -223,7 +238,7 @@ 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:
@@ -273,6 +288,9 @@ def egm():
                     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
                 cal_bd_np = (
                     bd_area_vec(
                         i_curt_samples,
@@ -286,7 +304,7 @@ def egm():
                         ground_angel,
                         rg_type,
                     )
-                    * thunder_density(i_curt_samples)
+                    * thunder_density(i_curt_samples, td, ip_a, ip_b)
                 )
                 calculus = np.sum(cal_bd_np[:-1] + cal_bd_np[1:]) / 2 * d_curt
                 # for i_curt in i_curt_samples[:-1]:

metadata.yml

@@ -1,4 +1,4 @@
-Version: 1.1.0
+Version: 1.2.0
 #CompanyName: My Imaginary Company
 #FileDescription: Simple App
 #InternalName: Simple App

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