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;   # Коэффициенты максимального тока, пускового и максимального момента, момент инерции двигателя [кг·м²]

# Расчёт параметров схемы замещения
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");

gr(aspectratio=:auto, xlims=:auto, ylims=:auto, size = (900,400))
I_graph = plot(t, Is; label = "I(t)", title = "Ток статора", ylabel = "I [A]", xlabel = "t [с]")
n_graph = plot(t, n; label = "n(t)", title = "Частота вращения ротора", ylabel = "n [у.е.]", xlabel = "t [с]")
plot(I_graph, n_graph)

gr(aspectratio=:auto, xlims=:auto, ylims=:auto, size = (900, 400))
plot(Is, n.*nₙ; label = "n(I)", title = "Электромеханическая харк-ка", ylabel = "n [об/мин]", xlabel = "I [А]")