#=
Параметризация радиатора О123.
Эта ячейка повторяет расчет параметров, которые вводятся в блок `Радиатор`:
масса, удельная теплоемкость, вектор перепадов температур и вектор отводимой мощности.
Степенная зависимость используется для естественной конвекции; линейная модель оставлена
как контрольная оценка по постоянному тепловому сопротивлению.
=#

using Printf

# O123的护照参数
P_nominal = 120.0      # 功率损耗，W
Rth_natural = 0.71     # 自然冷却时的热阻，°C/W
Rth_forced = 0.212     # 强制冷却时的热阻，°C/W
mass = 2.0             # 散热器重量，公斤
cp_aluminum = 900.0    # 铝合金比热容，J/(kg*K)

println("散热器o123的主要参数:")
println("Масса: $mass кг")
println("Удельная теплоемкость: $cp_aluminum Дж/(кг*К)")
println("Номинальная мощность: $P_nominal Вт")
println("Тепловое сопротивление, естественное охлаждение: $Rth_natural °C/Вт")
println("Тепловое сопротивление, принудительное охлаждение: $Rth_forced °C/Вт")

# 模型的动态部分：热容量设置加热的惯性。
heatsink_heat_capacity = mass * cp_aluminum

println()
println("===动态特性===")
println("散热器重量: ", mass, " 公斤")
println("比热容: ", cp_aluminum, " Дж/(кг*К)")
println("散热器容量: ", heatsink_heat_capacity, " Дж/К")

# 来自额定功率和Rth的自然对流的参考点。
delta_T_nominal = P_nominal * Rth_natural

println()
println("===参考点的计算===")
println("额定功率P= ", P_nominal, " 星期二")
println("热阻Rth= ", Rth_natural, " °C/Вт")
println("温差ΔT=P*Rth= ", round(delta_T_nominal, digits=2), " K")
println("参考点：（ΔT= ", round(delta_T_nominal, digits=2), " K, Q = ", P_nominal, " 星期二)")

# 对于自然对流，我们使用幂律Q=k*ΔT^n。
n = 1.25
k = P_nominal / (delta_T_nominal^n)

println()
println("===功率依赖参数===")
println("N=的指数 ", n)
println("系数k= ", round(k, digits=4))
println("公式：Q= ", round(k, digits=4), " * ΔT^", n)

# 向量ΔT对应于模型中辐射器的表格参数化。
delta_T_vector = [10.0, 30.0, 50.0, 70.0, 90.0]
Q_vector = [k * (dT^n) for dT in delta_T_vector]

println()
println("===用于模型参数化的向量===")
println("温差矢量ΔT,K:")
println(delta_T_vector)
println("输出功率矢量Q,W:")
println(round.(Q_vector, digits=1))
println()
println("数值表:")
println("＜ΔT,K＜Q,W＜")
println("├───────┼─────────┤")
for (i, dT) in enumerate(delta_T_vector)
    @printf("│ %6.1f │ %7.1f │\n", dT, Q_vector[i])
end

# 验证：在参考点，功率模型应返回P_nominal。
Q_check = k * (delta_T_nominal^n)
error = abs(Q_check - P_nominal)

println()
println("===检查计算===")
println("在ΔT= ", round(delta_T_nominal, digits=2), " K:")
println("  估计功率Q= ", round(Q_check, digits=2), " 星期二")
println("  额定功率= ", P_nominal, " 星期二")
println("  错误率= ", round(error, digits=4), " 星期二")

println()
println("="^60)
println("模型的汇总参数")
println("="^60)
println("参数化方法：表格依赖性Q（ΔT）")
println("对流类型：自然")
println("向量ΔT,K:[", join(round.(delta_T_vector, digits=1), ", "), "]")
println("向量Q,W:[", join(round.(Q_vector, digits=1), ", "), "]")
println("散热器重量: ", mass, " 公斤")
println("比热容: ", cp_aluminum, " Дж/(кг*К)")

# 线性模型仅用于与非线性自然对流进行比较。
Q_linear = [dT / Rth_natural for dT in delta_T_vector]

println()
println("===线性模型进行比较===")
println("热阻Rth= ", Rth_natural, " °C/Вт")
println("线性模型的向量Q，W:")
println(round.(Q_linear, digits=1))
println()
println("△ΔT,K△q幂律,W△Q线性,W△")
println("├───────┼─────────────────┼────────────────┤")
for (i, dT) in enumerate(delta_T_vector)
    @printf("│ %6.1f │ %15.1f │ %14.1f │\n", dT, Q_vector[i], Q_linear[i])
