From 77951ae54afa6c31aa99b5c84de5d090bf9b6522 Mon Sep 17 00:00:00 2001
From: facat <facat@facat.com>
Date: Tue, 21 Sep 2021 01:35:42 +0800
Subject: [PATCH] =?UTF-8?q?=E5=BE=AA=E7=8E=AF=E5=90=91=E9=87=8F=E5=8C=96?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 core.py | 62 ++++++++++++++++++++++++++++-----------------------------
 main.py | 36 +++++++++++++++++++++------------
 2 files changed, 53 insertions(+), 45 deletions(-)

diff --git a/core.py b/core.py
index 9996628..d011366 100644
--- a/core.py
+++ b/core.py
@@ -83,7 +83,7 @@ def thunder_density(i):  # l雷电流幅值密度函数
 
 
 def angel_density(angle):  # 入射角密度函数 angle单位是弧度
-    r = 0.75 * (math.cos(angle - math.pi / 2) ** 3)
+    r = 0.75 * (np.cos(angle - math.pi / 2) ** 3)
     return r
 
 
@@ -126,7 +126,7 @@ def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph):  # 暴露弧的角度
     return np.array([theta1, theta2])
 
 
-def distance_point_line(point_x, point_y, line_x, line_y, k):
+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
 
@@ -139,15 +139,8 @@ def func_calculus_pw(theta, max_w):
         abc = 123
         pass
     w_samples, d_w = np.linspace(0, max_w, int(max_w / w_fineness), retstep=True)
-    for cal_w in w_samples[:-1]:
-        r_pw += (
-            (
-                abs(angel_density(cal_w)) * math.sin(theta - (cal_w - math.pi / 2))
-                + abs(angel_density(cal_w + d_w))
-                * math.sin(theta - (cal_w + d_w - math.pi / 2))
-            )
-            / 2
-        ) * d_w
+    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
 
 
@@ -175,18 +168,18 @@ def calculus_bd(theta, rc, rs, rg, dgc, h_cav, h_gav):  # 对θ进行积分
             abc = 123
             tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k)
         # TODO to be removed
-        global 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
+        # 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
@@ -215,15 +208,20 @@ def bd_area(i_curt, u_ph, dgc, h_gav, h_cav):  # 暴露弧的投影面积
     theta_sample, d_theta = np.linspace(
         theta1, theta2, int((theta2 - theta1) / theta_fineness), retstep=True
     )
-    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
-        )
+    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, dgc, h_cav, h_gav)
+    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
 
 
diff --git a/main.py b/main.py
index b00ab9b..5524e29 100644
--- a/main.py
+++ b/main.py
@@ -1,3 +1,5 @@
+import numpy as np
+
 from core import *
 import timeit
 
@@ -104,19 +106,27 @@ def egm():
         i_curt_samples, d_curt = np.linspace(
             i_min, i_max, curt_segment_n + 1, retstep=True
         )
-        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
-            )
+        bd_area_vec = np.vectorize(bd_area)
+        cal_bd_np = bd_area_vec(
+            i_curt_samples, u_ph, dgc, h_gav, h_cav
+        ) * thunder_density(i_curt_samples)
+        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
         n_sf = (
             2 * ng / 10 * calculus
         )  # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13