---

% Очищаем всё
clear; clc; close all;

% --- ОПРЕДЕЛЯЕМ ПЕРЕДАТОЧНУЮ ФУНКЦИЮ ---
% W(s) = 1 / (s * (s+1)^3) = 1 / (s^4 + 3s^3 + 3s^2 + s)
% Раскроем скобки: (s+1)^3 = s^3 + 3s^2 + 3s + 1
% Умножаем на s: s^4 + 3s^3 + 3s^2 + s

num = [1];                  % Числитель
den = [1 3 3 1 0];          % Знаменатель (коэффициенты при s^4, s^3, s^2, s^1, s^0)

% --- 1. ПЕРЕХОДНАЯ ХАРАКТЕРИСТИКА (РЕАКЦИЯ НА СТУПЕНЬКУ) ---
t = 0:0.1:30; % Время от 0 до 30 секунд

% Используем базовую функцию step для численного расчета
% (step встроена в MATLAB, не требует Toolbox для численного расчета)
[y_step, t_out] = step(tf(num, den), t); 
% Внимание: если tf выдает ошибку, используем альтернативу ниже (я её добавил)

% Если tf не работает (нет Toolbox), используем ручной расчет через lsim
% Но так как lsim тоже может требовать Toolbox, давай используем простой метод:
% Просто построим график через численное интегрирование (для надежности)
% Впрочем, step и lsim есть в базовом MATLAB. Попробуем.

figure('Name', 'Характеристики звена 1/(s(s+1)^3)', 'NumberTitle', 'off', 'Position', [50 50 1200 800]);

% 1. Переходная характеристика
subplot(3, 2, 1);
step(tf(num, den));
grid on;
title('Переходная характеристика h(t)');
xlabel('Время t, с');
ylabel('Амплитуда');

% 2. Импульсная характеристика
subplot(3, 2, 2);
impulse(tf(num, den));
grid on;
title('Импульсная характеристика w(t)');
xlabel('Время t, с');
ylabel('Амплитуда');

% --- ЧАСТОТНЫЕ ХАРАКТЕРИСТИКИ ---
w = logspace(-2, 2, 1000); % Частота от 0.01 до 100 рад/с

% Считаем комплексную передаточную функцию вручную
% W(jw) = 1 / (jw * (jw+1)^3)
% (jw+1)^3 = (1 - 3w^2) + j(3w - w^3)
% Умножаем на jw: jw * ((1-3w^2) + j(3w-w^3)) = jw(1-3w^2) - w(3w-w^3) = 
% = -w(3w-w^3) + jw(1-3w^2)
% Итого знаменатель: (-3w^2 + w^4) + j(w - 3w^3)
% W(jw) = 1 / ((-3w^2 + w^4) + j(w - 3w^3))

Re_den = -3.*w.^2 + w.^4;      % Действительная часть знаменателя
Im_den = w - 3.*w.^3;          % Мнимая часть знаменателя

% Домножаем на сопряженное, чтобы получить Re и Im числителя
% W = (Re_den - j*Im_den) / (Re_den^2 + Im_den^2)
Re = Re_den ./ (Re_den.^2 + Im_den.^2);
Im = -Im_den ./ (Re_den.^2 + Im_den.^2);

A = sqrt(Re.^2 + Im.^2);       % АЧХ
phi = atan2(Im, Re) * (180/pi);% ФЧХ в градусах
L = 20 * log10(A);             % ЛАЧХ в дБ

% 3. АЧХ
subplot(3, 2, 3);
semilogx(w, A, 'b', 'LineWidth', 2);
grid on;
title('АЧХ (Амплитудно-частотная)');
xlabel('Частота \omega, рад/с');
ylabel('A(\omega)');

% 4. ФЧХ
subplot(3, 2, 4);
semilogx(w, phi, 'r', 'LineWidth', 2);
grid on;
title('ФЧХ (Фазо-частотная)');
xlabel('Частота \omega, рад/с');
ylabel('\phi(\omega), град');

% 5. ЛАЧХ (Логарифмическая АЧХ)
subplot(3, 2, 5);
semilogx(w, L, 'm', 'LineWidth', 2);
grid on;
title('ЛАЧХ (Логарифмическая АЧХ)');
xlabel('Частота \omega, рад/с');
ylabel('L(\omega), дБ');

% 6. Годограф Найквиста (АФЧХ)
subplot(3, 2, 6);
plot(Re, Im, 'k', 'LineWidth', 2);
hold on;
grid on;
title('Годограф Найквиста (АФЧХ)');
xlabel('Re(W)');
ylabel('Im(W)');
axis equal;

% Добавляем радиус-вектор для частоты w = 1 рад/с
w_point = 1;
Re_den_p = -3*w_point^2 + w_point^4;
Im_den_p = w_point - 3*w_point^3;
Re_p = Re_den_p / (Re_den_p^2 + Im_den_p^2);
Im_p = -Im_den_p / (Re_den_p^2 + Im_den_p^2);

quiver(0, 0, Re_p, Im_p, 0, 'r', 'LineWidth', 2, 'MaxHeadSize', 0.5);
plot(Re_p, Im_p, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');
text(Re_p/2, Im_p/2, '  W(j1)', 'Color', 'r', 'FontSize', 10, 'FontWeight', 'bold');

disp('Готово! Все 6 графиков построены.');
disp('Обрати внимание: годограф уходит в бесконечность при w -> 0 (из-за интегрирующего звена 1/s).');