egm/core.py

import math
import ezdxf
import numpy as np

gCAD = None
gMSP = None
gCount = 1


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, h_gav, h_cav, dgc, 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, h_cav)
        msp.add_circle((0, h_gav), rs, dxfattribs={"color": color})
        msp.add_line((0, 0), (0, h_gav))  # 地线
        msp.add_circle((dgc, h_cav), rc, dxfattribs={"color": color})
        msp.add_line((dgc, 0), (dgc, h_cav))  # 导线
        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(radius, rg, center_x, center_y):
    distance = distance_point_line(center_x, center_y, 0, rg, 0)  # 捕雷线到暴露圆中点的距离
    if distance > radius:
        return []
    else:
        r = (radius ** 2 - (rg - center_y) ** 2) ** 0.5 + center_x
        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 * (np.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)
    # r=14.7*(i**0.42)
    return r


def rg_fun(i_curt, h_cav):
    if h_cav < 40:
        rg = (3.6 + 1.7 ** math.log(43 - h_cav)) * (i_curt ** 0.65)
    else:
        rg = 5.5 * (i_curt ** 0.65)
    return rg


def intersection_angle(dgc, h_gav, h_cav, i_curt, u_ph):  # 暴露弧的角度
    rs = rs_fun(i_curt)
    rc = rc_fun(i_curt, u_ph)
    rg = rg_fun(i_curt, h_cav)
    circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)  # 两圆的交点
    circle_line_intersection = solve_circle_line_intersection(
        rc, rg, dgc, h_cav
    )  # 暴露圆和补雷线的交点
    np_circle_intersection = np.array(circle_intersection)
    # TODO to be removed
    if not circle_intersection:
        abc = 123
    theta2_line = np_circle_intersection - np.array([dgc, h_cav])
    theta2 = math.atan(theta2_line[1] / theta2_line[0])
    np_circle_line_intersection = np.array(circle_line_intersection)
    theta1_line = np_circle_line_intersection - np.array([dgc, h_cav])
    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
    r_pw = 0
    # TODO to be removed
    if int(max_w / w_fineness) < 0:
        abc = 123
        pass
    w_samples, d_w = np.linspace(0, max_w, int(max_w / w_fineness), retstep=True)
    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, dgc, h_cav, h_gav):  # 对θ进行积分
    max_w = 0
    # 求暴露弧上一点的切线
    line_x = math.cos(theta) * rc + dgc
    line_y = math.sin(theta) * rc + h_cav
    k = math.tan(theta + math.pi / 2)  # 入射角
    # 求保护弧到直线的距离，判断是否相交
    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)
        # TODO to be removed
        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
        # TODO to be removed
        if max_w < 0:
            abc = 123
            tangent_line_k(line_x, line_y, 0, h_gav, rs, init_k=k)
        # 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
    r = rc / math.cos(theta) * func_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), 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, 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


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
            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
                )
                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
    # 把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):  # 地闪密度
    return 0.023 * (td ** 1.3)