clear all;
%% Parametry
n1 = 7;
n2 = 9;
n3 = 13;
k = n3/n1;
m1 = n1;
m2 = 4*n2;
tau = 2*n2 + 3;
%% Obiekt
obiekt = tf([k], [1, m2, m1, 1], 'ioDelay', tau)
%% Charakterystyka impulsowa
figure(1)
[impulsY, impulsX] = impulse(obiekt);   %charakterystyka
odcietaMaxGlob = impulsX(find(impulsY == max(impulsY), 1)); %odcięta maks.  globalnego (kopfmuller)
rzednaMaxGlob = impulsY(find(impulsY == max(impulsY), 1)); %rzędna maks. globalnego (kopfmuller)

impulse(obiekt); hold on;   %wyświetlenie charakterystyki
plot([odcietaMaxGlob, odcietaMaxGlob], [0, rzednaMaxGlob], 'r'); hold off; %naniesienie maksimum
legend('Charakterystyka skokowa', ['Maksimum globalne = ', num2str(rzednaMaxGlob)]); %legenda
%% Charakterystyka skokowa
figure(2)
[stepY, ~] = step(obiekt, [odcietaMaxGlob-1, odcietaMaxGlob, odcietaMaxGlob+1]); % wybieramy pkt z charakterystyki
rzednaPPrzegiecia = stepY(2);  % rzędna punktu przegięcia
styczna1 = (stepY(3)-stepY(1))/2;  % współczynnik kierunkowy stycznej
[stepY, stepX] = step(obiekt); % pełna charakterystyka
wartUstalonaY = stepY(end); % wartość ustalona
x0 = stepX(find(abs(stepX*styczna1-rzednaPPrzegiecia)<0.01, 1)); % odniesienie przesuniecia
styczna2 = styczna1*(stepX-(odcietaMaxGlob-x0)); % przesunieta styczna
tauZ = stX(find(abs(styczY) < 0.01, 1)); % Opóźnienie czasowe zast.
tStSt = stX(find(abs(styczY-ustY)<0.02, 1));
TZ = tStSt-tauZ; % Stała czasowa zastępcza
kZ = ustY/1;