---

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

% --- ПАРАМЕТРЫ СИСТЕМЫ ---
% W(s) = 1 / (s * (s+1)^3) = 1 / (s^4 + 3s^3 + 3s^2 + s)
num = [1];
den = [1 3 3 1 0];

% --- ВРЕМЕННЫЕ ХАРАКТЕРИСТИКИ (вручную через численное интегрирование) ---
dt = 0.05;              % Шаг времени
t = 0:dt:40;            % Время от 0 до 40 секунд
N = length(t);

% Переходная характеристика h(t) - реакция на ступеньку
% Решаем дифур: y'''' + 3y''' + 3y'' + y' = 1 (при x=1(t))
% Вводим переменные состояния: x1 = y, x2 = y', x3 = y'', x4 = y'''
% x1' = x2
% x2' = x3
% x3' = x4
% x4' = -3x4 - 3x3 - 3x2 - x1 + 1

h = zeros(1, N);
x = [0; 0; 0; 0]; % Начальные условия
for i = 1:N-1
    h(i) = x(1);
    dx1 = x(2);
    dx2 = x(3);
    dx3 = x(4);
    dx4 = -3*x(4) - 3*x(3) - 3*x(2) - x(1) + 1;
    x = x + dt * [dx1; dx2; dx3; dx4];
end
h(N) = x(1);

% Импульсная характеристика w(t) - реакция на дельта-функцию
% Решаем дифур: y'''' + 3y''' + 3y'' + y' = delta(t)
% Начальные условия: x4(0) = 1 (единичный импульс в старшую производную)
w_imp = zeros(1, N);
x = [0; 0; 0; 1]; % Начальные условия для импульса
for i = 1:N-1
    w_imp(i) = x(1);
    dx1 = x(2);
    dx2 = x(3);
    dx3 = x(4);
    dx4 = -3*x(4) - 3*x(3) - 3*x(2) - x(1);
    x = x + dt * [dx1; dx2; dx3; dx4];
end
w_imp(N) = x(1);

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

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

Re_den = -3.*w.^2 + w.^4;
Im_den = w - 3.*w.^3;

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);              % ЛАЧХ в дБ

% --- СТРОИМ ВСЕ ГРАФИКИ ---
figure('Name', 'Характеристики звена 1/(s(s+1)^3)', 'NumberTitle', 'off', 'Position', [50 50 1200 800]);

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

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

% 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)');

% Добавляем радиус-вектор для частоты 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 графиков построены.');