end

string(
    "Cth=", round(heatsink_heat_capacity; digits=1),
    " J/K; delta_T_nom=", round(delta_T_nominal; digits=2),
    " K; Q_vector=", join(round.(Q_vector; digits=1), ","),
)

#=
Здесь собраны только аналитические оценки. Их цель - быстро проверить,
что паспортные данные и параметры блоков согласованы до запуска модели.
Сами формулы вынесены в текст выше, а ниже оставлены численные подстановки.
=#

# 护照和模型参数的计算
using Printf

f = 50.0
T = 1 / f
U_rms = 600.0
U_peak = sqrt(2) * U_rms

ITAVM = 200.0
ITRMS = 315.0
UT0 = 1.05
rT = 1.17e-3
P_TAV = UT0 * ITAVM + rT * ITRMS^2

Rth_jc = 0.080
Rth_ch = 0.020
Rth_ha_nat = 0.710
Rth_ha_forced = 0.212
Rth_ja_nat = Rth_jc + Rth_ch + Rth_ha_nat
Rth_ja_forced = Rth_jc + Rth_ch + Rth_ha_forced

T_ambient = 25.0
Tj_est_nat = T_ambient + P_TAV * Rth_ja_nat
Tj_est_forced = T_ambient + P_TAV * Rth_ja_forced

R_load = 6.0
L_load = 5e-3
X_L = 2*pi*f*L_load
Z_load = sqrt(R_load^2 + X_L^2)
tau_load = L_load / R_load

Rs = 100.0
Cs = 1e-4
tau_snubber = Rs * Cs

alpha_min = 90.0 - 82.0
alpha_max = 90.0 + 82.0
delay_min = alpha_min * T / 360
delay_max = alpha_max * T / 360

Vgate_source = 15.0
VGT = 3.5
IGT = 0.20
R2 = 7.0
Rg = 10.0
Ig_est = (Vgate_source - VGT) / (R2 + Rg)

rows = [
    ("Период сети T", T, "с"),
    ("Амплитуда сети Vm", U_peak, "В"),
    ("Потери тиристора при ITAVM/ITRMS", P_TAV, "Вт"),
    ("Rthja естественное", Rth_ja_nat, "К/Вт"),
    ("Rthja принудительное", Rth_ja_forced, "К/Вт"),
    ("Оценка Tj при естественном охлаждении", Tj_est_nat, "°C"),
    ("Оценка Tj при принудительном охлаждении", Tj_est_forced, "°C"),
    ("Индуктивное сопротивление нагрузки XL", X_L, "Ом"),
    ("Модуль импеданса нагрузки |Z|", Z_load, "Ом"),
    ("Постоянная времени нагрузки L/R", tau_load, "с"),
    ("Постоянная времени RC-цепи", tau_snubber, "с"),
    ("Минимальная задержка СИФУ", delay_min, "с"),
    ("Максимальная задержка СИФУ", delay_max, "с"),
    ("Оценка тока управления gate", Ig_est, "А"),
    ("Запас к 3*IGT", Ig_est / (3*IGT), "раз"),
]

for (name, value, unit) in rows
    @printf("%-46s = %12.6g %s\n", name, value, unit)
end

string("U_peak=", round(U_peak; digits=2), " V; P_TAV=", round(P_TAV; digits=2), " W; Ig/(3IGT)=", round(Ig_est/(3IGT); digits=3))

#=
Запуск модели.
Эта ячейка выполняет три задачи: готовит окружение,
запускает существующую модель Engee и приводит результаты к массивам.
=#

include("$(@__DIR__)/useful_functions.jl")

# 在上传之前，我们关闭同名模型，以免混合不同启动的结果。

try_close_model(model_name)
model_path = resolve_existing_file(model_file)
model = engee.load(model_path; force=true)
engee.open(model)

# 我们只记得状态栏中报告的初始设置。

original_stop_time = safe_model_param(model, "StopTime", "")
original_fixed_step = safe_model_param(model, "FixedStep", "")
original_solver = safe_model_param(model, "SolverName", "")

status_message = ""

# 当前会话中计算的时间限制。 不保存模型文件。
safe_set_model_param!(model, "StopTime", string(simulation_stop_time))

