m = "torsional_oscillations"
function start_model_engee(m)
    try
        engee.close(m, force=true) # закрытие модели если открыта
    catch err # в случае, если нет модели, которую нужно закрыть
        # модель не была открыта
    end
    
    try
        m = engee.load("$((@__DIR__))/$m.engee") # загрузка модели
        engee.run(m) # запуск модели
    catch err # в случае, если модель не загружена
        m = engee.load("$((@__DIR__))/$m.engee") # загрузка модели
        engee.run(m) # запуск модели
    end
end

import Pkg; 
Pkg.add(["DSP", "Peaks"])

using FFTW, DSP

function analyze_natural_frequencies(time_data, speed_data)
    # Подготовка данных для FFT
    dt = time_data[2] - time_data[1]
    fs = 1/dt  # Частота дискретизации
    
    # Применение оконной функции для лучшего разрешения
    windowed_data = speed_data .* hanning(length(speed_data))
    
    # Вычисление FFT
    fft_result = fft(windowed_data)
    
    # Частотная ось
    frequencies = fftfreq(length(fft_result), fs)
    
    # Амплитудный спектр в dB
    magnitude_db = 20 * log10.(abs.(fft_result) .+ 1e-10)
    
    return frequencies, magnitude_db
end

using Peaks

function find_natural_frequencies(frequencies, magnitude_db, min_height=20)
    # Поиск пиков в положительной части спектра  
    pos_freq_idx = frequencies .> 0
    pos_freqs = frequencies[pos_freq_idx]
    pos_magnitudes = magnitude_db[pos_freq_idx]
    
    # Поиск локальных максимумов
    peaks, _ = findmaxima(pos_magnitudes)
    
    # Фильтрация по высоте
    high_peaks = peaks[pos_magnitudes[peaks] .> min_height]
    
    natural_freqs = pos_freqs[high_peaks]
    amplitudes = pos_magnitudes[high_peaks]
    
    return natural_freqs, amplitudes
end

start_model_engee(m);

# Получение угловых скоростей всех инерций
omega1 = collect(simout["torsional_oscillations/FFT по скорости.1"]);

# Выделение времени и данных для анализа
time_data = omega1[:,1];
speed_data = omega1[:,2];  # Анализируем первую инерцию

# FFT анализ для определения собственных частот
frequencies, magnitude_db = analyze_natural_frequencies(time_data, speed_data);
natural_freqs, amplitudes = find_natural_frequencies(frequencies, magnitude_db);

# Вывод найденных собственных частот
println("Обнаруженные собственные частоты:")
for (i, (freq, amp)) in enumerate(zip(natural_freqs, amplitudes))
    println("$(i)-я мода: $(round(freq, digits=2)) Гц ($(round(amp .-30, digits=2)) dBm)")
end

using Plots
using Printf

function create_torsional_analysis_plots(time_data, speed_data, frequencies, magnitude_db)
    gr()
    
    # Преобразуем все данные в вещественные числа
    real_frequencies = real.(frequencies)
    real_magnitude_db = real.(magnitude_db)
    real_time_data = real.(time_data)
    real_speed_data = real.(speed_data)
    
    # График временного отклика
    p1 = plot(real_time_data, real_speed_data, 
              title="Крутильные колебания после импульса",
              xlabel="Время, с", 
              ylabel="Угловая скорость, рад/с",
              linewidth=2,
              linecolor=:blue)
    
    # График FFT спектра
    pos_idx = (real_frequencies .> 0) .& (real_frequencies .< 300)
    p2 = plot(real_frequencies[pos_idx], real_magnitude_db[pos_idx] .- 30,
              title="FFT анализ собственных частот", 
              xlabel="Частота, Гц",
              ylabel="Амплитуда, dBm",
              xlim=(0, 200),
              linewidth=2,
              linecolor=:red)
    
    # Выделение найденных собственных частот
    experimental_freqs = natural_freqs
    experimental_amps = amplitudes .- 30
    
    scatter!(p2, experimental_freqs, experimental_amps, 
            markersize=8, 
            markercolor=:green,
            markerstroke=:black,
            label="Собственные частоты")
    
    # Подписи частот
    for (i, (freq, amp)) in enumerate(zip(experimental_freqs, experimental_amps))
        annotate!(p2, freq+5, amp+3, text(@sprintf("%.1f Гц", freq), 8, :black))
    end
    
    # Объединение графиков
    combined_plot = plot(p1, p2, layout=(2,1), size=(600,800))
    
    return combined_plot
end

# Создание графиков
analysis_plots = create_torsional_analysis_plots(time_data, speed_data, frequencies, magnitude_db)

n = "torsional_oscillations_sin"

omega = 10
start_model_engee(n)
w = collect(simout["$n/Скорость.1"]);
torque = collect(simout["$n/Момент"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="Время, с", ylabel="Скорость \nвращения, рад/с", color="red")
p2 = plot(time_torque, torque_data, xlabel="Время, с", ylabel="Крутящий \nмомент, Н*м")
plot(p1,p2, layout=(2,1))

omega = 20
start_model_engee(n)
w = collect(simout["$n/Скорость.1"]);
torque = collect(simout["$n/Момент"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="Время, с", ylabel="Скорость \nвращения, рад/с", color="red")
p2 = plot(time_torque, torque_data, xlabel="Время, с", ylabel="Крутящий \nмомент, Н*м")
plot(p1,p2, layout=(2,1))

omega = 30
start_model_engee(n)
w = collect(simout["$n/Скорость.1"]);
torque = collect(simout["$n/Момент"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="Время, с", ylabel="Скорость \nвращения, рад/с", color="red")
p2 = plot(time_torque, torque_data, xlabel="Время, с", ylabel="Крутящий \nмомент, Н*м")
plot(p1,p2, layout=(2,1))

omega = natural_freqs[1] #Первая собственная частота
start_model_engee(n)
w = collect(simout["$n/Скорость.1"]);
torque = collect(simout["$n/Момент"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="Время, с", ylabel="Скорость \nвращения, рад/с", color="red")
p2 = plot(time_torque, torque_data, xlabel="Время, с", ylabel="Крутящий \nмомент, Н*м")
plot(p1,p2, layout=(2,1))

omega = 60
start_model_engee(n)
w = collect(simout["$n/Скорость.1"]);
torque = collect(simout["$n/Момент"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="Время, с", ylabel="Скорость \nвращения, рад/с", color="red")
p2 = plot(time_torque, torque_data, xlabel="Время, с", ylabel="Крутящий \nмомент, Н*м")
plot(p1,p2, layout=(2,1))