Pₙ = 7.5e3; Uₙ = 380; f₁ = 50;  # Механическая мощность на валу [Вт], линейное напряжение [В], частота [Гц]
nₙ = 1455; p = 2; η = 0.87;     # Номинальная скорость [об/мин], пары полюсов, КПД
cosϕ = 0.83; Iₙ = 15.8;        # Номинальный коэффициент мощности и линейный ток [А]
ki = 7; mₚ = 2.3; mₘ = 2.3; J = 0.02; # Пусковой ток/момент, максимальный момент, инерция [кг·м²]

# Расчёт параметров схемы замещения
Sₙ = Pₙ/(η*cosϕ); # Номинальная полная электрическая мощность [ВА]
U₁ = Uₙ/sqrt(3);                           # Номинальное фазное напряжение статора [В]
n₀ = 60*f₁/p;                             # Синхронная частота вращения [об/мин]
sₙ = (n₀-nₙ)/n₀; sₖ = sₙ*(mₘ+sqrt(mₘ^2-1)); # Номинальное и критическое скольжение
ω₀ = 2π*f₁/p; ωₙ = π*nₙ/30;                # Синхронная и номинальная углова скорость ротора [рад/с]
Mₙ = Pₙ/ωₙ; Mₘ = Mₙ*mₘ; Mₚ = Mₙ*mₚ;           # Номинальный, максимальный и пусковой момент [Н·м]
pₘ = 0.05*Pₙ;                              # Механические потери (5% от Pₙ) [Вт]

R₂ = 1/3*(Pₙ+pₘ)/(Iₙ^2*(1-sₙ)/sₙ);              # Приведённое активное сопротивление роторной обмотки [Ом]

C = 1.02;                   # Конструктивный коэффициент (предварительный)
    
R₁ = U₁*cosϕ*(1-η)/Iₙ - C^2*R₂ - pₘ/(3*Iₙ^2);  # Активное сопротивление статорной обмотки [Ом]

L₁σ = L₂σ = U₁/(4π*f₁*(1+C^2)*ki*Iₙ);          # Индуктивности рассеивания статора и приведённое ротора [Гн]
L₁ = L₂ = U₁/(2π*f₁*Iₙ*sqrt(1-cosϕ^2)-2/3*(2π*f₁*Mₘ*sₙ)/(p*U₁*sₖ)); # Индуктивности статора и приведённое ротора [Гн]
Lₘ = L₁ - L₁σ;                                 # Индуктивность магнитной цепи [Гн]

for i in 1:5

    C = 1 + L₁σ/Lₘ                             # Конструктивный коэффициент (уточнённый)
    
    R₁ = U₁*cosϕ*(1-η)/Iₙ - C^2*R₂ - pₘ/(3*Iₙ^2);  # Активное сопротивление статорной обмотки [Ом]

    L₁σ = L₂σ = U₁/(4π*f₁*(1+C^2)*ki*Iₙ);          # Индуктивности рассеивания статора и приведённое ротора [Гн]
    L₁ = L₂ = U₁/(2π*f₁*Iₙ*sqrt(1-cosϕ^2)-2/3*(2π*f₁*Mₘ*sₙ)/(p*U₁*sₖ)); # Индуктивности статора и приведённое ротора [Гн]
    Lₘ = L₁ - L₁σ;                                 # Индуктивность магнитной цепи [Гн]

    println("Итерация №$i: Уточнённый конструктивный коэффициент C равен $C")
end

    X₁σ = 2π*f₁*L₁σ;    X₂σ = 2π*f₁*L₂σ        # Реактивные сопротивления рассеивания статора и приведённое ротора [Ом]
    Xₘ = 2π*f₁*Lₘ;

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

simout = get_sim_results("im_parametrization.engee", example_path) #запускаем моделирование

t = get_sim_data(simout, "Ток статора", 1, "time");
Is = get_sim_data(simout, "Ток статора", 1, "value");
n = get_sim_data(simout, "Обороты (у.е)", 1, "value");

I_base = Sₙ/(3*Uₙ)              # База фазного тока для Δ-параметризации блока
I_A = Is .* I_base              # Terminal current из измерителя: о.е. -> А
n_rpm = n .* n₀                 # Rotor velocity: о.е. -> об/мин

gr(aspectratio=:auto, xlims=:auto, ylims=:auto, size=(900,400))
I_graph = plot(t, I_A; label="I(t)", title="Ток статора", ylabel="I [А]", xlabel="t [с]")
hline!(I_graph, [Iₙ]; linestyle=:dash, color=:red, label="In = $(Iₙ) А")
n_graph = plot(t, n_rpm; label="n(t)", title="Частота вращения ротора", ylabel="n [об/мин]", xlabel="t [с]")
hline!(n_graph, [nₙ]; linestyle=:dash, color=:green, label="nn = $(nₙ) об/мин")
plot(I_graph, n_graph)

# Электромеханическая характеристика в области номинального режима.
# Важно: S_base соответствует параметру S_rated модели, а Is — выходу Terminal current в о.е.
S_base = Pₙ/(η*cosϕ)
I_base = S_base/(3*Uₙ)
I_A = Is .* I_base
n_rpm = n .* n₀

# Установившаяся точка: усреднение за последние 0,05 с.
steady = t .>= (t[end] - 0.05)
I_ss = sum(I_A[steady])/sum(steady)
n_ss = sum(n_rpm[steady])/sum(steady)

# 1. Левый (основной) график
b = plot(I_A, n_rpm;
    label="n(I)", 
    title="Электромеханическая характеристика",
    xlabel="Ток статора, А", 
    ylabel="Скорость, об/мин",
    legend=:outerbottom)

x = [10, 10, 25, 25, 10]
y = [0.9 * n₀, 1.01 * n₀, 1.01 * n₀, 0.9 * n₀, 0.9 * n₀]
plot!(b, x, y; label=:none, linewidth=2)
# 2. Правый график (крупный план рабочей точки)
g = plot(I_A, n_rpm; 
    label = "Рабочая точка",
    title = "Увеличенная область",
    xlabel = "Ток статора, А",
    ylabel = "Скорость, об/мин",
    xlims=(10, 25), 
    ylims=(0.9 * n₀, 1.01 * n₀),
    legend=:outerbottom, legend_columns=1)

# Добавляем линии и точку на ПРАВЫЙ график (g)
vline!(g, [Iₙ]; linestyle=:dash, linewidth=2, color=:red,
    label="Номинальный ток: $(Iₙ) А")

hline!(g, [nₙ]; linestyle=:dash, linewidth=2, color=:green,
    label="Номинальная скорость: $(nₙ) об/мин")

scatter!(g, [I_ss], [n_ss]; color=:black, markersize=5,
    label="Модель: $(round(I_ss; digits=3)) А, $(round(n_ss; digits=3)) об/мин")

plot(b, g, layout = (1, 2), size = (1000, 500))