from core import *
import timeit


def egm():
    # TODO to be removed
    cccCount = 0
    avr_n_sf = 0  # 考虑电压的影响
    voltage_n = 3  # 工作电压分成多少份来计算
    ng = func_ng(20)
    h_whole = 140  # 杆塔全高
    insulator_c_len = 6.8  # 串子绝缘长度
    string_c_len = 9.2
    string_g_len = 0.5
    dgc = -0.9  # 导地线水平距离
    vertical_dgc = 2.7  # 导地线挂点垂直距离
    h_g_avr_sag = 11.67 * 2 / 3
    h_c_avr_sag = (14.43 - 11.67) * 2 / 3
    h_gav = h_whole - string_g_len - h_g_avr_sag  # 地线对地平均高
    h_cav = h_gav - string_c_len - vertical_dgc - h_c_avr_sag  # 导线对地平均高
    shield_angle = math.atan(dgc / (vertical_dgc + string_c_len)) * 180 / math.pi
    print(f"保护角{shield_angle:.3f}°")
    for u_bar in range(voltage_n):
        u_ph = (
            math.sqrt(2) * 750 * math.cos(2 * math.pi / voltage_n * u_bar) / 1.732
        )  # 运行相电压
        # 迭代法计算最大电流
        i_max = 0
        i_min = min_i(insulator_c_len, u_ph / 1.732)
        _min_i = i_min  # 尝试的最小电流
        _max_i = 200  # 尝试的最大电流
        # TODO remove it
        cad = Draw()
        cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2)
        for i_bar in np.linspace(_min_i, _max_i, int((_max_i - _min_i) / 0.1)):  # 雷电流
            # print(f"尝试计算电流为{i_bar:.2f}")
            rs = rs_fun(i_bar)
            rc = rc_fun(i_bar, u_ph)
            rg = rg_fun(i_bar, h_cav)
            #######
            cccCount += 1
            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)
            #######
            circle_intersection = solve_circle_intersection(rs, rc, h_gav, h_cav, dgc)
            if not circle_intersection:  # if circle_intersection is []
                # print("保护弧和暴露弧无交点，检查设置参数。程序退出。")
                continue
            circle_rc_line_intersection = solve_circle_line_intersection(
                rc, rg, dgc, h_cav
            )
            if not circle_rc_line_intersection:
                continue
            min_distance_intersection = (
                np.sum(
                    (
                        np.array(circle_intersection)
                        - np.array(circle_rc_line_intersection)
                    )
                    ** 2
                )
                ** 0.5
            )  # 计算两圆交点和地面直线交点的最小距离
            i_max = i_bar
            if min_distance_intersection < 0.1:
                break
            if circle_intersection[1] < circle_rc_line_intersection[1]:
                circle_rs_line_intersection = solve_circle_line_intersection(
                    rs, rg, 0, h_gav
                )
                # 判断与保护弧的交点是否在暴露弧外面
                distance = (
                    np.sum(
                        (np.array(circle_rs_line_intersection) - np.array([dgc, h_cav]))
                        ** 2
                    )
                    ** 0.5
                )
                if distance > rc:
                    print("暴露弧已经完全被屏蔽")
                    break
        cad.draw(i_min, u_ph, h_gav, h_cav, dgc, 2)
        cad.draw(i_max, u_ph, h_gav, h_cav, dgc, 6)
        cad.save()
        # 判断是否导线已经被完全保护
        if abs(i_max - _max_i) < 1e-5:
            print("无法找到最大电流，可能是杆塔较高。")
            print(f"最大电流设置为自然界最大电流{i_max}kA")
        print(f"最大电流为{i_max:.2f}")
        print(f"最小电流为{i_min:.2f}")
        curt_fineness = 0.1  # 电流积分细度
        if i_min > i_max or abs(i_min - i_max) < curt_fineness:
            print("最大电流小于最小电流，没有暴露弧，程序结束。")
            return
        # 开始积分
        curt_segment_n = int((i_max - i_min) / curt_fineness)  # 分成多少份
        calculus = 0
        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
            )
        n_sf = (
            2 * ng / 10 * calculus
        )  # 跳闸率 利用Q╱GDW 11452-2015 架空输电线路防雷导则的公式 Ng=0.023*Td^(1.3) 20天雷暴日地闪密度为1.13
        avr_n_sf += n_sf / voltage_n
        print(f"工作电压为{u_ph:.2f}kV时,跳闸率是{n_sf:.6}")
    print(f"跳闸率是{avr_n_sf:.6}")


def speed():
    a = 0
    for bar in range(100000000):
        a += bar


if __name__ == "__main__":
    run_time = timeit.timeit("egm()", globals=globals(), number=1)
    print(f"运行时间:{run_time:.2f}s")
    print("Finished.")