clear all;
%% Parametry
n1 = 7; %Czerski
n2 = 9; %Tyczyński
n3 = 13; %Giergielewicz
k = n3/n1;
m1 = n1;
m2 = 4*n2;
tau = 2*n2 + 3;
%% Obiekt
obiekt = tf([k], [1, m2, m1, 1], 'ioDelay', tau);
%% Impulse
[impulsY, impulsX] = impulse(obiekt);
MAXY = impulsY(find(impulsY == max(impulsY), 1));
pPX = impulsX(find(impulsY == MAXY), 1); %MAXY i MAXX to maksimum globalne
% figure(1)
% impulse(obiekt); hold on;
% plot(MAXX,MAXY,'rx');
% legend('Char. impulsowa', 'Max globalne')
%% Step
figure(2)
[styczY, ~] = step(obiekt, [MAXX-1, MAXX, MAXX+1]); %wyznaczamy punkty do wyliczenia stycznej
punktprzegieciaY = styczY(2);
styczna_a = (styczY(3)-styczY(1))/2; %współczynnik kierunkowy
[stepY, stepX] = step(obiekt);
przesuniecie = stepX(find(abs(stepX*styczna_a-punktprzegieciaY)<0.01, 1));
styczna_y = styczna_a*(stepX-(MAXX));
step(obiekt); hold on;
plot(stepX, styczna_y, 'r'); hold on;
plot(MAXX, punktprzegieciaY, 'gx')

legend('Charakterystyka skokowa', 'Styczna do pp', 'Punkt przegiecia');