Uc = sqrt(2)*573  # Амплитудное фазное напряжение сети (В)
fc = 50.0         # Частота сети (Гц)
Rc = 0.0528e-3; Xc = 0.468e-3; Lc = Xc/(2pi*50)  # Параметры кабеля: R (Ом), X (Ом), L (Гн)
Rt1 = 0.05e-3; Xt1 = 0.35e-3; Lt1 = Xt1/(2pi*50)  # Параметры трансформатора - первичная обмотка
Rt2 = 0.25e-3; Xt2 = 3.15e-3; Lt2 = Xt2/(2pi*50)  # Параметры трансформатора - вторичная обмотка

println("Параметры системы: Rc = $Rc Ом, Lc = $Lc Гн")
println("Параметры Трансформатора ВН: Rt1 = $Rt1 Ом, Lt1 = $Lt1 Гн")
println("Параметры Трансформатора НН: Rt2 = $Rt2 Ом, Lt2 = $Lt2 Гн")

iзаж = 0.5e4; Imax = 1e5; # Ток зажигания и максимальный ток, А
uзаж = 310.5; uп = uзаж/1.15; # Напряжение зажигания и напряжение постоянной величины, В
uт = (Imax+iзаж)/Imax * uп; # Вспомогательный параметр
T1 = 1e-2; T2 = 2e-2; # Тепловые постоянные времени дуги, с
Rs1 = 25.87*1e-3; # Сопротивление дуги в зарядном прмежутке, Ом

# Установка необходимых пакетов (если они ещё не установлены)
# JLD2 - для работы с файлами данных в формате HDF5/Julia
# FFTW - библиотека для быстрого преобразования Фурье (FFT)
# LinearAlgebra - стандартный модуль Julia для линейной алгебры
import Pkg; Pkg.add(["JLD2", "FFTW", "LinearAlgebra"])

# Подключение установленных пакетов
using JLD2, FFTW, LinearAlgebra

example_path = @__DIR__; # Получаем абсолютный путь к директории, содержащей текущий скрипт
cd(example_path); # Переходим в директорию примера
include("useful_functions.jl"); # Подключаем скрипт Julia со вспомогательными функциями

# @markdown **Определите источник данных для анализа:**
источник = "модель .engee" # @param ["модель .engee","файл .jld2"]
simout = get_sim_results(источник, example_path)

t = get_sim_data(simout, "Iarc_pu", 1, "time");
IarcA_pu = get_sim_data(simout, "Iarc_pu", 1, "value");
VarcA_pu = get_sim_data(simout, "Varc_pu", 1, "value");
IarcA = get_sim_data(simout, "Iarc", 1, "value");
VarcA = get_sim_data(simout, "Varc", 1, "value");
VcA = get_sim_data(simout, "Vc", 1, "value");

T = last(t);
ST = T/(size(t)[1]-1);
fs = Int(round(1/ST));

gr()
plot(t, [IarcA_pu, VarcA_pu];
     label = ["Ток, о.е." "Напряжение, о.е."],
     title = "Ток и напряжение ДСП-100НЗА в фазе А",
     xlabel = "Время, с")

plot(IarcA_pu, VarcA_pu;
    label = :none, title = "ВАХ ДСП-100НЗА",
    xlabel = "Ток Ia, о.е.", ylabel = "Напряжение Va, о.е.")

gr()
int = 2290:4301
plot(t[int],  VarcA[int];
    label = "Напряжение, В",
    title = "Период напряжения ДСП-100НЗА в фазе А",
    xlabel = "Время, с")

signal = VarcA[int]; t = t[int];

# Анализируем гармоники
results = analyze_harmonics(signal, fs, fc)

n, A, ph = get_harmonics_data(results)

# Выводим результаты
println("Гармоника | Амплитуда | Фаза (град)")
println("----------------------------------")
for (n, A, ph) in zip(n, A, ph)
    round(A, digits=3)!=0.0 ? println(lpad(n, 9), " | ", lpad(round(A, digits=3), 9), " | ", round((ph)*180/pi, digits=3)) : nothing
end

harm, harm2, harm3 = get_sequences(results, n);

gr()
plot( yaxis = :log, xlims=[0,55],ylims=[1e-1,1e3], title = "Амплитуды гармоник 1-55 напряжения в фазе А", xlabel="Порядок гармоник n", ylabel = "Логарифм амплитуды гармоники lg(An)")
plot!(harm[1], harm[2];
        markershape = :square, markersize=3, label = "Прямая последовательность")
plot!(harm2[1], harm2[2];
        markershape = :circle, markersize=3, label = "Обратная последовательность")
plot!(harm3[1], harm3[2];
        markershape = :ltriangle, markersize=5, label = "Тройная последовательность")
scatter!([harm[1][1]], [harm[2][1]];
        markershape = :star, markersize=7, label = "Базовая гармоника")

gr()
y = range(-π, π, length=9)
plot(xlims=[0,54], yticks = (y, ["$(round(yi, digits=2)) рад" for yi in y]), title = "Фазы гармоник 1-55 напряжения в фазе А", xlabel="Порядок гармоник n", ylabel = "Фаза гармоники φn")
plot!(harm[1], harm[3];
        markershape = :square, markersize=3, label = "Прямая последовательность")
plot!(harm2[1], harm2[3];
        markershape = :circle, markersize=3, label = "Обратная последовательность")
plot!(harm3[1], harm3[3];
        markershape = :ltriangle, markersize=5, label = "Тройная последовательность")
scatter!([harm[1][1]], [harm[3][1]];
        markershape = :star, markersize=7, label = "Базовая гармоника")

THD_V = norm(vcat(harm[2][2:end], harm2[2], harm3[2]))/harm[2][1]*100;
THD_V0 = norm(harm[2][2:end])/harm[2][1]*100;
THD_V2 = norm(harm2[2])/harm[2][1]*100;
THD_V3 = norm(harm3[2])/harm[2][1]*100;

println("Суммарный коэффициент гармонических искажений напряжения: ", round(THD_V, digits=2), "%.\nИз него, по последовательностям:")
println("- прямая: ", round(THD_V0, digits=2), "%,")
println("- обратная: ", round(THD_V2, digits=2), "%,")
println("- тройная: ", round(THD_V3, digits=2), "%.")

results_I = analyze_harmonics(IarcA, fs, fc);
n_I, A_I, ph_I = get_harmonics_data(results_I);
THD_I = norm(A_I[2:end])/A_I[1]*100

println("Суммарный коэффициент гармонических искажений тока: ", round(THD_I, digits=2), "%.")