# 当前会话中的信号的临时记录。 不保存原始模型文件。
log_ports = [
    "Генератор синусоиды/1",
    "Датчик напряжения/1",
    "Датчик тока/1",
    "Смещение/1", # 加热单元的温度从K转换到°C后
]

for port in log_ports
    safe_set_log(port)
end

# 我们在配置日志记录后运行模型;AnyMath结果被读取为一组命名信号。

results = engee.run(model)

# 我们提取三个关键信号：打开角度，负载电压和负载电流。

(angle_xy, angle_key) = get_signal(results, ["alpha", "Генератор синусоиды.1", "Генератор синусоиды.main_out", "Генератор синусоиды"])
(vload_xy, vload_key) = get_signal(results, ["Vload", "Датчик напряжения.Vload", "Датчик напряжения.1", "Датчик напряжения.v_out"])
(iload_xy, iload_key) = get_signal(results, ["Iload", "Датчик тока.Iload", "Датчик тока.1", "Датчик тока.i_out"])

# 在可视化之前，我们将瞬时信号细化到1000点。

t_alpha, alpha = thin(angle_xy[1], angle_xy[2]; max_points=max_plot_points)
t_v, vload = thin(vload_xy[1], vload_xy[2]; max_points=max_plot_points)
t_i, iload = thin(iload_xy[1], iload_xy[2]; max_points=max_plot_points)

# 我们计算一个网络周期的平均值,并将结果细化.

v_avg_full = moving_average_by_time(vload_xy[1], vload_xy[2], network_period)
i_avg_full = moving_average_by_time(iload_xy[1], iload_xy[2], network_period)
t_vavg, v_avg = thin(vload_xy[1], v_avg_full; max_points=max_plot_points)
t_iavg, i_avg = thin(iload_xy[1], i_avg_full; max_points=max_plot_points)

try_close_model(model_name)

# 该单元显示相位控制后的瞬时波形。
p1 = plot(t_v, [vload, iload]; label=["Vload" "Iload"], xlabel="t，与", ylabel="U，B;I，A", title="负载电压和电流", legend=:topright, size=(800, 400))
period = 600:800
p2 = plot(t_v[period], [vload[period], iload[period]]; label=["Vload" "Iload"], xlabel="t，与", ylabel="U，B;I，A", title="负载电压和电流", legend=:topright, size=(800, 400))
plot(p1, p2; layout=(2, 1), size=(800, 400))

# 在此，将周期的平均值与缓慢变化的打开角进行比较。

p3 = plot(t_vavg, [v_avg, alpha]; label=["平均(Vload),20毫秒窗口" "alpha"], xlabel="t，与", ylabel="Uavg，In;alpha，deg", title="平均电压", legend=:topright)
p4 = plot(t_iavg, [i_avg, alpha]; label=["平均（Iload），20毫秒窗口" "alpha"], xlabel="t，与", ylabel="Iavg，A;alpha，deg", title="平均电流", legend=:topright)
plot(p3, p4; layout=(2, 1), size=(800, 400))

# Alpha坐标中的图表显示了负载调整特性。

plot(alpha[40:end], [v_avg[40:end], i_avg[40:end]]; label=["Vload" "Iload"], xlabel="阿尔法，冰雹", ylabel="Uavg，B;Iavg，A", title="调整特性", legend=:topright, size=(800, 400))

# 该单元提取并绘制晶闸管外壳/热节点的温度。
# 名称列表已保持扩展，因为AnyMath可以通过行的名称返回信号。,
# 块名称或自动生成的端口名称。
(temp_xy, temp_key) = get_signal(results, [
    "temp_VT_C",
    "tempC",
    "Смещение.temp_C",
    "Смещение.1",
    "Смещение.main_out",
    "Датчик абсолютной температуры.1",
])

t_temp, temp_case = thin(temp_xy[1], temp_xy[2]; max_points=max_plot_points)
temp_rise = temp_case .- first(temp_case)

p_temp1 = plot(
    t_temp,
    temp_case;
    label=temp_key,
    xlabel="t，与",
    ylabel="Tcase, °C",
    title="晶闸管外壳的温度",
    legend=:topright,
)

p_temp2 = plot(
    t_temp,
    temp_rise;
    label="ΔTcase",
    xlabel="t，与",
    ylabel="ΔT, °C",
    title="壳体相对于初始值的过热",
    legend=:topright,
)

plot(p_temp1, p_temp2; layout=(2, 1), size=(800, 400))