egm/core.py

import math
import ezdxf
import numpy as np
from typing import List

gCAD = None
gMSP = None
gCount = 1


class Parameter:
    h_g_sag: float  # 地线弧垂
    h_c_sag: float  # 导线弧垂
    voltage_n: int  # 工作电压分成多少份来计算
    td: int  # 雷暴日
    insulator_c_len: float  # 串子绝缘长度
    string_c_len: float
    string_g_len: 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()


def rg_line_function_factory(_rg, ground_angel):  # 返回一个地面捕雷线的直线方程
    y_d = _rg / math.cos(ground_angel)  # y轴上的截距
    # 利用公式y-y0=k(x-x0) 得到直线公式
    y0 = y_d
    x0 = 0
    k = math.tan(math.pi - ground_angel)

    def f(x):
        return y0 + k * (x - x0)

    return f


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,
        rs_x,
        rs_y,
        rc_x,
        rc_y,
        rg_x,
        rg_y,
        rg_type,
        ground_angel,
        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, rc_y, u_ph, typ=rg_type)
        msp.add_circle((rs_x, rs_y), rs, dxfattribs={"color": color})
        msp.add_line((0, 0), (rs_x, rs_y))  # 地线
        msp.add_circle((rc_x, rc_y), rc, dxfattribs={"color": color + 2})
        msp.add_line((rc_x, 0), (rc_x, rc_y))  # 导线
        msp.add_line((rs_x, rs_y), (rc_x, rc_y))
        # 角度线
        circle_intersection = solve_circle_intersection(rs, rc, rs_x, rs_y, rc_x, rc_y)
        msp.add_line(
            (rc_x, rc_y), circle_intersection, dxfattribs={"color": color}
        )  # 地线
        if rg_type == "g":
            ground_angel_func = rg_line_function_factory(rg, ground_angel)
            msp.add_line(
                (0, ground_angel_func(0)),
                (2000, ground_angel_func(2000)),
                dxfattribs={"color": color},
            )
            circle_line_section = solve_circle_line_intersection(
                rc, rc_x, rc_y, ground_angel_func
            )
            if not circle_line_section:
                pass
            else:
                msp.add_line(
                    (rc_x, rc_y), circle_line_section, dxfattribs={"color": color}
                )  # 导线和圆的交点
        if rg_type == "c":
            msp.add_circle((rg_x, rg_y), rg, dxfattribs={"color": color})
            rg_rc_intersection = solve_circle_intersection(
                rg, rc, rg_x, rg_y, rc_x, rc_y
            )
            if rg_rc_intersection:
                msp.add_line(
                    (rc_x, rc_y), rg_rc_intersection, dxfattribs={"color": color}
                )  # 圆和圆的交点
        # 计算圆交点

        # msp.add_line((dgc, h_cav), circle_intersection)  # 导线

    def save(self):
        doc = self._doc
        doc.saveas("egm.dxf")

    def save_as(self, file_name):
        doc = self._doc
        doc.saveas(file_name)


# 圆交点
def solve_circle_intersection(
    radius1,
    radius2,
    center_x1,
    center_y1,
    center_x2,
    center_y2,
):
    # 用牛顿法求解
    x = radius2 + center_x2  # 初始值
    y = radius2 + center_y2  # 初始值
    # TODO 考虑出现2个解的情况
    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_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, center_x, center_y, ground_surface_func
):  # 返回交点的x和y坐标
    x0 = 0
    y0 = ground_surface_func(x0)
    x1 = 1
    y1 = ground_surface_func(x1)
    k = (y1 - y0) / (x1 - x0)
    distance = distance_point_line(center_x, center_y, x0, y0, k)  # 捕雷线到暴露圆中点的距离
    if distance > radius:
        return []
    else:
        # r = (radius ** 2 - (rg - center_y) ** 2) ** 0.5 + center_x
        a = center_x
        b = center_y
        c = y0
        d = x0
        bb = -2 * a + 2 * c * k - 2 * d * (k**2) - 2 * b * k
        aa = 1 + k**2
        rr = radius
        cc = (
            a**2
            + c**2
            - 2 * c * k * d
            + (k**2) * (d**2)
            - 2 * b * (c - k * d)
            + b**2
            - rr**2
        )
        _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
        assert abs(equ) < 1e-5
        return [_x, _y]


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


def thunder_density(i, td, ip_a, ip_b):  # 雷电流幅值密度函数
    # td = para.td
    r = None
    # 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单位是弧度
    r = 0.75 * abs((np.cos(angle - math.pi / 2) ** 3))  # 新版大手册公式3-275
    return r


def rs_fun(i):
    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) ** 1.125)  # 新版大手册公式3-272
    return r


# typ 如果是g，代表捕雷线公式，c代表暴露弧公式
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
        else:
            rg = 5.5 * (i_curt**0.65)  # 新版大手册公式3-273
    elif typ == "c":  # 此时返回的是圆半径
        rg = rc_fun(i_curt, u_ph)
    return rg


def intersection_angle(
    rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_angel, rg_type
):  # 暴露弧的角度
    rs = rs_fun(i_curt)
    rc = rc_fun(i_curt, u_ph)
    rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
    circle_intersection = solve_circle_intersection(
        rs, rc, rs_x, rs_y, rc_x, rc_y
    )  # 两圆的交点
    circle_line_or_rg_intersection = None
    rg_line_func = rg_line_function_factory(rg, ground_angel)
    if rg_type == "g":
        circle_line_or_rg_intersection = solve_circle_line_intersection(
            rc, rc_x, rc_y, rg_line_func
        )  # 暴露圆和补雷线的交点
    if rg_type == "c":
        circle_line_or_rg_intersection = solve_circle_intersection(
            rg, rc, rg_x, rg_y, rc_x, rc_y
        )  # 两圆的交点
    # TODO 应该是不存在落到地面线以下的情况，先把以下注释掉

    # if circle_line_or_rg_intersection:
    #     (
    #         circle_line_or_rg_intersection_x,
    #         circle_line_or_rg_intersection_y,
    #     ) = circle_line_or_rg_intersection
    #     if (
    #         ground_surface(rg, circle_line_or_rg_intersection_x)
    #         > circle_line_or_rg_intersection_y
    #     ):  # 交点在地面线以下，就可以不积分
    #         # 找到暴露弧和地面线的交点
    #         circle_line_or_rg_intersection = circle_ground_surface_intersection(
    #             rc, rc_x, rc_y, ground_surface
    #         )
    theta1 = None
    np_circle_intersection = np.array(circle_intersection)
    theta2_line = np_circle_intersection - np.array([rc_x, rc_y])
    theta2 = math.atan(theta2_line[1] / theta2_line[0])
    np_circle_line_or_rg_intersection = np.array(circle_line_or_rg_intersection)
    if not circle_line_or_rg_intersection:
        if rc_y - rc > rg_line_func(rc_x):  # rg在rc下面
            # 捕捉线太低了，对高塔无保护，θ_1从-90°开始计算，即从与地面垂直的角度开始就已经暴露了。
            theta1 = -math.pi / 2
    else:
        theta1_line = np_circle_line_or_rg_intersection - np.array([rc_x, rc_y])
        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
    segments = int(max_w / w_fineness)
    if segments < 2:  # 最大最小太小，没有可以积分的
        return 0
    w_samples, d_w = np.linspace(0, max_w, segments, retstep=True)
    # 童中宇 750KV信洛Ⅰ线雷电防护性能研究  公式 3-10
    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, rc_x, rc_y, rs_x, rs_y):  # 对θ进行积分
    max_w = 0
    # 求暴露弧上一点的切线
    line_x = math.cos(theta) * rc + rc_x
    line_y = math.sin(theta) * rc + rc_y
    k = math.tan(theta + math.pi / 2)  # 入射角
    # 求保护弧到直线的距离，判断是否相交
    d_to_rs = distance_point_line(rs_x, rs_y, line_x, line_y, k)
    if d_to_rs < rs:  # 相交
        # 要用过这一点到保护弧的切线
        new_k = tangent_line_k(line_x, line_y, rs_x, rs_y, 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
        # 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
    # 童中宇 750KV信洛Ⅰ线雷电防护性能研究  公式 3-10
    r = rc / math.cos(theta) * func_calculus_pw(theta, max_w)
    return r


def bd_area(
    i_curt, u_ph, rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, ground_angel, rg_type
):  # 暴露弧的投影面积
    theta1, theta2 = intersection_angle(
        rc_x, rc_y, rs_x, rs_y, rg_x, rg_y, i_curt, u_ph, ground_angel, rg_type
    )  # θ角度
    theta_fineness = 0.01
    rc = rc_fun(i_curt, u_ph)
    rs = rs_fun(i_curt)
    rg = rg_fun(i_curt, rc_y, u_ph, typ=rg_type)
    theta_segments = int((theta2 - theta1) / theta_fineness)
    if theta_segments < 2:
        return 0
    theta_sample, d_theta = np.linspace(theta1, theta2, theta_segments, 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, rc_x, rc_y, rs_x, rs_y)
    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
            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)

                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
                # 解决数值稳定性
                if bar == max_iteration - 1:
                    if abs(math.atan(k)) * 180 / math.pi > 89:
                        k_candidate[ind] = k
    # 把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):  # 地闪密度，通过雷暴日计算
    if para.ng > 0:
        r = para.ng
    else:
        r = 0.023 * (td**1.3)
    return r


# 圆和地面线的交点，只去正x轴上的。
def circle_ground_surface_intersection(radius, center_x, center_y, ground_surface):
    # 最笨的办法，一个个去试
    x_series = np.linspace(0, radius, int(radius / 0.001)) + center_x
    part_to_be_squared = (
        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
    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]
    r_y = ground_surface(r_x)
    return r_x, r_y


# u_ph是相电压
# insulator_c_len绝缘子闪络距离
def arc_possibility(rated_voltage, insulator_c_len):  # 建弧率
    # 50064 中附录给的公式
    # 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